--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/UMFPACK/UMFPACK/Demo/umf4hb64.f	Fri Feb 25 19:55:28 2005 +0000
@@ -0,0 +1,353 @@
+c=======================================================================
+c== umf4hb =============================================================
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c UMFPACK Version 4.4, Copyright (c) 2005 by Timothy A. Davis.  CISE
+c Dept, Univ. of Florida.  All Rights Reserved.  See ../Doc/License for
+c License.  web: http://www.cise.ufl.edu/research/sparse/umfpack
+c-----------------------------------------------------------------------
+
+c umf4hb64:
+c       read a sparse matrix in the Harwell/Boeing format, factorizes
+c       it, and solves Ax=b.  Also saves and loads the factors to/from a
+c       file.  Saving to a file is not required, it's just here to
+c       demonstrate how to use this feature of UMFPACK.  This program
+c       only works on square RUA-type matrices.
+c
+c       This is HIGHLY non-portable.  It may not work with your C and
+c       FORTRAN compilers.  See umf4_f77wrapper.c for more details.
+c
+c usage (for example):
+c
+c       in a Unix shell:
+c       umf4hb64 < HB/arc130.rua
+
+        integer*8
+     $          nzmax, nmax
+        parameter (nzmax = 5000000, nmax = 160000)
+        integer*8
+     $          Ap (nmax), Ai (nzmax), n, nz, totcrd, ptrcrd, i, j, p,
+     $          indcrd, valcrd, rhscrd, ncol, nrow, nrhs, nzrhs, nel,
+     $          numeric, symbolic, status, sys, filenum
+
+        character title*72, key*30, type*3, ptrfmt*16,
+     $          indfmt*16, valfmt*20, rhsfmt*20
+        double precision Ax (nzmax), x (nmax), b (nmax), aij, xj,
+     $          r (nmax), control (20), info (90)
+        character rhstyp*3
+
+c       ----------------------------------------------------------------
+c       read the Harwell/Boeing matrix
+c       ----------------------------------------------------------------
+
+        read (5, 10, err = 998)
+     $          title, key,
+     $          totcrd, ptrcrd, indcrd, valcrd, rhscrd,
+     $          type, nrow, ncol, nz, nel,
+     $          ptrfmt, indfmt, valfmt, rhsfmt
+        if (rhscrd .gt. 0) then
+c          new Harwell/Boeing format:
+           read (5, 20, err = 998) rhstyp, nrhs, nzrhs
+           endif
+10      format (a72, a8 / 5i14 / a3, 11x, 4i14 / 2a16, 2a20)
+20      format (a3, 11x, 2i14)
+
+        print *, 'Matrix key: ', key
+
+        n = nrow
+        if (type .ne. 'RUA' .or. nrow .ne. ncol) then
+           print *, 'Error: can only handle square RUA matrices'
+           stop
+        endif
+        if (n .ge. nmax .or. nz .gt. nzmax) then
+           print *, ' Matrix too big!'
+           stop
+        endif
+
+c       read the matrix (1-based)
+        read (5, ptrfmt, err = 998) (Ap (p), p = 1, ncol+1)
+        read (5, indfmt, err = 998) (Ai (p), p = 1, nz)
+        read (5, valfmt, err = 998) (Ax (p), p = 1, nz)
+
+c       ----------------------------------------------------------------
+c       create the right-hand-side, assume x (i) = 1 + i/n
+c       ----------------------------------------------------------------
+
+        do 30 i = 1,n
+            b (i) = 0
+30      continue
+c       b = A*x
+        do 50 j = 1,n
+            xj = j
+            xj = 1 + xj / n
+            do 40 p = Ap (j), Ap (j+1)-1
+                i = Ai (p)
+                aij = Ax (p)
+                b (i) = b (i) + aij * xj
+40          continue
+50      continue
+
+c       ----------------------------------------------------------------
+c       convert from 1-based to 0-based
+c       ----------------------------------------------------------------
+
+        do 60 j = 1, n+1
+            Ap (j) = Ap (j) - 1
+60      continue
+        do 70 p = 1, nz
+            Ai (p) = Ai (p) - 1
+70      continue
+
+c       ----------------------------------------------------------------
+c       factor the matrix and save to a file
+c       ----------------------------------------------------------------
+
+c       set default parameters
+        call umf4def (control)
+
+c       print control parameters.  set control (1) to 1 to print
+c       error messages only
+        control (1) = 2
+        call umf4pcon (control)
+
+c       pre-order and symbolic analysis
+        call umf4sym (n, n, Ap, Ai, Ax, symbolic, control, info)
+
+c       print statistics computed so far
+c       call umf4pinf (control, info) could also be done.
+        print 80, info (1), info (16),
+     $      (info (21) * info (4)) / 2**20,
+     $      (info (22) * info (4)) / 2**20,
+     $      info (23), info (24), info (25)
+80      format ('symbolic analysis:',/,
+     $      '   status:  ', f5.0, /,
+     $      '   time:    ', e10.2, ' (sec)'/,
+     $      '   estimates (upper bound) for numeric LU:', /,
+     $      '   size of LU:    ', f10.2, ' (MB)', /,
+     $      '   memory needed: ', f10.2, ' (MB)', /,
+     $      '   flop count:    ', e10.2, /
+     $      '   nnz (L):       ', f10.0, /
+     $      '   nnz (U):       ', f10.0)
+
+c       check umf4sym error condition
+        if (info (1) .lt. 0) then
+            print *, 'Error occurred in umf4sym: ', info (1)
+            stop
+        endif
+
+c       numeric factorization
+        call umf4num (Ap, Ai, Ax, symbolic, numeric, control, info)
+
+c       print statistics for the numeric factorization
+c       call umf4pinf (control, info) could also be done.
+        print 90, info (1), info (66),
+     $      (info (41) * info (4)) / 2**20,
+     $      (info (42) * info (4)) / 2**20,
+     $      info (43), info (44), info (45)
+90      format ('numeric factorization:',/,
+     $      '   status:  ', f5.0, /,
+     $      '   time:    ', e10.2, /,
+     $      '   actual numeric LU statistics:', /,
+     $      '   size of LU:    ', f10.2, ' (MB)', /,
+     $      '   memory needed: ', f10.2, ' (MB)', /,
+     $      '   flop count:    ', e10.2, /
+     $      '   nnz (L):       ', f10.0, /
+     $      '   nnz (U):       ', f10.0)
+
+c       check umf4num error condition
+        if (info (1) .lt. 0) then
+            print *, 'Error occurred in umf4num: ', info (1)
+            stop
+        endif
+
+c       save the symbolic analysis to the file s0.umf
+c       note that this is not needed until another matrix is
+c       factorized, below.
+	filenum = 0
+        call umf4ssym (symbolic, filenum, status)
+        if (status .lt. 0) then
+            print *, 'Error occurred in umf4ssym: ', status
+            stop
+        endif
+
+c       save the LU factors to the file n0.umf
+        call umf4snum (numeric, filenum, status)
+        if (status .lt. 0) then
+            print *, 'Error occurred in umf4snum: ', status
+            stop
+        endif
+
+c       free the symbolic analysis
+        call umf4fsym (symbolic)
+
+c       free the numeric factorization
+        call umf4fnum (numeric)
+
+c       No LU factors (symbolic or numeric) are in memory at this point.
+
+c       ----------------------------------------------------------------
+c       load the LU factors back in, and solve the system
+c       ----------------------------------------------------------------
+
+c       At this point the program could terminate and load the LU
+C       factors (numeric) from the n0.umf file, and solve the
+c       system (see below).  Note that the symbolic object is not
+c       required.
+
+c       load the numeric factorization back in (filename: n0.umf)
+        call umf4lnum (numeric, filenum, status)
+        if (status .lt. 0) then
+            print *, 'Error occurred in umf4lnum: ', status
+            stop
+        endif
+
+c       solve Ax=b, without iterative refinement
+        sys = 0
+        call umf4sol (sys, x, b, numeric, control, info)
+        if (info (1) .lt. 0) then
+            print *, 'Error occurred in umf4sol: ', info (1)
+            stop
+        endif
+
+c       free the numeric factorization
+        call umf4fnum (numeric)
+
+c       No LU factors (symbolic or numeric) are in memory at this point.
+
+c       print final statistics
+        call umf4pinf (control, info)
+
+c       print the residual.  x (i) should be 1 + i/n
+        call resid (n, nz, Ap, Ai, Ax, x, b, r)
+
+c       ----------------------------------------------------------------
+c       load the symbolic analysis back in, and factorize a new matrix
+c       ----------------------------------------------------------------
+
+c       Again, the program could terminate here, recreate the matrix,
+c       and refactorize.  Note that umf4sym is not called.
+
+c       load the symbolic factorization back in (filename: s0.umf)
+        call umf4lsym (symbolic, filenum, status)
+        if (status .lt. 0) then
+            print *, 'Error occurred in umf4lsym: ', status
+            stop
+        endif
+
+c       arbitrarily change the values of the matrix but not the pattern
+        do 100 p = 1, nz
+            Ax (p) = Ax (p) + 3.14159 / 100.0
+100     continue
+
+c       numeric factorization of the modified matrix
+        call umf4num (Ap, Ai, Ax, symbolic, numeric, control, info)
+        if (info (1) .lt. 0) then
+            print *, 'Error occurred in umf4num: ', info (1)
+            stop
+        endif
+
+c       free the symbolic analysis
+        call umf4fsym (symbolic)
+
+c       create a new right-hand-side, assume x (i) = 7 - i/n
+        do 110 i = 1,n
+            b (i) = 0
+110     continue
+c       b = A*x, with the modified matrix A (note that A is now 0-based)
+        do 130 j = 1,n
+            xj = j
+            xj = 7 - xj / n
+            do 120 p = Ap (j) + 1, Ap (j+1)
+                i = Ai (p) + 1
+                aij = Ax (p)
+                b (i) = b (i) + aij * xj
+120         continue
+130     continue
+
+c       ----------------------------------------------------------------
+c       solve Ax=b, with iterative refinement
+c       ----------------------------------------------------------------
+
+        sys = 0
+        call umf4solr (sys, Ap, Ai, Ax, x, b, numeric, control, info)
+        if (info (1) .lt. 0) then
+            print *, 'Error occurred in umf4solr: ', info (1)
+            stop
+        endif
+
+c       print the residual.  x (i) should be 7 - i/n
+        call resid (n, nz, Ap, Ai, Ax, x, b, r)
+
+c       ----------------------------------------------------------------
+c       solve Ax=b, without iterative refinement, broken into steps
+c       ----------------------------------------------------------------
+
+c       the factorization is PAQ=LU, PRAQ=LU, or P(R\A)Q=LU.
+
+c       x = R*b (or x=R\b, or x=b, as appropriate)
+        call umf4scal (x, b, numeric, status)
+        if (status .lt. 0) then
+            print *, 'Error occurred in umf4scal: ', status
+            stop
+        endif
+
+c       solve P'Lr=x for r (using r as workspace)
+        sys = 3
+        call umf4sol (sys, r, x, numeric, control, info)
+        if (info (1) .lt. 0) then
+            print *, 'Error occurred in umf4sol: ', info (1)
+            stop
+        endif
+
+c       solve UQ'x=r for x
+        sys = 9
+        call umf4sol (sys, x, r, numeric, control, info)
+        if (info (1) .lt. 0) then
+            print *, 'Error occurred in umf4sol: ', info (1)
+            stop
+        endif
+
+c       free the numeric factorization
+        call umf4fnum (numeric)
+
+c       print the residual.  x (i) should be 7 - i/n
+        call resid (n, nz, Ap, Ai, Ax, x, b, r)
+
+        stop
+998     print *, 'Read error: Harwell/Boeing matrix'
+        stop
+        end
+
+c=======================================================================
+c== resid ==============================================================
+c=======================================================================
+
+c Compute the residual, r = Ax-b, its max-norm, and print the max-norm
+C Note that A is zero-based.
+
+        subroutine resid (n, nz, Ap, Ai, Ax, x, b, r)
+        integer*8
+     $      n, nz, Ap (n+1), Ai (n), j, i, p
+        double precision Ax (nz), x (n), b (n), r (n), rmax, aij
+
+        do 10 i = 1, n
+            r (i) = -b (i)
+10      continue
+
+        do 30 j = 1,n
+            do 20 p = Ap (j) + 1, Ap (j+1)
+                i = Ai (p) + 1
+                aij = Ax (p)
+                r (i) = r (i) + aij * x (j)
+20          continue
+30      continue
+
+        rmax = 0
+        do 40 i = 1, n
+            rmax = max (rmax, r (i))
+40      continue
+
+        print *, 'norm (A*x-b): ', rmax
+        return
+        end