最速下降法求特征值matlab程序

最基本的特征值问题分为三类：
1、标准的线性特征值问题：
Ax=λx,A∈Cn∗n
2、普遍的线性特征值问题：
Ax=λBx,A、B∈Cn∗n
3、普遍的艾米特正定线性特征值问题：
Ax=λBx,A、B∈Cn∗n
A∗=A,B∗=B>0∈Cn∗n

下面给出用最速下降法求特征值的matlab程序示例：
**

第一个程序：

**

% Demo routine for the Steepest Decent methods for computing the smallest %   eigenvalue of Hermitian pencil A-\lambda B with B positive definite.clear; close all;load c6h6fe1n2500.mat;%   A matrix pencil arising from solving Korn-Sham equation %       for C6H6 by SCF.%   For more details: n=size(A,1);nrmAB=[norm(A,1),norm(B,1)];          tol=1e-8;x0=randn(n,1);% SD[lam, y, res, hsty, info]=SDgS(A, B, x0, [], nrmAB, [], tol);figure;semilogy(hsty(:,2), '+r' );xlabel('iteration')ylabel('relative residual norms')title('SD')% LOCG[lam, y, res, hsty, info]=LOCGgS(A, B, x0, [], nrmAB, [], tol);figure;semilogy(hsty(:,2), '+r' );xlabel('iteration')ylabel('relative residual norms')title('LOCG')

载入的数据的形式如下：

SDgS.m和LOCGgS.m的以及其他调用的子代码如下，输出在代码注释中解释很清楚了：

function [lam, y, res, hsty, info]=SDgS(A, B, x, BU, nrmAB, shift, tol)%%  Steepest decent method to compute the kth smallest eigenvalue of A-lamda B,%  where A and B are Hermitian, and B is positive definite, k=k0+1 and k0 is the number %  of columns in BU=B*U.%%  Deflation is done by shifting away known eigenvalues. Basically SD on%%      A+shift*B*U*(B*U)'=A+shift*BU*BU', B%% %%---------------------------------------------------------------%%  Input%%        A     n-by-n Hermitian matrix%        B     n-by-n Hermitian matrix and positive definite%        x     n-vector, initial guess to the kth eigenvector%        BU    n-by-k0, B*U, where U's columns span an approximate invariant subspace %              (accurate enough, however) associated with 1st to k0th smallest %              eigenvalues. Assume U'*B*U=I.%              Possibly k0=0, i.e., BU is empty array.%        nrmAB 2-by-1, nrmAB(1): estimate of ||A||%                      nrmAB(2): estimate of ||B||%        shift real for shifting away known eigenvalues%        tol   tolerance for testing convergence%%  Output%%        lam  converged eigenvalue%        y     corresponding eigenvector%        res   residual error  ||A y - lam B*y||%        hsty  *-by-2 %              hsty(:,1) -- eigenvalue approximations%              hsty(:,2) -- normalized residuals ||A yi - lam_i B*yi||/(||A||*||y_i||+|lam_i|*||B||*||y_i||)%        info  number of iterations% %---------------------------------------------------------------%n=max(size(A));maxitn=round(0.2*n); %%%%%%%%%%%%%%%%%%%%%%%%%0.2k0=size(BU,2); opts.Line=1;Bx=B*x; xBx=real(x'*Bx); nrmx=sqrt(xBx);x=(1/nrmx)*x; Bx=(1/nrmx)*Bx;if k0==0,   Ax=A*x; rho=x'*Ax;    r=Ax-rho*Bx;      err=norm(r)/( nrmAB(1)+abs(rho)*nrmAB(2) ); itn=0;   hsty=[rho err];   while err > tol & itn < maxitn,           Br=B*r; rBr=real(r'*Br); nrmr=sqrt(rBr); r=(1/nrmr)*r;       rhos=[rho r'*(A*r)];       opts.rhos=rhos;       [y, t, lam]=LNSRCHg(Ax,Bx,x,r,opts);       x=y;        Bx=B*x; xBx=real(x'*Bx); nrmx=sqrt(xBx);       x=(1/nrmx)*x; Bx=(1/nrmx)*Bx;       Ax=A*x; rho=x'*Ax;        r=Ax-rho*Bx;        err=norm(r)/( nrmAB(1)+abs(rho)*nrmAB(2) );       hsty=[hsty; rho err];       itn = itn+1;   endelse  % k0>0     Ax=A*x+shift*(BU*(BU'*x)); rho=x'*Ax;    r=Ax-rho*Bx;      err=norm(r)/( nrmAB(1)+abs(rho)*nrmAB(2) ); itn=0;   hsty=[rho err];   while err > tol & itn < maxitn,           Br=B*r; rBr=real(r'*Br); nrmr=sqrt(rBr); r=(1/nrmr)*r;       Ar=A*r+shift*(BU*(BU'*r));       rhos=[rho r'*(Ar)];       opts.rhos=rhos;       [y, t, lam]=LNSRCHg(Ax,Bx,x,r,opts);       x=y;        Bx=B*x; xBx=real(x'*Bx); nrmx=sqrt(xBx);       x=(1/nrmx)*x; Bx=(1/nrmx)*Bx;       Ax=A*x+shift*(BU*(BU'*x)); rho=x'*Ax;        r=Ax-rho*Bx;        err=norm(r)/( nrmAB(1)+abs(rho)*nrmAB(2) );       hsty=[hsty; rho err];       itn = itn+1;   endendinfo = itn; y=x; lam=rho; res=nrmr;function [lam, y, res, hsty, info]=LOCGgS(A, B, x, BU, nrmAB, shift, tol)%%  Locally Optimal Conjugate Gradient (LOCG) method to compute the kth smallest eigenvalue of A-lamda B,%  where A and B are Hermitian, and B is positive definite, k=k0+1 and k0 is the number %  of columns in U.%%  Deflation is done by shifting away known eigenvalues. Basically LOBCG on%%      A+shift*B*U*(B*U)'=A+shift*BU*BU', B%%---------------------------------------------------------------%%  Input%%        A     n-by-n Hermitian matrix%        B     n-by-n Hermitian matrix and positive definite%        x     n-vector, initial guess to the kth eigenvector%              associated with (k0+1)st to (k0+k)th smallest eigenvalues%        BU    n-by-k0, B*U, where U's columns span an approximate invariant subspace %              (accurate enough, however) associated with 1st to k0th smallest %              eigenvalues. Assume U'*B*U=I.%              Possibly k0=0, i.e., U is empty array.%        nrmAB 2-by-1, nrmAB(1): estimate of ||A||%                      nrmAB(2): estimate of ||B||%        shift real for shifting away known eigenvalues%        tol   tolerance for testing convergence%%  Output%%        lam   computed eigenvalues%        Y     corresponding eigenvectors%        res   residual error  ||A y - lam B y||%        hsty  *-by-2 %              hsty(:,1) -- eigenvalue approximations%              hsty(:,2) -- normalized residuals ||A yi - lam_i B*yi||/(||A||*||y_i||+|lam_i|*||B||*||y_i||)%        info  struct%              info.itn number of iterations%% %---------------------------------------------------------------%n=max(size(A)); k0=size(BU,2);maxitn=min(round(0.1*n),400);%0.1Bx=B*x; xBx=real(x'*Bx); nrmx=sqrt(xBx);x0=(1/nrmx)*x;  Bx=(1/nrmx)*Bx;if k0==0,   Ax=A*x0; rho=real(x0'*Ax);   r0=Ax-rho*Bx;    nrmr0=norm(r0);   err=nrmr0/(nrmAB(1)+abs(rho)*nrmAB(2));    hsty=[rho err];   itn=0;   Q=r0-x0*(Bx'*r0);  % now Q'*B*x0=0;   BQ=B*Q; pBp=real(Q'*(BQ)); nrmp=sqrt(pBp);  Q=(1/nrmp)*Q;%   disp('Here 1');%   disp(norm([x0 Q]'*B*[x0 Q]-eye(2),1));   AQ=A*Q; tmp=x0'*AQ;   Rho1=[rho     tmp;       conj(tmp) real(Q'*AQ)];   [v,rho]=mineig(Rho1,1);   x1=[x0,Q]*v; % x1'*B*x1 approx 1      y=Q;    Bx=B*x1; xBx=real(x1'*Bx); nrmx=sqrt(xBx);   x1=(1/nrmx)*x1;  Bx=(1/nrmx)*Bx;   Ax=[Ax,AQ]*(v/nrmx); % A*x1;   %Ax=A*x1;   r=Ax-rho*Bx;    while err > tol & itn < maxitn,      Q=MGSg(B,[y, r], x1);%      disp('Here 2')%      disp(norm([x1 Q]'*B*[x1 Q]-eye(3),1));      AQ=A*Q; tmp=x1'*AQ;      Rho1=[rho     tmp;           tmp'   Q'*AQ];      [v,rho]=mineig(Rho1,1);      x0=x1;      x1=[x0,Q]*v; % x1'*B*x1 approx 1         y=Q*v(2:3);       Bx=B*x1; xBx=real(x1'*Bx); nrmx=sqrt(xBx);      x1=(1/nrmx)*x1;  Bx=(1/nrmx)*Bx;      Ax=[Ax,AQ]*(v/nrmx); % A*x1;      %Ax=A*x1;      r=Ax-rho*Bx;       nrmr=norm(r);      err=nrmr/(nrmAB(1)+abs(rho)*nrmAB(2));      hsty=[hsty; rho err];      itn = itn+1;   endelse  % k0>0   Ax=A*x0+shift*(BU*(BU'*x0));    rho=real(x0'*Ax);   r0=Ax-rho*Bx;    nrmr0=norm(r0);   err=nrmr0/(nrmAB(1)+abs(rho)*nrmAB(2));    hsty=[rho err];   itn=0;   Q=r0-x0*(Bx'*r0);  % now Q'*B*x0=0;   BQ=B*Q; pBp=real(Q'*(BQ)); nrmp=sqrt(pBp);  Q=(1/nrmp)*Q;%   disp('Here 1');%   disp(norm([x0 Q]'*B*[x0 Q]-eye(2),1));   AQ=A*Q+shift*(BU*(BU'*Q)); tmp=x0'*AQ;   Rho1=[rho     tmp;       conj(tmp) real(Q'*AQ)];   [v,rho]=mineig(Rho1,1);   x1=[x0,Q]*v; % x1'*B*x1 approx 1      y=Q;    Bx=B*x1; xBx=real(x1'*Bx); nrmx=sqrt(xBx);   x1=(1/nrmx)*x1;  Bx=(1/nrmx)*Bx;   Ax=[Ax,AQ]*(v/nrmx); % A*x1+shift*(BU*(BU'*x1));   %Ax=A*x1+shift*(BU*(BU'*x1));   r=Ax-rho*Bx;    while err > tol & itn < maxitn,      Q=MGSg(B,[y, r], x1);%      disp('Here 2')%      disp(norm([x1 Q]'*B*[x1 Q]-eye(3),1));      AQ=A*Q+shift*(BU*(BU'*Q)); tmp=x1'*AQ;      Rho1=[rho     tmp;           tmp'   Q'*AQ];      [v,rho]=mineig(Rho1,1);      x0=x1;      x1=[x0,Q]*v; % x1'*B*x1 approx 1         y=Q*v(2:3);       Bx=B*x1; xBx=real(x1'*Bx); nrmx=sqrt(xBx);      x1=(1/nrmx)*x1;  Bx=(1/nrmx)*Bx;      Ax=[Ax,AQ]*(v/nrmx); % A*x1+shift*(BU*(BU'*x1));      %Ax=A*x1+shift*(BU*(BU'*x1));      r=Ax-rho*Bx;       nrmr=norm(r);      err=nrmr/(nrmAB(1)+abs(rho)*nrmAB(2));      hsty=[hsty; rho err];      itn = itn+1;   endendinfo = itn; y=x1; lam=rho; res=err;其他子代码： function  V=MGSg(B,X,U)%%---------------------------------------------------------------%  Modified Gram-Schmit on X with selective re-orthogonalization%  with respect to B-inner product%---------------------------------------------------------------%%  Input:%%      B      (n-by-n) Hermitian and positive definite%      X      (n-by-k)%      U      (n-by-k0) U'*B*U=I already. Optional argument%             Possibly k0=0, i.e., U is empty array%%  Output:%%      V      n-by-nV  B-Orthogonalized vectors from columns of X%%   %---------------------------------------------------------------[n,k]=size(X);rtol_re=1.0e-4; % relative tolerance to perform reorthogonalization%V = zeros(n,k); nV=1;if nargin==2,   nrm2 = sqrt(real(X(:,1)'*B*X(:,1)));    if nrm2>0,      V(:,nV)=X(:,1)/nrm2; nV=nV+1;   end   for j=2:k,     vh = X(:,j); nrm2 = sqrt(real(vh'*B*vh));     tolorth=rtol_re*nrm2;     %  by MGS     for i=1:nV-1,           vh = vh - V(:,i)*( V(:,i)'*B*vh );     end     nrm2=sqrt(real(vh'*B*vh));     %  >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>     %  Perform re-orthogonalization once by MGS when deemed necessary.     if nrm2 <= tolorth        for i=1:nV-1,           vh = vh - V(:,i)*( V(:,i)'*B*vh );        end        nrm2=sqrt(real(vh'*B*vh));     end     if nrm2>0        V(:,nV)=vh/nrm2; %R(j,j)=nrm2;        nV=nV+1;     end    endelseif nargin==3,   k0=size(U,2);   vh = X(:,1); nrm2 = sqrt(real(vh'*B*vh));   tolorth=rtol_re*nrm2;   % by MGS   for i=1:k0       vh = vh - U(:,i)*( ( U(:,i) )'*B* vh );   end   nrm2=sqrt(real(vh'*B*vh));   %  >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   %  Perform re-orthogonalization once by MGS when deemed necessary.   if nrm2 < tolorth          for i=1:k0          vh = vh - U(:,i)*( ( U(:,i) )'*B* vh );      end      nrm2=sqrt(real(vh'*B*vh));    end   if nrm2 > 0,      V(:,nV)=vh/nrm2; nV=nV+1;   end    for j=2:k,     vh = X(:,j); nrm2 = sqrt(real(vh'*B*vh));     tolorth=rtol_re*nrm2;     %  by MGS     for i=1:k0         vh = vh - U(:,i)*( ( U(:,i) )'*B* vh );     end     for i=1:nV-1,           vh = vh - V(:,i)*( V(:,i)'*B*vh );     end     nrm2=sqrt(real(vh'*B*vh));     %  >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>     %  Perform re-orthogonalization once by MGS when deemed necessary.     if nrm2 <= tolorth        for i=1:k0            vh = vh - U(:,i)*( ( U(:,i) )'*B* vh );        end        for i=1:nV-1,           vh = vh - V(:,i)*( V(:,i)'*B*vh );        end        nrm2=sqrt(real(vh'*B*vh));     end     if nrm2>0        V(:,nV)=vh/nrm2; %R(j,j)=nrm2;        nV=nV+1;     end   endendV=V(:,1:nV-1);function [y, t, lam]=LNSRCHg(Ax,Bx,x,p,opts)%%  It perform line search to solve%%            inf_t rho(x+t*p),%%  where%%       rho(x)=x' A x/(x' B x), Rayleigh quotient of a Hermitian matrix pencil A-lamda B%       with B positive definite%%%---------------------------------------------------------------%%  Input%%        Ax     n-by-1; it is A*x%        Bx     n-by-1; it is B*x%        x      n-vector, x'*B*x=1%        p      n-vector, p'*B*p=1%        opts   struct%               opts.rhos = [rho(x)  rho(p)]   %               opts.Line = 0,  through diff(rho(x+t*p),t)=0%                         = 1,  through solving [x,p]'*(A-lamda B)*[x,p]%%  Output%%        y     x+t*p, in case t=infty, y=p instead%        t     optimal argument%              = 0    inf_t rho(x+t*p)=rho(x):      no improvement%        lam  new approximate eigenvalue%%AA=[x'; p']*A*[x p];%if opts.rhos  == 0%   rhos=[x'*Ax  p'*A*p]; % because x'*B*x=p'*B*p=1 upon entry!%else   rhos=opts.rhos;%endtmp=Ax'*p;AA=[rhos(1)  tmp;    conj(tmp) rhos(2)];%BB=[x'; p']*B*[x p];tmp=Bx'*p;BB=[ 1    tmp;    conj(tmp) 1];%disp(cond(BB))% solve eig(AA,BB)tmp=BB(1,2)*BB(2,1); a=1.0-tmp;   if abs(a) <= 100*eps*(1.0+abs(tmp))     % BB is singular => x and p parallel since B is definite => no improvement   t = 0; y=x;     lam=rhos(1);   returnelse   b=-AA(1,1)-AA(2,2)+AA(1,2)*BB(2,1)+BB(1,2)*AA(2,1);   c=AA(1,1)*AA(2,2)-AA(1,2)*AA(2,1);   tmp = b*b-4*a*c;   % smallest eigebvalue lam is   if a>0      if b >= 0,         lam = (-b-sqrt(abs(tmp)))/(2*a);      else         lam = 2*c/(-b+sqrt(abs(tmp)));      end   else  % a>0      if b <= 0,         lam = (-b+sqrt(abs(tmp)))/(2*a);      else         lam = -2*c/(b+sqrt(abs(tmp)));      end   endendif opts.Line==0 % through diff(rho(x+t*p),t)=0   % a=AA(2,2)*(BB(1,2)+BB(2,1))-(AA(1,2)+AA(2,1))*BB(2,2);   tmp1 = AA(2,2)*(BB(1,2)+BB(2,1)); tmp2 = (AA(1,2)+AA(2,1))*BB(2,2);   a=tmp1-tmp2;   if abs(a) <= 4*eps*(abs(tmp1)+abs(tmp2))      a = 0;   end   % b=2*( AA(2,2)*BB(1,1)-AA(1,1)*BB(2,2) );   tmp3 = AA(2,2); tmp4 = AA(1,1);   b = 2*(tmp3-tmp4);   if abs(b) <= 8*eps*(abs(tmp3)+abs(tmp4))      b = 0;   end   % c=(AA(1,2)+AA(2,1))*BB(1,1)-AA(1,1)*(BB(1,2)+BB(2,1));   tmp5 = (AA(1,2)+AA(2,1)); tmp6 = AA(1,1)*(BB(1,2)+BB(2,1));   c = tmp5 - tmp6;   if abs(c) <= 4*eps*(abs(tmp5)+abs(tmp6))      c = 0;   end   if a ~= 0,      tmp = b*b-4*a*c;      if b <= 0,         t = (-b+sqrt(abs(tmp)))/(2*a);      else         t = -2*c/(b+sqrt(abs(tmp)));      end      y=x+t*p;   else  % a == 0      if  b == 0 % a = b == 0, too          y = x; t = 0;      % no improvement      elseif b < 0   % a = 0, b < 0          t=1/0; y=p;      else % a = 0, b>0         t=-c/b; y= x+t*p;      end   endelse % opts.Line==1, through solving [x,p]'*(A-lamda B)*[x,p]   y=[-(AA(1,2)-lam*BB(1,2)); AA(1,1)-lam];   % corresponding eigenvector y=[x,p]*y   if abs(y(1))>0.0      t=y(2)/y(1); y=x+t*p;   else      t=1/0.0; y=p;   endendfunction [V,D]=mineig(A,k)%% This function return k smallest eigenvalues of a Hermitian matrix A%%%---------------------------------------------------------------%%  Input%%        A     n-by-n; Hermitian%        k     integer: k<=n%              number of smallest eigenvalues asked%%  Output%%        V     n-by-k, eigenvector matrix with associated eigenvalues in D%        D     k-by-1, eigenvalues D(1) <= ... <= D(k)%%---------------------------------------------------------------n=size(A,1);if k>n   disp('mineig: too many eigenvalues asked');   returnend%[V,D]=eig(A);[V,D]=schur(A);D=real(diag(D)); [D,idx]=sort(D,'ascend');D=D(1:k); V=V(:,idx(1:k));err=norm(V'*V-eye(k),1);if err>=1000*n*eps   disp(strcat('mineig: returned eigenvectors are not orthonormal; err=',num2str(err)));   %[n norm(A'-A,1)]end

**

第二个程序：

**

% Demo routine for the Preconditioned Steepest Descent methods for %   computing the smallest eigenvalue of Hermitian pencil A-\lambda B %   with B positive definite.clear; close all;A = mmread('bcsstk11.mtx');%mmread以系数矩阵的方式读取数据B = mmread('bcsstm11.mtx');%   A matrix pencil bcsstk11 and bcsstm11 from Matrix Market:n=size(A,1);nrmAB=[norm(A,1),norm(B,1)]; %1范数就是绝对值按列求和取最大tol=1e-8;x0=randn(n,1);%产生一个随机列向量% Steepest Descent [lam, y, res, hsty, info]=SDgS(A, B, x0, [], nrmAB, [], tol);% function [lam, y, res, hsty, info]=SDgS(A, B, x, BU, nrmAB, shift, tol)% Steepest decent method to compute the kth smallest eigenvalue of A-lamda B,%    where A and B are Hermitian, and B is positive definite, k=k0+1 and k0 is the number %    of columns in BU=B*U.%    Deflation is done by shifting away known eigenvalues. Basically SD on%        A+shift*B*U*(B*U)'=A+shift*BU*BU', B%   Copyright by Song Lu, 7/14/2016%  ---------------------------------------------------------------%    Input%          A     n-by-n Hermitian matrix%          B     n-by-n Hermitian matrix and positive definite%          x     n-vector, initial guess to the kth eigenvector%          BU    n-by-k0, B*U, where U's columns span an approximate invariant subspace %                (accurate enough, however) associated with 1st to k0th smallest %                eigenvalues. Assume U'*B*U=I.%                Possibly k0=0, i.e., BU is empty array.%          nrmAB 2-by-1, nrmAB(1): estimate of ||A||%                        nrmAB(2): estimate of ||B||%          shift real for shifting away known eigenvalues%          tol   tolerance for testing convergence%  %    Output%  %          lam  converged eigenvalue%          y     corresponding eigenvector%          res   residual error  ||A y - lam B*y||%          hsty  *-by-2 %                hsty(:,1) -- eigenvalue approximations 每次迭代的特征值的变化%                hsty(:,2) -- normalized residuals ||A yi - lam_i B*yi||/(||A||*||y_i||+|lam_i|*||B||*||y_i||)%          info  number of iterations%   %  ---------------------------------------------------------------figure;semilogy(hsty(:,2), '+r' );%因为数据都集中在底部，直接plot不美观，我们改变一下y%轴的显示方式，即把数据按log10(y)等距显示，只是在显示上。即以10^m次方为一个刻度显示xlabel('iteration')ylabel('relative residual norms')title('SD')% Preconditioned SDopts.precond=1;opts.precond_one=1;opts.precond_sgm=0;[lam, y, res, hsty, info]=PSDgS(A, B, x0, [], nrmAB, [], tol, opts);figure;semilogy(hsty(:,2), '+r' );xlabel('iteration')ylabel('relative residual norms')title('PSD')% Locally optimal CG[lam, y, res, hsty, info]=LOCGgS(A, B, x0, [], nrmAB, [], tol);figure;semilogy(hsty(:,2), '+r' );xlabel('iteration')ylabel('relative residual norms')title('LOCG')% Preconditioned LOPCGopts.precond=1;opts.precond_one=1;opts.precond_sgm=0;[lam, y, res, hsty, info]=LOPCGgS(A, B, x0, [], nrmAB, [], tol, opts);figure;semilogy(hsty(:,2), '+r' );xlabel('iteration')ylabel('relative residual norms')title('LOPCG')

调用的子函数：function  [A,rows,cols,entries,rep,field,symm] = mmread(filename)%% function  [A] = mmread(filename)%% function  [A,rows,cols,entries,rep,field,symm] = mmread(filename)%%      Reads the contents of the Matrix Market file 'filename'%      into the matrix 'A'.  'A' will be either sparse or full,%      depending on the Matrix Market format indicated by%      'coordinate' (coordinate sparse storage), or%      'array' (dense array storage).  The data will be duplicated%      as appropriate if symmetry is indicated in the header.%%      Optionally, size information about the matrix can be %      obtained by using the return values rows, cols, and%      entries, where entries is the number of nonzero entries%      in the final matrix. Type information can also be retrieved%      using the optional return values rep (representation), field,%      and symm (symmetry).%mmfile = fopen(filename,'r');if ( mmfile == -1 ) disp(filename); error('File not found');end;header = fgets(mmfile);if (header == -1 )  error('Empty file.')end% NOTE: If using a version of Matlab for which strtok is not%       defined, substitute 'gettok' for 'strtok' in the %       following lines, and download gettok.m from the%       Matrix Market site.    [head0,header]   = strtok(header);  % see note above[head1,header]   = strtok(header);[rep,header]     = strtok(header);[field,header]   = strtok(header);[symm,header]    = strtok(header);head1 = lower(head1);rep   = lower(rep);field = lower(field);symm  = lower(symm);if ( length(symm) == 0 )   disp(['Not enough words in header line of file ',filename])    disp('Recognized format: ')   disp('%%MatrixMarket matrix representation field symmetry')   error('Check header line.')endif ( ~ strcmp(head0,'%%MatrixMarket') )   error('Not a valid MatrixMarket header.')endif (  ~ strcmp(head1,'matrix') )   disp(['This seems to be a MatrixMarket ',head1,' file.']);   disp('This function only knows how to read MatrixMarket matrix files.');   disp('  ');   error('  ');end% Read through comments, ignoring themcommentline = fgets(mmfile);while length(commentline) > 0 & commentline(1) == '%',  commentline = fgets(mmfile);end% Read size information, then branch according to% sparse or dense formatif ( strcmp(rep,'coordinate')) %  read matrix given in sparse                               %  coordinate matrix format  [sizeinfo,count] = sscanf(commentline,'%d%d%d');  while ( count == 0 )     commentline =  fgets(mmfile);     if (commentline == -1 )       error('End-of-file reached before size information was found.')     end     [sizeinfo,count] = sscanf(commentline,'%d%d%d');     if ( count > 0 & count ~= 3 )       error('Invalid size specification line.')     end  end  rows = sizeinfo(1);  cols = sizeinfo(2);  entries = sizeinfo(3);  if  ( strcmp(field,'real') )               % real valued entries:    [T,count] = fscanf(mmfile,'%f',3);    T = [T; fscanf(mmfile,'%f')];    if ( size(T) ~= 3*entries )       message = ...       str2mat('Data file does not contain expected amount of data.',...               'Check that number of data lines matches nonzero count.');       disp(message);       error('Invalid data.');    end    T = reshape(T,3,entries)';    A = sparse(T(:,1), T(:,2), T(:,3), rows , cols);  elseif   ( strcmp(field,'complex'))            % complex valued entries:    T = fscanf(mmfile,'%f',4);    T = [T; fscanf(mmfile,'%f')];    if ( size(T) ~= 4*entries )       message = ...       str2mat('Data file does not contain expected amount of data.',...               'Check that number of data lines matches nonzero count.');       disp(message);       error('Invalid data.');    end    T = reshape(T,4,entries)';    A = sparse(T(:,1), T(:,2), T(:,3) + T(:,4)*sqrt(-1), rows , cols);  elseif  ( strcmp(field,'pattern'))    % pattern matrix (no values given):    T = fscanf(mmfile,'%f',2);    T = [T; fscanf(mmfile,'%f')];    if ( size(T) ~= 2*entries )       message = ...       str2mat('Data file does not contain expected amount of data.',...               'Check that number of data lines matches nonzero count.');       disp(message);       error('Invalid data.');    end    T = reshape(T,2,entries)';    A = sparse(T(:,1), T(:,2), ones(entries,1) , rows , cols);  endelseif ( strcmp(rep,'array') ) %  read matrix given in dense                                %  array (column major) format  [sizeinfo,count] = sscanf(commentline,'%d%d');  while ( count == 0 )     commentline =  fgets(mmfile);     if (commentline == -1 )       error('End-of-file reached before size information was found.')     end     [sizeinfo,count] = sscanf(commentline,'%d%d');     if ( count > 0 & count ~= 2 )       error('Invalid size specification line.')     end  end  rows = sizeinfo(1);  cols = sizeinfo(2);  entries = rows*cols;  if  ( strcmp(field,'real') )               % real valued entries:    A = fscanf(mmfile,'%f',1);    A = [A; fscanf(mmfile,'%f')];    if ( strcmp(symm,'symmetric') | strcmp(symm,'hermitian') | strcmp(symm,'skew-symmetric') )       for j=1:cols-1,        currenti = j*rows;        A = [A(1:currenti); zeros(j,1);A(currenti+1:length(A))];      end    elseif ( ~ strcmp(symm,'general') )      disp('Unrecognized symmetry')      disp(symm)      disp('Recognized choices:')      disp('   symmetric')      disp('   hermitian')      disp('   skew-symmetric')      disp('   general')      error('Check symmetry specification in header.');    end      A = reshape(A,rows,cols);  elseif  ( strcmp(field,'complex'))         % complx valued entries:    tmpr = fscanf(mmfile,'%f',1);    tmpi = fscanf(mmfile,'%f',1);    A  = tmpr+tmpi*i;    for j=1:entries-1      tmpr = fscanf(mmfile,'%f',1);      tmpi = fscanf(mmfile,'%f',1);      A  = [A; tmpr + tmpi*i];    end    if ( strcmp(symm,'symmetric') | strcmp(symm,'hermitian') | strcmp(symm,'skew-symmetric') )       for j=1:cols-1,        currenti = j*rows;        A = [A(1:currenti); zeros(j,1);A(currenti+1:length(A))];      end    elseif ( ~ strcmp(symm,'general') )      disp('Unrecognized symmetry')      disp(symm)      disp('Recognized choices:')      disp('   symmetric')      disp('   hermitian')      disp('   skew-symmetric')      disp('   general')      error('Check symmetry specification in header.');    end    A = reshape(A,rows,cols);  elseif  ( strcmp(field,'pattern'))    % pattern (makes no sense for dense)   disp('Matrix type:',field)   error('Pattern matrix type invalid for array storage format.');  else                                 % Unknown matrix type   disp('Matrix type:',field)   error('Invalid matrix type specification. Check header against MM documentation.');  endend%% If symmetric, skew-symmetric or Hermitian, duplicate lower% triangular part and modify entries as appropriate:%if ( strcmp(symm,'symmetric') )   A = A + A.' - diag(diag(A));   entries = nnz(A);elseif ( strcmp(symm,'hermitian') )   A = A + A' - diag(diag(A));   entries = nnz(A);elseif ( strcmp(symm,'skew-symmetric') )   A = A - A';   entries = nnz(A);endfclose(mmfile);% Done.function [lam, y, res, hsty, info]=LOPCGgS(A, B, x, BU, nrmAB, shift, tol, opts)%%  Locally Optimal Preconditioned Conjugate Gradient (LOCG) method to compute the kth smallest %  eigenvalue of A-lamda B, where A and B are Hermitian, and B is positive definite, k=k0+1 and k0 %  is the number of columns in U.%%  Deflation is done by shifting away known eigenvalues. Basically LOBCG on%%      A+shift*B*U*(B*U)'=A+shift*BU*BU', B%%---------------------------------------------------------------%%  Input%%        A     n-by-n Hermitian matrix%        B     n-by-n Hermitian matrix and positive definite%        x     n-vector, initial guess to the kth eigenvector%              associated with (k0+1)st to (k0+k)th smallest eigenvalues%        BU    n-by-k0, B*U, where U's columns span an approximate invariant subspace %              (accurate enough, however) associated with 1st to k0th smallest %              eigenvalues. Assume U'*B*U=I.%              Possibly k0=0, i.e., U is empty array.%        nrmAB 2-by-1, nrmAB(1): estimate of ||A||%                      nrmAB(2): estimate of ||B||%        shift real for shifting away known eigenvalues%        tol   tolerance for testing convergence%        opts  options ...%              opts.precond=1,  use 1st type preconditioner K=(A-sgm B)^{-1}%                   1) opts.precond_one=1: direct solver%                   2) opts.precond_one=2: A-sgm B is (effectively) SPD. K is implemented by CG   %                   3) opts.precond_one=3: A-sgm B is indefinite. K is implemented by MINRES %              opts.precond=2,  use 2nd type preconditioner; suppose A-sgm B approx LDL'.%                   1) opts.precond_two=1: K=(LL')^{-1}%                   2) opts.precond_two=2: K=[LL'+shift*B*U*(B*U)]^{-1}%%              opts.precond_sgm  the shift parameter to build a preconditioner%%  Output%%        lam   computed eigenvalues%        Y     corresponding eigenvectors%        res   residual error  ||A y - lam B y||%        hsty  *-by-2 %              hsty(:,1) -- eigenvalue approximations%              hsty(:,2) -- normalized residuals ||A yi - lam_i B*yi||/(||A||*||y_i||+|lam_i|*||B||*||y_i||)%        info  struct%              info.itn number of iterations%% %---------------------------------------------------------------%n=max(size(A)); k0=size(BU,2);maxitn=min(round(0.1*n),400);sgm=opts.precond_sgm;C=A-sgm*B; if opts.precond==1   if opts.precond_one==1      [L,U,pp]=lu(C,'vector');   elseif opts.precond_one==2      CGopts.nitn=10;      CGopts.tol=1e-2;      CGopts.met=2;      CG_M=0;   elseif opts.precond_one==3      % TBI   endelseif opts.precond==2   iLUopts.thresh = 0;   iLUopts.udiag = 1;   iLUopts.milu = 1;   [L, U] = luinc(C, iLUopts);      %  Scale diagonals to get L.   d = diag (U);   d = sqrt(abs(d));   for i = 1:n      if (d(i) < 1e-2)         d(i) = 1.0;      end   end   L = L * spdiags(d, 0, n, n);end  Bx=B*x; xBx=real(x'*Bx); nrmx=sqrt(xBx);x0=(1/nrmx)*x;  Bx=(1/nrmx)*Bx;if k0==0,   if opts.precond==1      if opts.precond_one==1         % no action      elseif opts.precond_one==2         CGopts.update = 0;      elseif opts.precond_one==3             % TBI      end   elseif opts.precond==2      % no action   end       Ax=A*x0; rho=real(x0'*Ax);   r0=Ax-rho*Bx;    nrmr0=norm(r0);   err=nrmr0/(nrmAB(1)+abs(rho)*nrmAB(2));    hsty=[rho err];   itn=0;   if opts.precond==1      if opts.precond_one==1         % C*Kr=r0         Kr=U\(L\r0(pp));      elseif opts.precond_one==2         CGx0=zeros(n,1);         [Kr, error, iter, flag] = LinCG(C, CGx0, r0, CG_M, CGopts);       elseif opts.precond_one==3         % TBI      end   elseif opts.precond==2      % (L*L')*Kr=r0      Kr=(L')\(L\r0);   end %   p0=-Kr;   Q=Kr-x0*(Bx'*Kr);  % now Q'*B*x0=0;   BQ=B*Q; pBp=real(Q'*(BQ)); nrmp=sqrt(pBp);  Q=(1/nrmp)*Q;%   disp('Here 1');%   disp(norm([x0 Q]'*B*[x0 Q]-eye(2),1));   AQ=A*Q; tmp=x0'*AQ;   Rho=[rho     tmp;       conj(tmp) real(Q'*AQ)];   [v,rho]=mineig(Rho,1);   x1=[x0,Q]*v; % x1'*B*x1 approx 1      y=Q;    Bx=B*x1; xBx=real(x1'*Bx); nrmx=sqrt(xBx);   x1=(1/nrmx)*x1;  Bx=(1/nrmx)*Bx;   Ax=[Ax,AQ]*(v/nrmx); % A*x1;   %Ax=A*x1;   r=Ax-rho*Bx;    while err > tol & itn < maxitn,         if opts.precond==1         if opts.precond_one==1            % C*Kr=r            Kr=U\(L\r(pp));         elseif opts.precond_one==2            CGx0=zeros(n,1);            [Kr, error, iter, flag] = LinCG(C, CGx0, r, CG_M, CGopts);          elseif opts.precond_one==3            % TBI         end      elseif opts.precond==2         % (L*L')*Kr=r         Kr=(L')\(L\r);      end       Q=MGSg(B,[y, Kr], x1);%      disp('Here 2')%      disp(norm([x1 Q]'*B*[x1 Q]-eye(3),1));      AQ=A*Q; tmp=x1'*AQ;      Rho=[rho     tmp;           tmp'   Q'*AQ];      [v,rho]=mineig(Rho,1);      x0=x1;      x1=[x0,Q]*v; % x1'*B*x1 approx 1         y=Q*v(2:3);       Bx=B*x1; xBx=real(x1'*Bx); nrmx=sqrt(xBx);      x1=(1/nrmx)*x1;  Bx=(1/nrmx)*Bx;      Ax=[Ax,AQ]*(v/nrmx); % A*x1;      %Ax=A*x1;      r=Ax-rho*Bx;       nrmr=norm(r);      err=nrmr/(nrmAB(1)+abs(rho)*nrmAB(2));      hsty=[hsty; rho err];      itn = itn+1;   endelse  % k0>0   if opts.precond==1      if opts.precond_one==1         % [C+shift*V*V']^{-1} = C^{-1}-shift*C^{-1}*V*[I+shift*V'*C^{-1}*V]^{-1}*V'*C^{-1}, where V=BU         CiV=U\(L\BU(pp,:));         T=eye(k0)+shift*BU'*CiV; invT=inv(T);         CiViT=CiV*(shift*invT);      elseif opts.precond_one==2         CGopts.update = 1;         CGopts.shift =shift;         CGopts.V = BU;      elseif opts.precond_one==3             % TBI      end   elseif opts.precond==2      if opts.precond_two==2         % [C+shift*V*V']^{-1} = C^{-1}-shift*C^{-1}*V*[I+shift*V'*C^{-1}*V]^{-1}*V'*C^{-1}, where V=BU, C approx LL'         CiV=(L')\(L\BU);         T=eye(k0)+shift*BU'*CiV; invT=inv(T);         CiViT=CiV*(shift*invT);      end   end       Ax=A*x0+shift*(BU*(BU'*x0));    rho=real(x0'*Ax);   r0=Ax-rho*Bx;    nrmr0=norm(r0);   err=nrmr0/(nrmAB(1)+abs(rho)*nrmAB(2));    hsty=[rho err];   itn=0;   if opts.precond==1       if opts.precond_one==1          % [C+shift*BU*(BU)']*Kr=r0          Cinvr=U\(L\r0(pp));          Kr=Cinvr-CiViT*(CiV'*r0);       elseif opts.precond_one==2          CGx0=zeros(n,1);          [Kr, error, iter, flag] = LinCG(C, CGx0, r0, CG_M, CGopts);       elseif opts.precond_one==3          % TBI       end    elseif opts.precond==2       if opts.precond_two==1,          % (L*L')Kr=r0          Kr=(L')\(L\r0);       elseif opts.precond_two==2          % [LL'+shift*B*U*(B*U)]Kr=r0          Cinvr=(L')\(L\r0);          Kr=Cinvr-CiViT*(CiV'*r0);       end    end    Q=Kr-x0*(Bx'*Kr);  % now Q'*B*x0=0;   BQ=B*Q; pBp=real(Q'*(BQ)); nrmp=sqrt(pBp);  Q=(1/nrmp)*Q;%   disp('Here 1');%   disp(norm([x0 Q]'*B*[x0 Q]-eye(2),1));   AQ=A*Q+shift*(BU*(BU'*Q)); tmp=x0'*AQ;   Rho=[rho     tmp;       conj(tmp) real(Q'*AQ)];   [v,rho]=mineig(Rho,1);   x1=[x0,Q]*v; % x1'*B*x1 approx 1      y=Q;     Bx=B*x1; xBx=real(x1'*Bx); nrmx=sqrt(xBx);   x1=(1/nrmx)*x1;  Bx=(1/nrmx)*Bx;   Ax=[Ax,AQ]*(v/nrmx); % A*x1+shift*(BU*(BU'*x1));   %Ax=A*x1+shift*(BU*(BU'*x1));   r=Ax-rho*Bx;    while err > tol & itn < maxitn,         if opts.precond==1         if opts.precond_one==1            % [C+shift*BU*(BU)']*Kr=r            Cinvr=U\(L\r(pp));            Kr=Cinvr-CiViT*(CiV'*r);         elseif opts.precond_one==2            CGx0=zeros(n,1);            [Kr, error, iter, flag] = LinCG(C, CGx0, r, CG_M, CGopts);         elseif opts.precond_one==3            % TBI         end      elseif opts.precond==2         if opts.precond_two==1,            % (L*L')Kr=r            Kr=(L')\(L\r);         elseif opts.precond_two==2            % [LL'+shift*B*U*(B*U)]Kr=r            Cinvr=(L')\(L\r);            Kr=Cinvr-CiViT*(CiV'*r);         end      end          Q=MGSg(B,[y, Kr], x1);%      disp('Here 2')%      disp(norm([x1 Q]'*B*[x1 Q]-eye(3),1));      AQ=A*Q+shift*(BU*(BU'*Q)); tmp=x1'*AQ;      Rho=[rho     tmp;           tmp'   Q'*AQ];      [v,rho]=mineig(Rho,1);      x0=x1;      x1=[x0,Q]*v; % x1'*B*x1 approx 1         y=Q*v(2:3);       Bx=B*x1; xBx=real(x1'*Bx); nrmx=sqrt(xBx);      x1=(1/nrmx)*x1;  Bx=(1/nrmx)*Bx;      Ax=[Ax,AQ]*(v/nrmx); % A*x1+shift*(BU*(BU'*x1));      %Ax=A*x1+shift*(BU*(BU'*x1));      r=Ax-rho*Bx;       nrmr=norm(r);      err=nrmr/(nrmAB(1)+abs(rho)*nrmAB(2));      hsty=[hsty; rho err];      itn = itn+1;   endendinfo = itn; y=x1; lam=rho; res=err;function [lam, y, res, hsty, info]=PSDgS(A, B, x, BU, nrmAB, shift, tol, opts)%%  Preconditioned Steepest Decent Method to compute the kth smallest eigenvalue of A-lamda B,%  where A and B are Hermitian, and B is positive definite, k=k0+1 and k0 is the number %  of columns in BU=B*U.%%  Deflation is done by shifting away known eigenvalues. Basically SD on%%      A+shift*B*U*(B*U)'=A+shift*BU*BU', B%%%---------------------------------------------------------------%%  Input%%        A     n-by-n Hermitian matrix%        B     n-by-n Hermitian matrix and positive definite%        x     n-vector, initial guess to the kth eigenvector%        BU    n-by-k0, B*U, where U's columns span an approximate invariant subspace %              (accurate enough, however) associated with 1st to k0th smallest %              eigenvalues. Assume U'*B*U=I.%              Possibly k0=0, i.e., BU is empty array.%        nrmAB 2-by-1, nrmAB(1): estimate of ||A||%                      nrmAB(2): estimate of ||B||%        shift real for shifting away known eigenvalues%        tol   tolerance for testing convergence%        opts  options ...%              opts.precond=1,  use 1st type preconditioner K=(A-sgm B)^{-1}%                   1) opts.precond_one=1: direct solver%                   2) opts.precond_one=2: A-sgm B is (effectively) SPD. K is implemented by CG   %                   3) opts.precond_one=3: A-sgm B is indefinite. K is implemented by MINRES %              opts.precond=2,  use 2nd type preconditioner; suppose A-sgm B approx LDL'.%                   1) opts.precond_two=1: K=(LL')^{-1}%                   2) opts.precond_two=2: K=[LL'+shift*B*U*(B*U)]^{-1}%%              opts.precond_sgm  the shift parameter to build a preconditioner%%  Output%%        lam  converged eigenvalue%        y     corresponding eigenvector%        res   residual error  ||A y - lam B*y||%        hsty  *-by-2 %              hsty(:,1) -- eigenvalue approximations%              hsty(:,2) -- normalized residuals ||A yi - lam_i B*yi||/(||A||*||y_i||+|lam_i|*||B||*||y_i||)%        info  number of iterations% %---------------------------------------------------------------%n=max(size(A));maxitn=round(0.2*n); k0=size(BU,2);  LSopts.Line=1;   % Line Search opts howLS=0;         % =0 solve 2-by-2 eigenvalue problem                 % =1 use Line Search function LNSRCHg(...)sgm=opts.precond_sgm;C=A-sgm*B; if opts.precond==1   if opts.precond_one==1      [L,U,pp]=lu(C,'vector');   elseif opts.precond_one==2      CGopts.nitn=10;      CGopts.tol=1e-2;      CGopts.met=2;      CG_M=0;   elseif opts.precond_one==3      % TBI   endelseif opts.precond==2   iLUopts.thresh = 0;   iLUopts.udiag = 1;   iLUopts.milu = 1;   [L, U] = luinc(C, iLUopts);      %  Scale diagonals to get L.   d = diag (U);   d = sqrt(abs(d));   for i = 1:n      if (d(i) < 1e-2)         d(i) = 1.0;      end   end   L = L * spdiags(d, 0, n, n);end  Bx=B*x; xBx=real(x'*Bx); nrmx=sqrt(xBx);x=(1/nrmx)*x; Bx=(1/nrmx)*Bx;if k0==0,   if opts.precond==1      if opts.precond_one==1         % no action      elseif opts.precond_one==2         CGopts.update = 0;      elseif opts.precond_one==3             % TBI      end   elseif opts.precond==2      % no action   end       Ax=A*x; rho=real(x'*Ax);    r=Ax-rho*Bx;      err=norm(r)/( nrmAB(1)+abs(rho)*nrmAB(2) ); itn=0;   hsty=[rho err];   while err > tol & itn < maxitn,        if opts.precond==1          if opts.precond_one==1             % C*p=r             p=U\(L\r(pp));          elseif opts.precond_one==2             CGx0=zeros(n,1);             [p, error, iter, flag] = LinCG(C, CGx0, r, CG_M, CGopts);          elseif opts.precond_one==3             % TBI          end       elseif opts.precond==2          % (L*L')p=r          p=(L')\(L\r);       end          if howLS==1,          Bp=B*p; pBp=real(p'*Bp); nrmr=sqrt(pBp); p=(1/nrmr)*p;          rhos=[rho p'*(A*p)];          LSopts.rhos=rhos;          [y, t, lam]=LNSRCHg(Ax,Bx,x,p,LSopts);       else          q=p-x*(Bx'*p);  % now q'*B*x=0;          Bq=B*q; pBp=real(q'*(Bq)); nrmp=sqrt(pBp);  q=(1/nrmp)*q;          Aq=A*q; tmp=x'*Aq;          Rho1=[rho     tmp;              conj(tmp) real(q'*Aq)];         [v,lam]=mineig(Rho1,1);         y=[x,q]*v; % y'*B*y approx 1        end       x=y;       Bx=B*x; xBx=real(x'*Bx); nrmx=sqrt(xBx);       x=(1/nrmx)*x; Bx=(1/nrmx)*Bx;       Ax=A*x; rho=lam; % rho=x'*Ax;        r=Ax-rho*Bx;        %disp([rho, lam, rho-lam])       err=norm(r)/( nrmAB(1)+abs(rho)*nrmAB(2) );       hsty=[hsty; rho err];       itn = itn+1;   endelse  % k0>0   if opts.precond==1      if opts.precond_one==1         % [C+shift*V*V']^{-1} = C^{-1}-shift*C^{-1}*V*[I+shift*V'*C^{-1}*V]^{-1}*V'*C^{-1}, where V=BU         CiV=U\(L\BU(pp,:));         T=eye(k0)+shift*BU'*CiV; invT=inv(T);         CiViT=CiV*(shift*invT);      elseif opts.precond_one==2         CGopts.update = 1;         CGopts.shift =shift;         CGopts.V = BU;      elseif opts.precond_one==3             % TBI      end   elseif opts.precond==2      if opts.precond_two==2         % [C+shift*V*V']^{-1} = C^{-1}-shift*C^{-1}*V*[I+shift*V'*C^{-1}*V]^{-1}*V'*C^{-1}, where V=BU, C approx LL'         CiV=(L')\(L\BU);         T=eye(k0)+shift*BU'*CiV; invT=inv(T);         CiViT=CiV*(shift*invT);      end   end       Ax=A*x+shift*(BU*(BU'*x)); rho=real(x'*Ax);    r=Ax-rho*Bx;      err=norm(r)/( nrmAB(1)+abs(rho)*nrmAB(2) ); itn=0;   hsty=[rho err];   while err > tol & itn < maxitn,        if opts.precond==1          if opts.precond_one==1             % [C+shift*BU*(BU)']*p=r             Cinvr=U\(L\r(pp));             p=Cinvr-CiViT*(CiV'*r);          elseif opts.precond_one==2             CGx0=zeros(n,1);             [p, error, iter, flag] = LinCG(C, CGx0, r, CG_M, CGopts);          elseif opts.precond_one==3             % TBI          end       elseif opts.precond==2          if opts.precond_two==1,             % (L*L')p=r             p=(L')\(L\r);          elseif opts.precond_two==2             % [LL'+shift*B*U*(B*U)]p=r             Cinvr=(L')\(L\r);             p=Cinvr-CiViT*(CiV'*r);          end       end           if howLS==1,          % This option has difficulty reducing residual below abount 1e-12 on an example,          % But the other option doesn't. So I think the function LNSRCHg(...) is likely          % the reason.            Bp=B*p; pBp=real(p'*Bp); nrmr=sqrt(pBp); p=(1/nrmr)*p;          Ap=A*p+shift*(BU*(BU'*p));          rhos=[rho p'*(Ap)];          LSopts.rhos=rhos;          [y, t, lam]=LNSRCHg(Ax,Bx,x,p,LSopts);       else          q=p-x*(Bx'*p);  % now q'*B*x=0;          Bq=B*q; pBp=real(q'*(Bq)); nrmp=sqrt(pBp);  q=(1/nrmp)*q;          Aq=A*q+shift*(BU*(BU'*q)); tmp=x'*Aq;          Rho1=[rho     tmp;              conj(tmp) real(q'*Aq)];         [v,lam]=mineig(Rho1,1);         y=[x,q]*v; % y'*B*y approx 1        end%       Bp=B*p; pBp=real(p'*Bp); nrmr=sqrt(pBp); p=(1/nrmr)*p;%       Ap=A*p+shift*(BU*(BU'*p));%       rhos=[rho p'*(Ap)];%       LSopts.rhos=rhos;%       %       [y, t, lam]=LNSRCHg(Ax,Bx,x,p,LSopts);       x=y;               Bx=B*x; xBx=real(x'*Bx); nrmx=sqrt(xBx);       x=(1/nrmx)*x; Bx=(1/nrmx)*Bx;       Ax=A*x+shift*(BU*(BU'*x)); rho=lam; %rho=x'*Ax;        r=Ax-rho*Bx;        err=norm(r)/( nrmAB(1)+abs(rho)*nrmAB(2) );       hsty=[hsty; rho err];       itn = itn+1;   endendinfo = itn; y=x; lam=rho; res=err;