1 /* ========================================================================== */
2 /* === ldlsimple.c: a simple LDL main program ============================== */
3 /* ========================================================================== */
4
5 /* LDLSIMPLE: this is a very simple main program that illustrates the basic
6 * usage of the LDL routines. The output of this program is in ldlsimple.out.
7 * This program factorizes the matrix
8 *
9 * A =[ ...
10 * 1.7 0 0 0 0 0 0 0 .13 0
11 * 0 1. 0 0 .02 0 0 0 0 .01
12 * 0 0 1.5 0 0 0 0 0 0 0
13 * 0 0 0 1.1 0 0 0 0 0 0
14 * 0 .02 0 0 2.6 0 .16 .09 .52 .53
15 * 0 0 0 0 0 1.2 0 0 0 0
16 * 0 0 0 0 .16 0 1.3 0 0 .56
17 * 0 0 0 0 .09 0 0 1.6 .11 0
18 * .13 0 0 0 .52 0 0 .11 1.4 0
19 * 0 .01 0 0 .53 0 .56 0 0 3.1 ] ;
20 *
21 * and then solves a system Ax=b whose true solution is
22 * x = [.1 .2 .3 .4 .5 .6 .7 .8 .9 1]' ;
23 *
24 * Note that Li and Lx are statically allocated, with length 13. This is just
25 * enough to hold the factor L, but normally this size is not known until after
26 * ldl_symbolic has analyzed the matrix. The size of Li and Lx must be greater
27 * than or equal to lnz = Lp [N], which is 13 for this matrix L.
28 *
29 * Copyright (c) 2006 by Timothy A Davis, http://www.suitesparse.com.
30 * All Rights Reserved. See LDL/Doc/License.txt for the License.
31 */
32
33 #include <stdio.h>
34 #include "ldl.h"
35 #define N 10 /* A is 10-by-10 */
36 #define ANZ 19 /* # of nonzeros on diagonal and upper triangular part of A */
37 #define LNZ 13 /* # of nonzeros below the diagonal of L */
38
main(void)39 int main (void)
40 {
41 /* only the upper triangular part of A is required */
42 int Ap [N+1] = {0, 1, 2, 3, 4, 6, 7, 9, 11, 15, ANZ},
43 Ai [ANZ] = {0, 1, 2, 3, 1,4, 5, 4,6, 4,7, 0,4,7,8, 1,4,6,9 } ;
44 double Ax [ANZ] = {1.7, 1., 1.5, 1.1, .02,2.6, 1.2, .16,1.3, .09,1.6,
45 .13,.52,.11,1.4, .01,.53,.56,3.1},
46 b [N] = {.287, .22, .45, .44, 2.486, .72, 1.55, 1.424, 1.621, 3.759};
47 double Lx [LNZ], D [N], Y [N] ;
48 int Li [LNZ], Lp [N+1], Parent [N], Lnz [N], Flag [N], Pattern [N], d, i ;
49
50 /* factorize A into LDL' (P and Pinv not used) */
51 ldl_symbolic (N, Ap, Ai, Lp, Parent, Lnz, Flag, NULL, NULL) ;
52 printf ("Nonzeros in L, excluding diagonal: %d\n", Lp [N]) ;
53 d = ldl_numeric (N, Ap, Ai, Ax, Lp, Parent, Lnz, Li, Lx, D, Y, Pattern,
54 Flag, NULL, NULL) ;
55
56 if (d == N)
57 {
58 /* solve Ax=b, overwriting b with the solution x */
59 ldl_lsolve (N, b, Lp, Li, Lx) ;
60 ldl_dsolve (N, b, D) ;
61 ldl_ltsolve (N, b, Lp, Li, Lx) ;
62 for (i = 0 ; i < N ; i++) printf ("x [%d] = %g\n", i, b [i]) ;
63 }
64 else
65 {
66 printf ("ldl_numeric failed, D (%d,%d) is zero\n", d, d) ;
67 }
68 return (0) ;
69 }
70