1 /* ========================================================================== */
2 /* === ldlmain.c: LDL main program, for demo and testing ==================== */
3 /* ========================================================================== */
4
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
10 * ldlamd.out.
11 *
12 * The program reads in a set of matrices, in the Matrix/ subdirectory.
13 * The format of the files is as follows:
14 *
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
21 *
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.
25 *
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.
31 *
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
34 * command "ldlmain".
35 *
36 * Copyright (c) 2006 by Timothy A Davis, http://www.suitesparse.com.
37 * All Rights Reserved. See LDL/Doc/License.txt for the License.
38 */
39
40 #ifdef MATLAB_MEX_FILE
41 #ifndef LDL_LONG
42 #define LDL_LONG
43 #endif
44 #endif
45
46 #include <stdio.h>
47 #include <stdlib.h>
48 #include "ldl.h"
49
50 #define NMATRICES 30 /* number of test matrices in Matrix/ directory */
51 #define LEN 200 /* string length */
52
53 #ifdef USE_AMD /* get AMD include file, if using AMD */
54 #include "amd.h"
55 #define PROGRAM "ldlamd"
56 #else
57 #define PROGRAM "ldlmain"
58 #endif
59
60 #ifdef MATLAB_MEX_FILE
61 #include "mex.h"
62 #define EXIT_ERROR mexErrMsgTxt ("failure") ;
63 #define EXIT_OK
64 #else
65 #define EXIT_ERROR exit (EXIT_FAILURE) ;
66 #define EXIT_OK exit (EXIT_SUCCESS) ;
67 #endif
68
69 /* -------------------------------------------------------------------------- */
70 /* ALLOC_MEMORY: allocate a block of memory */
71 /* -------------------------------------------------------------------------- */
72
73 #define ALLOC_MEMORY(p,type,size) \
74 p = (type *) calloc ((((size) <= 0) ? 1 : (size)) , sizeof (type)) ; \
75 if (p == (type *) NULL) \
76 { \
77 printf (PROGRAM ": out of memory\n") ; \
78 EXIT_ERROR ; \
79 }
80
81 /* -------------------------------------------------------------------------- */
82 /* FREE_MEMORY: free a block of memory */
83 /* -------------------------------------------------------------------------- */
84
85 #define FREE_MEMORY(p,type) \
86 if (p != (type *) NULL) \
87 { \
88 free (p) ; \
89 p = (type *) NULL ; \
90 }
91
92 /* -------------------------------------------------------------------------- */
93 /* stand-alone main program, or MATLAB mexFunction */
94 /* -------------------------------------------------------------------------- */
95
96 #ifdef MATLAB_MEX_FILE
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])97 void mexFunction
98 (
99 int nargout,
100 mxArray *pargout[ ],
101 int nargin,
102 const mxArray *pargin[ ]
103 )
104 #else
105 int main (void)
106 #endif
107 {
108
109 /* ---------------------------------------------------------------------- */
110 /* local variables */
111 /* ---------------------------------------------------------------------- */
112
113 #ifdef USE_AMD
114 double Info [AMD_INFO] ;
115 #endif
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, ok ;
120 FILE *f ;
121 char s [LEN], filename [LEN] ;
122
123 /* ---------------------------------------------------------------------- */
124 /* check the error-checking routines with null matrices */
125 /* ---------------------------------------------------------------------- */
126
127 i = 1 ;
128 n = -1 ;
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))
133 {
134 printf (PROGRAM ": ldl error-checking routine failed\n") ;
135 EXIT_ERROR ;
136 }
137
138 /* ---------------------------------------------------------------------- */
139 /* read in a factorize a set of matrices */
140 /* ---------------------------------------------------------------------- */
141
142 for (matrix = 1 ; matrix <= NMATRICES ; matrix++)
143 {
144
145 /* ------------------------------------------------------------------ */
146 /* read in the matrix and the permutation */
147 /* ------------------------------------------------------------------ */
148
149 sprintf (filename, "../Matrix/A%02d", (int) matrix) ;
150 if ((f = fopen (filename, "r")) == (FILE *) NULL)
151 {
152 printf (PROGRAM ": could not open file: %s\n", filename) ;
153 EXIT_ERROR ;
154 }
155 printf ("\n\n--------------------------------------------------------");
156 printf ("\nInput file: %s\n", filename) ;
157 s [0] = 0 ;
158 ok = (fgets (s, LEN, f) != NULL) ;
159 printf ("%s", s) ;
160 printf ("--------------------------------------------------------\n\n");
161 n = 0 ;
162 if (ok)
163 {
164 ok = ok && (fscanf (f, LDL_ID " " LDL_ID, &n, &jumbled) == 2) ;
165 n = (n < 0) ? (0) : (n) ;
166 ALLOC_MEMORY (P, LDL_int, n) ;
167 ALLOC_MEMORY (Ap, LDL_int, n+1) ;
168 }
169 for (j = 0 ; ok && j <= n ; j++)
170 {
171 ok = ok && (fscanf (f, LDL_ID, &Ap [j]) == 1) ;
172 }
173 nz = 0 ;
174 if (ok)
175 {
176 nz = Ap [n] ;
177 ALLOC_MEMORY (Ai, LDL_int, nz) ;
178 ALLOC_MEMORY (Ax, double, nz) ;
179 }
180 for (p = 0 ; ok && p < nz ; p++)
181 {
182 ok = ok && (fscanf (f, LDL_ID , &Ai [p]) == 1) ;
183 }
184 for (p = 0 ; ok && p < nz ; p++)
185 {
186 ok = ok && (fscanf (f, "%lg", &Ax [p]) == 1) ;
187 }
188 for (j = 0 ; ok && j < n ; j++)
189 {
190 ok = ok && (fscanf (f, LDL_ID , &P [j]) == 1) ;
191 }
192 fclose (f) ;
193 if (!ok)
194 {
195 printf (PROGRAM ": error reading file: %s\n", filename) ;
196 EXIT_ERROR ;
197 }
198
199 /* ------------------------------------------------------------------ */
200 /* check the matrix A and the permutation P */
201 /* ------------------------------------------------------------------ */
202
203 ALLOC_MEMORY (Flag, LDL_int, n) ;
204
205 /* To test the error-checking routines, some of the input matrices
206 * are not valid. So this error is expected to occur. */
207 if (!LDL_valid_matrix (n, Ap, Ai) || !LDL_valid_perm (n, P, Flag))
208 {
209 printf (PROGRAM ": invalid matrix and/or permutation\n") ;
210 FREE_MEMORY (P, LDL_int) ;
211 FREE_MEMORY (Ap, LDL_int) ;
212 FREE_MEMORY (Ai, LDL_int) ;
213 FREE_MEMORY (Ax, double) ;
214 FREE_MEMORY (Flag, LDL_int) ;
215 continue ;
216 }
217
218 /* ------------------------------------------------------------------ */
219 /* get the AMD permutation, if available */
220 /* ------------------------------------------------------------------ */
221
222 #ifdef USE_AMD
223
224 /* recompute the permutation with AMD */
225 /* Assume that AMD produces a valid permutation P. */
226
227 #ifdef LDL_LONG
228
229 if (amd_l_order (n, Ap, Ai, P, (double *) NULL, Info) < AMD_OK)
230 {
231 printf (PROGRAM ": call to AMD failed\n") ;
232 EXIT_ERROR ;
233 }
234 amd_l_control ((double *) NULL) ;
235 amd_l_info (Info) ;
236
237 #else
238
239 if (amd_order (n, Ap, Ai, P, (double *) NULL, Info) < AMD_OK)
240 {
241 printf (PROGRAM ": call to AMD failed\n") ;
242 EXIT_ERROR ;
243 }
244 amd_control ((double *) NULL) ;
245 amd_info (Info) ;
246
247 #endif
248 #endif
249
250 /* ------------------------------------------------------------------ */
251 /* allocate workspace and the first part of LDL factorization */
252 /* ------------------------------------------------------------------ */
253
254 ALLOC_MEMORY (Pinv, LDL_int, n) ;
255 ALLOC_MEMORY (Y, double, n) ;
256 ALLOC_MEMORY (Pattern, LDL_int, n) ;
257 ALLOC_MEMORY (Lnz, LDL_int, n) ;
258 ALLOC_MEMORY (Lp, LDL_int, n+1) ;
259 ALLOC_MEMORY (Parent, LDL_int, n) ;
260 ALLOC_MEMORY (D, double, n) ;
261 ALLOC_MEMORY (B, double, n) ;
262 ALLOC_MEMORY (X, double, n) ;
263
264 /* ------------------------------------------------------------------ */
265 /* factorize twice, with and without permutation */
266 /* ------------------------------------------------------------------ */
267
268 for (trial = 1 ; trial <= 2 ; trial++)
269 {
270
271 if (trial == 1)
272 {
273 printf ("Factorize PAP'=LDL' and solve Ax=b\n") ;
274 Perm = P ;
275 PermInv = Pinv ;
276 }
277 else
278 {
279 printf ("Factorize A=LDL' and solve Ax=b\n") ;
280 Perm = (LDL_int *) NULL ;
281 PermInv = (LDL_int *) NULL ;
282 }
283
284 /* -------------------------------------------------------------- */
285 /* symbolic factorization to get Lp, Parent, Lnz, and Pinv */
286 /* -------------------------------------------------------------- */
287
288 LDL_symbolic (n, Ap, Ai, Lp, Parent, Lnz, Flag, Perm, PermInv) ;
289 lnz = Lp [n] ;
290
291 /* find # of nonzeros in L, and flop count for LDL_numeric */
292 flops = 0 ;
293 for (j = 0 ; j < n ; j++)
294 {
295 flops += ((double) Lnz [j]) * (Lnz [j] + 2) ;
296 }
297 printf ("Nz in L: "LDL_ID" Flop count: %g\n", lnz, flops) ;
298
299 /* -------------------------------------------------------------- */
300 /* allocate remainder of L, of size lnz */
301 /* -------------------------------------------------------------- */
302
303 ALLOC_MEMORY (Li, LDL_int, lnz) ;
304 ALLOC_MEMORY (Lx, double, lnz) ;
305
306 /* -------------------------------------------------------------- */
307 /* numeric factorization to get Li, Lx, and D */
308 /* -------------------------------------------------------------- */
309
310 d = LDL_numeric (n, Ap, Ai, Ax, Lp, Parent, Lnz, Li, Lx, D,
311 Y, Flag, Pattern, Perm, PermInv) ;
312
313 /* -------------------------------------------------------------- */
314 /* solve, or report singular case */
315 /* -------------------------------------------------------------- */
316
317 if (d != n)
318 {
319 printf ("Ax=b not solved since D("LDL_ID","LDL_ID") is zero.\n", d, d) ;
320 }
321 else
322 {
323 /* construct the right-hand-side, B */
324 for (i = 0 ; i < n ; i++)
325 {
326 B [i] = 1 + ((double) i) / 100 ;
327 }
328
329 /* solve Ax=b */
330 if (trial == 1)
331 {
332 /* the factorization is LDL' = PAP' */
333 LDL_perm (n, Y, B, P) ; /* y = Pb */
334 LDL_lsolve (n, Y, Lp, Li, Lx) ; /* y = L\y */
335 LDL_dsolve (n, Y, D) ; /* y = D\y */
336 LDL_ltsolve (n, Y, Lp, Li, Lx) ; /* y = L'\y */
337 LDL_permt (n, X, Y, P) ; /* x = P'y */
338 }
339 else
340 {
341 /* the factorization is LDL' = A */
342 for (i = 0 ; i < n ; i++) /* x = b */
343 {
344 X [i] = B [i] ;
345 }
346 LDL_lsolve (n, X, Lp, Li, Lx) ; /* x = L\x */
347 LDL_dsolve (n, X, D) ; /* x = D\x */
348 LDL_ltsolve (n, X, Lp, Li, Lx) ; /* x = L'\x */
349 }
350
351 /* compute the residual y = Ax-b */
352 /* note that this code can tolerate a jumbled matrix */
353 for (i = 0 ; i < n ; i++)
354 {
355 Y [i] = -B [i] ;
356 }
357 for (j = 0 ; j < n ; j++)
358 {
359 for (p = Ap [j] ; p < Ap [j+1] ; p++)
360 {
361 Y [Ai [p]] += Ax [p] * X [j] ;
362 }
363 }
364 /* rnorm = norm (y, inf) */
365 rnorm = 0 ;
366 for (i = 0 ; i < n ; i++)
367 {
368 r = (Y [i] > 0) ? (Y [i]) : (-Y [i]) ;
369 rnorm = (r > rnorm) ? (r) : (rnorm) ;
370 }
371 maxrnorm = (rnorm > maxrnorm) ? (rnorm) : (maxrnorm) ;
372 printf ("relative maxnorm of residual: %g\n", rnorm) ;
373 }
374
375 /* -------------------------------------------------------------- */
376 /* free the size-lnz part of L */
377 /* -------------------------------------------------------------- */
378
379 FREE_MEMORY (Li, LDL_int) ;
380 FREE_MEMORY (Lx, double) ;
381
382 }
383
384 /* free everything */
385 FREE_MEMORY (P, LDL_int) ;
386 FREE_MEMORY (Ap, LDL_int) ;
387 FREE_MEMORY (Ai, LDL_int) ;
388 FREE_MEMORY (Ax, double) ;
389 FREE_MEMORY (Pinv, LDL_int) ;
390 FREE_MEMORY (Y, double) ;
391 FREE_MEMORY (Flag, LDL_int) ;
392 FREE_MEMORY (Pattern, LDL_int) ;
393 FREE_MEMORY (Lnz, LDL_int) ;
394 FREE_MEMORY (Lp, LDL_int) ;
395 FREE_MEMORY (Parent, LDL_int) ;
396 FREE_MEMORY (D, double) ;
397 FREE_MEMORY (B, double) ;
398 FREE_MEMORY (X, double) ;
399 }
400
401 printf ("\nLargest residual during all tests: %g\n", maxrnorm) ;
402 if (maxrnorm < 1e-8)
403 {
404 printf ("\n" PROGRAM ": all tests passed\n") ;
405 EXIT_OK ;
406 }
407 else
408 {
409 printf ("\n" PROGRAM ": one more tests failed (residual too high)\n") ;
410 EXIT_ERROR ;
411 }
412 }
413