Mercurial > octave-nkf
diff liboctave/UMFPACK/AMD/Source/amd_preprocess.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/AMD/Source/amd_preprocess.c Fri Feb 25 19:55:28 2005 +0000 @@ -0,0 +1,225 @@ +/* ========================================================================= */ +/* === AMD_preprocess ====================================================== */ +/* ========================================================================= */ + +/* ------------------------------------------------------------------------- */ +/* AMD Version 1.1 (Jan. 21, 2004), Copyright (c) 2004 by Timothy A. Davis, */ +/* Patrick R. Amestoy, and Iain S. Duff. See ../README for License. */ +/* email: davis@cise.ufl.edu CISE Department, Univ. of Florida. */ +/* web: http://www.cise.ufl.edu/research/sparse/amd */ +/* ------------------------------------------------------------------------- */ + +/* Sorts, removes duplicate entries, and transposes from the nonzero pattern of + * a column-form matrix A, to obtain the matrix R. + * See amd.h for a complete description of AMD_preprocess + */ + +#include "amd_internal.h" + +GLOBAL Int AMD_preprocess /* returns AMD_OK if input is OK, AMD_INVALID + * if the matrix is invalid, or AMD_OUT_OF_MEMORY + * if out of memory for the 2n workspace. */ +( + Int n, /* input matrix: A is n-by-n */ + const Int Ap [ ], /* size n+1 */ + const Int Ai [ ], /* size nz = Ap [n] */ + + /* output matrix R: */ + Int Rp [ ], /* size n+1 */ + Int Ri [ ] /* size nz (or less, if duplicates present) */ +) +{ + /* --------------------------------------------------------------------- */ + /* local variables */ + /* --------------------------------------------------------------------- */ + + Int *Flag, *W ; + + /* --------------------------------------------------------------------- */ + /* check inputs (note: fewer restrictions than AMD_order) */ + /* --------------------------------------------------------------------- */ + + if (!AMD_preprocess_valid (n, Ap, Ai) || !Ri || !Rp) + { + return (AMD_INVALID) ; + } + + /* --------------------------------------------------------------------- */ + /* allocate workspace */ + /* --------------------------------------------------------------------- */ + + W = (Int *) ALLOCATE (MAX (n,1) * sizeof (Int)) ; + if (!W) + { + return (AMD_OUT_OF_MEMORY) ; + } + Flag = (Int *) ALLOCATE (MAX (n,1) * sizeof (Int)) ; + if (!Flag) + { + FREE (W) ; + return (AMD_OUT_OF_MEMORY) ; + } + + /* --------------------------------------------------------------------- */ + /* preprocess the matrix: sort, remove duplicates, and transpose */ + /* --------------------------------------------------------------------- */ + + AMD_wpreprocess (n, Ap, Ai, Rp, Ri, W, Flag) ; + + /* --------------------------------------------------------------------- */ + /* free the workspace */ + /* --------------------------------------------------------------------- */ + + FREE (W) ; + FREE (Flag) ; + return (AMD_OK) ; +} + + +/* ========================================================================= */ +/* === AMD_wpreprocess ===================================================== */ +/* ========================================================================= */ + +/* The AMD_wpreprocess routine is not user-callable. It does not check its + * input for errors or allocate workspace (that is done by the user-callable + * AMD_preprocess routine). It does handle the n=0 case. */ + +GLOBAL void AMD_wpreprocess +( + Int n, /* input matrix: A is n-by-n */ + const Int Ap [ ], /* size n+1 */ + const Int Ai [ ], /* size nz = Ap [n] */ + + /* output matrix R: */ + Int Rp [ ], /* size n+1 */ + Int Ri [ ], /* size nz (or less, if duplicates present) */ + + Int W [ ], /* workspace of size n */ + Int Flag [ ] /* workspace of size n */ +) +{ + + /* --------------------------------------------------------------------- */ + /* local variables */ + /* --------------------------------------------------------------------- */ + + Int i, j, p, p2 ; + + /* --------------------------------------------------------------------- */ + /* count the entries in each row of A (excluding duplicates) */ + /* --------------------------------------------------------------------- */ + + for (i = 0 ; i < n ; i++) + { + W [i] = 0 ; /* # of nonzeros in row i (excl duplicates) */ + Flag [i] = EMPTY ; /* Flag [i] = j if i appears in column j */ + } + for (j = 0 ; j < n ; j++) + { + p2 = Ap [j+1] ; + for (p = Ap [j] ; p < p2 ; p++) + { + i = Ai [p] ; + if (Flag [i] != j) + { + /* row index i has not yet appeared in column j */ + W [i]++ ; /* one more entry in row i */ + Flag [i] = j ; /* flag row index i as appearing in col j*/ + } + } + } + + /* --------------------------------------------------------------------- */ + /* compute the row pointers for R */ + /* --------------------------------------------------------------------- */ + + Rp [0] = 0 ; + for (i = 0 ; i < n ; i++) + { + Rp [i+1] = Rp [i] + W [i] ; + } + for (i = 0 ; i < n ; i++) + { + W [i] = Rp [i] ; + Flag [i] = EMPTY ; + } + + /* --------------------------------------------------------------------- */ + /* construct the row form matrix R */ + /* --------------------------------------------------------------------- */ + + /* R = row form of pattern of A */ + for (j = 0 ; j < n ; j++) + { + p2 = Ap [j+1] ; + for (p = Ap [j] ; p < p2 ; p++) + { + i = Ai [p] ; + if (Flag [i] != j) + { + /* row index i has not yet appeared in column j */ + Ri [W [i]++] = j ; /* put col j in row i */ + Flag [i] = j ; /* flag row index i as appearing in col j*/ + } + } + } + +#ifndef NDEBUG + for (j = 0 ; j < n ; j++) + { + ASSERT (W [j] == Rp [j+1]) ; + } + ASSERT (AMD_valid (n, n, Rp, Ri)) ; +#endif +} + + +/* ========================================================================= */ +/* === AMD_preprocess_valid ================================================ */ +/* ========================================================================= */ + +/* Not user-callable. Checks a matrix and returns TRUE if it is valid as input + * to AMD_wpreprocess, FALSE otherwise. */ + +GLOBAL Int AMD_preprocess_valid +( + Int n, + const Int Ap [ ], + const Int Ai [ ] +) +{ + Int i, j, p, nz ; + + if (n < 0 || !Ai || !Ap) + { + return (FALSE) ; + } + nz = Ap [n] ; + if (Ap [0] != 0 || nz < 0) + { + /* column pointers must start at Ap [0] = 0, and Ap [n] must be >= 0 */ + AMD_DEBUG0 (("column 0 pointer bad or nz < 0\n")) ; + return (FALSE) ; + } + for (j = 0 ; j < n ; j++) + { + if (Ap [j] > Ap [j+1]) + { + /* column pointers must be ascending */ + AMD_DEBUG0 (("column "ID" pointer bad\n", j)) ; + return (FALSE) ; + } + } + for (p = 0 ; p < nz ; p++) + { + i = Ai [p] ; + AMD_DEBUG3 (("row: "ID"\n", i)) ; + if (i < 0 || i >= n) + { + /* row index out of range */ + AMD_DEBUG0 (("index out of range, col "ID" row "ID"\n", j, i)) ; + return (FALSE) ; + } + } + return (TRUE) ; +}