Mercurial > octave-nkf
diff liboctave/UMFPACK/UMFPACK/Demo/umf4.c @ 5164:57077d0ddc8e
[project @ 2005-02-25 19:55:24 by jwe]
author | jwe |
---|---|
date | Fri, 25 Feb 2005 19:55:28 +0000 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/UMFPACK/UMFPACK/Demo/umf4.c Fri Feb 25 19:55:28 2005 +0000 @@ -0,0 +1,649 @@ +/* ========================================================================== */ +/* === 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) ; +}