//------------------------------------------------------------------------------ // LDL/Demo/ldlmain.c: demo program for LDL //------------------------------------------------------------------------------ // LDL, Copyright (c) 2005-2022 by Timothy A. Davis. All Rights Reserved. // SPDX-License-Identifier: LGPL-2.1+ //------------------------------------------------------------------------------ /* LDLMAIN: this main program has two purposes. It provides an example of how * to use the LDL routines, and it tests the package. The output of this * program is in ldlmain.out (without AMD) and ldlamd.out (with AMD). If you * did not download and install AMD, then the compilation of this program will * not work with USE_AMD defined. Compare your output with ldlmain.out and * ldlamd.out. * * The program reads in a set of matrices, in the Matrix/ subdirectory. * The format of the files is as follows: * * one line with the matrix description * one line with 2 integers: n jumbled * n+1 lines, containing the Ap array of size n+1 (column pointers) * nz lines, containing the Ai array of size nz (row indices) * nz lines, containing the Ax array of size nz (numerical values) * n lines, containing the P permutation array of size n * * The Ap, Ai, Ax, and P data structures are described in ldl.c. * The jumbled flag is 1 if the matrix could contain unsorted columns and/or * duplicate entries, and 0 otherwise. * * Once the matrix is read in, it is checked to see if it is valid. Some * matrices are invalid by design, to test the error-checking routines. If * valid, the matrix factorized twice (A and P*A*P'). A linear * system Ax=b is set up and solved, and the residual computed. * If any system is not solved accurately, this test will fail. */ #include #include #include "ldl.h" #define NMATRICES 30 /* number of test matrices in Matrix/ directory */ #define LEN 200 /* string length */ #ifdef USE_AMD /* get AMD include file, if using AMD */ #include "amd.h" #define PROGRAM "ldlamd" #if (AMD_VERSION < SUITESPARSE_VER_CODE(3,3)) #error "LDL:AMD @LDL_VERSION_MAJOR@.@LDL_VERSION_MINOR@.@LDL_VERSION_SUB@ requires AMD 3.3.1 or later" #endif #else #define PROGRAM "ldlmain" #endif #define EXIT_ERROR exit (EXIT_FAILURE) ; #define EXIT_OK exit (EXIT_SUCCESS) ; /* -------------------------------------------------------------------------- */ /* ALLOC_MEMORY: allocate a block of memory */ /* -------------------------------------------------------------------------- */ #define ALLOC_MEMORY(p,type,size) \ p = (type *) calloc ((((size) <= 0) ? 1 : (size)) , sizeof (type)) ; \ if (p == (type *) NULL) \ { \ printf (PROGRAM ": out of memory\n") ; \ EXIT_ERROR ; \ } /* -------------------------------------------------------------------------- */ /* FREE_MEMORY: free a block of memory */ /* -------------------------------------------------------------------------- */ #define FREE_MEMORY(p,type) \ if (p != (type *) NULL) \ { \ free (p) ; \ p = (type *) NULL ; \ } /* -------------------------------------------------------------------------- */ /* stand-alone main program */ /* -------------------------------------------------------------------------- */ int main (void) { /* ---------------------------------------------------------------------- */ /* local variables */ /* ---------------------------------------------------------------------- */ #ifdef USE_AMD double Info [AMD_INFO] ; #endif double r, rnorm, flops, maxrnorm = 0. ; double *Ax, *Lx, *B, *D, *X, *Y ; LDL_int matrix, *Ai, *Ap, *Li, *Lp, *P, *Pinv, *Perm, *PermInv, n, i, j, p, nz, *Flag, *Pattern, *Lnz, *Parent, trial, lnz, d, jumbled, ok ; FILE *f ; char s [LEN], filename [LEN] ; //-------------------------------------------------------------------------- // check the LDL version //-------------------------------------------------------------------------- printf ("LDL version %d.%d.%d, date: %s\n", LDL_MAIN_VERSION, LDL_SUB_VERSION, LDL_SUBSUB_VERSION, LDL_DATE) ; int version [3] ; ldl_version (version) ; if ((version [0] != LDL_MAIN_VERSION) || (version [1] != LDL_SUB_VERSION) || (version [2] != LDL_SUBSUB_VERSION)) { fprintf (stderr, "version in header does not match library\n") ; abort ( ) ; } /* ---------------------------------------------------------------------- */ /* check the error-checking routines with null matrices */ /* ---------------------------------------------------------------------- */ i = 1 ; n = -1 ; if (LDL_valid_perm (n, (LDL_int *) NULL, &i) || !LDL_valid_perm (0, (LDL_int *) NULL, &i) || LDL_valid_matrix (n, (LDL_int *) NULL, (LDL_int *) NULL) || LDL_valid_matrix (0, &i, &i)) { printf (PROGRAM ": ldl error-checking routine failed\n") ; EXIT_ERROR ; } /* ---------------------------------------------------------------------- */ /* read in a factorize a set of matrices */ /* ---------------------------------------------------------------------- */ for (matrix = 1 ; matrix <= NMATRICES ; matrix++) { /* ------------------------------------------------------------------ */ /* read in the matrix and the permutation */ /* ------------------------------------------------------------------ */ sprintf (filename, "Matrix/A%02d", (int) matrix) ; if ((f = fopen (filename, "r")) == (FILE *) NULL) { printf (PROGRAM ": could not open file: %s\n", filename) ; EXIT_ERROR ; } printf ("\n\n--------------------------------------------------------"); printf ("\nInput file: %s\n", filename) ; s [0] = 0 ; ok = (fgets (s, LEN, f) != NULL) ; printf ("%s", s) ; printf ("--------------------------------------------------------\n\n"); n = 0 ; if (ok) { ok = ok && (fscanf (f, LDL_ID " " LDL_ID, &n, &jumbled) == 2) ; n = (n < 0) ? (0) : (n) ; ALLOC_MEMORY (P, LDL_int, n) ; ALLOC_MEMORY (Ap, LDL_int, n+1) ; } for (j = 0 ; ok && j <= n ; j++) { ok = ok && (fscanf (f, LDL_ID, &Ap [j]) == 1) ; } nz = 0 ; if (ok) { nz = Ap [n] ; ALLOC_MEMORY (Ai, LDL_int, nz) ; ALLOC_MEMORY (Ax, double, nz) ; } for (p = 0 ; ok && p < nz ; p++) { ok = ok && (fscanf (f, LDL_ID , &Ai [p]) == 1) ; } for (p = 0 ; ok && p < nz ; p++) { ok = ok && (fscanf (f, "%lg", &Ax [p]) == 1) ; } for (j = 0 ; ok && j < n ; j++) { ok = ok && (fscanf (f, LDL_ID , &P [j]) == 1) ; } fclose (f) ; if (!ok) { printf (PROGRAM ": error reading file: %s\n", filename) ; EXIT_ERROR ; } /* ------------------------------------------------------------------ */ /* check the matrix A and the permutation P */ /* ------------------------------------------------------------------ */ ALLOC_MEMORY (Flag, LDL_int, n) ; /* To test the error-checking routines, some of the input matrices * are not valid. So this error is expected to occur. */ if (!LDL_valid_matrix (n, Ap, Ai) || !LDL_valid_perm (n, P, Flag)) { printf (PROGRAM ": invalid matrix and/or permutation\n") ; FREE_MEMORY (P, LDL_int) ; FREE_MEMORY (Ap, LDL_int) ; FREE_MEMORY (Ai, LDL_int) ; FREE_MEMORY (Ax, double) ; FREE_MEMORY (Flag, LDL_int) ; continue ; } /* ------------------------------------------------------------------ */ /* get the AMD permutation, if available */ /* ------------------------------------------------------------------ */ #ifdef USE_AMD /* recompute the permutation with AMD */ /* Assume that AMD produces a valid permutation P. */ #ifdef LDL_LONG if (amd_l_order (n, Ap, Ai, P, (double *) NULL, Info) < AMD_OK) { printf (PROGRAM ": call to AMD failed\n") ; EXIT_ERROR ; } amd_l_control ((double *) NULL) ; amd_l_info (Info) ; #else if (amd_order (n, Ap, Ai, P, (double *) NULL, Info) < AMD_OK) { printf (PROGRAM ": call to AMD failed\n") ; EXIT_ERROR ; } amd_control ((double *) NULL) ; amd_info (Info) ; #endif #endif /* ------------------------------------------------------------------ */ /* allocate workspace and the first part of LDL factorization */ /* ------------------------------------------------------------------ */ ALLOC_MEMORY (Pinv, LDL_int, n) ; ALLOC_MEMORY (Y, double, n) ; ALLOC_MEMORY (Pattern, LDL_int, n) ; ALLOC_MEMORY (Lnz, LDL_int, n) ; ALLOC_MEMORY (Lp, LDL_int, n+1) ; ALLOC_MEMORY (Parent, LDL_int, n) ; ALLOC_MEMORY (D, double, n) ; ALLOC_MEMORY (B, double, n) ; ALLOC_MEMORY (X, double, n) ; /* ------------------------------------------------------------------ */ /* factorize twice, with and without permutation */ /* ------------------------------------------------------------------ */ for (trial = 1 ; trial <= 2 ; trial++) { if (trial == 1) { printf ("Factorize PAP'=LDL' and solve Ax=b\n") ; Perm = P ; PermInv = Pinv ; } else { printf ("Factorize A=LDL' and solve Ax=b\n") ; Perm = (LDL_int *) NULL ; PermInv = (LDL_int *) NULL ; } /* -------------------------------------------------------------- */ /* symbolic factorization to get Lp, Parent, Lnz, and Pinv */ /* -------------------------------------------------------------- */ LDL_symbolic (n, Ap, Ai, Lp, Parent, Lnz, Flag, Perm, PermInv) ; lnz = Lp [n] ; /* find # of nonzeros in L, and flop count for LDL_numeric */ flops = 0 ; for (j = 0 ; j < n ; j++) { flops += ((double) Lnz [j]) * (Lnz [j] + 2) ; } printf ("Nz in L: "LDL_ID" Flop count: %g\n", lnz, flops) ; /* -------------------------------------------------------------- */ /* allocate remainder of L, of size lnz */ /* -------------------------------------------------------------- */ ALLOC_MEMORY (Li, LDL_int, lnz) ; ALLOC_MEMORY (Lx, double, lnz) ; /* -------------------------------------------------------------- */ /* numeric factorization to get Li, Lx, and D */ /* -------------------------------------------------------------- */ d = LDL_numeric (n, Ap, Ai, Ax, Lp, Parent, Lnz, Li, Lx, D, Y, Flag, Pattern, Perm, PermInv) ; /* -------------------------------------------------------------- */ /* solve, or report singular case */ /* -------------------------------------------------------------- */ if (d != n) { printf ("Ax=b not solved since D("LDL_ID","LDL_ID") is zero.\n", d, d) ; } else { /* construct the right-hand-side, B */ for (i = 0 ; i < n ; i++) { B [i] = 1 + ((double) i) / 100 ; } /* solve Ax=b */ if (trial == 1) { /* the factorization is LDL' = PAP' */ LDL_perm (n, Y, B, P) ; /* y = Pb */ LDL_lsolve (n, Y, Lp, Li, Lx) ; /* y = L\y */ LDL_dsolve (n, Y, D) ; /* y = D\y */ LDL_ltsolve (n, Y, Lp, Li, Lx) ; /* y = L'\y */ LDL_permt (n, X, Y, P) ; /* x = P'y */ } else { /* the factorization is LDL' = A */ for (i = 0 ; i < n ; i++) /* x = b */ { X [i] = B [i] ; } LDL_lsolve (n, X, Lp, Li, Lx) ; /* x = L\x */ LDL_dsolve (n, X, D) ; /* x = D\x */ LDL_ltsolve (n, X, Lp, Li, Lx) ; /* x = L'\x */ } /* compute the residual y = Ax-b */ /* note that this code can tolerate a jumbled matrix */ for (i = 0 ; i < n ; i++) { Y [i] = -B [i] ; } for (j = 0 ; j < n ; j++) { for (p = Ap [j] ; p < Ap [j+1] ; p++) { Y [Ai [p]] += Ax [p] * X [j] ; } } /* rnorm = norm (y, inf) */ rnorm = 0 ; for (i = 0 ; i < n ; i++) { r = (Y [i] > 0) ? (Y [i]) : (-Y [i]) ; rnorm = (r > rnorm) ? (r) : (rnorm) ; } maxrnorm = (rnorm > maxrnorm) ? (rnorm) : (maxrnorm) ; printf ("relative maxnorm of residual: %g\n", rnorm) ; } /* -------------------------------------------------------------- */ /* free the size-lnz part of L */ /* -------------------------------------------------------------- */ FREE_MEMORY (Li, LDL_int) ; FREE_MEMORY (Lx, double) ; } /* free everything */ FREE_MEMORY (P, LDL_int) ; FREE_MEMORY (Ap, LDL_int) ; FREE_MEMORY (Ai, LDL_int) ; FREE_MEMORY (Ax, double) ; FREE_MEMORY (Pinv, LDL_int) ; FREE_MEMORY (Y, double) ; FREE_MEMORY (Flag, LDL_int) ; FREE_MEMORY (Pattern, LDL_int) ; FREE_MEMORY (Lnz, LDL_int) ; FREE_MEMORY (Lp, LDL_int) ; FREE_MEMORY (Parent, LDL_int) ; FREE_MEMORY (D, double) ; FREE_MEMORY (B, double) ; FREE_MEMORY (X, double) ; } printf ("\nLargest residual during all tests: %g\n", maxrnorm) ; if (maxrnorm < 1e-8) { printf ("\n" PROGRAM ": all tests passed\n") ; EXIT_OK ; } else { printf ("\n" PROGRAM ": one more tests failed (residual too high)\n") ; EXIT_ERROR ; } }