C      ALGORITHM 695 , COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 17, NO. 3, SEPTEMBER, 1991, PP. 306-312.
c
c
c     Driver for  modified cholesky factorization algorithm.
c
c
      integer n,ndim
      double precision A(100,100)
      double precision Atwo(100,100)
      double precision g(100),E(100)
      double precision maxadd
      integer P(100)
      double precision y(100),b(100),x(100)
      double precision eps,tau1,tau2,sum
      integer ptran(100)
      integer z
      double precision high,low,temp
      logical errfnd
      ndim=100


C     mcheps subroutine computes machine precision,
C     The following line may be replaced by an assignment to eps
C     of  the correct machine precision constant for your machine.

      call mcheps(eps)

C     Tolerances used by modchl:
C     tau1 is used in determining when to switch to phase 2 and
C     tau2 is used in determining the amount to add to the diagonal
C     of the final 2X2 submatrix.

      tau1 = eps ** (1./3.)
      tau2 = eps ** (1./3.)

C     Initial seed for random number generator used to generate test
C     matrices
      z = 1000

C     high and low are the ranges of the eigenvalues for the test matrix
C     to be generated.

      high = 1.0
      low = -1.0

C     The first test problem will have dimension n=4, so that the entire
C     problem can be printed out.

      n = 4

      print *,'TEST PROBLEM #1'
      print *,'Test Matrix of size',n
      print *,'with eigenvalues within the range of ',low,' to ',high

      call mkmat(ndim,n,z,A,high,low,Atwo,g)

      print *,' '
      print *,'Original 4X4 matrix'

      do 25 i=1,n
  25     print 30,(A(i,j),j=1,n)
  30  format (4f20.8)

C     save the original matrix
      do 40 i=1,n
         do 40 j=1,n
 40         Atwo(i,j)=A(i,j)

      call modchl(ndim,n,A,g,eps,tau1,tau2,P,E)

      print *,' '
      print *,'Matrix after factorization with l in the lower triangle'

      do 50 i=1,n
  50     print 30,(A(i,j),j=1,n)

      print *,' '
      print *, 'Iteration  Permutation   Amt added to Aii'
      do 75 i=1,n
  75     print 80,i,P(i),e(i)
  80  format (i2,10x,i2,10x,f12.8)

      maxadd = E(n)

      print *,' '
      print *,'Maximum amount added to the diagonal is',maxadd

C
C     Generate  b for solve (using x(i)=i*2.0)
C
      do 90 i=1,n
         ptran(P(i))=i
 90   continue
      do 100 i=1,n
         Atwo(i,i)=Atwo(i,i)+E(ptran(i))
 100  continue
      do 120 i=1,n
         sum=0.0
         do 110 j=1,n
            sum= sum + Atwo(i,j)*(j*2.0)
 110     continue
         b(i)=sum
 120  continue

C
C     Call solve with A and P from modchl,
C     b is the input vector and x is the result vector
C        s.t.  solve computes (A+E)x=b.
C
      call solve(ndim,n,A,P,x,b,g,y)
C
C     compare the  correct answer to the solution found by solve
C
      errfnd = .false.
      do 130 i=1,n
         if (abs(x(i)-(i*2.0)) .gt. tau1) then
            if (.not. errfnd) then
               errfnd = .true.
               print *,' '
               print *,'Errors in Solve:'
               print *,'element #   correct answer   x (from solve)'
            end if
            temp = i * 2.0
            print 135,i,temp,x(i)
         end if
 130  continue
 135  format  (i3,10x,f8.4,7x,f8.4)


C     The next  test problem has size n=50
C     with eigenvalue range [-10000,-1].

      n = 50
      high = -1.0
      low = -10000.0

      print *,' '
      print *,'TEST PROBLEM #2'
      print *,'Test Matrix of size',n
      print *,'with eigenvalues in the range of ',low,' to ',high


      call mkmat(ndim,n,z,A,high,low,Atwo,g)

C     save the original matrix
      do 140 i=1,n
         do 140 j=1,n
 140         Atwo(i,j)=A(i,j)

      call modchl(ndim,n,A,g,eps,tau1,tau2,P,E)

      maxadd = E(n)

      print *,'Maximum amount added to the diagonal is',maxadd

C
C     Generate  b
C
      do 145 i=1,n
         b(i)=10.0 * i
 145  continue

C
C     Call solve
C
      call solve(ndim,n,A,P,x,b,g,y)
C
C     check the solution found by solve
C
      errfnd = .false.
      do 150 i=1,n
         ptran(P(i))=i
 150   continue
      do 160 i=1,n
         Atwo(i,i)=Atwo(i,i)+E(ptran(i))
 160  continue
      do 180 i=1,n
         sum=0.0
         do 170 j=1,n
            sum= sum + Atwo(i,j)*(x(j))
 170     continue
         if (abs(b(i)-sum) .gt. tau1) then
            if (.not. errfnd) then
               errfnd = .true.
               print *,' '
               print *,'Errors in Solve:'
               print *,'element #    b(i)         b(i)=(A+E)x'
               print *,'            (from input)  (x from solve)'
            end if
            print 135,i,b(i),sum
         end if
 180  continue

      stop
      end
c***********************************************************************
c       mcheps
c***********************************************************************
      subroutine mcheps(eps)
*
      double precision eps
*
      double precision temp
*
      temp = 1.0
*
 20   continue
         temp = temp / 2.0
      if ((1.0 + temp) .ne. 1.0) goto 20
*
      eps = temp * 2.0
*
      return
      end
*
C**********************************************************************
C
C     subroutine name : mkmat
C
C     purpose : Create an n dimensional matrix with eigenvalues
C             in the range of low to high by
C             forming the product q1*q2*q3*d*(q1*q2*q3)transpose,
C             where each qi is a Householder matrix, & each
C             diag element of d is in the desired eigenvalue
C             range.
C
C  input:   ndim  - largest dimension of matrix that will be used
C
C           n     - dimension of matrix
C
C           z     - initial seed for random number generator
C
C           high  - upper bound for eigenvalues of generated matrix
C
C           low   - lower bound for eigenvalues of generated matrix
C
C           q     - n*n work matrix
C
C           d     - n*1 work vector
C
C  output:  A     - generated matrix with eigenvalues in the range
C                   of low to high.
C
C           z     - next integer in the sequence generated by random
C                   number generator.
C
C
C**********************************************************************
*
      subroutine mkmat(ndim,n,z,A,high,low,q,v)
      integer n,ndim,z
      double precision A(ndim,n)
      double precision high,low
      double precision  q(ndim,n),v(n)
*
      integer i,j
      double precision r,drand,rand
      intrinsic abs
*
C     Make an orthonormal matrix A
      call mkorth(A,v,n,ndim,z)
*
*
C     Make an orthonormal matrix q
      call mkorth(q,v,n,ndim,z)
*
*
C     A = A * q
      call matmul(A,q,v,n,ndim)
*
*
C     Make an orthonormal matrix q
      call mkorth(q,v,n,ndim,z)
*
C     A = A * q
      call matmul(A,q,v,n,ndim)
*
*
C
C     q = A transpose
C
      do 20 i=1,n
         do 10 j=1,n
            q(i,j) = A(j,i)
  10     continue
  20  continue
*
*
C     Make a random vector v in the range of low to high
*
      r  = abs(high - low)
*
      do 30 i=1,n
         drand = rand(z)
         v(i) = low + r*drand
  30  continue
*
*
C     Make the first  diag element negative if range is big
*
      if ((high .gt. 100).and. (low .lt. 0)) then
         drand = rand(z)
         v(1) = -1.0 + drand
      end if
*
*
C     multiply  A = A * v
C     where v represents a diagonal matrix stored in a vector
*
      do 50 i=1,n
         do 40 j=1,n
            A(i,j) = A(i,j) * v(j)
  40     continue
  50  continue
*
C     A = A * q
      call matmul(A,q,v,n,ndim)
*
*
      return
      end
*
      subroutine mkorth(q,w,n,ndim,z)
C
C     purpose : make a Househoulder matrix by randomly generating
C               values in the range [-1,1] for an n dimension vector; w,
C               then computing the matrix
C               Q = I - (2/(2norm(w)**2)*w*(wtranspose)).
C     Input :   n,ndim,
C               z - seed for random number generator
C     Output :  w - an n*1 vector with values in the range [-1,1],
C               Q - a Househoulder matrix.
*
*
      integer n,ndim,z
      double precision q(ndim,n),w(n)
*
      double precision drand,norm2,rand
      integer j,k
*
*
      do 10 j=1,n
         drand  = rand(z)
         w(j) = -1.0 + (2.0 * drand)
  10  continue
*
      norm2 = 0.0
      do 20 j=1,n
            norm2 = norm2 + (w(j)**2)
  20  continue
      norm2 = 2.0/norm2
*
      do 40 j=1,n
         do 30 k=1,n
            q(j,k) = norm2 * w(j) * w(k)
            if (j .eq. k) then
               q(j,k) = 1.0 - q(j,k)
            else
                q(j,k) = 0.0 - q(j,k)
            end if
 30      continue
 40   continue
*
      return
      end
*
*
      subroutine matmul(a,b,v,n,ndim)
*
      double precision a(ndim,n), b(ndim,n), v(n)
      integer n,ndim
*
      double precision res
      integer i,j,k
*
      do 40 i=1,n
         do 20 j=1,n
             res = 0.0
             do 10 k=1,n
                res = res + a(i,k)*b(k,j)
  10         continue
             v(j) = res
  20     continue
         do 30 j=1,n
             a(i,j)=v(j)
  30     continue
  40  continue
*
      return
      end
*
*
C***********************************************************************
C     random number generator from:
C
C           Shrage, L.:A More Portable Random Number Generator,
C           ACM Trans. Math. Software, 5: 132-138(1979).
C
C     purpose: generates double precision numbers in the range of 0->1.
C
C     input:  ix  - initial value of seed
C     output: ix  - next integer in the random sequence
C                 in the range of 0 to 2^31 - 1.
C             rand - double precision number in the range of 0 -> 1.
C***********************************************************************
      double precision function rand(ix)
      integer a,p,ix,b15,b16,xhi,xalo,leftlo,fhi,k
      data a/16807/,b15/32768/,b16/65536/,p/2147483647/
*
      xhi=ix/b16
      xalo=(ix-xhi*b16)*a
      leftlo=xalo/b16
      fhi=xhi*a+leftlo
      k=fhi/b15
      ix=(((xalo-leftlo*b16)-p)+(fhi-k*b15)*b16)+k
      if (ix .lt. 0) ix=ix+p
      rand = dfloat(ix)*4.656612875e-10
      return
      end
C*********************************************************************
C
C       subroutine name: modchl
C
C       authors :  Elizabeth Eskow and Robert B. Schnabel
C
C       date    : December, 1988
C
C       purpose : perform a modified cholesky factorization
C                 of the form (Ptranspose)AP  + E = L(Ltranspose),
C       where L is stored in the lower triangle of the
C       original matrix A.
C       The factorization has 2 phases:
C        phase 1: Pivot on the maximum diagonal element.
C            Check that the normal cholesky update
C            would result in a positive diagonal
C            at the current iteration, and
C            if so, do the normal cholesky update,
C            otherwise switch to phase 2.
C        phase 2: Pivot on the minimum of the negatives
C            of the lower gerschgorin bound
C            estimates.
C            Compute the amount to add to the
C            pivot element and add this
C            to the pivot element.
C            Do the cholesky update.
C            Update the estimates of the
C            gerschgorin bounds.
C
C       input   : ndim    - largest dimension of matrix that
C                           will be used
C
C                 n       - dimension of matrix A
C
C                 A       - n*n symmetric matrix (only lower triangular
C            portion of A, including the main diagonal, is used)
C
C                 g       - n*1 work array
C
C                 mcheps - machine precision
C
C                tau1    - tolerance used for determining when to switch
C                          phase 2
C
C                tau2    - tolerance used for determining the maximum
C                          condition number of the final 2X2 submatrix.
C
C
C       output  : L     - stored in the matrix A (in lower triangular
C                           portion of A, including the main diagonal)
C
C                 P     - a record of how the rows and columns
C                         of the matrix were permuted while
C                         performing the decomposition
C
C                 E     - n*1 array, the ith element is the
C                         amount added to the diagonal of A
C                         at the ith iteration
C
C
C***********************************************************************
      subroutine modchl(ndim,n,A,g,mcheps,tau1,tau2,P,E)
*
      integer n,ndim
      double precision A(ndim,n),g(n),mcheps,tau1,tau2
      integer P(n)
      double precision E(n)
*
C
C  j        - current iteration number
C  iming    - index of the row with the min. of the
C           neg. lower Gersch. bounds
C  imaxd    - index of the row with the maximum diag.
C           element
C  i,itemp,jpl,k  - temporary integer variables
C  delta    - amount to add to Ajj at the jth iteration
C  gamma    - the maximum diagonal element of the original
C           matrix A.
C  normj    - the 1 norm of A(colj), rows j+1 --> n.
C  ming     - the minimum of the neg. lower Gersch. bounds
C  maxd     - the maximum diagonal element
C  taugam - tau1 * gamma
C  phase1      - logical, true if in phase1, otherwise false
C  delta1,temp,jdmin,tdmin,tempjj - temporary double precision vars.
C
*
      integer j,iming,i,imaxd,itemp,jp1,k
      double precision delta,gamma
      double precision normj, ming,maxd
      double precision delta1,temp,jdmin,tdmin,taugam,tempjj
      logical phase1
      intrinsic abs, max, sqrt, min
*
      call init(n, ndim, A, phase1, delta, P, g, E,
     *         ming,tau1,gamma,taugam)
C
C     check for n=1
C
      if (n.eq.1) then
         delta = (tau2 * abs(A(1,1))) - A(1,1)
         if (delta .gt. 0) E(1) = delta
         if (A(1,1) .eq. 0) E(1) = tau2
         A(1,1)=sqrt(A(1,1)+E(1))
      endif
C
      do 200 j = 1, n-1
C
C        PHASE 1
C
         if ( phase1 ) then
C
C           Find index of maximum diagonal element A(i,i) where i>=j
C
            maxd = A(j,j)
            imaxd = j
            do 20 i = j+1, n
               if (maxd .lt. A(i,i)) then
                  maxd = A(i,i)
                  imaxd = i
               end if
 20         continue
*
C
C           Pivot to the top the row and column with the max diag
C
            if (imaxd .ne. j) then
C
C              Swap row j with row of max diag
C
               do 30 i = 1, j-1
                  temp = A(j,i)
                  A(j,i) = A(imaxd,i)
                  A(imaxd,i) = temp
 30            continue
C
C              Swap colj and row maxdiag between j and maxdiag
C
               do 35 i = j+1,imaxd-1
                  temp = A(i,j)
                  A(i,j) = A(imaxd,i)
                  A(imaxd,i) = temp
 35            continue
C
C              Swap column j with column of max diag
C
               do 40 i = imaxd+1, n
                  temp = A(i,j)
                  A(i,j) = A(i,imaxd)
                  A(i,imaxd) = temp
 40            continue
C
C              Swap diag elements
C
               temp = A(j,j)
               A(j,j) = A(imaxd,imaxd)
               A(imaxd,imaxd) = temp
C
C              Swap elements of the permutation vector
C
               itemp = P(j)
               P(j) = P(imaxd)
               P(imaxd) = itemp
*
            end if
*
*
C           Check to see whether the normal cholesky update for this
C           iteration would result in a positive diagonal,
C           and if not then switch to phase 2.
*
            jp1 = j+1
            tempjj=A(j,j)
*
            if (tempjj.gt.0) then
*
               jdmin=A(jp1,jp1)
               do 60 i = jp1, n
                  temp = A(i,j) * A(i,j) / tempjj
                  tdmin = A(i,i) - temp
                  jdmin = min(jdmin, tdmin)
 60            continue
*
               if (jdmin .lt. taugam) phase1 = .false.
*
            else
*
               phase1 = .false.
*
            end if
*
            if (phase1) then
C
C              Do the normal cholesky update if still in phase 1
C
               A(j,j) = sqrt(A(j,j))
               tempjj = A(j,j)
               do 70 i = jp1, n
                  A(i,j) = A(i,j) / tempjj
 70            continue
               do 80 i=jp1,n
                  temp=A(i,j)
                  do 75 k = jp1, i
                     A(i,k) = A(i,k) - (temp * A(k,j))
 75               continue
 80            continue
*
               if (j .eq. n-1) A(n,n)=sqrt(A(n,n))
*
            else
*
C
C              Calculate the negatives of the lower gerschgorin bounds
C
               call gersch(ndim,n,A,j,g)
*
            end if
*
         end if
*
*
C
C        PHASE 2
C
         if (.not. phase1) then
*
            if (j .ne. n-1) then
C
C              Find the minimum negative gershgorin bound
C
*
               iming=j
               ming = g(j)
               do 90 i = j+1,n
                  if (ming .gt. g(i)) then
                     ming = g(i)
                     iming = i
                  end if
 90            continue
*
C
C               Pivot to the top the row and column with the
C               minimum negative gerschgorin bound
C
                if (iming .ne. j) then
C
C                  Swap row j with row of min gersch bound
C
                   do 100 i = 1, j-1
                      temp = A(j,i)
                       A(j,i) = A(iming,i)
                       A(iming,i) = temp
 100               continue
C
C                  Swap colj with row iming from j to iming
C
                   do 105 i = j+1,iming-1
                      temp = A(i,j)
                      A(i,j) = A(iming,i)
                      A(iming,i) = temp
 105              continue
C
C                 Swap column j with column of min gersch bound
C
                  do 110 i = iming+1, n
                     temp = A(i,j)
                     A(i,j) = A(i,iming)
                     A(i,iming) = temp
 110              continue
C
C                 Swap diagonal elements
C
                  temp = A(j,j)
                  A(j,j) = A(iming,iming)
                  A(iming,iming) = temp
C
C                 Swap elements of the permutation vector
C
                  itemp = P(j)
                  P(j) = P(iming)
                  P(iming) = itemp
C
C                 Swap elements of the negative gerschgorin bounds vecto
C
                  temp = g(j)
                  g(j) = g(iming)
                  g(iming) = temp
*
               end if
C
C              Calculate delta and add to the diagonal.
C              delta=max{0,-A(j,j) + max{normj,taugam},delta_previous}
C              where normj=sum of |A(i,j)|,for i=1,n,
C              delta_previous is the delta computed at the previous iter
C              and taugam is tau1*gamma.
C
*
               normj = 0.0
               do 140 i = j+1, n
                  normj = normj + abs(A(i,j))
 140           continue
*
               temp = max(normj,taugam)
               delta1 = temp - A(j,j)
               temp = 0.0
               delta1 = max(temp, delta1)
               delta = max(delta1,delta)
               E(j) =  delta
               A(j,j) = A(j,j) + E(j)
C
C              Update the gerschgorin bound estimates
C              (note: g(i) is the negative of the
C               Gerschgorin lower bound.)
C
               if (A(j,j) .ne. normj) then
                  temp = (normj/A(j,j)) - 1.0
*
                  do 150 i = j+1, n
                     g(i) = g(i) + abs(A(i,j)) * temp
 150              continue
*
               end if
C
C              Do the cholesky update
C
               A(j,j) = sqrt(A(j,j))
               tempjj = A(j,j)
               do 160 i = j+1, n
                  A(i,j) = A(i,j) / tempjj
 160           continue
               do 180 i = j+1, n
                  temp = A(i,j)
                  do 170 k = j+1, i
                     A(i,k) = A(i,k) - (temp * A(k,j))
 170              continue
 180           continue
*
            else
*
               call fin2x2(ndim, n, A, E, j, tau2, delta,gamma)
*
            end if
*
         end if
*
 200   continue
*
      return
      end
C***********************************************************************
C       subroutine name : init
C
C       purpose : set up for start of cholesky factorization
C
C       input : n, ndim, A, tau1
C
C       output : phase1    - boolean value set to true if in phase one,
C             otherwise false.
C      delta     - amount to add to Ajj at iteration j
C      P,g,E - described above in modchl
C      ming      - the minimum negative gerschgorin bound
C      gamma     - the maximum diagonal element of A
C      taugam  - tau1 * gamma
C
C***********************************************************************
      subroutine init(n,ndim,A,phase1,delta,P,g,E,ming,
     *                tau1,gamma,taugam)
*
      integer n,ndim
      double precision A(ndim,n)
      logical phase1
      double precision delta,g(n),E(n)
      integer P(n)
      double precision ming,tau1,gamma,taugam
      intrinsic abs, max
*
*
      phase1 = .true.
      delta = 0.0
      ming = 0.0
      do 10 i=1,n
         P(i)=i
         g(i)= 0.0
         E(i) = 0.0
 10   continue
*
C
C     Find the maximum magnitude of the diagonal elements.
C     If any diagonal element is negative, then phase1 is false.
C
      gamma = 0.0
      do 20 i=1,n
         gamma=max(gamma,abs(A(i,i)))
         if (A(i,i) .lt. 0.0) phase1 = .false.
 20   continue
*
      taugam = tau1 * gamma
*
C
C     If not in phase1, then calculate the initial gerschgorin bounds
C     needed for the start of phase2.
C
      if ( .not.phase1) call gersch(ndim,n,A,1,g)
*
      return
      end
C***********************************************************************
C
C       subroutine name : gersch
C
C       purpose : Calculate the negative of the gerschgorin bounds
C                 called once at the start of phase II.
C
C       input   : ndim, n, A, j
C
C       output  : g - an n vector containing the negatives of the
C           Gerschgorin bounds.
C
C***********************************************************************
      subroutine gersch(ndim, n, A, j, g)
*
      integer ndim, n, j
      double precision A(ndim,n), g(n)
*
      integer i, k
      double precision offrow
      intrinsic abs
*
      do 30 i = j, n
         offrow = 0.0
         do 10 k = j, i-1
            offrow = offrow + abs(A(i,k))
 10      continue
         do 20 k = i+1, n
            offrow = offrow + abs(A(k,i))
 20      continue
            g(i) = offrow - A(i,i)
 30   continue
*
      return
      end
C***********************************************************************
C
C  subroutine name : fin2x2
C
C  purpose : Handles final 2X2 submatrix in Phase II.
C            Finds eigenvalues of final 2 by 2 submatrix,
C            calculates the amount to add to the diagonal,
C            adds to the final 2 diagonal elements,
C            and does the final update.
C
C  input : ndim, n, A, E, j, tau2,
C          delta - amount added to the diagonal in the
C                  previous iteration
C
C  output : A - matrix with complete L factor in the lower triangle,
C           E - n*1 vector containing the amount added to the diagonal
C               at each iteration,
C           delta - amount added to diagonal elements n-1 and n.
C
C***********************************************************************
      subroutine fin2x2(ndim, n, A, E, j, tau2, delta,gamma)
*
      integer ndim, n, j
      double precision A(ndim,n), E(n), tau2, delta,gamma
*
      double precision t1, t2, t3,lmbd1,lmbd2,lmbdhi,lmbdlo
      double precision delta1, temp
      intrinsic sqrt, max, min
*
C
C     Find eigenvalues of final 2 by 2 submatrix
C
      t1 = A(n-1,n-1) + A(n,n)
      t2 = A(n-1,n-1) - A(n,n)
      t3 = sqrt(t2*t2 + 4.0*A(n,n-1)*A(n,n-1))
      lmbd1 = (t1 - t3)/2.
      lmbd2 = (t1 + t3)/2.
      lmbdhi = max(lmbd1,lmbd2)
      lmbdlo = min(lmbd1,lmbd2)
C
C     Find delta such that:
C     1.  the l2 condition number of the final
C     2X2 submatrix + delta*I <= tau2
C     2. delta >= previous delta,
C     3. lmbdlo + delta >= tau2 * gamma,
C     where lmbdlo is the smallest eigenvalue of the final
C     2X2 submatrix
C
*
      delta1=(lmbdhi-lmbdlo)/(1.0-tau2)
      delta1= max(delta1,gamma)
      delta1= tau2 * delta1 - lmbdlo
      temp = 0.0
      delta = max(delta, temp)
      delta = max(delta, delta1)
*
      if (delta .gt. 0.0) then
         A(n-1,n-1) = A(n-1,n-1) + delta
         A(n,n) = A(n,n) + delta
         E(n-1) = delta
         E(n) = delta
      end if
C
C     Final update
C
      A(n-1,n-1) = sqrt(A(n-1,n-1))
      A(n,n-1) = A(n,n-1)/A(n-1,n-1)
      A(n,n) = A(n,n) - (A(n,n-1)**2)
      A(n,n) = sqrt(A(n,n))
*
      return
      end
C**********************************************************************
C
C     subroutine name : solve
C
C     purpose :  solves (LLtranspose)P(x)=P(b), where L is stored in
C                the lower triangle of A, P is the record of the
C                permutations  performed in forming L,b is the rhs vecto
C                and x is the result.
C                L is the result of the modified cholesky factorization
C                which computes Ptranspose(A+E)P=LLtranspose.
C
C  input:   ndim  - largest dimension of matrix that will be used
C
C           n     - dimension of matrix
C
C           A     - n*n array (only lower triangle portion,
C                   including the diagonal, is used).
C
C           P     - n*1 integer vector that contains a record
C                   of permutations performed in forming L.
C                   i.e. each Pi was initialized to i, and if
C                   row and columns i & j were interchanged when
C                   computing L, the values of Pi and Pj were swapped.
C
C           b     - n*1 double precision vector which is the
C                   right hand side of the system to be solved.
C
C           y,z   - n*1 double precision  work vectors.
C
C  output:  x     - n*1 double precision result vector
C
C
C**********************************************************************
*
      subroutine solve(ndim,n,A,P,x,b,y,z)
      integer n,ndim
      double precision A(ndim,n)
      double precision x(n),b(n)
      integer P(n)
      double precision y(n),z(n)
*
      integer i,j
      double precision sum
C
C     Solve Ly=Pb
C
      y(1) = b(P(1))/A(1,1)
      do 20 i=2,n
         sum=0.0
         do 10 j=1,i-1
            sum=sum + (A(i,j)*y(j))
 10      continue
         y(i) = (b(P(i)) - sum) / A(i,i)
 20   continue
C
C     Solve Ltranspose z = y
C
      z(n) = y(n)/A(n,n)
      do 40 i=n-1,1,-1
         sum=0.0
         do 30 j=i+1,n
            sum = sum + (A(j,i) * z(j))
 30      continue
         z(i) = (y(i) - sum) / A(i,i)
 40   continue
C
C     x = Ptranspose z
C
      do 50 i=1,n
         x(P(i))=z(i)
 50   continue
*
      return
      end
*
*
*
*
TEST PROBLEM #1
Test Matrix of size  4
with eigenvalues within the range of    -1.0000000000000 to     1.000000
*
Original 4X4 matrix
          0.35711021         -0.10302945          0.02737268         -0.
         -0.10302945          0.25254612          0.07358379         -0.
          0.02737268          0.07358379          0.23396662         -0.
         -0.04594879         -0.38451624         -0.28782367          0.
*
Matrix after factorization with l in the lower triangle
          0.59758699         -0.10302945          0.02737268         -0.
         -0.07689054          0.82587804          0.07358379         -0.
          0.04580534         -0.34424172          0.49639272         -0.
         -0.17240912         -0.48163633         -0.16986202          0.
*
Iteration  Permutation   Amt added to Aii
 1           1            0.00000000
 2           4            0.13303961
 3           3            0.13303961
 4           2            0.13303961
*
Maximum amount added to the diagonal is   0.13303960618874
*
TEST PROBLEM #2
Test Matrix of size  50
with eigenvalues in the range of   -10000.0000000000 to    -1.0000000000
Maximum amount added to the diagonal is    11499.231418878
