[isdcgsresult,n,b,x,nelt,ia,ja,a,isym,matvec,msolve,itol,tol,itmax,iter,err,ierr,iunit,r,r0,p,q,u,v1,v2,rwork,iwork,ak,bk,bnrm,solnrm]=isdcgs(n,b,x,nelt,ia,ja,a,isym,matvec,msolve,itol,tol,itmax,iter,err,ierr,iunit,r,r0,p,q,u,v1,v2,rwork,iwork,ak,bk,bnrm,

function [isdcgsresult,n,b,x,nelt,ia,ja,a,isym,matvec,msolve,itol,tol,itmax,iter,err,ierr,iunit,r,r0,p,q,u,v1,v2,rwork,iwork,ak,bk,bnrm,solnrm]=isdcgs(n,b,x,nelt,ia,ja,a,isym,matvec,msolve,itol,tol,itmax,iter,err,ierr,iunit,r,r0,p,q,u,v1,v2,rwork,iwork,ak,bk,bnrm,solnrm);
isdcgsresult=[];
persistent i ; 

;
%***BEGIN PROLOGUE  ISDCGS
%***SUBSIDIARY
%***PURPOSE  Preconditioned BiConjugate Gradient Squared Stop Test.
%            This routine calculates the stop test for the BiConjugate
%            Gradient Squared 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 (ISSCGS-S, ISDCGS-D)
%***KEYWORDS  ITERATIVE PRECONDITION, NON-SYMMETRIC LINEAR SYSTEM, 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 [phone]
%             [e-mail]
%***DESCRIPTION
%
% *Usage:
%     INTEGER  N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX, ITER
%     INTEGER  IERR, IUNIT, IWORK(USER DEFINED)
%     doubleprecision B(N), X(N), A(N), TOL, ERR, R(N), R0(N), P(N)
%     doubleprecision Q(N), U(N), V1(N), V2(N)
%     doubleprecision RWORK(USER DEFINED), AK, BK, BNRM, SOLNRM
%     EXTERNAL MATVEC, MSOLVE
%
%     IF( ISDCGS(N, B, X, NELT, IA, JA, A, ISYM, MATVEC, MSOLVE, ITOL,
%    $     TOL, ITMAX, ITER, ERR, IERR, IUNIT, R, R0, P, Q, U, V1,
%    $     V2, RWORK, IWORK, AK, BK, BNRM, SOLNRM) .NE. 0 )
%    $     THEN ITERATION DONE
%
% *Arguments:
% N      :IN       Integer
%         Order of the Matrix.
% B      :IN       doubleprecision B(N).
%         Right-hand side vector.
% X      :INOUT    doubleprecision 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       doubleprecision A(NELT).
%         These arrays contain the matrix data structure for A.
%         It could take any form.  See 'Description' in SLAP routine
%         DCGS 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.
% MATVEC :EXT      External.
%         Name of a routine which  performs the matrix vector multiply
%         operation  Y = A*X  given A and X.  The  name of  the MATVEC
%         routine must  be declared external  in the  calling program.
%         The calling sequence of MATVEC is:
%             CALL MATVEC( N, X, Y, NELT, IA, JA, A, ISYM )
%         Where N is the number of unknowns, Y is the product A*X upon
%         return,  X is an input  vector.  NELT, IA,  JA, A, and  ISYM
%         define the SLAP matrix data structure.
% 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 of 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 define the SLAP matrix data structure.
%         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.
% 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.
%         This routine must calculate the residual from R = A*X - B.
%         This is unnatural and hence expensive for this type of iter-
%         ative method.  ITOL=2 is *STRONGLY* recommended.
%         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 time a vector is the pre-
%         conditioning step.  This is the *NATURAL* stopping for this
%         iterative method and is *STRONGLY* recommended.
%         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.
%         Current iteration count.  (Must be zero on first call.)
%         ITMAX iterations.
% ERR    :OUT      doubleprecision.
%         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       doubleprecision R(N).
%         The residual r = b - Ax.
% R0     :WORK     doubleprecision R0(N).
% P      :DUMMY    doubleprecision P(N).
% Q      :DUMMY    doubleprecision Q(N).
% U      :DUMMY    doubleprecision U(N).
% V1     :DUMMY    doubleprecision V1(N).
%         doubleprecision arrays used for workspace.
% V2     :WORK     doubleprecision V2(N).
%         If ITOL==1 then V2 is used to hold A * X - B on every call.
%         If ITOL==2 then V2 is used to hold M-inv * B on the first
%         call.
%         If ITOL==11 then V2 is used to X - SOLN.
% RWORK  :WORK     doubleprecision RWORK(USER DEFINED).
%         doubleprecision 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       doubleprecision.
%         Current iterate BiConjugate Gradient iteration parameter.
% BK     :IN       doubleprecision.
%         Current iterate BiConjugate Gradient iteration parameter.
% BNRM   :INOUT    doubleprecision.
%         Norm of the right hand side.  Type of norm depends on ITOL.
%         Calculated only on the first call.
% SOLNRM :INOUT    doubleprecision.
%         2-Norm of the truemlv solution, SOLN.  Only computed and used
%         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  DCGS
%***ROUTINES CALLED  D1MACH, DNRM2
%***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)
%   891003  Removed C***REFER TO line, per MKS.
%   910411  Prologue converted to Version 4.0 format.  (BAB)
%   910502  Removed MATVEC and MSOLVE from ROUTINES CALLED list.  (FNF)
%   910506  Made subsidiary to DCGS.  (FNF)
%   920407  COMMON BLOCK renamed DSLBLK.  (WRB)
%   920511  Added complete declaration section.  (WRB)
%   920930  Corrected to not print AK,BK when ITER=0.  (FNF)
%   921026  Changed 1.0E10 to D1MACH(2) and corrected D to E in
%           output format.  (FNF)
%   921113  Corrected C***CATEGORY line.  (FNF)
%***end PROLOGUE  ISDCGS
%     .. Scalar Arguments ..
%     .. Array Arguments ..
rwork_shape=size(rwork);rwork=reshape(rwork,1,[]);
iwork_shape=size(iwork);iwork=reshape(iwork,1,[]);
%     .. subroutine Arguments ..
%     .. Arrays in Common ..
global dslblk_1; if isempty(dslblk_1), dslblk_1=zeros(1,1); end;
%     .. Local Scalars ..
if isempty(i), i=0; end;
%     .. External Functions ..
%     .. Common blocks ..
% common :: ;
%% common /dslblk/ soln;
%% common /dslblk/ dslblk_1;
%***FIRST EXECUTABLE STATEMENT  ISDCGS
isdcgsresult = 0;
%
if( itol==1 )
%         err = ||Residual||/||RightHandSide|| (2-Norms).
if( iter==0 )
[ bnrm ,n,b]=dnrm2(n,b,1);
end;
[n,x,v2,nelt,ia,ja,a,isym]=matvec(n,x,v2,nelt,ia,ja,a,isym);
for i = 1 : n;
v2(i) = v2(i) - b(i);
end; i = fix(n+1);
err = dnrm2(n,v2,1)./bnrm;
elseif( itol==2 ) ;
%                  -1              -1
%         err = ||M  Residual||/||M  RightHandSide|| (2-Norms).
if( iter==0 )
[n,b,v2,nelt,ia,ja,a,isym,rwork,iwork]=msolve(n,b,v2,nelt,ia,ja,a,isym,rwork,iwork);
[bnrm ,n,v2]=dnrm2(n,v2,1);
end;
err = dnrm2(n,r,1)./bnrm;
elseif( itol==11 ) ;
%         err = ||x-TrueSolution||/||TrueSolution|| (2-Norms).
if( iter==0 )
[ solnrm ,n,dslblk_1]=dnrm2(n,dslblk_1,1);
end;
for i = 1 : n;
v2(i) = x(i) - dslblk_1(i);
end; i = fix(n+1);
err = dnrm2(n,v2,1)./solnrm;
else;
%
%         If we get here ITOL is not one of the acceptable values.
[err ]=d1mach(2);
ierr = 3;
end;
%
%         Print the error and Coefficients AK, BK on each step,
%         if desired.
if( iunit~=0 )
if( iter==0 )
writef(iunit,[' Preconditioned BiConjugate Gradient Squared for ','N, ITOL = ','%5i','%5i', '\n ' ,' ITER','   Error Estimate','            Alpha','             Beta' ' \n'], n , itol);
%format (' Preconditioned BiConjugate Gradient Squared for ','N, ITOL = ',i5,i5,/' ITER','   Error Estimate','            Alpha','             Beta');
writef(iunit,[repmat(' ',1,1),'%4i',repmat(' ',1,1),'%16.7f',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',repmat(' ',1,1),'%16.7f' ' \n'], iter , err , ak , bk);
end;
end;
if( err<=tol )
isdcgsresult = 1;
end;
%
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(3)), assignin('caller','FUntemp',x); evalin('caller',[inputname(3),'=FUntemp;']); end
if csnil&&~isempty(inputname(24)), assignin('caller','FUntemp',v2); evalin('caller',[inputname(24),'=FUntemp;']); end
if csnil&&~isempty(inputname(23)), assignin('caller','FUntemp',v1); evalin('caller',[inputname(23),'=FUntemp;']); end
if csnil&&~isempty(inputname(22)), assignin('caller','FUntemp',u); evalin('caller',[inputname(22),'=FUntemp;']); end
if csnil&&~isempty(inputname(12)), assignin('caller','FUntemp',tol); evalin('caller',[inputname(12),'=FUntemp;']); end
if csnil&&~isempty(inputname(30)), assignin('caller','FUntemp',solnrm); evalin('caller',[inputname(30),'=FUntemp;']); end
if csnil&&~isempty(inputname(25)), assignin('caller','FUntemp',rwork); evalin('caller',[inputname(25),'=FUntemp;']); end
if csnil&&~isempty(inputname(19)), assignin('caller','FUntemp',r0); evalin('caller',[inputname(19),'=FUntemp;']); end
if csnil&&~isempty(inputname(18)), assignin('caller','FUntemp',r); evalin('caller',[inputname(18),'=FUntemp;']); end
if csnil&&~isempty(inputname(21)), assignin('caller','FUntemp',q); evalin('caller',[inputname(21),'=FUntemp;']); end
if csnil&&~isempty(inputname(20)), assignin('caller','FUntemp',p); evalin('caller',[inputname(20),'=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(10)), assignin('caller','FUntemp',msolve); evalin('caller',[inputname(10),'=FUntemp;']); end
if csnil&&~isempty(inputname(9)), assignin('caller','FUntemp',matvec); 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(29)), assignin('caller','FUntemp',bnrm); evalin('caller',[inputname(29),'=FUntemp;']); end
if csnil&&~isempty(inputname(28)), assignin('caller','FUntemp',bk); 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(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,d16.7,1x,d16.7,1x,d16.7);
%------------- LAST LINE OF ISDCGS FOLLOWS ----------------------------
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(3)), assignin('caller','FUntemp',x); evalin('caller',[inputname(3),'=FUntemp;']); end
if csnil&&~isempty(inputname(24)), assignin('caller','FUntemp',v2); evalin('caller',[inputname(24),'=FUntemp;']); end
if csnil&&~isempty(inputname(23)), assignin('caller','FUntemp',v1); evalin('caller',[inputname(23),'=FUntemp;']); end
if csnil&&~isempty(inputname(22)), assignin('caller','FUntemp',u); evalin('caller',[inputname(22),'=FUntemp;']); end
if csnil&&~isempty(inputname(12)), assignin('caller','FUntemp',tol); evalin('caller',[inputname(12),'=FUntemp;']); end
if csnil&&~isempty(inputname(30)), assignin('caller','FUntemp',solnrm); evalin('caller',[inputname(30),'=FUntemp;']); end
if csnil&&~isempty(inputname(25)), assignin('caller','FUntemp',rwork); evalin('caller',[inputname(25),'=FUntemp;']); end
if csnil&&~isempty(inputname(19)), assignin('caller','FUntemp',r0); evalin('caller',[inputname(19),'=FUntemp;']); end
if csnil&&~isempty(inputname(18)), assignin('caller','FUntemp',r); evalin('caller',[inputname(18),'=FUntemp;']); end
if csnil&&~isempty(inputname(21)), assignin('caller','FUntemp',q); evalin('caller',[inputname(21),'=FUntemp;']); end
if csnil&&~isempty(inputname(20)), assignin('caller','FUntemp',p); evalin('caller',[inputname(20),'=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(10)), assignin('caller','FUntemp',msolve); evalin('caller',[inputname(10),'=FUntemp;']); end
if csnil&&~isempty(inputname(9)), assignin('caller','FUntemp',matvec); 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(29)), assignin('caller','FUntemp',bnrm); evalin('caller',[inputname(29),'=FUntemp;']); end
if csnil&&~isempty(inputname(28)), assignin('caller','FUntemp',bk); 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(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 ISDGMR