## "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
```