c=======================================================================
c== umf4zhb ============================================================
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 umf4zhb:
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 CUA-type matrices.
c
c       This is HIGHLY non-portable.  It may not work with your C and
c       FORTRAN compilers.  See umf4z_f77wrapper.c for more details.
c
c usage (for example):
c
c       in a Unix shell:
c       umf4zhb < HB/arc130.cua

        integer
     $          nzmax, nmax
        parameter (nzmax = 5000000, nmax = 160000)
        integer
     $          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),
     $          control (20), info (90)
	complex*16 AA (nzmax), XX (nmax), BB (nmax), r (nmax), aij, xj
	double precision Az (nmax), xz (nmax), bz (nmax), xi, xr
        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. 'CUA' .or. nrow .ne. ncol) then
           print *, 'Error: can only handle square CUA 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) (AA (p), p = 1, nz)

	do 15 p = 1, nz
	    Ax (p) = dble (AA (p))
	    Az (p) = imag (AA (p))
15	continue

c       ----------------------------------------------------------------
c       create the right-hand-side, assume
c	x (i) = (1 + i/n), (n + i/100)
c       ----------------------------------------------------------------

        do 30 i = 1,n
            BB (i) = 0
30      continue
c       b = A*x
        do 50 j = 1,n
            xr = j
	    xi = n
	    xi = xi + xr/100
	    xr = 1 + xr / n
            xj = dcmplx (xr, xi)
            do 40 p = Ap (j), Ap (j+1)-1
                i = Ai (p)
                aij = AA (p)
                BB (i) = BB (i) + aij * xj
40          continue
50      continue
        do 32 i = 1,n
            b  (i) = dble (BB (i))
            bz (i) = imag (BB (i))
32      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 umf4zdef (control)

c       print control parameters.  set control (1) to 1 to print
c       error messages only
        control (1) = 2
        call umf4zpcon (control)

c       pre-order and symbolic analysis
        call umf4zsym (n, n, Ap, Ai, Ax, Az, symbolic, control, info)

c       print statistics computed so far
c       call umf4zpinf (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 umf4zsym error condition
        if (info (1) .lt. 0) then
            print *, 'Error occurred in umf4zsym: ', info (1)
            stop
        endif

c       numeric factorization
        call umf4znum (Ap, Ai, Ax, Az, symbolic, numeric, control, info)

c       print statistics for the numeric factorization
c       call umf4zpinf (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 umf4znum error condition
        if (info (1) .lt. 0) then
            print *, 'Error occurred in umf4znum: ', info (1)
            stop
        endif

c       save the symbolic analysis to the file s42.umf
c       note that this is not needed until another matrix is
c       factorized, below.
	filenum = 42
        call umf4zssym (symbolic, filenum, status)
        if (status .lt. 0) then
            print *, 'Error occurred in umf4zssym: ', status
            stop
        endif

c       save the LU factors to the file n0.umf
        call umf4zsnum (numeric, filenum, status)
        if (status .lt. 0) then
            print *, 'Error occurred in umf4zsnum: ', status
            stop
        endif

c       free the symbolic analysis
        call umf4zfsym (symbolic)

c       free the numeric factorization
        call umf4zfnum (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 umf4zlnum (numeric, filenum, status)
        if (status .lt. 0) then
            print *, 'Error occurred in umf4zlnum: ', status
            stop
        endif

c       solve Ax=b, without iterative refinement
        sys = 0
        call umf4zsol (sys, x, xz, b, bz, numeric, control, info)
        if (info (1) .lt. 0) then
            print *, 'Error occurred in umf4zsol: ', info (1)
            stop
        endif
        do 33 i = 1,n
            XX (i) = dcmplx (x (i), xz (i))
33      continue

c       free the numeric factorization
        call umf4zfnum (numeric)

c       No LU factors (symbolic or numeric) are in memory at this point.

c       print final statistics
        call umf4zpinf (control, info)

c       print the residual.  x (i) should be 1 + i/n
        call resid (n, nz, Ap, Ai, AA, XX, BB, 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, A, x, b, r)
        integer
     $      n, nz, Ap (n+1), Ai (n), j, i, p
        complex*16 A (nz), x (n), b (n), r (n), aij
	double precision rmax

        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 = A (p)
                r (i) = r (i) + aij * x (j)
20          continue
30      continue

        rmax = 0
        do 40 i = 1, n
            rmax = max (rmax, abs (r (i)))
40      continue

        print *, 'norm (A*x-b): ', rmax
        return
        end