/* This is a modified version of the ldl.c file released by Timothy A. Davis */ /* in the LDL package and carrying the copyright shown below. The */ /* modifications are to replace scratch arrays passed as arguments by */ /* dynamically allocated arrays. */ /* Douglas Bates (Nov., 2004) */ /* ========================================================================== */ /* === ldl.c: sparse LDL' factorization and solve package =================== */ /* ========================================================================== */ /* LDL: a simple set of routines for sparse LDL' factorization. These routines * are not terrifically fast (they do not use dense matrix kernels), but the * code is very short. The purpose is to illustrate the algorithms in a very * concise manner, primarily for educational purposes. Although the code is * very concise, this package is slightly faster than the built-in sparse * Cholesky factorization in MATLAB 6.5 (chol), when using the same input * permutation. * * The routines compute the LDL' factorization of a real sparse symmetric * matrix A (or PAP' if a permutation P is supplied), and solve upper * and lower triangular systems with the resulting L and D factors. If A is * positive definite then the factorization will be accurate. A can be * indefinite (with negative values on the diagonal D), but in this case no * guarantee of accuracy is provided, since no numeric pivoting is performed. * * The n-by-n sparse matrix A is in compressed-column form. The nonzero values * in column j are stored in Ax [Ap [j] ... Ap [j+1]-1], with corresponding row * indices in Ai [Ap [j] ... Ap [j+1]-1]. Ap [0] = 0 is required, and thus * nz = Ap [n] is the number of nonzeros in A. Ap is an int array of size n+1. * The int array Ai and the double array Ax are of size nz. This data structure * is identical to the one used by MATLAB, except for the following * generalizations. The row indices in each column of A need not be in any * particular order, although they must be in the range 0 to n-1. Duplicate * entries can be present; any duplicates are summed. That is, if row index i * appears twice in a column j, then the value of A (i,j) is the sum of the two * entries. The data structure used here for the input matrix A is more * flexible than MATLAB's, which requires sorted columns with no duplicate * entries. * * Only the diagonal and upper triangular part of A (or PAP' if a permutation * P is provided) is accessed. The lower triangular parts of the matrix A or * PAP' can be present, but they are ignored. * * The optional input permutation is provided as an array P of length n. If * P [k] = j, the row and column j of A is the kth row and column of PAP'. * If P is present then the factorization is LDL' = PAP' or L*D*L' = A(P,P) in * 0-based MATLAB notation. If P is not present (a null pointer) then no * permutation is performed, and the factorization is LDL' = A. * * The lower triangular matrix L is stored in the same compressed-column * form (an int Lp array of size n+1, an int Li array of size Lp [n], and a * double array Lx of the same size as Li). It has a unit diagonal, which is * not stored. The row indices in each column of L are always returned in * ascending order, with no duplicate entries. This format is compatible with * MATLAB, except that it would be more convenient for MATLAB to include the * unit diagonal of L. Doing so here would add additional complexity to the * code, and is thus omitted in the interest of keeping this code short and * readable. * * The elimination tree is held in the Parent [0..n-1] array. It is normally * not required by the user, but it is required by R_ldl_numeric. The diagonal * matrix D is held as an array D [0..n-1] of size n. * * -------------------- * C-callable routines: * -------------------- * * R_ldl_symbolic: Given the pattern of A, computes the Lp and Parent * arrays required by R_ldl_numeric. Takes time proportional to the * number of nonzeros in L. Computes the inverse Pinv of P if P is * provided. Also returns Lnz, the count of nonzeros in each column * of L below the diagonal (this is not required by R_ldl_numeric). * R_ldl_numeric: Given the pattern and numerical values of A, the Lp * array, the Parent array, and P and Pinv if applicable, computes * the pattern and numerical values of L and D. * R_ldl_lsolve: Solves Lx=b for a dense vector b. * R_ldl_dsolve: Solves Dx=b for a dense vector b. * R_ldl_ltsolve: Solves L'x=b for a dense vector b. * R_ldl_perm: Computes x=Pb for a dense vector b. * R-ldl_permt: Computes x=P'b for a dense vector b. * R_ldl_valid_perm: checks the validity of a permutation vector * R_ldl_valid_matrix: checks the validity of the sparse matrix A * * ---------------------------- * Limitations of this package: * ---------------------------- * * In the interest of keeping this code simple and readable, * R_ldl_symbolic and R_ldl_numeric assume their inputs are valid. * You can check your own inputs prior to calling these routines with * the R_ldl_valid_perm and R_ldl_valid_matrix routines. Except for * the two R_ldl_valid_* routines, no routine checks to see if the * array arguments are present (non-NULL). Like all C routines, no * routine can determine if the arrays are long enough and don't * overlap. * * The R_ldl_numeric does check the numerical factorization, however. * It returns n if the factorization is successful. If D (k,k) is * zero, then k is returned, and L is only partially computed. * * No pivoting to control fill-in is performed, which is often * critical for obtaining good performance. I recommend that you * compute the permutation P using AMD or SYMAMD (approximate minimum * degree ordering routines), or an appropriate graph-partitioning * based ordering. See the ldldemo.m routine for an example in * MATLAB, and the ldlmain.c stand-alone C program for examples of how * to find P. Routines for manipulating compressed-column matrices * are available in UMFPACK. AMD, SYMAMD, UMFPACK, and this LDL * package are all available at * http://www.cise.ufl.edu/research/sparse. * * ------------------------- * Possible simplifications: * ------------------------- * * These routines could be made even simpler with a few additional * assumptions. If no input permutation were performed, the caller * would have to permute the matrix first, but the computation of * Pinv, and the use of P and Pinv could be removed. If only the * diagonal and upper triangular part of A or PAP' are present, then * the tests in the "if (i < k)" statement in R_ldl_symbolic and "if * (i <= k)" in R_ldl_numeric, are always true, and could be removed * (i can equal k in R_ldl_symbolic, but then the body of the if * statement would correctly do no work since Flag [k] == k). If we * could assume that no duplicate entries are present, then the * statement Y [i] += Ax [p] could be replaced with Y [i] = Ax [p] in * R_ldl_numeric. * * -------------------------- * Description of the method: * -------------------------- * * LDL computes the symbolic factorization by finding the pattern of L * one row at a time. It does this based on the following theory. * Consider a sparse system Lx=b, where L, x, and b, are all sparse, * and where L comes from a Cholesky (or LDL') factorization. The * elimination tree (etree) of L is defined as follows. The parent of * node j is the smallest k > j such that L (k,j) is nonzero. Node j * has no parent if column j of L is completely zero below the * diagonal (j is a root of the etree in this case). The nonzero * pattern of x is the union of the paths from each node i to the * root, for each nonzero b (i). To compute the numerical solution to * Lx=b, we can traverse the columns of L corresponding to nonzero * values of x. This traversal does not need to be done in the order * 0 to n-1. It can be done in any "topological" order, such that x * (i) is computed before x (j) if i is a descendant of j in the * elimination tree. * * The row-form of the LDL' factorization is shown in the MATLAB * function ldlrow.m in this LDL package. Note that row k of L is * found via a sparse triangular solve of L (1:k-1, 1:k-1) \ A (1:k-1, * k), to use 1-based MATLAB notation. Thus, we can start with the * nonzero pattern of the kth column of A (above the diagonal), follow * the paths up to the root of the etree of the (k-1)-by-(k-1) leading * submatrix of L, and obtain the pattern of the kth row of L. Note * that we only need the leading (k-1)-by-(k-1) submatrix of L to do * this. The elimination tree can be constructed as we go. * * The symbolic factorization does the same thing, except that it * discards the pattern of L as it is computed. It simply counts the * number of nonzeros in each column of L and then constructs the Lp * index array when it's done. The symbolic factorization does not * need to do this in topological order. Compare R_ldl_symbolic with * the first part of R_ldl_numeric, and note that the while (len > 0) * loop is not present in R_ldl_symbolic. * * LDL Version 1.0 (Dec. 31, 2003), Copyright (c) 2003 by Timothy A * Davis, University of Florida. All Rights Reserved. Developed * while on sabbatical at Stanford University and Lawrence Berkeley * National Laboratory. Refer to the README file for the License. * Available at http://www.cise.ufl.edu/research/sparse. */ #include "R_ldl.h" /* ========================================================================== */ /* === R_ldl_symbolic ======================================================= */ /* ========================================================================== */ /** * The input to this routine is a sparse matrix A, stored in column * form, and an optional permutation P. The output is the elimination * tree and the number of nonzeros in each column of L. Parent [i] = * k if k is the parent of i in the tree. The Parent array is * required by R_ldl_numeric. Lnz [k] gives the number of nonzeros in * the kth column of L, excluding the diagonal. * * If P is NULL, then it is ignored. The factorization will be LDL' = * A. Pinv is not computed. In this case, neither P nor Pinv are * required by R_ldl_numeric. * * If P is not NULL, then it is assumed to be a valid permutation. If * row and column j of A is the kth pivot, the P [k] = j. The * factorization will be LDL' = PAP', or A (p,p) in MATLAB notation. * The inverse permutation Pinv is computed, where Pinv [j] = k if P * [k] = j. In this case, both P and Pinv are required as inputs to * R_ldl_numeric. * * The floating-point operation count of the subsequent call to * R_ldl_numeric is not returned, but could be computed after * R_ldl_symbolic is done. It is the sum of (Lnz [k]) * (Lnz [k] + 2) * for k = 0 to n-1. * * @param n A and L are n-by-n, where n >= 0 * @param Ap column pointers of size n+1 * @param Ai row indices of size nz=Ap[n] * @param Lp column pointers of size n+1 * @param Parent elimination tree of size n * @param P optional permutation vector of size n [use (int *) NULL for none] * @param Pinv optional inverse permutation [not used if P is NULL] */ void R_ldl_symbolic(int n, const int Ap[], const int Ai[], int Lp[], int Parent[], const int P[], int Pinv[]) { int i, k, p, kk, p2; int *Flag = Calloc(n, int); int *Lnz = Calloc(n, int); if (P) { /* If P is present then compute Pinv, the inverse of P */ for (k = 0 ; k < n ; k++) { Pinv [P [k]] = k ; } } for (k = 0 ; k < n ; k++) { /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */ Parent [k] = -1 ; /* parent of k is not yet known */ Flag [k] = k ; /* mark node k as visited */ Lnz [k] = 0 ; /* count of nonzeros in column k of L */ kk = (P) ? (P [k]) : (k) ; /* kth original, or permuted, column */ p2 = Ap [kk+1] ; for (p = Ap [kk] ; p < p2 ; p++) { /* A (i,k) is nonzero (original or permuted A) */ i = (Pinv) ? (Pinv [Ai [p]]) : (Ai [p]) ; if (i < k) { /* follow path from i to root of etree, stop at flagged node */ for ( ; Flag [i] != k ; i = Parent [i]) { /* find parent of i if not yet determined */ if (Parent [i] == -1) { Parent [i] = k ; } Lnz [i]++ ; /* L (k,i) is nonzero */ Flag [i] = k ; /* mark i as visited */ } } } } /* construct Lp index array from Lnz column counts */ Lp [0] = 0 ; for (k = 0 ; k < n ; k++) { Lp [k+1] = Lp [k] + Lnz [k] ; } Free(Flag); Free(Lnz); } /** * Given a sparse matrix A (the arguments n, Ap, Ai, and Ax) and its * symbolic analysis (Lp and Parent, and optionally P and Pinv), * compute the numeric LDL' factorization of A or PAP'. The outputs * of this routine are arguments Li, Lx, and D. * * @param n A and L are n-by-n, where n >= 0 * @param Ap column pointer array of size n+1 * @param Ai row index array of size nz=Ap[n] (upper triangle only) * @param Ax array of non-zero matrix elements of size nz=Ap[n] * @param Lp column pointer array of size n+1 * @param Parent elimination tree of size n * @param Li row index array of size lnz=Lp[n] * @param Lx non-zero off-diagonal elements of L (size lnz=Lp[n]) * @param D vector of diagonal elements (size n) * @param P optional permutation vector of size n [use (int *) NULL for none] * @param Pinv optional inverse permutation [use (int *) NULL for none] * * @return n if successful, k if D (k,k) is zero */ int R_ldl_numeric(int n, const int Ap[], const int Ai[], const double Ax[], const int Lp[], const int Parent[], int Li[], double Lx[], double D[], const int P[], const int Pinv[]) { double yi, l_ki ; int i, k, p, kk, p2, len, top ; int *Lnz = Calloc(n, int), *Pattern = Calloc(n, int), *Flag = Calloc(n, int); double *Y = Calloc(n, double); for (k = 0 ; k < n ; k++) { /* compute nonzero Pattern of kth row of L, in topological order */ Y [k] = 0.0 ; /* Y (0:k) is now all zero */ top = n ; /* stack for pattern is empty */ Flag [k] = k ; /* mark node k as visited */ Lnz [k] = 0 ; /* count of nonzeros in column k of L */ kk = (P) ? (P [k]) : (k) ; /* kth original, or permuted, column */ p2 = Ap [kk+1] ; for (p = Ap [kk] ; p < p2 ; p++) { i = (Pinv) ? (Pinv [Ai [p]]) : (Ai [p]) ; /* get A(i,k) */ if (i <= k) { Y [i] += Ax [p] ; /* scatter A(i,k) into Y (sum duplicates) */ /* follow path from i to root of etree, stop at flagged node */ for (len = 0 ; Flag [i] != k ; i = Parent [i]) { Pattern [len++] = i ; /* L (k,i) is nonzero */ Flag [i] = k ; /* mark i as visited */ } while (len > 0) /* push path on top of stack */ { Pattern [--top] = Pattern [--len] ; } } } /* Pattern [top ... n-1] now contains nonzero pattern of L (:,k) */ /* compute numerical values kth row of L (a sparse triangular solve) */ D [k] = Y [k] ; /* get D (k,k) and clear Y (k) */ Y [k] = 0.0 ; for ( ; top < n ; top++) { i = Pattern [top] ; yi = Y [i] ; /* get and clear Y (i) */ Y [i] = 0.0 ; p2 = Lp [i] + Lnz [i] ; for (p = Lp [i] ; p < p2 ; p++) { Y [Li [p]] -= Lx [p] * yi ; } l_ki = yi / D [i] ; /* the nonzero entry L (k,i) */ D [k] -= l_ki * yi ; Li [p] = k ; /* store L(k,k )in column form of L */ Lx [p] = l_ki ; Lnz [i]++ ; /* increment count of nonzeros in col i */ } if (D [k] == 0.0) { Free(Y); Free(Pattern); Free(Flag); Free(Lnz); return (k) ; /* failure, D (k,k) is zero */ } } Free(Y); Free(Pattern); Free(Flag); Free(Lnz); return (n) ; /* success, diagonal of D is all nonzero */ } /** * solve Lx=b * * @param n L is n-by-n, where n >= 0 * @param X size n. right-hand-side on input, soln. on output * @param Lp column pointer array of size n+1 * @param Li row index array of size lnz=Lp[n] * @param Lx non-zero off-diagonal elements (size lnz=Lp[n]) */ void R_ldl_lsolve (int n, double X[], const int Lp[], const int Li[], const double Lx[]) { int j, p, p2 ; for (j = 0 ; j < n ; j++) { p2 = Lp [j+1] ; for (p = Lp [j] ; p < p2 ; p++) { X [Li [p]] -= Lx [p] * X [j] ; } } } /** * solve Dx=b * * @param n L is n-by-n, where n >= 0 * @param X size n. right-hand-side on input, soln. on output * @param D diagonal elements of size n */ void R_ldl_dsolve(int n, double X[], const double D[]) { int j ; for (j = 0 ; j < n ; j++) { X [j] /= D [j] ; } } /** * solve L'x=b * * @param n L is n-by-n, where n >= 0 * @param X size n. right-hand-side on input, soln. on output * @param Lp column pointer array of size n+1 * @param Li row index array of size lnz=Lp[n] * @param Lx non-zero off-diagonal elements (size lnz=Lp[n]) */ void R_ldl_ltsolve (int n, double X[], const int Lp[], const int Li[], const double Lx[]) { int j, p, p2 ; for (j = n-1 ; j >= 0 ; j--) { p2 = Lp [j+1] ; for (p = Lp [j] ; p < p2 ; p++) { X [j] -= Lx [p] * X [Li [p]] ; } } } /** * permute a vector, x=Pb * * @param n size of X, B, and P * @param X output of size n * @param B input of size n * @param P input permutation array of size n */ void R_ldl_perm(int n, double X[], const double B[], const int P[]) { int j ; for (j = 0 ; j < n ; j++) { X [j] = B [P [j]] ; } } /** * permute a vector, x=P'b * * @param n size of X, B, and P * @param X output of size n * @param B input of size n * @param P input permutation array of size n */ void R_ldl_permt(int n, double X[], const double B[], const int P[]) { int j ; for (j = 0 ; j < n ; j++) { X [P [j]] = B [j] ; } } /** * Check if a permutation vector is valid * * @param n size of permutation * @param P input of size n, a permutation of 0:n-1 * * @return 1 if valid, otherwise 0 */ int R_ldl_valid_perm (int n, const int P[]) { int j, k ; int *Flag = (int *) R_alloc(n, sizeof(int)); if (n < 0 || !Flag) { return (0) ; /* n must be >= 0, and Flag must be present */ } if (!P) { return (1) ; /* If NULL, P is assumed to be the identity perm. */ } for (j = 0 ; j < n ; j++) { Flag [j] = 0 ; /* clear the Flag array */ } for (k = 0 ; k < n ; k++) { j = P [k] ; if (j < 0 || j >= n || Flag [j] != 0) { return (0) ; /* P is not valid */ } Flag [j] = 1 ; } return (1) ; /* P is valid */ } /** * This routine checks to see if a sparse matrix A is valid for input * to R_ldl_symbolic and R_ldl_numeric. It returns 1 if the matrix is * valid, 0 otherwise. A is in sparse column form. The numerical * values in column j are stored in Ax [Ap [j] ... Ap [j+1]-1], with * row indices in Ai [Ap [j] ... Ap [j+1]-1]. The Ax array is not * checked. * * @param n A is n by n (n >= 0) * @param Ap column pointer array of size n+1 * @param Ai row index array of size nz=Ap[n] * * @return 1 if valid sparse matrix, otherwise 0 */ int R_ldl_valid_matrix (int n, const int Ap[], const int Ai[]) { int j, p ; if (n < 0 || !Ap || !Ai || Ap [0] != 0) { return (0) ; /* n must be >= 0, and Ap and Ai must be present */ } for (j = 0 ; j < n ; j++) { if (Ap [j] > Ap [j+1]) { return (0) ; /* Ap must be monotonically nondecreasing */ } } for (p = 0 ; p < Ap [n] ; p++) { if (Ai [p] < 0 || Ai [p] >= n) { return (0) ; /* row indices must be in the range 0 to n-1 */ } } return (1) ; /* matrix is valid */ }

Generated by Doxygen 1.6.0 Back to index