!  File:      gsfct.double.f
!  Date:      12/31/01  (downloaded fall 2001)
!  Author:    C. Paciorek - taken off the web
!  Purpose:   Cholesky factorization of pos def matrix stored in sparse format
!             double precision version
!  Compile:   g77 -c -ffree-form
!  Usage:     
!  Arguments: 
!  Requires:  
!  Reference: sparspak.f90 from the web -- George & Liu 1981, see below

subroutine gsfct (n,ixlnz,lnz,xnzsub,nzsub,diag,link,first,temp,eps,war)
!
!*******************************************************************************
!
!! GS_FACTOR performs the symmetric factorization for a general sparse
!  system, stored in the compressed subscript data format.
!
!  Code taken from sparspak.f90 on the web (CJP 9/4/01) and amended to use
!  Schervish' generalized Choleski factor (Lockwood et al.)
!
!  This is the double precision version
!
!  Reference:
!
!    Alan George and J W Liu,
!    Computer Solution of Large Sparse Positive Definite Systems,
!    Prentice Hall, 1981.
!
!  Parameters:
!
!    Input, integer N, the order of the matrix.
!
!  ixlnz-index vector for lnz.  ixlnz(i) points to the
!  start of nonzeros in column i of factor l.
!  (xnzsub,nzsub)-the compressed subscript data
!  structure for factor l.
!
!  updated parameters-
!  lnz-on input,contains nonzeros of a,and on
!  return,the nonzeros of l.
!  diag-the diagonal of l overwrites that of a.
!
!  working parameters-
!  link-at step j,the list in
!  link(j),link(link(j)),.
!  consists of those columns that will modify
!  the column l(*,j).
!  first-temporary vector to point to the first
!  nonzero in each column that will be used
!  next for modification.
!
  integer n
!  CJP changed reals to double precision 9/26/01
  double precision diag(n)
  double precision diagj
  integer first(n)
  integer i
  integer ii
  integer istop
  integer istrt
  integer isub
  integer ixlnz(n+1)
  integer j
  integer k
  integer kfirst
  integer link(n)
  double precision ljk
  double precision lnz(*)
  integer newk
  integer nzsub(*)
  double precision temp(n)
  integer xnzsub(*)
  double precision eps
  integer war
  integer writeout
  integer numcol
!
!  Initialize working vectors.  
!
!  CJP changed these to for loops instead of single statements
      do 100 i= 1, n
  link(i) = 0
  temp(i) = 0.0
 100  continue
!
!  Compute column l(*,j) for j = 1,....,N.
!
  numcol=0
  do j = 1, n
!
!  For each column l(*,k) that affects l(*,j).
!
    diagj = 0.0
    newk = link(j)
    writeout=0
   20   continue

    k = newk
    if ( k == 0 ) then
      goto 40
    end if
    newk = link(k)
!
!  Outer product modification of L(*,J) by L(*,K) starting at FIRST(K) of L(*,K).
!
    kfirst = first(k)
    ljk = lnz(kfirst)
    diagj = diagj + ljk*ljk
      if(writeout==1) then
         write(*,*) diagj, ljk
         end if
    istrt = kfirst+1
    istop = ixlnz(k+1)-1

    if ( istop < istrt ) then
      goto 20
    end if
!
!  Before modification, update vectors first,
!  and link for future modification steps.
!
    first(k) = istrt
    i = xnzsub(k) + ( kfirst - ixlnz(k) ) + 1
    isub = nzsub(i)
    link(k) = link(isub)
    link(isub) = k
!
!  The actual mod is saved in vector temp.
!
    do ii = istrt, istop
      isub = nzsub(i)
      temp(isub) = temp(isub) + lnz(ii) * ljk
      if(writeout==1) then
         write( *,* ) temp(isub), lnz(ii), ljk
         end if
      i = i + 1
    end do

    goto 20
!
!  Apply the modifications accumulated in temp to column L(*,J).
!
   40   continue

    diagj = diag(j) - diagj

!    if ( diagj > eps ) then
!       write ( *, * ) 'GSFCT - matrix pos def in col: ',j,diagj
!       end if
    if ( diagj <= eps ) then

       if( war == 1 ) then
          war=j
!          write ( *, * ) ' '
!          write ( *, * ) 'GSFCT - mat not pos def in column: ', j,diagj
!          write ( *, * ) '  Zero or negative diagonal entry!'
!          write ( *, * ) '  DIAG(J) = ', diagj
!          write ( *, * ) '  for diagonal J = ', j
!          write ( *, * ) ' '
!          write ( *, * ) '  The matrix is not positive definite.'
      end if

       numcol=numcol+1
       diagj=0;             ! this block of code added by CP 9/4/01 for gen chol
       diag(j)=0;
       istrt = ixlnz(j)
       istop = ixlnz(j+1) - 1
       if ( istop >= istrt ) then
          first(j)=istrt
          i=xnzsub(j)
          isub=nzsub(i)
          link(j)=link(isub)
          link(isub)=j
          do ii = istrt, istop
             isub=nzsub(i)
             lnz(ii)=0
             temp(isub)=0.0
             i=i+1
          end do

       end if
       goto 60
              ! through here added by CP
    end if

    diagj = sqrt ( diagj )
    diag(j) = diagj
    istrt = ixlnz(j)
    istop = ixlnz(j+1) - 1

    if ( istop >= istrt ) then

      first(j) = istrt
      i = xnzsub(j)
      isub = nzsub(i)
      link(j) = link(isub)
      link(isub) = j

      do ii = istrt, istop
        isub = nzsub(i)
        lnz(ii) = ( lnz(ii) - temp(isub) ) / diagj
        if(writeout==1) then
           write(*,*) lnz(ii), temp(isub), diagj
           end if
        temp(isub) = 0.0
        i = i+1
      end do

    end if

 60   end do

      if(war>=1) then
         if(numcol>0) then
            write(*,*) 'zeroed out columns:',numcol
         end if
      end if

  return
end






