/* ========================================================================== */
/* === umf4 ================================================================= */
/* ========================================================================== */
/* -------------------------------------------------------------------------- */
/* UMFPACK Version 4.4, Copyright (c) 2005 by Timothy A. Davis. CISE Dept, */
/* Univ. of Florida. All Rights Reserved. See ../Doc/License for License. */
/* web: http://www.cise.ufl.edu/research/sparse/umfpack */
/* -------------------------------------------------------------------------- */
/* Demo program for UMFPACK. Reads in a triplet-form matrix in the
* directory tmp/A, whose size and # of nonzeros are in the file tmp/Asize.
* Then calls UMFPACK to analyze, factor, and solve the system.
*
* Syntax:
*
* umf4 default "auto" strategy, 1-norm row scaling
* umf4 a default "auto" strategy, 1-norm row scaling
* umf4 u unsymmetric strategy, 1-norm row scaling
* umf4 s symmetric strategy, 1-norm row scaling
* umf4 2 2-by-2 strategy, maxnorm row scaling
* umf4 A default "auto" strategy, maxnorm row scaling
* umf4 U unsymmetric strategy, maxnorm row scaling
* umf4 S symmetric strategy, maxnorm row scaling
* umf4 T 2-by-2 strategy , maxnorm row scaling
*
* To test a matrix in the Harwell/Boeing format, do the following:
*
* readhb < HB/arc130.rua > tmp/A
* readhb_size < HB/arc130.rua > tmp/Asize
* umf4
*
* The above options do not drop any nonzero entry in L or U. To compute an
* incomplete factorization, you can add a second argument to give the drop
* tolerance. Entries less than or equal to the drop tolerance are then
* removed from L and U during factorization, unless dropping those entries
* does not save any memory space. For example:
*
* umf4 a 1e-6 default "auto" strategy, 1-norm row scaling,
* drop tolerance of 1e-6.
*
* Note that adding a drop tolerance can lead to an apparent (but not real)
* increase in peak memory usage. This is illustrated in the arc130.rua
* matrix. With a drop tolerance, garbage collection happens to be avoided
* for this matrix. During garbage collection, both internal and external
* fragmentation in the memory space is removed. Peak memory usage includes
* all internal memory fragmentation, even though this can be removed via
* garbage collection.
*
* Control parameters can also be set in the optional tmp/control.umf4 file.
* The right-hand-side can be provided in the optional tmp/b file. The solution
* is written to tmp/x, and the output statistics are written to tmp/info.umf4.
*
* After the matrix is factorized, solved, and the LU factors deallocated,
* this program then test the AMD ordering routine. This call to AMD is NOT
* part of the UMFPACK analysis, factorize, or solve steps. It is just a
* separate test of the AMD ordering routine. If the matrix is unsymmetric,
* AMD orders the pattern of A+A'.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "umfpack.h"
#include "amd.h"
#define SMAX 256
#define ABS(x) ((x) >= 0 ? (x) : -(x))
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
#define XTRUE(i,n) (1.0 + ((double) i) / ((double) n))
#ifndef FALSE
#define FALSE 0
#endif
#ifndef TRUE
#define TRUE 1
#endif
/* -------------------------------------------------------------------------- */
/* err: compute the relative error, ||x-xtrue||/||xtrue|| */
/* -------------------------------------------------------------------------- */
static double err
(
int n,
double x [ ]
)
{
int i ;
double enorm, e, abse, absxtrue, xnorm ;
enorm = 0 ;
xnorm = 0 ;
for (i = 0 ; i < n ; i++)
{
if (isnan (x [i]))
{
enorm = x [i] ;
break ;
}
e = x [i] - XTRUE (i,n) ;
abse = ABS (e) ;
enorm = MAX (enorm, abse) ;
}
for (i = 0 ; i < n ; i++)
{
/* XTRUE is positive, but do this in case XTRUE is redefined */
absxtrue = ABS (XTRUE (i,n)) ;
xnorm = MAX (xnorm, absxtrue) ;
}
if (xnorm == 0)
{
xnorm = 1 ;
}
return (enorm / xnorm) ;
}
/* -------------------------------------------------------------------------- */
/* resid: compute the relative residual, ||Ax-b||/||b|| or ||A'x-b||/||b|| */
/* -------------------------------------------------------------------------- */
static double resid
(
int n,
int Ap [ ],
int Ai [ ],
double Ax [ ],
double x [ ],
double r [ ],
double b [ ],
int transpose
)
{
int i, j, p ;
double rnorm, absr, absb, bnorm ;
for (i = 0 ; i < n ; i++)
{
r [i] = 0 ;
}
if (transpose)
{
for (j = 0 ; j < n ; j++)
{
for (p = Ap [j] ; p < Ap [j+1] ; p++)
{
i = Ai [p] ;
r [j] += Ax [p] * x [i] ;
}
}
}
else
{
for (j = 0 ; j < n ; j++)
{
for (p = Ap [j] ; p < Ap [j+1] ; p++)
{
i = Ai [p] ;
r [i] += Ax [p] * x [j] ;
}
}
}
for (i = 0 ; i < n ; i++)
{
r [i] -= b [i] ;
}
rnorm = 0. ;
bnorm = 0. ;
for (i = 0 ; i < n ; i++)
{
if (isnan (r [i]))
{
rnorm = r [i] ;
break ;
}
absr = ABS (r [i]) ;
rnorm = MAX (rnorm, absr) ;
}
for (i = 0 ; i < n ; i++)
{
if (isnan (b [i]))
{
bnorm = b [i] ;
break ;
}
absb = ABS (b [i]) ;
bnorm = MAX (bnorm, absb) ;
}
if (bnorm == 0)
{
bnorm = 1 ;
}
return (rnorm / bnorm) ;
}
/* -------------------------------------------------------------------------- */
/* Atimesx: compute y = A*x or A'*x, where x (i) = 1 + i/n */
/* -------------------------------------------------------------------------- */
static void Atimesx
(
int n,
int Ap [ ],
int Ai [ ],
double Ax [ ],
double y [ ],
int transpose
)
{
int i, j, p ;
for (i = 0 ; i < n ; i++)
{
y [i] = 0 ;
}
if (transpose)
{
for (j = 0 ; j < n ; j++)
{
for (p = Ap [j] ; p < Ap [j+1] ; p++)
{
i = Ai [p] ;
y [j] += Ax [p] * XTRUE (i,n) ;
}
}
}
else
{
for (j = 0 ; j < n ; j++)
{
for (p = Ap [j] ; p < Ap [j+1] ; p++)
{
i = Ai [p] ;
y [i] += Ax [p] * XTRUE (j,n) ;
}
}
}
}
/* -------------------------------------------------------------------------- */
/* main program */
/* -------------------------------------------------------------------------- */
int main (int argc, char **argv)
{
int i, j, k, n, nz, *Ap, *Ai, *Ti, *Tj, status, *Pamd, nrow, ncol, rhs ;
double *Ax, *b, *x, Control [UMFPACK_CONTROL], Info [UMFPACK_INFO], aij,
*Tx, *r, amd_Control [AMD_CONTROL], amd_Info [AMD_INFO], tamd [2],
stats [2], droptol ;
void *Symbolic, *Numeric ;
FILE *f, *f2 ;
char s [SMAX] ;
/* ---------------------------------------------------------------------- */
/* set controls */
/* ---------------------------------------------------------------------- */
printf ("\n===========================================================\n"
"=== UMFPACK v4.4 ==========================================\n"
"===========================================================\n") ;
umfpack_di_defaults (Control) ;
Control [UMFPACK_PRL] = 3 ;
Control [UMFPACK_BLOCK_SIZE] = 32 ;
f = fopen ("tmp/control.umf4", "r") ;
if (f != (FILE *) NULL)
{
printf ("Reading control file tmp/control.umf4\n") ;
for (i = 0 ; i < UMFPACK_CONTROL ; i++)
{
fscanf (f, "%lg\n", & Control [i]) ;
}
fclose (f) ;
}
if (argc > 1)
{
char *s = argv [1] ;
/* get the strategy */
if (s [0] == 'u')
{
Control [UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC ;
}
else if (s [0] == 'a')
{
Control [UMFPACK_STRATEGY] = UMFPACK_STRATEGY_AUTO ;
}
else if (s [0] == 's')
{
Control [UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC ;
}
else if (s [0] == '2')
{
Control [UMFPACK_STRATEGY] = UMFPACK_STRATEGY_2BY2 ;
}
else if (s [0] == 'U')
{
Control [UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC ;
Control [UMFPACK_SCALE] = UMFPACK_SCALE_MAX ;
}
else if (s [0] == 'A')
{
Control [UMFPACK_STRATEGY] = UMFPACK_STRATEGY_AUTO ;
Control [UMFPACK_SCALE] = UMFPACK_SCALE_MAX ;
}
else if (s [0] == 'S')
{
Control [UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC ;
Control [UMFPACK_SCALE] = UMFPACK_SCALE_MAX ;
}
else if (s [0] == 'T')
{
Control [UMFPACK_STRATEGY] = UMFPACK_STRATEGY_2BY2 ;
Control [UMFPACK_SCALE] = UMFPACK_SCALE_MAX ;
}
else
{
printf ("unrecognized strategy: %s\n", argv [1]) ;
}
if (s [1] == 'n')
{
/* no aggressive absorption */
Control [UMFPACK_AGGRESSIVE] = FALSE ;
}
}
if (argc > 2)
{
/* get the drop tolerance */
sscanf (argv [2], "%lg", &droptol) ;
printf ("droptol %g\n", droptol) ;
Control [UMFPACK_DROPTOL] = droptol ;
}
umfpack_di_report_control (Control) ;
/* ---------------------------------------------------------------------- */
/* open the matrix file (tmp/A) */
/* ---------------------------------------------------------------------- */
printf ("File: tmp/A\n") ;
f = fopen ("tmp/A", "r") ;
if (!f)
{
printf ("Unable to open file\n") ;
exit (1) ;
}
/* ---------------------------------------------------------------------- */
/* get n and nz */
/* ---------------------------------------------------------------------- */
printf ("File: tmp/Asize\n") ;
f2 = fopen ("tmp/Asize", "r") ;
if (f2)
{
fscanf (f2, "%d %d %d\n", &nrow, &ncol, &nz) ;
fclose (f2) ;
}
else
{
nrow = 1 ;
ncol = 1 ;
}
nz = 0 ;
while (fgets (s, SMAX, f) != (char *) NULL)
{
sscanf (s, "%d %d %lg", &i, &j, &aij) ;
#ifdef ZERO_BASED
/* matrix is zero based */
i++ ;
j++ ;
#endif
nrow = MAX (nrow, i) ;
ncol = MAX (ncol, j) ;
nz++ ;
}
fclose (f) ;
n = MAX (nrow, ncol) ;
printf ("n %d nrow %d ncol %d nz %d\n", n, nrow, ncol, nz) ;
/* ---------------------------------------------------------------------- */
/* allocate space for the input triplet form */
/* ---------------------------------------------------------------------- */
Ti = (int *) malloc (nz * sizeof (int)) ;
Tj = (int *) malloc (nz * sizeof (int)) ;
Tx = (double *) malloc (nz * sizeof (double)) ;
if (!Ti || !Tj || !Tx)
{
printf ("out of memory for input matrix\n") ;
exit (1) ;
}
/* ---------------------------------------------------------------------- */
/* read in the triplet form */
/* ---------------------------------------------------------------------- */
f2 = fopen ("tmp/A", "r") ;
if (!f2)
{
printf ("Unable to open file\n") ;
exit (1) ;
}
k = 0 ;
while (fgets (s, SMAX, f2) != (char *) NULL)
{
sscanf (s, "%d %d %lg", &i, &j, &aij) ;
#ifndef ZERO_BASED
i-- ; /* convert to 0-based */
j-- ;
#endif
if (k >= nz)
{
printf ("Error! Matrix size is wrong\n") ;
exit (1) ;
}
Ti [k] = i ;
Tj [k] = j ;
Tx [k] = aij ;
k++ ;
}
fclose (f2) ;
(void) umfpack_di_report_triplet (nrow, ncol, nz, Ti, Tj, Tx, Control) ;
/* ---------------------------------------------------------------------- */
/* convert to column form */
/* ---------------------------------------------------------------------- */
/* convert to column form */
Ap = (int *) malloc ((n+1) * sizeof (int)) ;
Ai = (int *) malloc (nz * sizeof (int)) ;
Ax = (double *) malloc (nz * sizeof (double)) ;
b = (double *) malloc (n * sizeof (double)) ;
r = (double *) malloc (n * sizeof (double)) ;
x = (double *) malloc (n * sizeof (double)) ;
if (!Ap || !Ai || !Ax || !b || !r)
{
printf ("out of memory") ;
exit (1) ;
}
umfpack_tic (stats) ;
status = umfpack_di_triplet_to_col (nrow, ncol, nz, Ti, Tj, Tx, Ap, Ai, Ax,
(int *) NULL) ;
umfpack_toc (stats) ;
printf ("triplet-to-col time: wall %g cpu %g\n", stats [0], stats [1]) ;
if (status != UMFPACK_OK)
{
umfpack_di_report_status (Control, status) ;
printf ("umfpack_di_triplet_to_col failed") ;
exit (1) ;
}
/* print the column-form of A */
(void) umfpack_di_report_matrix (nrow, ncol, Ap, Ai, Ax, 1, Control) ;
/* b = A * xtrue */
rhs = FALSE ;
if (nrow == ncol)
{
f = fopen ("tmp/b", "r") ;
if (f != (FILE *) NULL)
{
printf ("Reading tmp/b\n") ;
rhs = TRUE ;
for (i = 0 ; i < n ; i++)
{
fscanf (f, "%lg\n", &b [i]) ;
}
fclose (f) ;
}
else
{
Atimesx (n, Ap, Ai, Ax, b, FALSE) ;
}
}
/* ---------------------------------------------------------------------- */
/* free the triplet form */
/* ---------------------------------------------------------------------- */
free (Ti) ;
free (Tj) ;
free (Tx) ;
/* ---------------------------------------------------------------------- */
/* symbolic factorization */
/* ---------------------------------------------------------------------- */
status = umfpack_di_symbolic (nrow, ncol, Ap, Ai, Ax, &Symbolic,
Control, Info) ;
umfpack_di_report_info (Control, Info) ;
if (status != UMFPACK_OK)
{
umfpack_di_report_status (Control, status) ;
printf ("umfpack_di_symbolic failed") ;
exit (1) ;
}
/* print the symbolic factorization */
(void) umfpack_di_report_symbolic (Symbolic, Control) ;
/* ---------------------------------------------------------------------- */
/* numeric factorization */
/* ---------------------------------------------------------------------- */
status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, Control, Info);
if (status < UMFPACK_OK)
{
umfpack_di_report_info (Control, Info) ;
umfpack_di_report_status (Control, status) ;
fprintf (stderr, "umfpack_di_numeric failed: %d\n", status) ;
printf ("umfpack_di_numeric failed\n") ;
exit (1) ;
}
/* print the numeric factorization */
(void) umfpack_di_report_numeric (Numeric, Control) ;
/* ---------------------------------------------------------------------- */
/* solve Ax=b */
/* ---------------------------------------------------------------------- */
if (nrow == ncol && status == UMFPACK_OK)
{
status = umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, x, b, Numeric,
Control, Info) ;
umfpack_di_report_info (Control, Info) ;
umfpack_di_report_status (Control, status) ;
if (status < UMFPACK_OK)
{
printf ("umfpack_di_solve failed\n") ;
exit (1) ;
}
(void) umfpack_di_report_vector (n, x, Control) ;
printf ("relative maxnorm of residual, ||Ax-b||/||b||: %g\n",
resid (n, Ap, Ai, Ax, x, r, b, FALSE)) ;
if (!rhs)
{
printf ("relative maxnorm of error, ||x-xtrue||/||xtrue||: %g\n\n",
err (n, x)) ;
}
f = fopen ("tmp/x", "w") ;
if (f != (FILE *) NULL)
{
printf ("Writing tmp/x\n") ;
for (i = 0 ; i < n ; i++)
{
fprintf (f, "%30.20e\n", x [i]) ;
}
fclose (f) ;
}
else
{
printf ("Unable to write output x!\n") ;
exit (1) ;
}
f = fopen ("tmp/info.umf4", "w") ;
if (f != (FILE *) NULL)
{
printf ("Writing tmp/info.umf4\n") ;
for (i = 0 ; i < UMFPACK_INFO ; i++)
{
fprintf (f, "%30.20e\n", Info [i]) ;
}
fclose (f) ;
}
else
{
printf ("Unable to write output info!\n") ;
exit (1) ;
}
}
else
{
/* don't solve, just report the results */
umfpack_di_report_info (Control, Info) ;
umfpack_di_report_status (Control, status) ;
}
/* ---------------------------------------------------------------------- */
/* free the Symbolic and Numeric factorization */
/* ---------------------------------------------------------------------- */
umfpack_di_free_symbolic (&Symbolic) ;
umfpack_di_free_numeric (&Numeric) ;
printf ("umf4 done, strategy: %g\n", Control [UMFPACK_STRATEGY]) ;
/* ---------------------------------------------------------------------- */
/* test just AMD ordering (not part of UMFPACK, but a separate test) */
/* ---------------------------------------------------------------------- */
/* first make the matrix square */
if (ncol < n)
{
for (j = ncol+1 ; j <= n ; j++)
{
Ap [j] = Ap [ncol] ;
}
}
printf (
"\n\n===========================================================\n"
"=== AMD ===================================================\n"
"===========================================================\n") ;
printf ("\n\n------- Now trying the AMD ordering. This not part of\n"
"the UMFPACK analysis or factorization, above, but a separate\n"
"test of just the AMD ordering routine.\n") ;
Pamd = (int *) malloc (n * sizeof (int)) ;
if (!Pamd)
{
printf ("out of memory\n") ;
exit (1) ;
}
amd_defaults (amd_Control) ;
amd_control (amd_Control) ;
umfpack_tic (tamd) ;
status = amd_order (n, Ap, Ai, Pamd, amd_Control, amd_Info) ;
umfpack_toc (tamd) ;
printf ("AMD ordering time: cpu %10.2f wall %10.2f\n",
tamd [1], tamd [0]) ;
if (status != AMD_OK)
{
printf ("amd failed: %d\n", status) ;
exit (1) ;
}
amd_info (amd_Info) ;
free (Pamd) ;
printf ("AMD test done\n") ;
return (0) ;
}
syntax highlighted by Code2HTML, v. 0.9.1