function[SDL]=gama(AL,SJ,r)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   SDL=gama(AL,SJ,r)                      %
%   MATLAB function finds the matrix SDL   %
%                                          %
%   SDL= [SD(s1)al1,...,SD(sL)alL]         % 
%                                          %
%   SDL dimensions are  m(r+1)xL           %
%   where SD(s)=[I,...,I.s^r]' -- m(r+1)xm %
%   degree r is a integer                  %
%   AL = [al1,...,alL] is a given matrix   %
%   AL dimensions are  mxL                 %
%   SJ is a vector containing the inter-   %
%   polation points sj (distinct)          %
%   SJ = [s1,s2,...,sL]                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   used function: 0                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   S. Pejchova, November 10th 1992        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          %   dimensions check             %
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[m,L]=size(AL);
[rSJ,cSJ]=size(SJ);

if (L~=cSJ)
   error('gama: Inconsistent dimensions of input matrices');
end

I=eye(m);
    for j=1:L
        SD=I;
        for i=1:r

          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          % Calculation SD=[SD(s1),...,SD(sL)]    %
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

          SD=[SD;I*SJ(j)^i];
        end
        if j==1
          SDL=SD*AL(:,1);
        else
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          % Create SDL=[SD(s1)al1,...,SD(sL)alL]  %
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

          SDL=[SDL,SD*AL(:,j)];
        end
    end