my_pgm.asv

求解线性系统的Krylov方法的工具箱

function [x,iter,flg,inner]=my_pgm(A,b,tol,restart,maxit,M1,M2,p_type,x0)
% [x,iter,flg]=my_ipgm(A,b,tol,restart,maxit,M1,M2,p_type,x0)
% this is the implementation of preconditioned GMRES method
% m is the number of restart, 
% p is the preconditioning type, p=0: implicit preconditioner; p=1 explicit
% preconditioner
% flg=0 success; flg=1 fail;

% Developed by: Plum_Liliang UESTC China
% Date        : 2006-06-08

% checking input arguments and assign default value
[m,n]=size(A);
if m~=n
    error('The coefficient matrix must be square!')
end 
[m_b,n_b]=size(b);
if n_b~=1
    error('b must be a vector!')
end

if m~=m_b
    error('The right-hand vecotor must have the same length of A!')
end

if nargin<2
    error('Not enough input arguments!');
end
if nargin<3 | tol==[]
    tol=1e-6;
end
if nargin<4 | restart==[]
    restart=min(10); % the maximum outer iteration is assigned to 10 defaultly
end
if nargin<5 | maxit==[]
    maxit=min(20,n);
end

% if 'b' is zero then the solution is zero
norm_b=norm(b);
if norm_b==0
    x=0;
    iter=0;
    flg=0;
    inner=0;
    x=0;
    return;
end

% Check the preconditioners
if nargin>=6 & ~isempty(M1)
    existM1=1;
    if ~isequal(size(M1),[m,n])
        error('The preconditioners M1 should match the size of A!')
    end
else
    existM1=0;
end
if nargin>=7 & ~isempty(M2)
    existM2=1;
    if ~isequal(size(M2),[m,n])
        error('The preconditioners M2 should match the size of A!')
    end
else
    existM2=0;
end

if nargin<8 | p_type==[]
    p_type=0;
end

% check the initial guess
if nargin<9 | x0==[]
    x0=zeros(n,1); % the initial guess is defaultly 0-vector;
else
    [m_x,n_x]=size(x0);
    if m_x~=n | n_x~=1
        error ('The initial guess should be has the same size with the right-hand vector!');
    end
end

if nargin>9
    error('Too many input arguments!')
end
% the main process of GMRES
V=zeros(n,restart+1);
E1=zeros(n+1,1);E1(1,1)=1;
H=zeros(restart+1,restart+1);
flg=1;
for j=1:maxit
    % set up of the inner iteration 
    r=b-A*x0;
    % preconditioning
    if p_type==0
        if existM1==1
            r=M1\r;
        end
        if existM2==2
            r=M2\r;
        end
    elseif p_type==1
        if existM1==1
            r=M1*r;
        end
        if existM2==2
            r=M2*r;
        end
    end
    norm_r=norm(r,2);
    V(:,1)=r/norm_r;
    s_r=norm_r*E1;
    for i=1:restart
        w=A*V(:,i);
         % preconditioning
         if p_type==0
             if existM1==1
                 w=M1\w;
             end
             if existM2==1
                 w=M2\w;
             end
         elseif p_type==1
             if existM1==1
                 w=M1*w;
             end
             if existM2==1
                 w=M2*w;
             end
         end
         %construct the orthogonal vectors
         for k=1:i
             H(k,i)=w'*V(:,k);
             w=w-H(k,i)*V(:,k);
         end
         H(i+1,i)=norm(w,2);
         V(:,i+1)=w/H(i+1,i);
         % premultiply H with 
         for k=1:i-1
             h_temp1=c(k)*H(k,i)-s(k)*H(k+1,i);
             h_temp2=s(k)*H(k,i)+c(k)*H(k+1,i);
             H(k,i)=h_temp1; H(k+1,i)=h_temp2;
         end
         temp=sqrt(H(i,i)^2+H(i+1,i)^2);
         c(i)=H(i,i)/temp;
         s(i)=-H(i+1,i)/temp;
         H(i,i)=temp;H(i+1,i)=0;
         s_r_temp=c(i)*s_r(i)-s(i)*s_r(i+1);
         s_r(i+1)=s(i)*s_r(i)+c(i)*s_r(i+1);
         s_r(i)=s_r_temp;
         % converge
         if abs(s_r(i+1))<tol
             flg=0;
             inner=i;
             iter=j;
             break;
         end
     end % end  i=1:restart
     if flg==1
         y=H(1:restart,1:restart)\s_r(1:restart,1);
         x0=V(:,1:restart)*y;
     elseif flg==0
         y=H(1:inner,1:inner)\s_r(1:inner,1);
         x0=V(:,1:inner)*y;
         break;
     end
 end % j=1:maxit
 x=x0;
 if flg==1
     inner=restart;
     iter=maxit;
 end