my_ipgmres.asv

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

function [x,iter,flg,inner,operations]=my_ipgmres(A,b,tol,restart,maxit,M,x0,epsilon,scale,p_type)
% IPGMRES 
% [x,iter,flg,inner]=my_ipgmres(A,b,tol,restart,maxit,M,x0,epsilon,scale,p_type)
% this is the implementation of preconditioned GMRES method
% restart is the number of restart, 
% epsilon and scale controle the inexact precondition 
% p_type is the kind of the solver for precondition 1-CG; 2-MINRES; 3-GMRES
% flg=0 success; flg=1 fail;

% Developed by: Plum_Liliang UESTC China
% Date        : 2006-07-01

% 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) | isempty(tol)
    tol=1e-6;
end
if nargin<4 | isempty(restart)
    restart=min(10); % the maximum outer iteration is assigned to 10 defaultly
end
if nargin<5 | isempty(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(M)
    existM=1;
    if ~isequal(size(M),[m,n])
        error('The preconditioners M should match the size of A!')
    end
else
    existM=0;
    error('The preconditioners M should be assigned for this method!')
end
% check the initial guess
if nargin<7 | isempty(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<8 | isempty(epsilon)
    epsilon=0.9;
end
if nargin<9 | isempty(scale)
    scale=0.1;
end
if nargin<10 | isempty(p_type)
    p_type=3;
end
if nargin>10
    error('Too many input arguments!')
end
V=zeros(n,restart+1);
Z=V;
E1=zeros(n+1,1);E1(1,1)=1;
H=zeros(restart+1,restart+1);
flg=1;
operations=0;
op_p=0; op_ind=1;
for j=1:maxit
    epsilon=epsilon*scale;
    % set up of the inner iteration 
    r=b-A*x0;
    norm_r=norm(r,2);
    s_r=norm_r*E1;
    V(:,1)=r/norm_r;
    for i=1:restart
        % preconditioning
        if existM
            switch p_type % the initial guess is the right-hand vector
                case 1 % M is spd and CG method is applied.[x,iter,flg,inner,operations]=my_pgm(A,b,tol,restart,maxit,M1,M2,p_type,x0)
                    [Z(:,i),iter,flg,op_p] = my_pcg(M,V(:,i),epsilon,maxit*5,[],[],[],V(:,i));
                    if flg~=0
                        warning('Precondition is not convergent!');
                    end
                case 2
                    [Z(:,i),flg,relres,iter_p,resvec,resveccg] = minres(M,V(:,i),epsilon,maxit*5,[],[],V(:,i));
                    if flg~=0
                        warning('Precondition is not convergent!');                
                    end
                case 3
                    [Z(:,i),iter,flg,inner,op_p]=my_pgm(M,V(:,i),epsilon,20,10,[],[],[],V(:,i));
                    if flg~=0
                        warning('Precondition is not convergent!');
                    end
            end % switch
        end  % if exitsM
        w=A*Z(:,i);
        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);
        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;
         operations=op_p+operations+(2*i+2)*n+nz;
         %%%%%%%%%%%%%%%%%%%---record--%%%%%%%%%%
         
         % 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=x0+Z(:,1:restart)*y;
     elseif flg==0
         y=H(1:inner,1:inner)\s_r(1:inner,1);
         x0=x0+Z(:,1:inner)*y;
         break;
     end
 end % j=1:maxit
 x=x0;
 if flg==1
     inner=restart;
     iter=maxit;
 end