'**************************************************************
'*       Program to test Cholesky Method for symmetric        *
'*               and positive definite matrices               *
'*                                                            *
'*                      Basic version by J-P Moreau, Paris    *
'*                              (www.jpmoreau.fr)             *
'* ---------------------------------------------------------- *
'* Reference:                                                 *
'*                                                            *
'*   "Numerical Algorithms with C, By Gisela Engeln-Muellges  *
'*    and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].     *
'* ---------------------------------------------------------- *
'* SAMPLE RUN:                                                *
'*                                                            *
'* Input file: matsym.dat                                     *
'*                                                            *
'* 4                                                          *
'* 5.0 -1.0 -1.0 -1.0  1.0                                    *
'*      5.0 -1.0 -1.0  1.0                                    *
'*           5.0 -1.0  1.0                                    *
'*                5.0  1.0                                    *
'*                                                            *
'* Output file: matsym.lst                                    *
'*                                                            *
'* ---------------------------------------------------------- *
'*      Solving a linear system by Cholesky Method            *
'* ---------------------------------------------------------- *
'*  The matrix of left hand coefficients must be a symmetric  *
'*  positive definite matrix (else error code = 2)            *
'*                                                            *
'*  Size of system: 4                                         *
'*                                                            *
'*  Input Linear System:                                      *
'*                                                            *
'*   5.000000   -1.000000   -1.000000   -1.000000    1.000000 *
'*                                                            *
'*  -1.000000    5.000000   -1.000000   -1.000000    1.000000 *
'*                                                            *
'*  -1.000000   -1.000000    5.000000   -1.000000    1.000000 *
'*                                                            *
'*  -1.000000   -1.000000   -1.000000    5.000000    1.000000 *
'*                                                            *
'*  Solution:                                                 *
'*                                                            *
'*   0.500000    0.500000    0.500000    0.500000             *
'* ---------------------------------------------------------- *
'*                                                            *
'**************************************************************
defint i-n
defdbl a-h,o-z

' Only diagonal and upper-diagonal elements are read

' main variables

' a, org : matrix (n,n)
' b, x   : vector (n)
' n, i, j, mode, irc : INTEGER
' finput$, foutput$, s$ : STRINGS


'begin main program

  EPSQUAD = 1e-12#
  ZERO = 0#

  cls
  print
  input " Input data file name (without *.dat): ", s$

  finput$  = s$ + ".dat"
  foutput$ = s$ + ".lst"

' open input, output files

  open finput$  for input  as #1
  open foutput$ for output as #2

' print header to output file

  print #2,"----------------------------------------------------------"
  print #2,"     Solving a linear system by Cholesky Method          "
  print #2,"----------------------------------------------------------"
  print #2," The matrix of left hand coefficients must be a symmetric "
  print #2," positive definite matrix (else error code = 2)           "
  print #2,

' read size n of system from input file
  input #1, n

  if n < 1 then
    print " n must be > 0 !"
    irc = 1
    end
  end if

  dim a(n,n), org(n,n), b(n), x(n)

  for i = 0 to n-1
    'read left hand coefficients of line i
    for j = i to n-1
      input #1, tmp
      org(i,j) = tmp
      org(j,i) = tmp    'symmetric term
    next j
    'read right hand coefficient of line i
    input #1, b(i)
  next i
  close #1

  'print data to output file

  print #2," Size of system: "; n
  print #2,
  print #2," Input Linear System:"
  print #2,

  f$="#####.######"

  for i = 0 to n-1
    for j = 0 to n-1
      a(i,j) = org(i,j)
      print #2,using f$; a(i,j);
    next j
    print #2,using f$; b(i)
    print #2,
  next i

  mode = 0    'Factor matrix a and solve linear system

  'call Cholesky procedure
  gosub 1000

  'print results to output file
  if irc = 0 then
    print #2," Solution:"
    print #2,
    for i=0 to n-1
      print #2,using f$; x(i);
    next i
  else
    print #2," ERROR choly: code = "; irc
  end if

  print #2,
  print #2,"----------------------------------------------------------"
  close #2

  print
  print " Results in file "; foutput$

END 'of main program


1000 'Subroutine Choly
'======================================================================
'*                                                                    *
'*  The cholesky subroutine solves linear systems :  a * x = b        *
'*  for positive definite symmetric n x n matrices mat using the      *
'*  Cholesky method.                                                  *
'*                                                                    *
'*  a must be symmetric and positive definite, or the method will     *
'*  fail. i.e. for all nonzero vectors y we must have y' * a * y > 0  *
'*                                                                    *
'*  cholesky uses only the lower triangle of matrix a.                *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Application:                                                     *
'*   ============                                                     *
'*                                                                    *
'*      Solve linear systems with positive definite system matrices   *
'*      efficeiently.                                                 *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Control parameter:                                               *
'*   ==================                                               *
'*      mode     integer                                              *
'*       = 0     Factor matrix a and solve linear system              *
'*       = 1     Compute factorization only                           *
'*       = 2     Solve system only; for this call the factorization   *
'*               must be stored in matrix a from chodec. Used to      *
'*               solve for many right hand sides.                     *
'*                                                                    *
'*   Input parameters:                                                *
'*   ================                                                 *
'*      n        integer (n > 0)                                      *
'*               Dimension of a, size of b and x                      *
'*      a        REAL  matrix (n,n)                                   *
'*                 mode = 0, 1: Matrix of the linear system           *
'*                 mode = 2   : Cholesky factor                       *
'*      b        REAL  vector (n)       (for mode = 0, 2)             *
'*               Right hand side of system of equations               *
'*                                                                    *
'*   Output parameters:                                               *
'*   ==================                                               *
'*      a        REAL matrix (n,n)      (for mode = 0, 1)             *
'*               Cholesky decomposition of input matrix mat           *
'*      x        REAL vector (n)        (for mode = 0, 2)             *
'*               solution vector                                      *
'*                                                                    *
'*   Return value rc:                                                 *
'*   ===============                                                  *
'*      = 0      all ok                                               *
'*      = 1      n < 1 or false input parameter                       *
'*      = 2      Matrix is not positive definite                      *
'*      = 3      wrong value for mode                                 *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Subroutines in use:                                              *
'*   ==================                                               *
'*      2000: chodec: Factorization                                   *
'*      3000: chosol: Solving system                                  *
'*                                                                    *
'======================================================================

  irc=3

  if n < 1 then
    irc=1
    return
  end if

  if mode=0 then  'Factor matrix and solve system
    'chodec (n, a, irc)
    gosub 2000
    if irc = 0 then
      'chosol(n, a, b, x, irc)
      gosub 3000
    end if
  end if

  if mode=1 then   'Factor only
    'chodec (n, a, irc)
    gosub 2000
  end if

  if mode=2 then   'Ssolve only
    'chosol (n, a, b, x, irc)
    gosub 3000
  end if

Return


2000 'Procedure chodec
'======================================================================
'*                                                                    *
'*  chodec decomposes the symmetric positive definite matrix mat.     *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Input parameters:                                                *
'*   ================                                                 *
'*      n        integer (n > 0)                                      *
'*               Dimension of matrix a                                *
'*      a        REAL matrix (n,n)                                    *
'*               Matrix of left hand coefficients                     *
'*                                                                    *
'*   Output parameter:                                                *
'*   ================                                                 *
'*      a        REAL matrix (n,n)                                    *
'*               Cholesky decomposition in the lower triangle         *
'*                                                                    *
'*   Return value rc:                                                 *
'*   ===============                                                  *
'*      = 0      all ok                                               *
'*      = 1      n < 1                                                *
'*      = 2      Matrix not  positive definite                        *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Functions in use:   SQR (square root in Basic)                   *
'*   ================                                                 *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Constants in use:   EPSQUAD (small number)                       *
'*   ================                                                 *
'*                                                                    *
'======================================================================
  if n < 1 then
    irc=1                                      ' n < 1  error !
    return
  end if

  if a(0,0) < EPSQUAD  then                    ' matrix a not positive definite
    irc=2
    return
  end if
  a(0,0) = SQR(a(0,0))
  for j = 1 to n-1
    a(j,0) = a(j,0) / a(0,0)
  next j

  for i = 1 to n-1
    sum = a(i,i)
    for j = 0 to i-1
      sum = sum - a(i,j)^2
    next j

    if sum < EPSQUAD then                      ' matrix a not positive definite
      irc=2
      return
    end if
    a(i,i) = SQR(sum)
    for j = i + 1 to n-1
      sum = a(j,i)
      for k = 0 to i-1
        sum = sum - a(i,k) * a(j,k)
      next k
      a(j,i) = sum / a(i,i)
    next j
  next i

  irc = 0
return


3000 'Procedure chosol
'======================================================================
'*                                                                    *
'*  chosol finds the solution x of the linear system B' *  B * x = b  *
'*  for a lower triangular nonsingular matrix B as supplied in chodec.*
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Input parameters:                                                *
'*   ================                                                 *
'*      n        integer  (n > 0)                                     *
'*               Dimension of matrix a, size of vectors b and x       *
'*      a        REAL  matrix (n,n)                                   *
'*               lower triangular matrix as supplied by  chodec       *
'*      b        REAL  vector b(n)                                    *
'*               Right hand side                                      *
'*                                                                    *
'*   Output parameter:                                                *
'*   =================                                                *
'*      x        REAL  vector x(n)                                    *
'*               solution vector                                      *
'*                                                                    *
'*   Return value rc:                                                 *
'*   ===============                                                  *
'*      = 0      all ok                                               *
'*      = 1      improper lwer triangular matrix or  n < 1            *
'*                                                                    *
'======================================================================
  if n < 1 then   ' n < 1 error !
    irc=1
    return
  end if

  if a(0,0) = ZERO then                  ' improper factor matrix
    irc=1
    return
  end if

  x(0) = b(0) / a(0,0)                   ' update right hand side
  for k = 1 to n-1
    sum = ZERO
    for j = 0 to k-1
      sum = sum + a(k,j) * x(j)
    next j
    if a(k,k) = ZERO then
      irc=1
      return
    end if
    x(k) = (b(k) - sum) / a(k,k)
  next k

  x(n-1) = x(n-1) / a(n-1,n-1)            ' back substitution
  for k=n-2 to 0 step -1
    sum = ZERO
    for j = k + 1 to n-1
      sum = sum + a(j,k) * x(j)
    next j
    x(k) = (x(k) - sum) / a(k,k)
  next k

  irc = 0
Return

'end of file tcholy1.bas