| [issomnresult,n,b,x,nelt,ia,ja,a,isym,msolve,nsave,itol,tol,itmax,iter,err,ierr,iunit,r,z,p,ap,emap,dz,csav,rwork,iwork,ak,bnrm,solnrm]=issomn(n,b,x,nelt,ia,ja,a,isym,msolve,nsave,itol,tol,itmax,iter,err,ierr,iunit,r,z,p,ap,emap,dz,csav,rwork,iwork,ak,bnr |
function [issomnresult,n,b,x,nelt,ia,ja,a,isym,msolve,nsave,itol,tol,itmax,iter,err,ierr,iunit,r,z,p,ap,emap,dz,csav,rwork,iwork,ak,bnrm,solnrm]=issomn(n,b,x,nelt,ia,ja,a,isym,msolve,nsave,itol,tol,itmax,iter,err,ierr,iunit,r,z,p,ap,emap,dz,csav,rwork,iwork,ak,bnrm,solnrm);
issomnresult=[];
persistent i ;
;
%***BEGIN PROLOGUE ISSOMN
%***SUBSIDIARY
%***PURPOSE Preconditioned Orthomin Stop Test.
% This routine calculates the stop test for the Orthomin
% 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 SINGLE PRECISION (ISSOMN-S, ISDOMN-D)
%***KEYWORDS ITERATIVE PRECONDITION, NON-SYMMETRIC LINEAR SYSTEM,
% ORTHOMIN, SLAP, SPARSE, STOP TEST
%***AUTHOR Greenbaum, Anne, (Courant Institute)
% Seager, Mark K., (LLNL)
% Lawrence Livermore National Laboratory
% PO BOX 808, L-60
% Livermore, CA 94550 (510) 423-3141
% seager@llnl.gov
%***DESCRIPTION
%
% *Usage:
% INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, NSAVE, ITOL, ITMAX
% INTEGER ITER, IERR, IUNIT, IWORK(USER DEFINED)
% REAL B(N), X(N), A(NELT), TOL, ERR, R(N), Z(N)
% REAL P(N,0:NSAVE), AP(N,0:NSAVE), EMAP(N,0:NSAVE)
% REAL DZ(N), CSAV(NSAVE), RWORK(USER DEFINED), AK
% REAL BNRM, SOLNRM
% EXTERNAL MSOLVE
%
% IF( ISSOMN(N, B, X, NELT, IA, JA, A, ISYM, MSOLVE, NSAVE,
% $ ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, R, Z, P, AP,
% $ EMAP, DZ, CSAV, RWORK, IWORK, AK, BNRM, SOLNRM)
% $ .NE.0 ) THEN ITERATION CONVERGED
%
% *Arguments:
% N :IN Integer.
% Order of the matrix.
% B :IN Real B(N).
% Right-hand side vector.
% X :IN Real X(N).
% On input X is your initial guess for solution vector.
% On output X is the final approximate solution.
% NELT :IN Integer.
% Number of Non-Zeros stored in A.
% IA :IN Integer IA(NELT).
% JA :IN Integer JA(NELT).
% A :IN Real A(NELT).
% These arrays should hold the matrix A in either the SLAP
% Triad format or the SLAP Column format. See 'Description'
% in the SSDOMN or SSLUOM prologue.
% 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 real 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.
% NSAVE :IN Integer.
% Number of direction vectors to savemlv and orthogonalize against.
% ITOL :IN Integer.
% Flag to indicate type of convergence criterion.
% If ITOL=1, iteration stops when the 2-norm of the residual
% divided by the 2-norm of the right-hand side is less than TOL.
% If ITOL=2, iteration stops when the 2-norm of M-inv times the
% residual divided by the 2-norm of M-inv times the right hand
% side is less than TOL, where M-inv is the inverse of the
% diagonal of A.
% 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 /SSLBLK/ 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 /SSLBLK/ 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 Real.
% Convergence criterion, as described above.
% ITMAX :IN Integer.
% Maximum number of iterations.
% ITER :IN Integer.
% Current iteration count. (Must be zero on first call.)
% ERR :OUT Real.
% Error estimate of error in final approximate solution, as
% defined by ITOL.
% IERR :OUT Integer.
% Error flag. IERR is set to 3 if ITOL is not one of the
% acceptable values, see above.
% 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 :IN Real R(N).
% The residual R = B-AX.
% Z :WORK Real Z(N).
% P :IN Real P(N,0:NSAVE).
% Workspace used to hold the conjugate direction vector(s).
% AP :IN Real AP(N,0:NSAVE).
% Workspace used to hold the matrix A times the P vector(s).
% EMAP :IN Real EMAP(N,0:NSAVE).
% Workspace used to hold M-inv times the AP vector(s).
% DZ :WORK Real DZ(N).
% Workspace.
% CSAV :DUMMY Real CSAV(NSAVE)
% Reserved for future use.
% RWORK :WORK Real RWORK(USER DEFINED).
% Real array that can be used for workspace in MSOLVE.
% IWORK :WORK Integer IWORK(USER DEFINED).
% Integer array that can be used for workspace in MSOLVE.
% AK :IN Real.
% Current iterate Orthomin iteration parameter.
% BNRM :OUT Real.
% Current solution B-norm, if ITOL = 1 or 2.
% SOLNRM :OUT Real.
% truemlv solution norm, if ITOL = 11.
%
% *function Return Values:
% 0 : Error estimate (determined by ITOL) is *NOT* less than the
% specified tolerance, TOL. The iteration must continue.
% 1 : Error estimate (determined by ITOL) is less than the
% specified tolerance, TOL. The iteration can be considered
% complete.
%
% *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.
%
%***SEE ALSO SOMN, SSDOMN, SSLUOM
%***ROUTINES CALLED R1MACH, SNRM2
%***COMMON BLOCKS SSLBLK
%***REVISION HISTORY (YYMMDD)
% 871119 DATE WRITTEN
% 881213 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)
% 891003 Removed C***REFER TO line, per MKS.
% 910411 Prologue converted to Version 4.0 format. (BAB)
% 910502 Removed MSOLVE from ROUTINES CALLED list. (FNF)
% 910506 Made subsidiary to SOMN. (FNF)
% 920407 COMMON BLOCK renamed SSLBLK. (WRB)
% 920511 Added complete declaration section. (WRB)
% 920930 Corrected to not print AK when ITER=0. (FNF)
% 921026 Changed 1.0E10 to R1MACH(2). (FNF)
% 921113 Corrected C***CATEGORY line. (FNF)
%***end PROLOGUE ISSOMN
% .. Scalar Arguments ..
% .. Array Arguments ..
ap_orig=ap;ap_shape=[n,0:nsave];ap=reshape([ap_orig(1:min(prod(ap_shape),numel(ap_orig))),zeros(1,max(0,prod(ap_shape)-numel(ap_orig)))],ap_shape);
emap_orig=emap;emap_shape=[n,0:nsave];emap=reshape([emap_orig(1:min(prod(emap_shape),numel(emap_orig))),zeros(1,max(0,prod(emap_shape)-numel(emap_orig)))],emap_shape);
p_orig=p;p_shape=[n,0:nsave];p=reshape([p_orig(1:min(prod(p_shape),numel(p_orig))),zeros(1,max(0,prod(p_shape)-numel(p_orig)))],p_shape);
rwork_shape=size(rwork);rwork=reshape(rwork,1,[]);
iwork_shape=size(iwork);iwork=reshape(iwork,1,[]);
% .. subroutine Arguments ..
% .. Arrays in Common ..
global sslblk_1; if isempty(sslblk_1), sslblk_1=zeros(1,1); end;
% .. Local Scalars ..
if isempty(i), i=0; end;
% .. External Functions ..
% .. Common blocks ..
% common :: ;
%% common /sslblk/ soln;
%% common /sslblk/ sslblk_1;
%***FIRST EXECUTABLE STATEMENT ISSOMN
issomnresult = 0;
%
if( itol==1 )
% err = ||Residual||/||RightHandSide|| (2-Norms).
if( iter==0 )
[ bnrm ,n,b]=snrm2(n,b,1);
end;
err = snrm2(n,r,1)./bnrm;
elseif( itol==2 ) ;
% -1 -1
% err = ||M Residual||/||M RightHandSide|| (2-Norms).
if( iter==0 )
[n,b,dz,nelt,ia,ja,a,isym,rwork,iwork]=msolve(n,b,dz,nelt,ia,ja,a,isym,rwork,iwork);
[bnrm ,n,dz]=snrm2(n,dz,1);
end;
err = snrm2(n,z,1)./bnrm;
elseif( itol==11 ) ;
% err = ||x-TrueSolution||/||TrueSolution|| (2-Norms).
if( iter==0 )
[ solnrm ,n,sslblk_1]=snrm2(n,sslblk_1,1);
end;
for i = 1 : n;
dz(i) = x(i) - sslblk_1(i);
end; i = fix(n+1);
err = snrm2(n,dz,1)./solnrm;
else;
%
% If we get here ITOL is not one of the acceptable values.
[err ]=r1mach(2);
ierr = 3;
end;
%
if( iunit~=0 )
if( iter==0 )
writef(iunit,[' Preconditioned Orthomin(','%3i',') for ','N, ITOL = ','%5i','%5i', '\n ' ,' ITER',' Error Estimate',' Alpha' ' \n'], nsave , n , itol);
%format (' Preconditioned Orthomin(',i3,') for ','N, ITOL = ',i5,i5,/' ITER',' Error Estimate',' Alpha');
writef(iunit,[repmat(' ',1,1),'%4i',repmat(' ',1,1),'%16.7f',repmat(' ',1,1),'%16.7f' ' \n'], iter , err);
else;
writef(iunit,[repmat(' ',1,1),'%4i',repmat(' ',1,1),'%16.7f',repmat(' ',1,1),'%16.7f' ' \n'], iter , err , ak);
end;
end;
if( err<=tol )
issomnresult = 1;
end;
%
ap_orig(1:prod(ap_shape))=ap;ap=ap_orig;
emap_orig(1:prod(emap_shape))=emap;emap=emap_orig;
p_orig(1:prod(p_shape))=p;p=p_orig;
rwork_shape=zeros(rwork_shape);rwork_shape(:)=rwork(1:numel(rwork_shape));rwork=rwork_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_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(3)), assignin('caller','FUntemp',x); evalin('caller',[inputname(3),'=FUntemp;']); end
if csnil&&~isempty(inputname(12)), assignin('caller','FUntemp',tol); evalin('caller',[inputname(12),'=FUntemp;']); end
if csnil&&~isempty(inputname(29)), assignin('caller','FUntemp',solnrm); evalin('caller',[inputname(29),'=FUntemp;']); end
if csnil&&~isempty(inputname(25)), assignin('caller','FUntemp',rwork); evalin('caller',[inputname(25),'=FUntemp;']); end
if csnil&&~isempty(inputname(18)), assignin('caller','FUntemp',r); evalin('caller',[inputname(18),'=FUntemp;']); end
if csnil&&~isempty(inputname(20)), assignin('caller','FUntemp',p); evalin('caller',[inputname(20),'=FUntemp;']); end
if csnil&&~isempty(inputname(10)), assignin('caller','FUntemp',nsave); evalin('caller',[inputname(10),'=FUntemp;']); end
if csnil&&~isempty(inputname(4)), assignin('caller','FUntemp',nelt); evalin('caller',[inputname(4),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',n); evalin('caller',[inputname(1),'=FUntemp;']); end
if csnil&&~isempty(inputname(9)), assignin('caller','FUntemp',msolve); evalin('caller',[inputname(9),'=FUntemp;']); end
if csnil&&~isempty(inputname(6)), assignin('caller','FUntemp',ja); evalin('caller',[inputname(6),'=FUntemp;']); end
if csnil&&~isempty(inputname(26)), assignin('caller','FUntemp',iwork); evalin('caller',[inputname(26),'=FUntemp;']); end
if csnil&&~isempty(inputname(17)), assignin('caller','FUntemp',iunit); evalin('caller',[inputname(17),'=FUntemp;']); end
if csnil&&~isempty(inputname(11)), assignin('caller','FUntemp',itol); evalin('caller',[inputname(11),'=FUntemp;']); end
if csnil&&~isempty(inputname(13)), assignin('caller','FUntemp',itmax); evalin('caller',[inputname(13),'=FUntemp;']); end
if csnil&&~isempty(inputname(14)), assignin('caller','FUntemp',iter); evalin('caller',[inputname(14),'=FUntemp;']); end
if csnil&&~isempty(inputname(8)), assignin('caller','FUntemp',isym); evalin('caller',[inputname(8),'=FUntemp;']); end
if csnil&&~isempty(inputname(16)), assignin('caller','FUntemp',ierr); evalin('caller',[inputname(16),'=FUntemp;']); end
if csnil&&~isempty(inputname(5)), assignin('caller','FUntemp',ia); evalin('caller',[inputname(5),'=FUntemp;']); end
if csnil&&~isempty(inputname(15)), assignin('caller','FUntemp',err); evalin('caller',[inputname(15),'=FUntemp;']); end
if csnil&&~isempty(inputname(22)), assignin('caller','FUntemp',emap); evalin('caller',[inputname(22),'=FUntemp;']); end
if csnil&&~isempty(inputname(23)), assignin('caller','FUntemp',dz); evalin('caller',[inputname(23),'=FUntemp;']); end
if csnil&&~isempty(inputname(24)), assignin('caller','FUntemp',csav); evalin('caller',[inputname(24),'=FUntemp;']); end
if csnil&&~isempty(inputname(28)), assignin('caller','FUntemp',bnrm); evalin('caller',[inputname(28),'=FUntemp;']); end
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',b); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(21)), assignin('caller','FUntemp',ap); evalin('caller',[inputname(21),'=FUntemp;']); end
if csnil&&~isempty(inputname(27)), assignin('caller','FUntemp',ak); evalin('caller',[inputname(27),'=FUntemp;']); end
if csnil&&~isempty(inputname(7)), assignin('caller','FUntemp',a); evalin('caller',[inputname(7),'=FUntemp;']); end
return;
%format(1x,i4,1x,e16.7,1x,e16.7);
%------------- LAST LINE OF ISSOMN FOLLOWS ----------------------------
ap_orig(1:prod(ap_shape))=ap;ap=ap_orig;
emap_orig(1:prod(emap_shape))=emap;emap=emap_orig;
p_orig(1:prod(p_shape))=p;p=p_orig;
rwork_shape=zeros(rwork_shape);rwork_shape(:)=rwork(1:numel(rwork_shape));rwork=rwork_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_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(3)), assignin('caller','FUntemp',x); evalin('caller',[inputname(3),'=FUntemp;']); end
if csnil&&~isempty(inputname(12)), assignin('caller','FUntemp',tol); evalin('caller',[inputname(12),'=FUntemp;']); end
if csnil&&~isempty(inputname(29)), assignin('caller','FUntemp',solnrm); evalin('caller',[inputname(29),'=FUntemp;']); end
if csnil&&~isempty(inputname(25)), assignin('caller','FUntemp',rwork); evalin('caller',[inputname(25),'=FUntemp;']); end
if csnil&&~isempty(inputname(18)), assignin('caller','FUntemp',r); evalin('caller',[inputname(18),'=FUntemp;']); end
if csnil&&~isempty(inputname(20)), assignin('caller','FUntemp',p); evalin('caller',[inputname(20),'=FUntemp;']); end
if csnil&&~isempty(inputname(10)), assignin('caller','FUntemp',nsave); evalin('caller',[inputname(10),'=FUntemp;']); end
if csnil&&~isempty(inputname(4)), assignin('caller','FUntemp',nelt); evalin('caller',[inputname(4),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',n); evalin('caller',[inputname(1),'=FUntemp;']); end
if csnil&&~isempty(inputname(9)), assignin('caller','FUntemp',msolve); evalin('caller',[inputname(9),'=FUntemp;']); end
if csnil&&~isempty(inputname(6)), assignin('caller','FUntemp',ja); evalin('caller',[inputname(6),'=FUntemp;']); end
if csnil&&~isempty(inputname(26)), assignin('caller','FUntemp',iwork); evalin('caller',[inputname(26),'=FUntemp;']); end
if csnil&&~isempty(inputname(17)), assignin('caller','FUntemp',iunit); evalin('caller',[inputname(17),'=FUntemp;']); end
if csnil&&~isempty(inputname(11)), assignin('caller','FUntemp',itol); evalin('caller',[inputname(11),'=FUntemp;']); end
if csnil&&~isempty(inputname(13)), assignin('caller','FUntemp',itmax); evalin('caller',[inputname(13),'=FUntemp;']); end
if csnil&&~isempty(inputname(14)), assignin('caller','FUntemp',iter); evalin('caller',[inputname(14),'=FUntemp;']); end
if csnil&&~isempty(inputname(8)), assignin('caller','FUntemp',isym); evalin('caller',[inputname(8),'=FUntemp;']); end
if csnil&&~isempty(inputname(16)), assignin('caller','FUntemp',ierr); evalin('caller',[inputname(16),'=FUntemp;']); end
if csnil&&~isempty(inputname(5)), assignin('caller','FUntemp',ia); evalin('caller',[inputname(5),'=FUntemp;']); end
if csnil&&~isempty(inputname(15)), assignin('caller','FUntemp',err); evalin('caller',[inputname(15),'=FUntemp;']); end
if csnil&&~isempty(inputname(22)), assignin('caller','FUntemp',emap); evalin('caller',[inputname(22),'=FUntemp;']); end
if csnil&&~isempty(inputname(23)), assignin('caller','FUntemp',dz); evalin('caller',[inputname(23),'=FUntemp;']); end
if csnil&&~isempty(inputname(24)), assignin('caller','FUntemp',csav); evalin('caller',[inputname(24),'=FUntemp;']); end
if csnil&&~isempty(inputname(28)), assignin('caller','FUntemp',bnrm); evalin('caller',[inputname(28),'=FUntemp;']); end
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',b); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(21)), assignin('caller','FUntemp',ap); evalin('caller',[inputname(21),'=FUntemp;']); end
if csnil&&~isempty(inputname(27)), assignin('caller','FUntemp',ak); evalin('caller',[inputname(27),'=FUntemp;']); end
if csnil&&~isempty(inputname(7)), assignin('caller','FUntemp',a); evalin('caller',[inputname(7),'=FUntemp;']); end
end
%DECK ISWAP
|
|
