"Fossies" - the Fresh Open Source Software Archive

Member "cloc-1.82/tests/inputs/Lanczos.m" (3 May 2019, 3072 Bytes) of package /linux/privat/cloc-1.82.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Matlab source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file.

    1 function [phiKM, AscendingLambda] = Lanczos(K, M, sigma, Jmax)
    2 [rows,cols] = size(K);                                                       
    3 if (rows ~= cols)                                                       
    4   fprintf('Lanczos needs square matrices');                
    5 end                                                       
    6 Z       = K - sigma*M;                               % initialize some
    7 Q       = zeros(rows,Jmax+1);                        % variables
    8 T       = zeros(Jmax,Jmax);                          %
    9 rRand   = randn(rows,1);                             %
   10 rOld = rRand;
   11 betaOld = sqrt(rOld'*M*rOld);                        %
   12 for j = 1:Jmax,                                      %
   13   Q(:,j+1) = rOld/betaOld;                           %
   14   u = Z \ (M*Q(:,j+1) - Z*Q(:,j)*betaOld);           % D.S.Scott's formulation
   15   alpha = Q(:,j+1)'*M*u;                             % of the recurrence
   16   r     = u - alpha*Q(:,j+1);                        %
   17   for i=1:3,                                         %
   18     sum = zeros(rows,1);                             % repeat a full orhto-
   19     for k=2:j+1,                                     % gonalization three
   20       sum = sum + (Q(:,k)'*M*r)*Q(:,k);              % times to ensure
   21     end;                                             % high quality
   22     r = r - sum;                                     % solutions
   23   end;                                               %
   24   beta = sqrt(r'*M*r);                               %
   25   T(j,j)   = alpha;                                  %
   26   if (j ~= Jmax)                                     % augment [T] with new
   27     T(j+1,j) = beta;                                 % alpha_i, beta_i+1
   28     T(j,j+1) = beta;                                 %
   29   end;                                               %
   30   Jactual = j;                                       %
   31   if (abs(beta) < 1.0e-12)                           % singular beta; going
   32     break                                            % any more will introduce
   33   end                                                % spurious modes
   34   betaOld = beta;                                    %
   35   rOld    = r;                                       %
   36 end                                                  %
   37 [phiT,lambdaT] = eig(T(1:Jactual,1:Jactual));        % solve [T]{y} = L{y}
   38 lambdaKM   = zeros(Jactual,1);                       %
   39 for j = 1:Jactual,                                   % invert and shift the
   40   lambdaKM(j) = sigma + 1/lambdaT(j,j);              % eigenvalues to the
   41 end                                                  % user's domain
   42 [AscendingLambda, ordering] = sort(lambdaKM);        %
   43 phiKM      = zeros(rows,Jactual);                    % sort the eigenvalues
   44 UnOrdphiKM = zeros(rows,Jactual);                    % in ascending order
   45 UnOrdphiKM = Q(:,2:Jactual+1)*phiT;                  %
   46 for j = 1:Jactual,                                   % resequence the e-vectors
   47   phiKM(:,j) = UnOrdphiKM(:,ordering(j));            % to correspond to the
   48 end                                                  % e-values