最基本的特征值问题分为三类：
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.2
k0=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;
end
else % 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;
end
end
info = 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.1
Bx=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;
end
else % 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;
end
end
info = 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
end
elseif 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
end
end
V=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;
%end
tmp=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);
return
else
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
end
end
if 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
end
else % 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;
end
end
function [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');
return
end
%[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 SD
opts.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 LOPCG
opts.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.')
end
if ( ~ strcmp(head0,'%%MatrixMarket') )
error('Not a valid MatrixMarket header.')
end
if ( ~ 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 them
commentline = fgets(mmfile);
while length(commentline) > 0 & commentline(1) == '%',
commentline = fgets(mmfile);
end
% Read size information, then branch according to
% sparse or dense format
if ( 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);
end
elseif ( 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.');
end
end
%
% 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);
end
fclose(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
end
elseif 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;
end
else % 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;
end
end
info = 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
end
elseif 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;
end
else % 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;
end
end
info = itn; y=x; lam=rho; res=err;
扫一扫
2769
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
