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

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

;
%***BEGIN PROLOGUE  ISDBCG
%***SUBSIDIARY
%***PURPOSE  Preconditioned BiConjugate Gradient Stop Test.
%            This routine calculates the stop test for the BiConjugate
%            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  D2A4, D2B4
%***TYPE      doubleprecision (ISSBCG-S, ISDBCG-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), Z(N), P(N)
%     doubleprecision RR(N), ZZ(N), PP(N), DZ(N)
%     doubleprecision RWORK(USER DEFINED), AK, BK, BNRM, SOLNRM
%     EXTERNAL MSOLVE
%
%     IF( ISDBCG(N, B, X, NELT, IA, JA, A, ISYM, MSOLVE, ITOL, TOL,
%    $     ITMAX, ITER, ERR, IERR, IUNIT, R, Z, P, RR, ZZ, PP, 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      :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 the SLAP
%         routine DBCG 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 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.
%         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 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.
% Z      :WORK     doubleprecision Z(N).
% P      :DUMMY    doubleprecision P(N).
% RR     :DUMMY    doubleprecision RR(N).
% ZZ     :DUMMY    doubleprecision ZZ(N).
% PP     :DUMMY    doubleprecision PP(N).
%         doubleprecision arrays used for workspace.
% DZ     :WORK     doubleprecision DZ(N).
%         If ITOL==0 then DZ is used to hold M-inv * B on the first
%         call.  If ITOL==11 then DZ is used to hold X-SOLN.
% RWORK  :WORK     doubleprecision RWORK(USER DEFINED).
%         doubleprecision array that can be used for workspace in
%         MSOLVE and MTSOLV.
% IWORK  :WORK     Integer IWORK(USER DEFINED).
%         Integer array that can be used for workspace in MSOLVE
%         and MTSOLV.
% 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  DBCG
%***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 MSOLVE from ROUTINES CALLED list.  (FNF)
%   910506  Made subsidiary to DBCG.  (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  ISDBCG
%     .. 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  ISDBCG
isdbcgresult = 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,b,dz,nelt,ia,ja,a,isym,rwork,iwork]=msolve(n,b,dz,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 BiConjugate Gradient for N, ITOL = ','%5i','%5i', '\n ' ,' ITER','   Error Estimate','            Alpha','             Beta' ' \n'], n , itol);
%format (' Preconditioned BiConjugate 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 )
isdbcgresult = 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(21)), assignin('caller','FUntemp',zz); evalin('caller',[inputname(21),'=FUntemp;']); end
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(29)), assignin('caller','FUntemp',solnrm); evalin('caller',[inputname(29),'=FUntemp;']); end
if csnil&&~isempty(inputname(24)), assignin('caller','FUntemp',rwork); evalin('caller',[inputname(24),'=FUntemp;']); end
if csnil&&~isempty(inputname(20)), assignin('caller','FUntemp',rr); evalin('caller',[inputname(20),'=FUntemp;']); end
if csnil&&~isempty(inputname(17)), assignin('caller','FUntemp',r); evalin('caller',[inputname(17),'=FUntemp;']); end
if csnil&&~isempty(inputname(22)), assignin('caller','FUntemp',pp); evalin('caller',[inputname(22),'=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(25)), assignin('caller','FUntemp',iwork); evalin('caller',[inputname(25),'=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(23)), assignin('caller','FUntemp',dz); evalin('caller',[inputname(23),'=FUntemp;']); end
if csnil&&~isempty(inputname(28)), assignin('caller','FUntemp',bnrm); evalin('caller',[inputname(28),'=FUntemp;']); end
if csnil&&~isempty(inputname(27)), assignin('caller','FUntemp',bk); evalin('caller',[inputname(27),'=FUntemp;']); end
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',b); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(26)), assignin('caller','FUntemp',ak); evalin('caller',[inputname(26),'=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 ISDBCG 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(21)), assignin('caller','FUntemp',zz); evalin('caller',[inputname(21),'=FUntemp;']); end
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(29)), assignin('caller','FUntemp',solnrm); evalin('caller',[inputname(29),'=FUntemp;']); end
if csnil&&~isempty(inputname(24)), assignin('caller','FUntemp',rwork); evalin('caller',[inputname(24),'=FUntemp;']); end
if csnil&&~isempty(inputname(20)), assignin('caller','FUntemp',rr); evalin('caller',[inputname(20),'=FUntemp;']); end
if csnil&&~isempty(inputname(17)), assignin('caller','FUntemp',r); evalin('caller',[inputname(17),'=FUntemp;']); end
if csnil&&~isempty(inputname(22)), assignin('caller','FUntemp',pp); evalin('caller',[inputname(22),'=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(25)), assignin('caller','FUntemp',iwork); evalin('caller',[inputname(25),'=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(23)), assignin('caller','FUntemp',dz); evalin('caller',[inputname(23),'=FUntemp;']); end
if csnil&&~isempty(inputname(28)), assignin('caller','FUntemp',bnrm); evalin('caller',[inputname(28),'=FUntemp;']); end
if csnil&&~isempty(inputname(27)), assignin('caller','FUntemp',bk); evalin('caller',[inputname(27),'=FUntemp;']); end
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',b); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(26)), assignin('caller','FUntemp',ak); evalin('caller',[inputname(26),'=FUntemp;']); end
if csnil&&~isempty(inputname(7)), assignin('caller','FUntemp',a); evalin('caller',[inputname(7),'=FUntemp;']); end
end
%DECK ISDCG