| [isdgmrresult,n,b,x,xl,nelt,ia,ja,a,isym,msolve,nmsl,itol,tol,itmax,iter,err,iunit,r,z,dz,rwork,iwork,rnrm,bnrm,sb,sx,jscal,kmp,lgmr,maxl,maxlp1,v,q,snormw,prod,r0nrm,hes,jpre]=isdgmr(n,b,x,xl,nelt,ia,ja,a,isym,msolve,nmsl,itol,tol,itmax,iter,err,iunit,r, |
function [isdgmrresult,n,b,x,xl,nelt,ia,ja,a,isym,msolve,nmsl,itol,tol,itmax,iter,err,iunit,r,z,dz,rwork,iwork,rnrm,bnrm,sb,sx,jscal,kmp,lgmr,maxl,maxlp1,v,q,snormw,prod,r0nrm,hes,jpre]=isdgmr(n,b,x,xl,nelt,ia,ja,a,isym,msolve,nmsl,itol,tol,itmax,iter,err,iunit,r,z,dz,rwork,iwork,rnrm,bnrm,sb,sx,jscal,kmp,lgmr,maxl,maxlp1,v,q,snormw,prod,r0nrm,hes,jpre);
isdgmrresult=[];
persistent dxnrm fuzz i ielmax rat ratmax solnrm tem ;
;
%***BEGIN PROLOGUE ISDGMR
%***SUBSIDIARY
%***PURPOSE Generalized Minimum Residual Stop Test.
% This routine calculates the stop test for the Generalized
% Minimum RESidual (GMRES) iteration scheme. It returns a
% non-zero if the error estimate (the type of which is
% determined by ITOL) is less than the user specified
% tolerance TOL.
%***LIBRARY SLATEC (SLAP)
%***CATEGORY D2A4, D2B4
%***TYPE doubleprecision (ISSGMR-S, ISDGMR-D)
%***KEYWORDS GMRES, LINEAR SYSTEM, SLAP, SPARSE, STOP TEST
%***AUTHOR Brown, Peter, (LLNL), pnbrown@llnl.gov
% Hindmarsh, Alan, (LLNL), alanh@llnl.gov
% Seager, Mark K., (LLNL), seager@llnl.gov
% Lawrence Livermore National Laboratory
% PO Box 808, L-60
% Livermore, CA 94550 (510) 423-3141
%***DESCRIPTION
%
% *Usage:
% INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, NMSL, ITOL
% INTEGER ITMAX, ITER, IUNIT, IWORK(USER DEFINED), JSCAL
% INTEGER KMP, LGMR, MAXL, MAXLP1, JPRE
% doubleprecision B(N), X(N), XL(MAXL), A(NELT), TOL, ERR,
% $ R(N), Z(N), DZ(N), RWORK(USER DEFINED),
% $ RNRM, BNRM, SB(N), SX(N), V(N,MAXLP1),
% $ Q(2*MAXL), SNORMW, PROD, R0NRM,
% $ HES(MAXLP1,MAXL)
% EXTERNAL MSOLVE
%
% IF (ISDGMR(N, B, X, XL, NELT, IA, JA, A, ISYM, MSOLVE,
% $ NMSL, ITOL, TOL, ITMAX, ITER, ERR, IUNIT, R, Z, DZ,
% $ RWORK, IWORK, RNRM, BNRM, SB, SX, JSCAL,
% $ KMP, LGMR, MAXL, MAXLP1, V, Q, SNORMW, PROD, R0NRM,
% $ HES, JPRE) .NE. 0) THEN ITERATION DONE
%
% *Arguments:
% N :IN Integer.
% Order of the Matrix.
% B :IN doubleprecision B(N).
% Right-hand-side vector.
% X :IN doubleprecision X(N).
% Approximate solution vector as of the last restart.
% XL :OUT doubleprecision XL(N)
% An array of length N used to hold the approximate
% solution as of the current iteration. Only computed by
% this routine when ITOL=11.
% NELT :IN Integer.
% Number of Non-Zeros stored in A.
% IA :IN Integer IA(NELT).
% JA :IN Integer JA(NELT).
% A :IN doubleprecision A(NELT).
% These arrays contain the matrix data structure for A.
% It could take any form. See 'Description', in the DGMRES,
% DSLUGM and DSDGMR routines for more details.
% ISYM :IN Integer.
% Flag to indicate symmetric storage format.
% If ISYM=0, all non-zero entries of the matrix are stored.
% If ISYM=1, the matrix is symmetric, and only the upper
% or lower triangle of the matrix is stored.
% MSOLVE :EXT External.
% Name of a routine which solves a linear system Mz = r for z
% given r with the preconditioning matrix M (M is supplied via
% RWORK and IWORK arrays. The name of the MSOLVE routine must
% be declared external in the calling program. The calling
% sequence to MSOLVE is:
% CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK)
% Where N is the number of unknowns, R is the right-hand side
% vector and Z is the solution upon return. NELT, IA, JA, A and
% ISYM are defined as above. RWORK is a doubleprecision array
% that can be used to pass necessary preconditioning information
% and/or workspace to MSOLVE. IWORK is an integer work array
% for the same purpose as RWORK.
% NMSL :INOUT Integer.
% A counter for the number of calls to MSOLVE.
% ITOL :IN Integer.
% Flag to indicate the type of convergence criterion used.
% ITOL=0 Means the iteration stops when the test described
% below on the residual RL is satisfied. This is
% the 'Natural Stopping Criteria' for this routine.
% Other values of ITOL cause extra, otherwise
% unnecessary, computation per iteration and are
% therefore much less efficient.
% ITOL=1 Means the iteration stops when the first test
% described below on the residual RL is satisfied,
% and there is either right or no preconditioning
% being used.
% ITOL=2 Implies that the user is using left
% preconditioning, and the second stopping criterion
% below is used.
% ITOL=3 Means the iteration stops when the third test
% described below on Minv*Residual is satisfied, and
% there is either left or no preconditioning begin
% used.
% ITOL=11 is often useful for checking and comparing
% different routines. For this case, the user must
% supply the 'exact' solution or a very accurate
% approximation (one with an error much less than
% TOL) through a common block,
% COMMON /DSLBLK/ SOLN( )
% If ITOL=11, iteration stops when the 2-norm of the
% difference between the iterative approximation and
% the user-supplied solution divided by the 2-norm
% of the user-supplied solution is less than TOL.
% Note that this requires the user to set up the
% 'COMMON /DSLBLK/ SOLN(LENGTH)' in the calling
% routine. The routine with this declaration should
% be loaded before the stop test so that the correct
% length is used by the loader. This procedure is
% not standard Fortran and may not work correctly on
% your system (although it has worked on every
% system the authors have tried). If ITOL is not 11
% then this common block is indeed standard Fortran.
% TOL :IN doubleprecision.
% Convergence criterion, as described above.
% ITMAX :IN Integer.
% Maximum number of iterations.
% ITER :IN Integer.
% The iteration for which to check for convergence.
% ERR :OUT doubleprecision.
% Error estimate of error in final approximate solution, as
% defined by ITOL. Letting norm() denote the Euclidean
% norm, ERR is defined as follows..
%
% If ITOL=0, then ERR = norm(SB*(B-A*X(L)))/norm(SB*B),
% for right or no preconditioning, and
% ERR = norm(SB*(M-inverse)*(B-A*X(L)))/
% norm(SB*(M-inverse)*B),
% for left preconditioning.
% If ITOL=1, then ERR = norm(SB*(B-A*X(L)))/norm(SB*B),
% since right or no preconditioning
% being used.
% If ITOL=2, then ERR = norm(SB*(M-inverse)*(B-A*X(L)))/
% norm(SB*(M-inverse)*B),
% since left preconditioning is being
% used.
% If ITOL=3, then ERR = Max |(Minv*(B-A*X(L)),i)/x(i)|
% i=1,n
% If ITOL=11, then ERR = norm(SB*(X(L)-SOLN))/norm(SB*SOLN).
% IUNIT :IN Integer.
% Unit number on which to write the error at each iteration,
% if this is desired for monitoring convergence. If unit
% number is 0, no writing will occur.
% R :INOUT doubleprecision R(N).
% Work array used in calling routine. It contains
% information necessary to compute the residual RL = B-A*XL.
% Z :WORK doubleprecision Z(N).
% Workspace used to hold the pseudo-residual M z = r.
% DZ :WORK doubleprecision DZ(N).
% Workspace used to hold temporary vector(s).
% RWORK :WORK doubleprecision RWORK(USER DEFINED).
% doubleprecision array that can be used by MSOLVE.
% IWORK :WORK Integer IWORK(USER DEFINED).
% Integer array that can be used by MSOLVE.
% RNRM :IN doubleprecision.
% Norm of the current residual. Type of norm depends on ITOL.
% BNRM :IN doubleprecision.
% Norm of the right hand side. Type of norm depends on ITOL.
% SB :IN doubleprecision SB(N).
% Scaling vector for B.
% SX :IN doubleprecision SX(N).
% Scaling vector for X.
% JSCAL :IN Integer.
% Flag indicating if scaling arrays SB and SX are being
% used in the calling routine DPIGMR.
% JSCAL=0 means SB and SX are not used and the
% algorithm will perform as if all
% SB(i) = 1 and SX(i) = 1.
% JSCAL=1 means only SX is used, and the algorithm
% performs as if all SB(i) = 1.
% JSCAL=2 means only SB is used, and the algorithm
% performs as if all SX(i) = 1.
% JSCAL=3 means both SB and SX are used.
% KMP :IN Integer
% The number of previous vectors the new vector VNEW
% must be made orthogonal to. (KMP <= MAXL)
% LGMR :IN Integer
% The number of GMRES iterations performed on the current call
% to DPIGMR (i.e., # iterations since the last restart) and
% the current order of the upper Hessenberg
% matrix HES.
% MAXL :IN Integer
% The maximum allowable order of the matrix H.
% MAXLP1 :IN Integer
% MAXPL1 = MAXL + 1, used for dynamic dimensioning of HES.
% V :IN doubleprecision V(N,MAXLP1)
% The N by (LGMR+1) array containing the LGMR
% orthogonal vectors V(*,1) to V(*,LGMR).
% Q :IN doubleprecision Q(2*MAXL)
% A doubleprecision array of length 2*MAXL containing the
% components of the Givens rotations used in the QR
% decomposition of HES.
% SNORMW :IN doubleprecision
% A scalar containing the scaled norm of VNEW before it
% is renormalized in DPIGMR.
% PROD :IN doubleprecision
% The product s1*s2*...*sl = the product of the sines of the
% Givens rotations used in the QR factorization of the
% Hessenberg matrix HES.
% R0NRM :IN doubleprecision
% The scaled norm of initial residual R0.
% HES :IN doubleprecision HES(MAXLP1,MAXL)
% The upper triangular factor of the QR decomposition
% of the (LGMR+1) by LGMR upper Hessenberg matrix whose
% entries are the scaled inner-products of A*V(*,I)
% and V(*,K).
% JPRE :IN Integer
% Preconditioner type flag.
% (See description of IGWK(4) in DGMRES.)
%
% *Description
% When using the GMRES solver, the preferred value for ITOL
% is 0. This is due to the fact that when ITOL=0 the norm of
% the residual required in the stopping test is obtained for
% free, since this value is already calculated in the GMRES
% algorithm. The variable RNRM contains the appropriate
% norm, which is equal to norm(SB*(RL - A*XL)) when right or
% no preconditioning is being performed, and equal to
% norm(SB*Minv*(RL - A*XL)) when using left preconditioning.
% Here, norm() is the Euclidean norm. Nonzero values of ITOL
% require additional work to calculate the actual scaled
% residual or its scaled/preconditioned form, and/or the
% approximate solution XL. Hence, these values of ITOL will
% not be as efficient as ITOL=0.
%
% *Cautions:
% This routine will attempt to write to the Fortran logical output
% unit IUNIT, if IUNIT ~= 0. Thus, the user must make sure that
% this logical unit is attached to a file or terminal before calling
% this routine with a non-zero value for IUNIT. This routine does
% not check for the validity of a non-zero IUNIT unit number.
%
% This routine does not verify that ITOL has a valid value.
% The calling routine should make such a test before calling
% ISDGMR, as is done in DGMRES.
%
%***SEE ALSO DGMRES
%***ROUTINES CALLED D1MACH, DCOPY, DNRM2, DRLCAL, DSCAL, DXLCAL
%***COMMON BLOCKS DSLBLK
%***REVISION HISTORY (YYMMDD)
% 890404 DATE WRITTEN
% 890404 Previous REVISION DATE
% 890915 Made changes requested at July 1989 CML Meeting. (MKS)
% 890922 Numerous changes to prologue to make closer to SLATEC
% standard. (FNF)
% 890929 Numerous changes to reduce SP/DP differences. (FNF)
% 910411 Prologue converted to Version 4.0 format. (BAB)
% 910502 Corrected conversion errors, etc. (FNF)
% 910502 Removed MSOLVE from ROUTINES CALLED list. (FNF)
% 910506 Made subsidiary to DGMRES. (FNF)
% 920407 COMMON BLOCK renamed DSLBLK. (WRB)
% 920511 Added complete declaration section. (WRB)
% 921026 Corrected D to E in output format. (FNF)
% 921113 Corrected C***CATEGORY line. (FNF)
%***end PROLOGUE ISDGMR
% .. Scalar Arguments ..
% .. Array Arguments ..
a_shape=size(a);a=reshape(a,1,[]);
b_shape=size(b);b=reshape(b,1,[]);
dz_shape=size(dz);dz=reshape(dz,1,[]);
hes_orig=hes;hes_shape=[maxlp1,maxl];hes=reshape([hes_orig(1:min(prod(hes_shape),numel(hes_orig))),zeros(1,max(0,prod(hes_shape)-numel(hes_orig)))],hes_shape);
q_shape=size(q);q=reshape(q,1,[]);
r_shape=size(r);r=reshape(r,1,[]);
rwork_shape=size(rwork);rwork=reshape(rwork,1,[]);
sb_shape=size(sb);sb=reshape(sb,1,[]);
sx_shape=size(sx);sx=reshape(sx,1,[]);
v_shape=size(v);v=reshape([v(:).',zeros(1,ceil(numel(v)./prod([n])).*prod([n])-numel(v))],n,[]);
x_shape=size(x);x=reshape(x,1,[]);
xl_shape=size(xl);xl=reshape(xl,1,[]);
z_shape=size(z);z=reshape(z,1,[]);
ia_shape=size(ia);ia=reshape(ia,1,[]);
iwork_shape=size(iwork);iwork=reshape(iwork,1,[]);
ja_shape=size(ja);ja=reshape(ja,1,[]);
% .. subroutine Arguments ..
% .. Arrays in Common ..
global dslblk_1; if isempty(dslblk_1), dslblk_1=zeros(1,1); end;
% .. Local Scalars ..
if isempty(dxnrm), dxnrm=0; end;
if isempty(fuzz), fuzz=0; end;
if isempty(rat), rat=0; end;
if isempty(ratmax), ratmax=0; end;
if isempty(solnrm), solnrm=0; end;
if isempty(tem), tem=0; end;
if isempty(i), i=0; end;
if isempty(ielmax), ielmax=0; end;
% .. External Functions ..
% .. External Subroutines ..
% .. Intrinsic Functions ..
% .. Common blocks ..
% common :: ;
%% common /dslblk/ soln;
%% common /dslblk/ dslblk_1;
% .. Save statement ..
%***FIRST EXECUTABLE STATEMENT ISDGMR
isdgmrresult = 0;
%
% use input from DPIGMR to determine if stop conditions are met.
%
if( itol==0 )
err = rnrm./bnrm;
end;
if((itol>0) &&(itol<=3) )
%
% use DRLCAL to calculate the scaled residual vector.
% Store answer in R.
%
if( lgmr~=0 )
[n,kmp,lgmr,maxl,v,q,r,snormw,prod,r0nrm]=drlcal(n,kmp,lgmr,maxl,v,q,r,snormw,prod,r0nrm);
end;
if( itol<=2 )
% err = ||Residual||/||RightHandSide||(2-Norms).
err = dnrm2(n,r,1)./bnrm;
%
% Unscale R by R0NRM*PROD when KMP < MAXL.
%
if((kmp<maxl) &&(lgmr~=0) )
tem = 1.0d0./(r0nrm.*prod);
[n,tem,r]=dscal(n,tem,r,1);
end;
elseif( itol==3 ) ;
% err = Max |(Minv*Residual,i)/x(i)|
% When JPRE < 0, R already contains Minv*Residual.
if( jpre>0 )
[n,r,dz,nelt,ia,ja,a,isym,rwork,iwork]=msolve(n,r,dz,nelt,ia,ja,a,isym,rwork,iwork);
nmsl = fix(nmsl + 1);
end;
%
% Unscale R by R0NRM*PROD when KMP < MAXL.
%
if((kmp<maxl) &&(lgmr~=0) )
tem = 1.0d0./(r0nrm.*prod);
[n,tem,r]=dscal(n,tem,r,1);
end;
%
[fuzz ]=d1mach(1);
ielmax = 1;
ratmax = abs(dz(1))./max(abs(x(1)),fuzz);
for i = 2 : n;
rat = abs(dz(i))./max(abs(x(i)),fuzz);
if( rat>ratmax )
ielmax = fix(i);
ratmax = rat;
end;
end; i = fix(n+1);
err = ratmax;
if( ratmax<=tol )
isdgmrresult = 1;
end;
if( iunit>0 )
writef(iunit,[repmat(' ',1,1),' ITER = ','%5i',' IELMAX = ','%5i',' |R(IELMAX)/X(IELMAX)| = ','%12.5f' ' \n'], iter , ielmax , ratmax);
end;
%format (1X,' ITER = ',i5,' IELMAX = ',i5,' |R(IELMAX)/X(IELMAX)| = ',d12.5);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
dz_shape=zeros(dz_shape);dz_shape(:)=dz(1:numel(dz_shape));dz=dz_shape;
hes_orig(1:prod(hes_shape))=hes;hes=hes_orig;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
rwork_shape=zeros(rwork_shape);rwork_shape(:)=rwork(1:numel(rwork_shape));rwork=rwork_shape;
sb_shape=zeros(sb_shape);sb_shape(:)=sb(1:numel(sb_shape));sb=sb_shape;
sx_shape=zeros(sx_shape);sx_shape(:)=sx(1:numel(sx_shape));sx=sx_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
xl_shape=zeros(xl_shape);xl_shape(:)=xl(1:numel(xl_shape));xl=xl_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
ia_shape=zeros(ia_shape);ia_shape(:)=ia(1:numel(ia_shape));ia=ia_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
ja_shape=zeros(ja_shape);ja_shape(:)=ja(1:numel(ja_shape));ja=ja_shape;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(19)), assignin('caller','FUntemp',z); evalin('caller',[inputname(19),'=FUntemp;']); end
if csnil&&~isempty(inputname(4)), assignin('caller','FUntemp',xl); evalin('caller',[inputname(4),'=FUntemp;']); end
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',x); evalin('caller',[inputname(3),'=FUntemp;']); end
if csnil&&~isempty(inputname(32)), assignin('caller','FUntemp',v); evalin('caller',[inputname(32),'=FUntemp;']); end
if csnil&&~isempty(inputname(13)), assignin('caller','FUntemp',tol); evalin('caller',[inputname(13),'=FUntemp;']); end
if csnil&&~isempty(inputname(26)), assignin('caller','FUntemp',sx); evalin('caller',[inputname(26),'=FUntemp;']); end
if csnil&&~isempty(inputname(34)), assignin('caller','FUntemp',snormw); evalin('caller',[inputname(34),'=FUntemp;']); end
if csnil&&~isempty(inputname(25)), assignin('caller','FUntemp',sb); evalin('caller',[inputname(25),'=FUntemp;']); end
if csnil&&~isempty(inputname(21)), assignin('caller','FUntemp',rwork); evalin('caller',[inputname(21),'=FUntemp;']); end
if csnil&&~isempty(inputname(23)), assignin('caller','FUntemp',rnrm); evalin('caller',[inputname(23),'=FUntemp;']); end
if csnil&&~isempty(inputname(36)), assignin('caller','FUntemp',r0nrm); evalin('caller',[inputname(36),'=FUntemp;']); end
if csnil&&~isempty(inputname(18)), assignin('caller','FUntemp',r); evalin('caller',[inputname(18),'=FUntemp;']); end
if csnil&&~isempty(inputname(33)), assignin('caller','FUntemp',q); evalin('caller',[inputname(33),'=FUntemp;']); end
if csnil&&~isempty(inputname(35)), assignin('caller','FUntemp',prod); evalin('caller',[inputname(35),'=FUntemp;']); end
if csnil&&~isempty(inputname(11)), assignin('caller','FUntemp',nmsl); evalin('caller',[inputname(11),'=FUntemp;']); end
if csnil&&~isempty(inputname(5)), assignin('caller','FUntemp',nelt); evalin('caller',[inputname(5),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',n); evalin('caller',[inputname(1),'=FUntemp;']); end
if csnil&&~isempty(inputname(10)), assignin('caller','FUntemp',msolve); evalin('caller',[inputname(10),'=FUntemp;']); end
if csnil&&~isempty(inputname(31)), assignin('caller','FUntemp',maxlp1); evalin('caller',[inputname(31),'=FUntemp;']); end
if csnil&&~isempty(inputname(30)), assignin('caller','FUntemp',maxl); evalin('caller',[inputname(30),'=FUntemp;']); end
if csnil&&~isempty(inputname(29)), assignin('caller','FUntemp',lgmr); evalin('caller',[inputname(29),'=FUntemp;']); end
if csnil&&~isempty(inputname(28)), assignin('caller','FUntemp',kmp); evalin('caller',[inputname(28),'=FUntemp;']); end
if csnil&&~isempty(inputname(27)), assignin('caller','FUntemp',jscal); evalin('caller',[inputname(27),'=FUntemp;']); end
if csnil&&~isempty(inputname(38)), assignin('caller','FUntemp',jpre); evalin('caller',[inputname(38),'=FUntemp;']); end
if csnil&&~isempty(inputname(7)), assignin('caller','FUntemp',ja); evalin('caller',[inputname(7),'=FUntemp;']); end
if csnil&&~isempty(inputname(22)), assignin('caller','FUntemp',iwork); evalin('caller',[inputname(22),'=FUntemp;']); end
if csnil&&~isempty(inputname(17)), assignin('caller','FUntemp',iunit); evalin('caller',[inputname(17),'=FUntemp;']); end
if csnil&&~isempty(inputname(12)), assignin('caller','FUntemp',itol); evalin('caller',[inputname(12),'=FUntemp;']); end
if csnil&&~isempty(inputname(14)), assignin('caller','FUntemp',itmax); evalin('caller',[inputname(14),'=FUntemp;']); end
if csnil&&~isempty(inputname(15)), assignin('caller','FUntemp',iter); evalin('caller',[inputname(15),'=FUntemp;']); end
if csnil&&~isempty(inputname(9)), assignin('caller','FUntemp',isym); evalin('caller',[inputname(9),'=FUntemp;']); end
if csnil&&~isempty(inputname(6)), assignin('caller','FUntemp',ia); evalin('caller',[inputname(6),'=FUntemp;']); end
if csnil&&~isempty(inputname(37)), assignin('caller','FUntemp',hes); evalin('caller',[inputname(37),'=FUntemp;']); end
if csnil&&~isempty(inputname(16)), assignin('caller','FUntemp',err); evalin('caller',[inputname(16),'=FUntemp;']); end
if csnil&&~isempty(inputname(20)), assignin('caller','FUntemp',dz); evalin('caller',[inputname(20),'=FUntemp;']); end
if csnil&&~isempty(inputname(24)), assignin('caller','FUntemp',bnrm); evalin('caller',[inputname(24),'=FUntemp;']); end
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',b); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(8)), assignin('caller','FUntemp',a); evalin('caller',[inputname(8),'=FUntemp;']); end
return;
end;
end;
if( itol==11 )
%
% use DXLCAL to calculate the approximate solution XL.
%
if((lgmr~=0) &&(iter>0) )
xl_orig=xl; [n,lgmr,x,xl,dumvar5,hes,maxlp1,q,v,r0nrm,dz,sx,jscal,jpre,msolve,nmsl,rwork,iwork,nelt,ia,ja,a,isym]=dxlcal(n,lgmr,x,xl,xl,hes,maxlp1,q,v,r0nrm,dz,sx,jscal,jpre,msolve,nmsl,rwork,iwork,nelt,ia,ja,a,isym); xl(dumvar5~=xl_orig)=dumvar5(dumvar5~=xl_orig);
elseif( iter==0 ) ;
% Copy X to XL to check if initial guess is good enough.
[n,x,dumvar3,xl]=dcopy(n,x,1,xl,1);
else;
% Return since this is the first call to DPIGMR on a restart.
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
dz_shape=zeros(dz_shape);dz_shape(:)=dz(1:numel(dz_shape));dz=dz_shape;
hes_orig(1:prod(hes_shape))=hes;hes=hes_orig;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
rwork_shape=zeros(rwork_shape);rwork_shape(:)=rwork(1:numel(rwork_shape));rwork=rwork_shape;
sb_shape=zeros(sb_shape);sb_shape(:)=sb(1:numel(sb_shape));sb=sb_shape;
sx_shape=zeros(sx_shape);sx_shape(:)=sx(1:numel(sx_shape));sx=sx_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
xl_shape=zeros(xl_shape);xl_shape(:)=xl(1:numel(xl_shape));xl=xl_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
ia_shape=zeros(ia_shape);ia_shape(:)=ia(1:numel(ia_shape));ia=ia_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
ja_shape=zeros(ja_shape);ja_shape(:)=ja(1:numel(ja_shape));ja=ja_shape;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(19)), assignin('caller','FUntemp',z); evalin('caller',[inputname(19),'=FUntemp;']); end
if csnil&&~isempty(inputname(4)), assignin('caller','FUntemp',xl); evalin('caller',[inputname(4),'=FUntemp;']); end
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',x); evalin('caller',[inputname(3),'=FUntemp;']); end
if csnil&&~isempty(inputname(32)), assignin('caller','FUntemp',v); evalin('caller',[inputname(32),'=FUntemp;']); end
if csnil&&~isempty(inputname(13)), assignin('caller','FUntemp',tol); evalin('caller',[inputname(13),'=FUntemp;']); end
if csnil&&~isempty(inputname(26)), assignin('caller','FUntemp',sx); evalin('caller',[inputname(26),'=FUntemp;']); end
if csnil&&~isempty(inputname(34)), assignin('caller','FUntemp',snormw); evalin('caller',[inputname(34),'=FUntemp;']); end
if csnil&&~isempty(inputname(25)), assignin('caller','FUntemp',sb); evalin('caller',[inputname(25),'=FUntemp;']); end
if csnil&&~isempty(inputname(21)), assignin('caller','FUntemp',rwork); evalin('caller',[inputname(21),'=FUntemp;']); end
if csnil&&~isempty(inputname(23)), assignin('caller','FUntemp',rnrm); evalin('caller',[inputname(23),'=FUntemp;']); end
if csnil&&~isempty(inputname(36)), assignin('caller','FUntemp',r0nrm); evalin('caller',[inputname(36),'=FUntemp;']); end
if csnil&&~isempty(inputname(18)), assignin('caller','FUntemp',r); evalin('caller',[inputname(18),'=FUntemp;']); end
if csnil&&~isempty(inputname(33)), assignin('caller','FUntemp',q); evalin('caller',[inputname(33),'=FUntemp;']); end
if csnil&&~isempty(inputname(35)), assignin('caller','FUntemp',prod); evalin('caller',[inputname(35),'=FUntemp;']); end
if csnil&&~isempty(inputname(11)), assignin('caller','FUntemp',nmsl); evalin('caller',[inputname(11),'=FUntemp;']); end
if csnil&&~isempty(inputname(5)), assignin('caller','FUntemp',nelt); evalin('caller',[inputname(5),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',n); evalin('caller',[inputname(1),'=FUntemp;']); end
if csnil&&~isempty(inputname(10)), assignin('caller','FUntemp',msolve); evalin('caller',[inputname(10),'=FUntemp;']); end
if csnil&&~isempty(inputname(31)), assignin('caller','FUntemp',maxlp1); evalin('caller',[inputname(31),'=FUntemp;']); end
if csnil&&~isempty(inputname(30)), assignin('caller','FUntemp',maxl); evalin('caller',[inputname(30),'=FUntemp;']); end
if csnil&&~isempty(inputname(29)), assignin('caller','FUntemp',lgmr); evalin('caller',[inputname(29),'=FUntemp;']); end
if csnil&&~isempty(inputname(28)), assignin('caller','FUntemp',kmp); evalin('caller',[inputname(28),'=FUntemp;']); end
if csnil&&~isempty(inputname(27)), assignin('caller','FUntemp',jscal); evalin('caller',[inputname(27),'=FUntemp;']); end
if csnil&&~isempty(inputname(38)), assignin('caller','FUntemp',jpre); evalin('caller',[inputname(38),'=FUntemp;']); end
if csnil&&~isempty(inputname(7)), assignin('caller','FUntemp',ja); evalin('caller',[inputname(7),'=FUntemp;']); end
if csnil&&~isempty(inputname(22)), assignin('caller','FUntemp',iwork); evalin('caller',[inputname(22),'=FUntemp;']); end
if csnil&&~isempty(inputname(17)), assignin('caller','FUntemp',iunit); evalin('caller',[inputname(17),'=FUntemp;']); end
if csnil&&~isempty(inputname(12)), assignin('caller','FUntemp',itol); evalin('caller',[inputname(12),'=FUntemp;']); end
if csnil&&~isempty(inputname(14)), assignin('caller','FUntemp',itmax); evalin('caller',[inputname(14),'=FUntemp;']); end
if csnil&&~isempty(inputname(15)), assignin('caller','FUntemp',iter); evalin('caller',[inputname(15),'=FUntemp;']); end
if csnil&&~isempty(inputname(9)), assignin('caller','FUntemp',isym); evalin('caller',[inputname(9),'=FUntemp;']); end
if csnil&&~isempty(inputname(6)), assignin('caller','FUntemp',ia); evalin('caller',[inputname(6),'=FUntemp;']); end
if csnil&&~isempty(inputname(37)), assignin('caller','FUntemp',hes); evalin('caller',[inputname(37),'=FUntemp;']); end
if csnil&&~isempty(inputname(16)), assignin('caller','FUntemp',err); evalin('caller',[inputname(16),'=FUntemp;']); end
if csnil&&~isempty(inputname(20)), assignin('caller','FUntemp',dz); evalin('caller',[inputname(20),'=FUntemp;']); end
if csnil&&~isempty(inputname(24)), assignin('caller','FUntemp',bnrm); evalin('caller',[inputname(24),'=FUntemp;']); end
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',b); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(8)), assignin('caller','FUntemp',a); evalin('caller',[inputname(8),'=FUntemp;']); end
return;
end;
%
if((jscal==0) ||(jscal==2) )
% err = ||x-TrueSolution||/||TrueSolution||(2-Norms).
if( iter==0 )
[ solnrm ,n,dslblk_1]=dnrm2(n,dslblk_1,1);
end;
for i = 1 : n;
dz(i) = xl(i) - dslblk_1(i);
end; i = fix(n+1);
err = dnrm2(n,dz,1)./solnrm;
else;
if( iter==0 )
solnrm = 0;
for i = 1 : n;
solnrm = solnrm +(sx(i).*dslblk_1(i)).^2;
end; i = fix(n+1);
solnrm = sqrt(solnrm);
end;
dxnrm = 0;
for i = 1 : n;
dxnrm = dxnrm +(sx(i).*(xl(i)-dslblk_1(i))).^2;
end; i = fix(n+1);
dxnrm = sqrt(dxnrm);
% err = ||SX*(x-TrueSolution)||/||SX*TrueSolution|| (2-Norms).
err = dxnrm./solnrm;
end;
end;
%
if( iunit~=0 )
if( iter==0 )
writef(iunit,[' Generalized Minimum Residual(','%3i','%3i',') for ','N, ITOL = ','%5i','%5i', '\n ' ,' ITER',' Natural Err Est',' Error Estimate' ' \n'], n , itol , maxl , kmp);
%format (' Generalized Minimum Residual(',i3,i3,') for ','N, ITOL = ',i5,i5,/' ITER',' Natural Err Est',' Error Estimate');
end;
writef(iunit,[repmat(' ',1,1),'%4i',repmat(' ',1,1),'%16.7f',repmat(' ',1,1),'%16.7f' ' \n'], iter , rnrm./bnrm , err);
%format(1x,i4,1x,d16.7,1x,d16.7);
end;
if( err<=tol )
isdgmrresult = 1;
end;
%
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
dz_shape=zeros(dz_shape);dz_shape(:)=dz(1:numel(dz_shape));dz=dz_shape;
hes_orig(1:prod(hes_shape))=hes;hes=hes_orig;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
rwork_shape=zeros(rwork_shape);rwork_shape(:)=rwork(1:numel(rwork_shape));rwork=rwork_shape;
sb_shape=zeros(sb_shape);sb_shape(:)=sb(1:numel(sb_shape));sb=sb_shape;
sx_shape=zeros(sx_shape);sx_shape(:)=sx(1:numel(sx_shape));sx=sx_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
xl_shape=zeros(xl_shape);xl_shape(:)=xl(1:numel(xl_shape));xl=xl_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
ia_shape=zeros(ia_shape);ia_shape(:)=ia(1:numel(ia_shape));ia=ia_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
ja_shape=zeros(ja_shape);ja_shape(:)=ja(1:numel(ja_shape));ja=ja_shape;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(19)), assignin('caller','FUntemp',z); evalin('caller',[inputname(19),'=FUntemp;']); end
if csnil&&~isempty(inputname(4)), assignin('caller','FUntemp',xl); evalin('caller',[inputname(4),'=FUntemp;']); end
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',x); evalin('caller',[inputname(3),'=FUntemp;']); end
if csnil&&~isempty(inputname(32)), assignin('caller','FUntemp',v); evalin('caller',[inputname(32),'=FUntemp;']); end
if csnil&&~isempty(inputname(13)), assignin('caller','FUntemp',tol); evalin('caller',[inputname(13),'=FUntemp;']); end
if csnil&&~isempty(inputname(26)), assignin('caller','FUntemp',sx); evalin('caller',[inputname(26),'=FUntemp;']); end
if csnil&&~isempty(inputname(34)), assignin('caller','FUntemp',snormw); evalin('caller',[inputname(34),'=FUntemp;']); end
if csnil&&~isempty(inputname(25)), assignin('caller','FUntemp',sb); evalin('caller',[inputname(25),'=FUntemp;']); end
if csnil&&~isempty(inputname(21)), assignin('caller','FUntemp',rwork); evalin('caller',[inputname(21),'=FUntemp;']); end
if csnil&&~isempty(inputname(23)), assignin('caller','FUntemp',rnrm); evalin('caller',[inputname(23),'=FUntemp;']); end
if csnil&&~isempty(inputname(36)), assignin('caller','FUntemp',r0nrm); evalin('caller',[inputname(36),'=FUntemp;']); end
if csnil&&~isempty(inputname(18)), assignin('caller','FUntemp',r); evalin('caller',[inputname(18),'=FUntemp;']); end
if csnil&&~isempty(inputname(33)), assignin('caller','FUntemp',q); evalin('caller',[inputname(33),'=FUntemp;']); end
if csnil&&~isempty(inputname(35)), assignin('caller','FUntemp',prod); evalin('caller',[inputname(35),'=FUntemp;']); end
if csnil&&~isempty(inputname(11)), assignin('caller','FUntemp',nmsl); evalin('caller',[inputname(11),'=FUntemp;']); end
if csnil&&~isempty(inputname(5)), assignin('caller','FUntemp',nelt); evalin('caller',[inputname(5),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',n); evalin('caller',[inputname(1),'=FUntemp;']); end
if csnil&&~isempty(inputname(10)), assignin('caller','FUntemp',msolve); evalin('caller',[inputname(10),'=FUntemp;']); end
if csnil&&~isempty(inputname(31)), assignin('caller','FUntemp',maxlp1); evalin('caller',[inputname(31),'=FUntemp;']); end
if csnil&&~isempty(inputname(30)), assignin('caller','FUntemp',maxl); evalin('caller',[inputname(30),'=FUntemp;']); end
if csnil&&~isempty(inputname(29)), assignin('caller','FUntemp',lgmr); evalin('caller',[inputname(29),'=FUntemp;']); end
if csnil&&~isempty(inputname(28)), assignin('caller','FUntemp',kmp); evalin('caller',[inputname(28),'=FUntemp;']); end
if csnil&&~isempty(inputname(27)), assignin('caller','FUntemp',jscal); evalin('caller',[inputname(27),'=FUntemp;']); end
if csnil&&~isempty(inputname(38)), assignin('caller','FUntemp',jpre); evalin('caller',[inputname(38),'=FUntemp;']); end
if csnil&&~isempty(inputname(7)), assignin('caller','FUntemp',ja); evalin('caller',[inputname(7),'=FUntemp;']); end
if csnil&&~isempty(inputname(22)), assignin('caller','FUntemp',iwork); evalin('caller',[inputname(22),'=FUntemp;']); end
if csnil&&~isempty(inputname(17)), assignin('caller','FUntemp',iunit); evalin('caller',[inputname(17),'=FUntemp;']); end
if csnil&&~isempty(inputname(12)), assignin('caller','FUntemp',itol); evalin('caller',[inputname(12),'=FUntemp;']); end
if csnil&&~isempty(inputname(14)), assignin('caller','FUntemp',itmax); evalin('caller',[inputname(14),'=FUntemp;']); end
if csnil&&~isempty(inputname(15)), assignin('caller','FUntemp',iter); evalin('caller',[inputname(15),'=FUntemp;']); end
if csnil&&~isempty(inputname(9)), assignin('caller','FUntemp',isym); evalin('caller',[inputname(9),'=FUntemp;']); end
if csnil&&~isempty(inputname(6)), assignin('caller','FUntemp',ia); evalin('caller',[inputname(6),'=FUntemp;']); end
if csnil&&~isempty(inputname(37)), assignin('caller','FUntemp',hes); evalin('caller',[inputname(37),'=FUntemp;']); end
if csnil&&~isempty(inputname(16)), assignin('caller','FUntemp',err); evalin('caller',[inputname(16),'=FUntemp;']); end
if csnil&&~isempty(inputname(20)), assignin('caller','FUntemp',dz); evalin('caller',[inputname(20),'=FUntemp;']); end
if csnil&&~isempty(inputname(24)), assignin('caller','FUntemp',bnrm); evalin('caller',[inputname(24),'=FUntemp;']); end
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',b); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(8)), assignin('caller','FUntemp',a); evalin('caller',[inputname(8),'=FUntemp;']); end
return;
%------------- LAST LINE OF ISDGMR FOLLOWS ----------------------------
end
%DECK ISDIR
|
|
