[a,lda,n,v,itask,ind,work]=dpofs(a,lda,n,v,itask,ind,work);

function [a,lda,n,v,itask,ind,work]=dpofs(a,lda,n,v,itask,ind,work);
%***BEGIN PROLOGUE  DPOFS
%***PURPOSE  Solve a positive definite symmetric system of linear
%            equations.
%***LIBRARY   SLATEC
%***CATEGORY  D2B1B
%***TYPE      doubleprecision (SPOFS-S, DPOFS-D, CPOFS-C)
%***KEYWORDS  HERMITIAN, LINEAR EQUATIONS, POSITIVE DEFINITE, SYMMETRIC
%***AUTHOR  Voorhees, E. A., (LANL)
%***DESCRIPTION
%
%    subroutine DPOFS solves a  positive definite symmetric
%    NxN system of doubleprecision linear equations using
%    LINPACK subroutines DPOCO and DPOSL.  That is, if A is an
%    NxN doubleprecision positive definite symmetric matrix and if
%    X and B are doubleprecision N-vectors, then DPOFS solves
%    the equation
%
%                          A*X=B.
%
%    The matrix A is first factored into upper and lower tri-
%    angular matrices R and R-TRANPOSE.  These factors are used to
%    find the solution vector X.  An approximate condition number is
%    calculated to provide a rough estimate of the number of
%    digits of accuracy in the computed solution.
%
%    If the equation A*X=B is to be solved for more than one vector
%    B, the factoring of A does not need to be performed again and
%    the option only to solve (ITASK .GT. 1) will be faster for
%    the succeeding solutions.  In this case, the contents of A,
%    LDA, and N must not have been altered by the user following
%    factorization (ITASK=1).  IND will not be changed by DPOFS
%    in this case.
%
%  Argument Description ***
%
%    A      doubleprecision(LDA,N)
%             on entry, the doubly subscripted array with dimension
%               (LDA,N) which contains the coefficient matrix.  Only
%               the upper triangle, including the diagonal, of the
%               coefficient matrix need be entered and will subse-
%               quently be referenced and changed by the routine.
%             on return, A contains in its upper triangle an upper
%               triangular matrix R such that A = (R-TRANPOSE) * R .
%    LDA    INTEGER
%             the leading dimension of the array A.  LDA must be great-
%             er than or equal to N.  (terminal error message IND=-1)
%    N      INTEGER
%             the order of the matrix A.  N must be greater
%             than or equal to 1.  (terminal error message IND=-2)
%    V      doubleprecision(N)
%             on entry, the singly subscripted array(vector) of di-
%               mension N which contains the right hand side B of a
%               system of simultaneous linear equations  A*X=B.
%             on return, V contains the solution vector, X .
%    ITASK  INTEGER
%             If ITASK = 1, the matrix A is factored and then the
%               linear equation is solved.
%             If ITASK .GT. 1, the equation is solved using the existing
%               factored matrix A.
%             If ITASK .LT. 1, then terminal error message IND=-3 is
%               printed.
%    IND    INTEGER
%             GT. 0  IND is a rough estimate of the number of digits
%                     of accuracy in the solution, X.
%             LT. 0  See error message corresponding to IND below.
%    WORK   doubleprecision(N)
%             a singly subscripted array of dimension at least N.
%
%  Error Messages Printed ***
%
%    IND=-1  Terminal   N is greater than LDA.
%    IND=-2  Terminal   N is less than 1.
%    IND=-3  Terminal   ITASK is less than 1.
%    IND=-4  Terminal   The matrix A is computationally singular or
%                         is not positive definite.  A solution
%                         has not been computed.
%    IND=-10 Warning    The solution has no apparent significance.
%                         The solution may be inaccurate or the
%                         matrix A may be poorly scaled.
%
%               Note-  The above Terminal(*fatal*) Error Messages are
%                      designed to be handled by XERMSG in which
%                      LEVEL=1 (recoverable) and IFLAG=2 .  LEVEL=0
%                      for warning error messages from XERMSG.  Unless
%                      the user provides otherwise, an error message
%                      will be printed followed by an abort.
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  D1MACH, DPOCO, DPOSL, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   800514  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   890831  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%   900510  Convert XERRWV calls to XERMSG calls.  (RWC)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  DPOFS
%
persistent info rcond xern1 xern2 ; 

if isempty(info), info=0; end;
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([lda])).*prod([lda])-numel(a))],lda,[]);
v_shape=size(v);v=reshape(v,1,[]);
work_shape=size(work);work=reshape(work,1,[]);
if isempty(rcond), rcond=0; end;
if isempty(xern1), xern1=repmat(' ',1,8); end;
if isempty(xern2), xern2=repmat(' ',1,8); end;
%***FIRST EXECUTABLE STATEMENT  DPOFS
if( lda<n )
ind = -1;
xern1=sprintf(['%8i'], lda);
xern2=sprintf(['%8i'], n);
xermsg('SLATEC','DPOFS',['LDA = ',[xern1,[' IS LESS THAN N = ',xern2]]],-1,1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
%
if( n<=0 )
ind = -2;
xern1=sprintf(['%8i'], n);
xermsg('SLATEC','DPOFS',['N = ',[xern1,' IS LESS THAN 1']],-2,1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
%
if( itask<1 )
ind = -3;
xern1=sprintf(['%8i'], itask);
xermsg('SLATEC','DPOFS',['ITASK = ',[xern1,' IS LESS THAN 1']],-3,1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
%
if( itask==1 )
%
%        FACTOR MATRIX A INTO R
%
[a,lda,n,rcond,work,info]=dpoco(a,lda,n,rcond,work,info);
%
%        CHECK FOR POSITIVE DEFINITE MATRIX
%
if( info~=0 )
ind = -4;
xermsg('SLATEC','DPOFS','SINGULAR OR NOT POSITIVE DEFINITE - NO SOLUTION',-4,1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
%
%        COMPUTE IND (ESTIMATE OF NO. OF SIGNIFICANT DIGITS)
%        AND CHECK FOR IND GREATER THAN ZERO
%
ind = fix(-log10(d1mach(4)./rcond));
if( ind==0 )
ind = -10;
xermsg('SLATEC','DPOFS','SOLUTION MAY HAVE NO SIGNIFICANCE',-10,0);
end;
end;
%
%     SOLVE AFTER FACTORING
%
[a,lda,n,v]=dposl(a,lda,n,v);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
end
%DECK DPOLCF