1
/* ========================================================================== */
2
/* === ldlmain.c: LDL main program, for demo and testing ==================== */
3
/* ========================================================================== */
5
/* LDLMAIN: this main program has two purposes. It provides an example of how
6
* to use the LDL routines, and it tests the package. The output of this
7
* program is in ldlmain.out (without AMD) and ldlamd.out (with AMD). If you
8
* did not download and install AMD, then the compilation of this program will
9
* not work with USE_AMD defined. Compare your output with ldlmain.out and
12
* The program reads in a set of matrices, in the Matrix/ subdirectory.
13
* The format of the files is as follows:
15
* one line with the matrix description
16
* one line with 2 integers: n jumbled
17
* n+1 lines, containing the Ap array of size n+1 (column pointers)
18
* nz lines, containing the Ai array of size nz (row indices)
19
* nz lines, containing the Ax array of size nz (numerical values)
20
* n lines, containing the P permutation array of size n
22
* The Ap, Ai, Ax, and P data structures are described in ldl.c.
23
* The jumbled flag is 1 if the matrix could contain unsorted columns and/or
24
* duplicate entries, and 0 otherwise.
26
* Once the matrix is read in, it is checked to see if it is valid. Some
27
* matrices are invalid by design, to test the error-checking routines. If
28
* valid, the matrix factorized twice (A and P*A*P'). A linear
29
* system Ax=b is set up and solved, and the residual computed.
30
* If any system is not solved accurately, this test will fail.
32
* This program can also be compiled as a MATLAB mexFunction, with the command
33
* "mex ldlmain.c ldl.c". You can then run the program in MATLAB, with the
36
* LDL Copyright (c) by Timothy A Davis,
37
* University of Florida. All Rights Reserved. See README for the License.
40
#ifdef MATLAB_MEX_FILE
50
#define NMATRICES 30 /* number of test matrices in Matrix/ directory */
51
#define LEN 200 /* string length */
53
#ifdef USE_AMD /* get AMD include file, if using AMD */
55
#define PROGRAM "ldlamd"
57
#define PROGRAM "ldlmain"
60
#ifdef MATLAB_MEX_FILE
62
#define EXIT_ERROR mexErrMsgTxt ("failure") ;
65
#define EXIT_ERROR exit (EXIT_FAILURE) ;
66
#define EXIT_OK exit (EXIT_SUCCESS) ;
69
/* -------------------------------------------------------------------------- */
70
/* ALLOC_MEMORY: allocate a block of memory */
71
/* -------------------------------------------------------------------------- */
73
#define ALLOC_MEMORY(p,type,size) \
74
p = (type *) malloc ((((size) <= 0) ? 1 : (size)) * sizeof (type)) ; \
75
if (p == (type *) NULL) \
77
printf (PROGRAM ": out of memory\n") ; \
81
/* -------------------------------------------------------------------------- */
82
/* FREE_MEMORY: free a block of memory */
83
/* -------------------------------------------------------------------------- */
85
#define FREE_MEMORY(p,type) \
86
if (p != (type *) NULL) \
92
/* -------------------------------------------------------------------------- */
93
/* stand-alone main program, or MATLAB mexFunction */
94
/* -------------------------------------------------------------------------- */
96
#ifdef MATLAB_MEX_FILE
102
const mxArray *pargin[ ]
109
/* ---------------------------------------------------------------------- */
110
/* local variables */
111
/* ---------------------------------------------------------------------- */
114
double Info [AMD_INFO] ;
116
double r, rnorm, flops, maxrnorm = 0. ;
117
double *Ax, *Lx, *B, *D, *X, *Y ;
118
LDL_int matrix, *Ai, *Ap, *Li, *Lp, *P, *Pinv, *Perm, *PermInv, n, i, j, p,
119
nz, *Flag, *Pattern, *Lnz, *Parent, trial, lnz, d, jumbled ;
123
/* ---------------------------------------------------------------------- */
124
/* check the error-checking routines with null matrices */
125
/* ---------------------------------------------------------------------- */
129
if (LDL_valid_perm (n, (LDL_int *) NULL, &i)
130
|| !LDL_valid_perm (0, (LDL_int *) NULL, &i)
131
|| LDL_valid_matrix (n, (LDL_int *) NULL, (LDL_int *) NULL)
132
|| LDL_valid_matrix (0, &i, &i))
134
printf (PROGRAM ": ldl error-checking routine failed\n") ;
138
/* ---------------------------------------------------------------------- */
139
/* read in a factorize a set of matrices */
140
/* ---------------------------------------------------------------------- */
142
for (matrix = 1 ; matrix <= NMATRICES ; matrix++)
145
/* ------------------------------------------------------------------ */
146
/* read in the matrix and the permutation */
147
/* ------------------------------------------------------------------ */
149
sprintf (s, "../Matrix/A%02d", (int) matrix) ;
150
if ((f = fopen (s, "r")) == (FILE *) NULL)
152
printf (PROGRAM ": could not open file: %s\n", s) ;
156
printf ("\n\n--------------------------------------------------------");
157
printf ("\nInput matrix: %s", s) ;
158
printf ("--------------------------------------------------------\n\n");
159
fscanf (f, LDL_ID " " LDL_ID, &n, &jumbled) ;
160
n = (n < 0) ? (0) : (n) ;
161
ALLOC_MEMORY (P, LDL_int, n) ;
162
ALLOC_MEMORY (Ap, LDL_int, n+1) ;
163
for (j = 0 ; j <= n ; j++)
165
fscanf (f, LDL_ID, &Ap [j]) ;
168
ALLOC_MEMORY (Ai, LDL_int, nz) ;
169
ALLOC_MEMORY (Ax, double, nz) ;
170
for (p = 0 ; p < nz ; p++)
172
fscanf (f, LDL_ID , &Ai [p]) ;
174
for (p = 0 ; p < nz ; p++)
176
fscanf (f, "%lg", &Ax [p]) ;
178
for (j = 0 ; j < n ; j++)
180
fscanf (f, LDL_ID , &P [j]) ;
184
/* ------------------------------------------------------------------ */
185
/* check the matrix A and the permutation P */
186
/* ------------------------------------------------------------------ */
188
ALLOC_MEMORY (Flag, LDL_int, n) ;
190
/* To test the error-checking routines, some of the input matrices
191
* are not valid. So this error is expected to occur. */
192
if (!LDL_valid_matrix (n, Ap, Ai) || !LDL_valid_perm (n, P, Flag))
194
printf (PROGRAM ": invalid matrix and/or permutation\n") ;
195
FREE_MEMORY (P, LDL_int) ;
196
FREE_MEMORY (Ap, LDL_int) ;
197
FREE_MEMORY (Ai, LDL_int) ;
198
FREE_MEMORY (Ax, double) ;
199
FREE_MEMORY (Flag, LDL_int) ;
203
/* ------------------------------------------------------------------ */
204
/* get the AMD permutation, if available */
205
/* ------------------------------------------------------------------ */
209
/* recompute the permutation with AMD */
210
/* Assume that AMD produces a valid permutation P. */
214
if (amd_l_order (n, Ap, Ai, P, (double *) NULL, Info) < AMD_OK)
216
printf (PROGRAM ": call to AMD failed\n") ;
219
amd_l_control ((double *) NULL) ;
224
if (amd_order (n, Ap, Ai, P, (double *) NULL, Info) < AMD_OK)
226
printf (PROGRAM ": call to AMD failed\n") ;
229
amd_control ((double *) NULL) ;
235
/* ------------------------------------------------------------------ */
236
/* allocate workspace and the first part of LDL factorization */
237
/* ------------------------------------------------------------------ */
239
ALLOC_MEMORY (Pinv, LDL_int, n) ;
240
ALLOC_MEMORY (Y, double, n) ;
241
ALLOC_MEMORY (Pattern, LDL_int, n) ;
242
ALLOC_MEMORY (Lnz, LDL_int, n) ;
243
ALLOC_MEMORY (Lp, LDL_int, n+1) ;
244
ALLOC_MEMORY (Parent, LDL_int, n) ;
245
ALLOC_MEMORY (D, double, n) ;
246
ALLOC_MEMORY (B, double, n) ;
247
ALLOC_MEMORY (X, double, n) ;
249
/* ------------------------------------------------------------------ */
250
/* factorize twice, with and without permutation */
251
/* ------------------------------------------------------------------ */
253
for (trial = 1 ; trial <= 2 ; trial++)
258
printf ("Factorize PAP'=LDL' and solve Ax=b\n") ;
264
printf ("Factorize A=LDL' and solve Ax=b\n") ;
265
Perm = (LDL_int *) NULL ;
266
PermInv = (LDL_int *) NULL ;
269
/* -------------------------------------------------------------- */
270
/* symbolic factorization to get Lp, Parent, Lnz, and Pinv */
271
/* -------------------------------------------------------------- */
273
LDL_symbolic (n, Ap, Ai, Lp, Parent, Lnz, Flag, Perm, PermInv) ;
276
/* find # of nonzeros in L, and flop count for LDL_numeric */
278
for (j = 0 ; j < n ; j++)
280
flops += ((double) Lnz [j]) * (Lnz [j] + 2) ;
282
printf ("Nz in L: "LDL_ID" Flop count: %g\n", lnz, flops) ;
284
/* -------------------------------------------------------------- */
285
/* allocate remainder of L, of size lnz */
286
/* -------------------------------------------------------------- */
288
ALLOC_MEMORY (Li, LDL_int, lnz) ;
289
ALLOC_MEMORY (Lx, double, lnz) ;
291
/* -------------------------------------------------------------- */
292
/* numeric factorization to get Li, Lx, and D */
293
/* -------------------------------------------------------------- */
295
d = LDL_numeric (n, Ap, Ai, Ax, Lp, Parent, Lnz, Li, Lx, D,
296
Y, Flag, Pattern, Perm, PermInv) ;
298
/* -------------------------------------------------------------- */
299
/* solve, or report singular case */
300
/* -------------------------------------------------------------- */
304
printf ("Ax=b not solved since D("LDL_ID","LDL_ID") is zero.\n", d, d) ;
308
/* construct the right-hand-side, B */
309
for (i = 0 ; i < n ; i++)
311
B [i] = 1 + ((double) i) / 100 ;
317
/* the factorization is LDL' = PAP' */
318
LDL_perm (n, Y, B, P) ; /* y = Pb */
319
LDL_lsolve (n, Y, Lp, Li, Lx) ; /* y = L\y */
320
LDL_dsolve (n, Y, D) ; /* y = D\y */
321
LDL_ltsolve (n, Y, Lp, Li, Lx) ; /* y = L'\y */
322
LDL_permt (n, X, Y, P) ; /* x = P'y */
326
/* the factorization is LDL' = A */
327
for (i = 0 ; i < n ; i++) /* x = b */
331
LDL_lsolve (n, X, Lp, Li, Lx) ; /* x = L\x */
332
LDL_dsolve (n, X, D) ; /* x = D\x */
333
LDL_ltsolve (n, X, Lp, Li, Lx) ; /* x = L'\x */
336
/* compute the residual y = Ax-b */
337
/* note that this code can tolerate a jumbled matrix */
338
for (i = 0 ; i < n ; i++)
342
for (j = 0 ; j < n ; j++)
344
for (p = Ap [j] ; p < Ap [j+1] ; p++)
346
Y [Ai [p]] += Ax [p] * X [j] ;
349
/* rnorm = norm (y, inf) */
351
for (i = 0 ; i < n ; i++)
353
r = (Y [i] > 0) ? (Y [i]) : (-Y [i]) ;
354
rnorm = (r > rnorm) ? (r) : (rnorm) ;
356
maxrnorm = (rnorm > maxrnorm) ? (rnorm) : (maxrnorm) ;
357
printf ("relative maxnorm of residual: %g\n", rnorm) ;
360
/* -------------------------------------------------------------- */
361
/* free the size-lnz part of L */
362
/* -------------------------------------------------------------- */
364
FREE_MEMORY (Li, LDL_int) ;
365
FREE_MEMORY (Lx, double) ;
369
/* free everything */
370
FREE_MEMORY (P, LDL_int) ;
371
FREE_MEMORY (Ap, LDL_int) ;
372
FREE_MEMORY (Ai, LDL_int) ;
373
FREE_MEMORY (Ax, double) ;
374
FREE_MEMORY (Pinv, LDL_int) ;
375
FREE_MEMORY (Y, double) ;
376
FREE_MEMORY (Flag, LDL_int) ;
377
FREE_MEMORY (Pattern, LDL_int) ;
378
FREE_MEMORY (Lnz, LDL_int) ;
379
FREE_MEMORY (Lp, LDL_int) ;
380
FREE_MEMORY (Parent, LDL_int) ;
381
FREE_MEMORY (D, double) ;
382
FREE_MEMORY (B, double) ;
383
FREE_MEMORY (X, double) ;
386
printf ("\nLargest residual during all tests: %g\n", maxrnorm) ;
389
printf ("\n" PROGRAM ": all tests passed\n") ;
394
printf ("\n" PROGRAM ": one more tests failed (residual too high)\n") ;
added
removed
LDL/Demo/ldlmain.c
