[isdcgresult,n,b,x,nelt,ia,ja,a,isym,msolve,itol,tol,itmax,iter,err,ierr,iunit,r,z,p,dz,rwork,iwork,ak,bk,bnrm,solnrm]=isdcg(n,b,x,nelt,ia,ja,a,isym,msolve,itol,tol,itmax,iter,err,ierr,iunit,r,z,p,dz,rwork,iwork,ak,bk,bnrm,solnrm);

function [isdcgresult,n,b,x,nelt,ia,ja,a,isym,msolve,itol,tol,itmax,iter,err,ierr,iunit,r,z,p,dz,rwork,iwork,ak,bk,bnrm,solnrm]=isdcg(n,b,x,nelt,ia,ja,a,isym,msolve,itol,tol,itmax,iter,err,ierr,iunit,r,z,p,dz,rwork,iwork,ak,bk,bnrm,solnrm);
isdcgresult=[];
persistent i ; 

;
%***BEGIN PROLOGUE  ISDCG
%***SUBSIDIARY
%***PURPOSE  Preconditioned Conjugate Gradient Stop Test.
%            This routine calculates the stop test for the Conjugate
%            Gradient 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  D2B4
%***TYPE      doubleprecision (ISSCG-S, ISDCG-D)
%***KEYWORDS  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), Z(N)
%     doubleprecision P(N), DZ(N), RWORK(USER DEFINED), AK, BK
%     doubleprecision BNRM, SOLNRM
%     EXTERNAL MSOLVE
%
%     IF( ISDCG(N, B, X, NELT, IA, JA, A, ISYM, MSOLVE, ITOL, TOL,
%    $     ITMAX, ITER, ERR, IERR, IUNIT, R, Z, P, DZ, 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      :IN       doubleprecision X(N).
%         The current approximate solution vector.
% 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 should hold the matrix A in either the SLAP
%         Triad format or the SLAP Column format.  See 'Description'
%         in the DCG, DSDCG or DSICCG routines.
% 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.
% 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 /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.)
% ERR    :OUT      doubleprecision.
%         Error estimate of error in the X(N) 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.
% Z      :WORK     doubleprecision Z(N).
%         Workspace used to hold the pseudo-residual M Z = R.
% P      :IN       doubleprecision P(N).
%         The conjugate direction vector.
% 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.
% AK     :IN       doubleprecision.
% BK     :IN       doubleprecision.
%         Current conjugate gradient parameters alpha and beta.
% 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  DCG, DSDCG, DSICCG
%***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)
%   890921  Removed TeX from comments.  (FNF)
%   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 DCG.  (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)
%***end PROLOGUE  ISDCG
%     .. 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  ISDCG
isdcgresult = 0;
%
if( itol==1 )
%         err = ||Residual||/||RightHandSide|| (2-Norms).
if( iter==0 )
[ bnrm ,n,b]=dnrm2(n,b,1);
end;
err = dnrm2(n,r,1)./bnrm;
elseif( itol==2 ) ;
%                  -1              -1
%         err = ||M  Residual||/||M  RightHandSide|| (2-Norms).
if( iter==0 )
[n,r,z,nelt,ia,ja,a,isym,rwork,iwork]=msolve(n,r,z,nelt,ia,ja,a,isym,rwork,iwork);
[bnrm ,n,dz]=dnrm2(n,dz,1);
end;
err = dnrm2(n,z,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;
dz(i) = x(i) - dslblk_1(i);
end; i = fix(n+1);
err = dnrm2(n,dz,1)./solnrm;
else;
%
%         If we get here ITOL is not one of the acceptable values.
[err ]=d1mach(2);
ierr = 3;
end;
%
if( iunit~=0 )
if( iter==0 )
writef(iunit,[' Preconditioned Conjugate Gradient for ','N, ITOL = ','%5i','%5i', '\n ' ,' ITER','   Error Estimate','            Alpha','             Beta' ' \n'], n , itol);
%format (' Preconditioned Conjugate Gradient 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 )
isdcgresult = 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(18)), assignin('caller','FUntemp',z); evalin('caller',[inputname(18),'=FUntemp;']); end
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',x); evalin('caller',[inputname(3),'=FUntemp;']); end
if csnil&&~isempty(inputname(11)), assignin('caller','FUntemp',tol); evalin('caller',[inputname(11),'=FUntemp;']); end
if csnil&&~isempty(inputname(26)), assignin('caller','FUntemp',solnrm); evalin('caller',[inputname(26),'=FUntemp;']); end
if csnil&&~isempty(inputname(21)), assignin('caller','FUntemp',rwork); evalin('caller',[inputname(21),'=FUntemp;']); end
if csnil&&~isempty(inputname(17)), assignin('caller','FUntemp',r); evalin('caller',[inputname(17),'=FUntemp;']); end
if csnil&&~isempty(inputname(19)), assignin('caller','FUntemp',p); evalin('caller',[inputname(19),'=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(22)), assignin('caller','FUntemp',iwork); evalin('caller',[inputname(22),'=FUntemp;']); end
if csnil&&~isempty(inputname(16)), assignin('caller','FUntemp',iunit); evalin('caller',[inputname(16),'=FUntemp;']); end
if csnil&&~isempty(inputname(10)), assignin('caller','FUntemp',itol); evalin('caller',[inputname(10),'=FUntemp;']); end
if csnil&&~isempty(inputname(12)), assignin('caller','FUntemp',itmax); evalin('caller',[inputname(12),'=FUntemp;']); end
if csnil&&~isempty(inputname(13)), assignin('caller','FUntemp',iter); evalin('caller',[inputname(13),'=FUntemp;']); end
if csnil&&~isempty(inputname(8)), assignin('caller','FUntemp',isym); evalin('caller',[inputname(8),'=FUntemp;']); end
if csnil&&~isempty(inputname(15)), assignin('caller','FUntemp',ierr); evalin('caller',[inputname(15),'=FUntemp;']); end
if csnil&&~isempty(inputname(5)), assignin('caller','FUntemp',ia); evalin('caller',[inputname(5),'=FUntemp;']); end
if csnil&&~isempty(inputname(14)), assignin('caller','FUntemp',err); evalin('caller',[inputname(14),'=FUntemp;']); end
if csnil&&~isempty(inputname(20)), assignin('caller','FUntemp',dz); evalin('caller',[inputname(20),'=FUntemp;']); end
if csnil&&~isempty(inputname(25)), assignin('caller','FUntemp',bnrm); evalin('caller',[inputname(25),'=FUntemp;']); end
if csnil&&~isempty(inputname(24)), assignin('caller','FUntemp',bk); 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(23)), assignin('caller','FUntemp',ak); evalin('caller',[inputname(23),'=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 ISDCG 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(18)), assignin('caller','FUntemp',z); evalin('caller',[inputname(18),'=FUntemp;']); end
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',x); evalin('caller',[inputname(3),'=FUntemp;']); end
if csnil&&~isempty(inputname(11)), assignin('caller','FUntemp',tol); evalin('caller',[inputname(11),'=FUntemp;']); end
if csnil&&~isempty(inputname(26)), assignin('caller','FUntemp',solnrm); evalin('caller',[inputname(26),'=FUntemp;']); end
if csnil&&~isempty(inputname(21)), assignin('caller','FUntemp',rwork); evalin('caller',[inputname(21),'=FUntemp;']); end
if csnil&&~isempty(inputname(17)), assignin('caller','FUntemp',r); evalin('caller',[inputname(17),'=FUntemp;']); end
if csnil&&~isempty(inputname(19)), assignin('caller','FUntemp',p); evalin('caller',[inputname(19),'=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(22)), assignin('caller','FUntemp',iwork); evalin('caller',[inputname(22),'=FUntemp;']); end
if csnil&&~isempty(inputname(16)), assignin('caller','FUntemp',iunit); evalin('caller',[inputname(16),'=FUntemp;']); end
if csnil&&~isempty(inputname(10)), assignin('caller','FUntemp',itol); evalin('caller',[inputname(10),'=FUntemp;']); end
if csnil&&~isempty(inputname(12)), assignin('caller','FUntemp',itmax); evalin('caller',[inputname(12),'=FUntemp;']); end
if csnil&&~isempty(inputname(13)), assignin('caller','FUntemp',iter); evalin('caller',[inputname(13),'=FUntemp;']); end
if csnil&&~isempty(inputname(8)), assignin('caller','FUntemp',isym); evalin('caller',[inputname(8),'=FUntemp;']); end
if csnil&&~isempty(inputname(15)), assignin('caller','FUntemp',ierr); evalin('caller',[inputname(15),'=FUntemp;']); end
if csnil&&~isempty(inputname(5)), assignin('caller','FUntemp',ia); evalin('caller',[inputname(5),'=FUntemp;']); end
if csnil&&~isempty(inputname(14)), assignin('caller','FUntemp',err); evalin('caller',[inputname(14),'=FUntemp;']); end
if csnil&&~isempty(inputname(20)), assignin('caller','FUntemp',dz); evalin('caller',[inputname(20),'=FUntemp;']); end
if csnil&&~isempty(inputname(25)), assignin('caller','FUntemp',bnrm); evalin('caller',[inputname(25),'=FUntemp;']); end
if csnil&&~isempty(inputname(24)), assignin('caller','FUntemp',bk); 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(23)), assignin('caller','FUntemp',ak); evalin('caller',[inputname(23),'=FUntemp;']); end
if csnil&&~isempty(inputname(7)), assignin('caller','FUntemp',a); evalin('caller',[inputname(7),'=FUntemp;']); end
end
%DECK ISDCGN