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
