1 /* ========================================================================== */
2 /* === Tcov/solve =========================================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Tcov Module. Copyright (C) 2005-2006, Timothy A. Davis
7 * http://www.suitesparse.com
8 * -------------------------------------------------------------------------- */
9
10 /* Test CHOLMOD for solving various systems of linear equations. */
11
12 #include "cm.h"
13
14 #define NFTYPES 17
15 Int ll_types [NFTYPES] = { 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0 } ;
16 Int pk_types [NFTYPES] = { 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0 } ;
17 Int mn_types [NFTYPES] = { 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0 } ;
18 Int co_types [NFTYPES] = { 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1 } ;
19
20 #define NRHS 9
21
22 #ifndef NDEBUG
23 #ifndef EXTERN
24 #define EXTERN extern
25 #endif
26 EXTERN int cholmod_dump, cholmod_l_dump ;
27 #endif
28
29 /* ========================================================================== */
30 /* === test_solver ========================================================== */
31 /* ========================================================================== */
32
33 /* Test solve(A) with various control parameters */
34
test_solver(cholmod_sparse * A)35 double test_solver (cholmod_sparse *A)
36 {
37 double err, maxerr = 0 ;
38 int save ;
39
40 for (cm->postorder = 0 ; cm->postorder <= 1 ; cm->postorder++)
41 {
42
43 my_srand (42) ; /* RAND reset */
44
45 /* simplicial, no extra memory */
46 printf ("test_solver: simplicial, no extra memory\n") ;
47 cm->supernodal = CHOLMOD_SIMPLICIAL ;
48 cm->grow2 = 0 ;
49 err = solve (A) ;
50 MAXERR (maxerr, err, 1) ;
51 printf ("test_solver err: %6.2e\n", err) ;
52
53 /* simplicial, extra space in columns of L */
54 printf ("test_solver: simplicial, extra space in columns of L\n") ;
55 cm->grow2 = 5 ;
56 err = solve (A) ;
57 MAXERR (maxerr, err, 1) ;
58 printf ("test_solver err: %6.2e\n", err) ;
59
60 /* supernodal */
61 printf ("test_solver: supernodal\n") ;
62 cm->supernodal = CHOLMOD_SUPERNODAL ;
63 err = solve (A) ;
64 MAXERR (maxerr, err, 1) ;
65 printf ("test_solver err: %6.2e\n", err) ;
66
67 /* supernodal, without final resymbol */
68 printf ("test_solver: supernodal, without final resymbol\n") ;
69 cm->final_resymbol = FALSE ;
70 err = solve (A) ;
71 MAXERR (maxerr, err, 1) ;
72 printf ("test_solver err: %6.2e\n", err) ;
73
74 /* supernodal, with resymbol, final_super false */
75 printf ("test_solver: supernodal, with resymbol\n") ;
76 cm->supernodal = CHOLMOD_SUPERNODAL ;
77 cm->final_asis = FALSE ;
78 cm->final_resymbol = TRUE ;
79 cm->final_super = FALSE ;
80 err = solve (A) ;
81 MAXERR (maxerr, err, 1) ;
82 printf ("test_solver err: %6.2e\n", err) ;
83
84 /* supernodal, with resymbol, final_super tree */
85 printf ("test_solver: supernodal, with resymbol\n") ;
86 cm->supernodal = CHOLMOD_SUPERNODAL ;
87 cm->final_asis = FALSE ;
88 cm->final_resymbol = TRUE ;
89 cm->final_super = TRUE ;
90 err = solve (A) ;
91 MAXERR (maxerr, err, 1) ;
92 printf ("test_solver err: %6.2e\n", err) ;
93
94 /* simplicial LL' */
95 printf ("test_solver: simplicial LL', try NESDIS instead of METIS\n") ;
96 cm->supernodal = CHOLMOD_SIMPLICIAL ;
97 cm->final_ll = TRUE ;
98 cm->default_nesdis = TRUE ;
99 err = solve (A) ;
100 cm->default_nesdis = FALSE ;
101 MAXERR (maxerr, err, 1) ;
102 cm->final_ll = FALSE ;
103 printf ("test_solver err: %6.2e\n", err) ;
104
105 /* supernodal, without final resymbol, and no relaxed supernodes */
106 printf (
107 "test_solver: supernodal, without final resymbol, and no relaxed supernodes\n") ;
108 cm->supernodal = CHOLMOD_SUPERNODAL ;
109 cm->final_asis = TRUE ;
110 cm->nrelax [0] = 0 ;
111 cm->nrelax [1] = 0 ;
112 cm->nrelax [2] = 0 ;
113 cm->zrelax [0] = 0 ;
114 cm->zrelax [1] = 0 ;
115 cm->zrelax [2] = 0 ;
116 cm->grow0 = 1 ;
117 cm->grow1 = 1 ;
118 cm->grow2 = 0 ;
119 err = solve (A) ;
120 MAXERR (maxerr, err, 1) ;
121 printf ("test_solver err: %6.2e\n", err) ;
122
123 /* ------------------------------------------------------------------ */
124 /* restore defaults */
125 /* ------------------------------------------------------------------ */
126
127 cm->dbound = 0.0 ;
128
129 cm->grow0 = 1.2 ;
130 cm->grow1 = 1.2 ;
131 cm->grow2 = 5 ;
132
133 cm->final_asis = TRUE ;
134 cm->final_super = TRUE ;
135 cm->final_ll = FALSE ;
136 cm->final_pack = TRUE ;
137 cm->final_monotonic = TRUE ;
138 cm->final_resymbol = FALSE ;
139
140 cm->supernodal = CHOLMOD_AUTO ;
141 cm->nrelax [0] = 4 ;
142 cm->nrelax [1] = 16 ;
143 cm->nrelax [2] = 48 ;
144 cm->zrelax [0] = 0.8 ;
145 cm->zrelax [1] = 0.1 ;
146 cm->zrelax [2] = 0.05 ;
147
148 /* do not restore these defaults: */
149 /*
150 cm->maxrank = ...
151 cm->metis_memory = 2.0
152 cm->metis_nswitch = 3000
153 cm->metis_dswitch = 0.66
154 cm->print = 3
155 cm->precise = FALSE
156 */
157 }
158
159 progress (1, '.') ;
160 return (maxerr) ;
161 }
162
163
164 /* ========================================================================== */
165 /* === solve ================================================================ */
166 /* ========================================================================== */
167
168 /* solve Ax=b or AA'x=b, systems involving just L, D, or L', and update/downdate
169 * the system. Returns the worst-case residual. This routine keeps going if
170 * it runs out of memory (unless the error handler terminates it), because
171 * it is used both normally and in the memory tests.
172 */
173
solve(cholmod_sparse * A)174 double solve (cholmod_sparse *A)
175 {
176 double r, enorm, snorm, maxerr = 0, gnorm, anorm, bnorm, xnorm, norm, rcond;
177 cholmod_factor *L, *Lcopy, *L2 ;
178 cholmod_sparse *Lmat, *Lo, *D, *Up, *S, *LD, *LDL, *E, *I, *C, *CC, *Ct,
179 *Ssym, *Cperm, *C2, *S2, *H, *F, *Lo1, *Aboth, *Bset ;
180 cholmod_dense *X, *Cdense, *B, *Bcomplex, *Bzomplex, *Breal, *W, *R,
181 *A3, *C3, *E3 ;
182 cholmod_dense *B2, *B2complex, *B2zomplex, *B2real, *Ywork, *Ework ;
183 cholmod_sparse *AFt, *AF, *G, *RowK, *Bsparse, *Xsparse ;
184 double *Cx ;
185 double *B2x ;
186 Int *P, *cset, *fset, *Parent, *Post, *RowCount, *ColCount,
187 *First, *Level, *rcount, *ccount, *Lp, *Li ;
188 Int p, i, j, k, n, nrhs, save, save2, csize, rank, nrow, ncol, is_ll,
189 xtype, isreal, prefer_zomplex, Lxtype, xtype2, save3 ;
190 cm->blas_ok = TRUE ;
191
192 Ywork = NULL ;
193 Ework = NULL ;
194
195 if (cm->print > 1)
196 {
197 printf ("============================================== in solve:\n") ;
198 }
199
200 if (A == NULL)
201 {
202 ERROR (CHOLMOD_INVALID, "nothing to solve") ;
203 return (1) ;
204 }
205
206 /* ---------------------------------------------------------------------- */
207 /* construct right-hand-side (Ax=b if symmetric, AA'x=b otherwise) */
208 /* ---------------------------------------------------------------------- */
209
210 n = A->nrow ;
211 nrow = A->nrow ;
212 ncol = A->ncol ;
213 xtype = A->xtype ;
214 isreal = (xtype == CHOLMOD_REAL) ;
215 B = rhs (A, NRHS, n) ;
216 anorm = CHOLMOD(norm_sparse) (A, 1, cm) ;
217 save = cm->final_asis ;
218 cm->final_asis = TRUE ;
219
220 /* contents of these will be revised later */
221 Bzomplex = CHOLMOD(copy_dense) (B, cm) ;
222 Bcomplex = CHOLMOD(copy_dense) (Bzomplex, cm) ;
223 Breal = CHOLMOD(copy_dense) (Bzomplex, cm) ;
224
225 /* create another right-hand-side, slightly different */
226 B2 = CHOLMOD(copy_dense) (B, cm) ;
227 if (B2 != NULL)
228 {
229 B2x = (double *) B2->x ;
230 B2x [0] = 42 ;
231 }
232
233 /* ---------------------------------------------------------------------- */
234 /* analyze */
235 /* ---------------------------------------------------------------------- */
236
237 if (n < 100 && cm->nmethods == 1 && cm->method[0].ordering == CHOLMOD_GIVEN)
238 {
239 Int *UserPerm = prand (nrow) ; /* RAND */
240 L = CHOLMOD(analyze_p) (A, UserPerm, NULL, 0, cm) ;
241 OK (CHOLMOD(print_common) ("with UserPerm", cm)) ;
242 CHOLMOD(free) (nrow, sizeof (Int), UserPerm, cm) ;
243 }
244 else
245 {
246 L = CHOLMOD(analyze) (A, cm) ;
247 }
248
249 /* test rowadd on a symbolic factor */
250 if (isreal)
251 {
252 RowK = CHOLMOD(spzeros) (n, 1, 0, CHOLMOD_REAL, cm) ;
253 Lcopy = CHOLMOD(copy_factor) (L, cm) ;
254 if (n > 0)
255 {
256 CHOLMOD(rowadd) (0, RowK, Lcopy, cm) ;
257 CHOLMOD(check_factor) (Lcopy, cm) ;
258 CHOLMOD(print_factor) (Lcopy, "Lcopy, now numeric", cm) ;
259 }
260 CHOLMOD(free_sparse) (&RowK, cm) ;
261 CHOLMOD(free_factor) (&Lcopy, cm) ;
262 }
263
264 /* ---------------------------------------------------------------------- */
265 /* factorize */
266 /* ---------------------------------------------------------------------- */
267
268 CHOLMOD(factorize) (A, L, cm) ;
269
270 /* ---------------------------------------------------------------------- */
271 /* various solves */
272 /* ---------------------------------------------------------------------- */
273
274 if (B != NULL)
275 {
276 /* B->ncol = 1 ; */
277 prefer_zomplex = (B->xtype == CHOLMOD_ZOMPLEX) ;
278 }
279 else
280 {
281 prefer_zomplex = FALSE ;
282 }
283
284 for (k = 0 ; k <= 1 ; k++)
285 {
286
287 cm->prefer_zomplex = k ;
288
289 /* compute the residual, X complex if L or B not real */
290 X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
291 r = resid (A, X, B) ;
292 MAXERR (maxerr, r, 1) ;
293 CHOLMOD(free_dense) (&X, cm) ;
294
295 /* now test the same thing, but with cholmod_solve2 */
296 CHOLMOD(solve2) (CHOLMOD_A, L, B, NULL, &X, NULL, &Ywork, &Ework, cm) ;
297 r = resid (A, X, B) ;
298 MAXERR (maxerr, r, 1) ;
299 CHOLMOD(solve2) (CHOLMOD_A, L, B2, NULL, &X, NULL, &Ywork, &Ework, cm) ;
300 r = resid (A, X, B2) ;
301 MAXERR (maxerr, r, 1) ;
302
303 CHOLMOD(free_dense) (&X, cm) ;
304 CHOLMOD(free_dense) (&Ywork, cm) ;
305 CHOLMOD(free_dense) (&Ework, cm) ;
306
307 /* see what happens if X has the wrong size or type on input. This
308 * will free X on input, and reallocate it the right size and type */
309 X = CHOLMOD(ones) (1, 1, CHOLMOD_REAL, cm) ;
310 CHOLMOD(solve2) (CHOLMOD_A, L, B, NULL, &X, NULL, &Ywork, &Ework, cm) ;
311 r = resid (A, X, B) ;
312 MAXERR (maxerr, r, 1) ;
313 CHOLMOD(free_dense) (&Ywork, cm) ;
314
315 /* now X is OK, and so is E, but Y is null */
316 CHOLMOD(solve2) (CHOLMOD_A, L, B, NULL, &X, NULL, &Ywork, &Ework, cm) ;
317 r = resid (A, X, B) ;
318 MAXERR (maxerr, r, 1) ;
319 CHOLMOD(free_dense) (&X, cm) ;
320 CHOLMOD(free_dense) (&Ywork, cm) ;
321 CHOLMOD(free_dense) (&Ework, cm) ;
322
323 /* zomplex right-hand-side */
324 CHOLMOD(dense_xtype) (CHOLMOD_ZOMPLEX, Bzomplex, cm) ;
325 if (Bzomplex != NULL && B != NULL && B->xtype == CHOLMOD_REAL
326 && Bzomplex->xtype == CHOLMOD_ZOMPLEX)
327 {
328 /* add an arbitrary imaginary part */
329 double *Bz = Bzomplex->z ;
330 for (j = 0 ; j < NRHS ; j++)
331 {
332 for (i = 0 ; i < n ; i++)
333 {
334 Bz [i+j*n] = (double) (i+j*n) ;
335 }
336 }
337 }
338 X = CHOLMOD(solve) (CHOLMOD_A, L, Bzomplex, cm) ;
339 r = resid (A, X, Bzomplex) ;
340 MAXERR (maxerr, r, 1) ;
341 CHOLMOD(free_dense) (&X, cm) ;
342
343 B2zomplex = CHOLMOD(copy_dense) (Bzomplex, cm) ;
344 if (B2zomplex != NULL)
345 {
346 B2x = (double *) B2zomplex->x ;
347 B2x [0] = 99 ;
348 }
349
350 /* now test the same thing, but with cholmod_solve2 */
351 CHOLMOD(solve2) (CHOLMOD_A, L, Bzomplex, NULL, &X, NULL,
352 &Ywork, &Ework, cm) ;
353 r = resid (A, X, Bzomplex) ;
354 MAXERR (maxerr, r, 1) ;
355 CHOLMOD(solve2) (CHOLMOD_A, L, B2zomplex, NULL, &X, NULL,
356 &Ywork, &Ework, cm) ;
357 r = resid (A, X, B2zomplex) ;
358 MAXERR (maxerr, r, 1) ;
359 CHOLMOD(free_dense) (&X, cm) ;
360 CHOLMOD(free_dense) (&Ywork, cm) ;
361 CHOLMOD(free_dense) (&Ework, cm) ;
362
363 /* complex right-hand-side */
364 CHOLMOD(dense_xtype) (CHOLMOD_COMPLEX, Bcomplex, cm) ;
365 X = CHOLMOD(solve) (CHOLMOD_A, L, Bcomplex, cm) ;
366 r = resid (A, X, Bcomplex) ;
367 MAXERR (maxerr, r, 1) ;
368 CHOLMOD(free_dense) (&X, cm) ;
369
370 B2complex = CHOLMOD(copy_dense) (Bcomplex, cm) ;
371 if (B2complex != NULL)
372 {
373 B2x = (double *) B2complex->x ;
374 B2x [0] = 777 ;
375 }
376
377 /* now test the same thing, but with cholmod_solve2 */
378 CHOLMOD(solve2) (CHOLMOD_A, L, Bcomplex, NULL, &X, NULL,
379 &Ywork, &Ework, cm) ;
380 r = resid (A, X, Bcomplex) ;
381 MAXERR (maxerr, r, 1) ;
382 CHOLMOD(solve2) (CHOLMOD_A, L, B2complex, NULL, &X, NULL,
383 &Ywork, &Ework, cm) ;
384 r = resid (A, X, B2complex) ;
385 MAXERR (maxerr, r, 1) ;
386 CHOLMOD(free_dense) (&X, cm) ;
387 CHOLMOD(free_dense) (&Ywork, cm) ;
388 CHOLMOD(free_dense) (&Ework, cm) ;
389
390 /* real right-hand-side */
391 CHOLMOD(dense_xtype) (CHOLMOD_REAL, Breal, cm) ;
392 X = CHOLMOD(solve) (CHOLMOD_A, L, Breal, cm) ;
393 r = resid (A, X, Breal) ;
394 MAXERR (maxerr, r, 1) ;
395 CHOLMOD(free_dense) (&X, cm) ;
396
397 B2real = CHOLMOD(copy_dense) (Breal, cm) ;
398 if (B2real != NULL)
399 {
400 B2x = (double *) B2real->x ;
401 B2x [0] = 1234 ;
402 }
403
404 /* now test the same thing, but with cholmod_solve2 */
405 CHOLMOD(solve2) (CHOLMOD_A, L, Breal, NULL, &X, NULL,
406 &Ywork, &Ework, cm) ;
407 r = resid (A, X, Breal) ;
408 MAXERR (maxerr, r, 1) ;
409 CHOLMOD(solve2) (CHOLMOD_A, L, B2real, NULL, &X, NULL,
410 &Ywork, &Ework, cm) ;
411 r = resid (A, X, B2real) ;
412 MAXERR (maxerr, r, 1) ;
413 CHOLMOD(free_dense) (&X, cm) ;
414 CHOLMOD(free_dense) (&Ywork, cm) ;
415 CHOLMOD(free_dense) (&Ework, cm) ;
416
417 /* sparse solve of Ax=b, b real */
418 Bsparse = CHOLMOD(dense_to_sparse) (Breal, TRUE, cm) ;
419 Xsparse = CHOLMOD(spsolve) (CHOLMOD_A, L, Bsparse, cm) ;
420 r = resid_sparse (A, Xsparse, Bsparse) ;
421 MAXERR (maxerr, r, 1) ;
422 CHOLMOD(free_sparse) (&Bsparse, cm) ;
423 CHOLMOD(free_sparse) (&Xsparse, cm) ;
424
425 /* sparse solve of Ax=b, b complex */
426 Bsparse = CHOLMOD(dense_to_sparse) (Bcomplex, TRUE, cm) ;
427 Xsparse = CHOLMOD(spsolve) (CHOLMOD_A, L, Bsparse, cm) ;
428 r = resid_sparse (A, Xsparse, Bsparse) ;
429 MAXERR (maxerr, r, 1) ;
430 CHOLMOD(free_sparse) (&Bsparse, cm) ;
431 CHOLMOD(free_sparse) (&Xsparse, cm) ;
432
433 /* sparse solve of Ax=b, b zomplex */
434 Bsparse = CHOLMOD(dense_to_sparse) (Bzomplex, TRUE, cm) ;
435 Xsparse = CHOLMOD(spsolve) (CHOLMOD_A, L, Bsparse, cm) ;
436 r = resid_sparse (A, Xsparse, Bsparse) ;
437 MAXERR (maxerr, r, 1) ;
438 CHOLMOD(free_sparse) (&Bsparse, cm) ;
439 CHOLMOD(free_sparse) (&Xsparse, cm) ;
440
441 CHOLMOD(free_dense) (&B2zomplex, cm) ;
442 CHOLMOD(free_dense) (&B2complex, cm) ;
443 CHOLMOD(free_dense) (&B2real, cm) ;
444 }
445
446 cm->prefer_zomplex = FALSE ;
447
448 /* ---------------------------------------------------------------------- */
449 /* sparse solve to compute inv(A) */
450 /* ---------------------------------------------------------------------- */
451
452 CHOLMOD(print_sparse) (A, "A", cm) ;
453 CHOLMOD(print_factor) (L, "L", cm) ;
454 rcond = CHOLMOD(rcond) (L, cm) ;
455 if (cm->print > 1)
456 {
457 printf ("rcond: %g\n", rcond) ;
458 }
459
460 if (n < 100 && A->stype != 0)
461 {
462 /* solve A*C=I, so C should equal A inverse */
463 I = CHOLMOD(speye) (n, n, CHOLMOD_REAL, cm) ;
464 C = CHOLMOD(spsolve) (CHOLMOD_A, L, I, cm) ;
465 /* compute norm of A*C-I */
466 if (isreal && n > 10)
467 {
468 /* A and C are large and real */
469 E = CHOLMOD(ssmult) (A, C, 0, TRUE, FALSE, cm) ;
470 F = CHOLMOD(add) (E, I, minusone, one, TRUE, FALSE, cm) ;
471 r = CHOLMOD(norm_sparse) (F, 1, cm) ;
472 CHOLMOD(free_sparse) (&E, cm) ;
473 CHOLMOD(free_sparse) (&F, cm) ;
474 }
475 else
476 {
477 /* There is no complex ssmult or add, so use the BLAS.
478 * Also test sparse_to_dense for small symmetric matrices. */
479 A3 = CHOLMOD(sparse_to_dense) (A, cm) ;
480 C3 = CHOLMOD(sparse_to_dense) (C, cm) ;
481 xtype2 = isreal ? CHOLMOD_REAL : CHOLMOD_COMPLEX ;
482 CHOLMOD(dense_xtype) (xtype2, A3, cm) ;
483 CHOLMOD(dense_xtype) (xtype2, C3, cm) ;
484 E3 = CHOLMOD(eye) (n, n, xtype2, cm) ;
485 if (A3 != NULL && C3 != NULL && E3 != NULL)
486 {
487 /* E3 = A3*C3-I */
488 if (isreal)
489 {
490 BLAS_dgemm ("N", "N", n, n, n, one, A3->x, n, C3->x, n,
491 minusone, E3->x, n) ;
492 }
493 else
494 {
495 BLAS_zgemm ("N", "N", n, n, n, one, A3->x, n, C3->x, n,
496 minusone, E3->x, n) ;
497 }
498 OK (cm->blas_ok) ;
499 }
500 r = CHOLMOD(norm_dense) (E3, 1, cm) ;
501 CHOLMOD(free_dense) (&A3, cm) ;
502 CHOLMOD(free_dense) (&C3, cm) ;
503 CHOLMOD(free_dense) (&E3, cm) ;
504 }
505 MAXERR (maxerr, r, 1) ;
506 CHOLMOD(free_sparse) (&I, cm) ;
507 CHOLMOD(free_sparse) (&C, cm) ;
508 }
509
510 /* ---------------------------------------------------------------------- */
511 /* change complexity of L and solve again; test copy/change routines */
512 /* ---------------------------------------------------------------------- */
513
514 /* change to complex, otherwise leave as */
515 Lcopy = CHOLMOD(copy_factor) (L, cm) ;
516 CHOLMOD(factor_xtype) (CHOLMOD_COMPLEX, Lcopy, cm) ;
517 X = CHOLMOD(solve) (CHOLMOD_A, Lcopy, B, cm) ;
518 r = resid (A, X, B) ;
519 MAXERR (maxerr, r, 1) ;
520 CHOLMOD(free_dense) (&X, cm) ;
521 CHOLMOD(free_factor) (&Lcopy, cm) ;
522
523 /* change to zomplex LDL' */
524 Lxtype = (L == NULL) ? CHOLMOD_REAL : (L->xtype) ;
525 Lcopy = CHOLMOD(copy_factor) (L, cm) ;
526 CHOLMOD(check_factor) (L, cm) ;
527 CHOLMOD(print_factor) (Lcopy, "Lcopy", cm) ;
528 CHOLMOD(change_factor) (Lxtype, FALSE, FALSE, TRUE, TRUE, Lcopy, cm) ;
529 CHOLMOD(factor_xtype) (CHOLMOD_ZOMPLEX, Lcopy, cm) ;
530 X = CHOLMOD(solve) (CHOLMOD_A, Lcopy, B, cm) ;
531 r = resid (A, X, B) ;
532 MAXERR (maxerr, r, 1) ;
533 CHOLMOD(free_dense) (&X, cm) ;
534 CHOLMOD(free_factor) (&Lcopy, cm) ;
535
536 Lcopy = CHOLMOD(copy_factor) (L, cm) ;
537 CHOLMOD(change_factor) (Lxtype, TRUE, FALSE, FALSE, FALSE, Lcopy, cm) ;
538 CHOLMOD(check_factor) (L, cm) ;
539 CHOLMOD(print_factor) (Lcopy, "Lcopy LL unpacked", cm) ;
540 CHOLMOD(free_factor) (&Lcopy, cm) ;
541
542 CHOLMOD(free_factor) (&L, cm) ;
543 cm->final_asis = save ;
544
545 /* ---------------------------------------------------------------------- */
546 /* solve again, but use cm->final_asis as given */
547 /* ---------------------------------------------------------------------- */
548
549 if (n < 100 && cm->nmethods == 1 && cm->method[0].ordering == CHOLMOD_GIVEN)
550 {
551 Int *UserPerm = prand (nrow) ; /* RAND */
552 L = CHOLMOD(analyze_p) (A, UserPerm, NULL, 0, cm) ;
553 CHOLMOD(free) (nrow, sizeof (Int), UserPerm, cm) ;
554 }
555 else
556 {
557 L = CHOLMOD(analyze) (A, cm) ;
558 }
559
560 CHOLMOD(print_factor) (L, "Lsymbolic", cm) ;
561 CHOLMOD(factorize) (A, L, cm) ;
562
563 /* turn off memory tests [ */
564 save3 = my_tries ;
565 my_tries = -1 ;
566
567 CHOLMOD(print_factor) (L, "Lnumeric for solver tests", cm) ;
568 CHOLMOD(print_dense) (B, "B for solver tests", cm) ;
569 CHOLMOD(print_dense) (Breal, "Breal for solver tests", cm) ;
570 CHOLMOD(print_dense) (Bcomplex, "Bcomplex for solver tests", cm) ;
571 CHOLMOD(print_dense) (Bzomplex, "Bzomplex for solver tests", cm) ;
572
573 if (B != NULL && Breal != NULL && Bcomplex != NULL && Bzomplex != NULL)
574 {
575 for (nrhs = 1 ; nrhs <= NRHS ; nrhs++)
576 {
577 for (cm->prefer_zomplex = 0 ; cm->prefer_zomplex <= 1 ;
578 cm->prefer_zomplex++)
579 {
580 B->ncol = nrhs ;
581 Breal->ncol = nrhs ;
582 Bcomplex->ncol = nrhs ;
583 Bzomplex->ncol = nrhs ;
584
585 X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
586 r = resid (A, X, B) ;
587 CHOLMOD(free_dense) (&X, cm) ;
588 MAXERR (maxerr, r, 1) ;
589
590 X = CHOLMOD(solve) (CHOLMOD_A, L, Breal, cm) ;
591 r = resid (A, X, Breal) ;
592 CHOLMOD(free_dense) (&X, cm) ;
593 MAXERR (maxerr, r, 1) ;
594
595 X = CHOLMOD(solve) (CHOLMOD_A, L, Bcomplex, cm) ;
596 r = resid (A, X, Bcomplex) ;
597 CHOLMOD(free_dense) (&X, cm) ;
598 MAXERR (maxerr, r, 1) ;
599
600 X = CHOLMOD(solve) (CHOLMOD_A, L, Bzomplex, cm) ;
601 r = resid (A, X, Bzomplex) ;
602 CHOLMOD(free_dense) (&X, cm) ;
603 MAXERR (maxerr, r, 1) ;
604 }
605 }
606 }
607 cm->prefer_zomplex = FALSE ;
608
609 /* turn memory tests back on, where we left off ] */
610 my_tries = save3 ;
611
612 Lcopy = CHOLMOD(copy_factor) (L, cm) ;
613
614 /* ---------------------------------------------------------------------- */
615 /* convert L to LDL' packed or LL packed */
616 /* ---------------------------------------------------------------------- */
617
618 printf ("before change factor : %d\n", L ? L->is_super : -1) ;
619 is_ll = (L == NULL) ? FALSE : (L->is_ll) ;
620 Lxtype = (L == NULL) ? CHOLMOD_REAL : (L->xtype) ;
621 CHOLMOD(change_factor) (Lxtype, is_ll, FALSE, TRUE, TRUE, Lcopy, cm) ;
622 printf ("after change factor : %d\n", L ? L->is_super : -1) ;
623
624 /* ---------------------------------------------------------------------- */
625 /* extract L, D, and L' as matrices from Lcopy */
626 /* ---------------------------------------------------------------------- */
627
628 CHOLMOD(resymbol) (A, NULL, 0, TRUE, Lcopy, cm) ;
629
630 Lmat = CHOLMOD(factor_to_sparse) (Lcopy, cm) ;
631 CHOLMOD(check_sparse) (Lmat, cm) ;
632
633 I = CHOLMOD(speye) (n, n, CHOLMOD_REAL, cm) ;
634
635 Lo = NULL ;
636 D = NULL ;
637
638 Lxtype = (Lmat == NULL) ? CHOLMOD_REAL : (Lmat->xtype) ;
639
640 if (isreal)
641 {
642 /* use band and add */
643 if (!is_ll)
644 {
645 /* factorization is LDL' = Lo*D*Up */
646 Lo1 = CHOLMOD(band) (Lmat, -n, -1, 1, cm) ;
647 Lo = CHOLMOD(add) (Lo1, I, one, one, TRUE, TRUE, cm) ;
648 CHOLMOD(free_sparse) (&Lo1, cm) ;
649 D = CHOLMOD(band) (Lmat, 0, 0, 1, cm) ;
650 }
651 else
652 {
653 /* factorization is LL' = Lo*D*Up */
654 Lo = CHOLMOD(band) (Lmat, -n, 0, 1, cm) ;
655 D = CHOLMOD(speye) (n, n, Lxtype, cm) ;
656 }
657 }
658 else
659 {
660 /* band and add do not work for c/zomplex matrices*/
661 D = CHOLMOD(speye) (n, n, Lxtype, cm) ;
662 Lo = CHOLMOD(copy_sparse) (Lmat, cm) ;
663 if (!is_ll && D != NULL && Lo != NULL)
664 {
665 /* factorization is LDL' = Lo*D*Up */
666 double *Dx = D->x ;
667 double *Lx = Lo->x ;
668 Lp = Lo->p ;
669 for (k = 0 ; k < n ; k++)
670 {
671 p = Lp [k] ;
672 if (Lxtype == CHOLMOD_COMPLEX)
673 {
674 Dx [2*k] = Lx [2*p] ;
675 Lx [2*p] = 1 ;
676 }
677 else
678 {
679 Dx [k] = Lx [p] ;
680 Lx [p] = 1 ;
681 }
682 }
683 }
684 }
685
686 Up = CHOLMOD(transpose) (Lo, 2, cm) ;
687
688 /* ---------------------------------------------------------------------- */
689 /* compute 1-norm of (Lo*D*Up - PAP') or (Lo*D*Up - PAA'P') */
690 /* ---------------------------------------------------------------------- */
691
692 P = (L != NULL) ? (L->Perm) : NULL ;
693 S = NULL ;
694 G = NULL ;
695
696 if (isreal)
697 {
698
699 if (A->stype == 0)
700 {
701 /* G = A*A', try with fset = prand (ncol) */
702 fset = prand (ncol) ; /* RAND */
703 AFt = CHOLMOD(ptranspose) (A, 1, NULL, fset, ncol, cm) ;
704 AF = CHOLMOD(transpose) (AFt, 1, cm) ;
705 CHOLMOD(free) (ncol, sizeof (Int), fset, cm) ;
706 G = CHOLMOD(ssmult) (AF, AFt, 0, TRUE, TRUE, cm) ;
707
708 /* also try aat */
709 H = CHOLMOD(aat) (AF, NULL, 0, 1, cm) ;
710 E = CHOLMOD(add) (G, H, one, minusone, TRUE, FALSE, cm) ;
711 enorm = CHOLMOD(norm_sparse) (E, 0, cm) ;
712 gnorm = CHOLMOD(norm_sparse) (G, 0, cm) ;
713 MAXERR (maxerr, enorm, gnorm) ;
714 if (cm->print > 1)
715 {
716 printf ("enorm %g gnorm %g hnorm %g\n", enorm, gnorm,
717 CHOLMOD(norm_sparse) (H, 0, cm)) ;
718 }
719 if (gnorm > 0)
720 {
721 enorm /= gnorm ;
722 }
723 OK (enorm < 1e-8) ;
724 CHOLMOD(free_sparse) (&AFt, cm) ;
725 CHOLMOD(free_sparse) (&AF, cm) ;
726 CHOLMOD(free_sparse) (&E, cm) ;
727 CHOLMOD(free_sparse) (&H, cm) ;
728 }
729 else
730 {
731 /* G = A */
732 G = CHOLMOD(copy) (A, 0, 1, cm) ;
733 }
734
735 if (A->stype == 0)
736 {
737 /* S = PAA'P' */
738 S = CHOLMOD(submatrix) (G, P, n, P, n, TRUE, FALSE, cm) ;
739 }
740 else
741 {
742 /* S = PAP' */
743 Aboth = CHOLMOD(copy) (A, 0, 1, cm) ;
744 S = CHOLMOD(submatrix) (Aboth, P, n, P, n, TRUE, FALSE, cm) ;
745 CHOLMOD(free_sparse) (&Aboth, cm) ;
746 }
747
748 if (n < NSMALL)
749 {
750 /* only do this for small test matrices, since L*D*L' can have many
751 * nonzero entries */
752
753 /* E = L*D*L' - S */
754 LD = CHOLMOD(ssmult) (Lo, D, 0, TRUE, FALSE, cm) ;
755 LDL = CHOLMOD(ssmult) (LD, Up, 0, TRUE, FALSE, cm) ;
756 CHOLMOD(free_sparse) (&LD, cm) ;
757 E = CHOLMOD(add) (LDL, S, one, minusone, TRUE, FALSE, cm) ;
758 CHOLMOD(free_sparse) (&LDL, cm) ;
759
760 /* e = norm (E) / norm (S) */
761 enorm = CHOLMOD(norm_sparse) (E, 1, cm) ;
762 snorm = CHOLMOD(norm_sparse) (S, 0, cm) ;
763 MAXERR (maxerr, enorm, snorm) ;
764 CHOLMOD(free_sparse) (&E, cm) ;
765 }
766
767 /* check the row/col counts */
768 RowCount = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
769 ColCount = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
770 Parent = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
771 Post = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
772 First = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
773 Level = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
774 rcount = CHOLMOD(calloc) (n, sizeof (Int), cm) ;
775 ccount = CHOLMOD(calloc) (n, sizeof (Int), cm) ;
776
777 if (S != NULL && Lmat != NULL && RowCount != NULL && ColCount != NULL &&
778 Parent != NULL && Post != NULL && First != NULL && Level != NULL &&
779 rcount != NULL && ccount != NULL)
780 {
781 S->stype = 1 ;
782
783 CHOLMOD(etree) (S, Parent, cm) ;
784 CHOLMOD(print_parent) (Parent, n, "Parent", cm) ;
785 CHOLMOD(postorder) (Parent, n, NULL, Post, cm) ;
786 CHOLMOD(print_perm) (Post, n, n, "Post", cm) ;
787 S->stype = -1 ;
788 CHOLMOD(rowcolcounts) (S, NULL, 0, Parent, Post, RowCount, ColCount,
789 First, Level, cm) ;
790
791 Lp = Lmat->p ;
792 Li = Lmat->i ;
793 OK (Lmat->packed) ;
794 for (j = 0 ; j < n ; j++)
795 {
796 for (p = Lp [j] ; p < Lp [j+1] ; p++)
797 {
798 i = Li [p] ;
799 rcount [i]++ ;
800 ccount [j]++ ;
801 }
802 }
803 /* a singular matrix will only be partially factorized */
804 if (L->minor == (size_t) n)
805 {
806 for (j = 0 ; j < n ; j++)
807 {
808 OK (ccount [j] == ColCount [j]) ;
809 }
810 }
811 for (i = 0 ; i < (Int) (L->minor) ; i++)
812 {
813 OK (rcount [i] == RowCount [i]) ;
814 }
815 }
816
817 CHOLMOD(free) (n, sizeof (Int), RowCount, cm) ;
818 CHOLMOD(free) (n, sizeof (Int), ColCount, cm) ;
819 CHOLMOD(free) (n, sizeof (Int), Parent, cm) ;
820 CHOLMOD(free) (n, sizeof (Int), Post, cm) ;
821 CHOLMOD(free) (n, sizeof (Int), First, cm) ;
822 CHOLMOD(free) (n, sizeof (Int), Level, cm) ;
823 CHOLMOD(free) (n, sizeof (Int), rcount, cm) ;
824 CHOLMOD(free) (n, sizeof (Int), ccount, cm) ;
825
826 CHOLMOD(free_sparse) (&S, cm) ;
827 }
828
829 CHOLMOD(free_factor) (&Lcopy, cm) ;
830
831 /* ---------------------------------------------------------------------- */
832 /* solve other systems */
833 /* ---------------------------------------------------------------------- */
834
835 /* turn off memory tests [ */
836 save3 = my_tries ;
837 my_tries = -1 ;
838
839 for (nrhs = 1 ; nrhs <= 4 ; nrhs++) /* reduced here (6 to 4) */
840 {
841
842 if (B == NULL)
843 {
844 break ;
845 }
846
847 B->ncol = nrhs ;
848
849 /* solve Ax=b and compare with the sparse solve */
850 X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
851 r = resid (A, X, B) ;
852 MAXERR (maxerr, r, 1) ;
853 CHOLMOD(free_dense) (&X, cm) ;
854
855 /* solve LDL'x=b and compare with sparse solve */
856 X = CHOLMOD(solve) (CHOLMOD_LDLt, L, B, cm) ;
857 /* printf ("LDL'x=b %p %p %p\n", Lo, X, B) ; */
858 r = resid3 (Lo, D, Up, X, B) ;
859 MAXERR (maxerr, r, 1) ;
860 CHOLMOD(free_dense) (&X, cm) ;
861
862 /* solve LDx=b and compare with sparse solve */
863 X = CHOLMOD(solve) (CHOLMOD_LD, L, B, cm) ;
864 /* printf ("LDx=b %p %p %p\n", Lo, X, B) ; */
865 r = resid3 (Lo, D, NULL, X, B) ;
866 MAXERR (maxerr, r, 1) ;
867 CHOLMOD(free_dense) (&X, cm) ;
868
869 /* solve DL'x=b and compare with sparse solve */
870 X = CHOLMOD(solve) (CHOLMOD_DLt, L, B, cm) ;
871 /* printf ("DL'x=b %p %p %p\n", D, X, B) ; */
872 r = resid3 (D, Up, NULL, X, B) ;
873 MAXERR (maxerr, r, 1) ;
874 CHOLMOD(free_dense) (&X, cm) ;
875
876 /* solve Lx=b and compare with sparse solve */
877 X = CHOLMOD(solve) (CHOLMOD_L, L, B, cm) ;
878 /* printf ("Lx=b %p %p %p\n", Lo, X, B) ; */
879 r = resid3 (Lo, NULL, NULL, X, B) ;
880 MAXERR (maxerr, r, 1) ;
881 CHOLMOD(free_dense) (&X, cm) ;
882
883 /* solve L'x=b and compare with sparse solve */
884 X = CHOLMOD(solve) (CHOLMOD_Lt, L, B, cm) ;
885 /* printf ("L'x=b %p %p %p\n", Up, X, B) ; */
886 r = resid3 (Up, NULL, NULL, X, B) ;
887 MAXERR (maxerr, r, 1) ;
888 CHOLMOD(free_dense) (&X, cm) ;
889
890 /* solve Dx=b and compare with sparse solve */
891 X = CHOLMOD(solve) (CHOLMOD_D, L, B, cm) ;
892 /* printf ("Dx=b %p %p %p\n", D, X, B) ; */
893 r = resid3 (D, NULL, NULL, X, B) ;
894 MAXERR (maxerr, r, 1) ;
895 CHOLMOD(free_dense) (&X, cm) ;
896
897 save2 = cm->prefer_zomplex ;
898 for (k = 0 ; k <= 1 ; k++)
899 {
900 cm->prefer_zomplex = k ;
901
902 /* x=Pb and compare with sparse solve */
903 X = CHOLMOD(solve) (CHOLMOD_P, L, B, cm) ;
904 r = pnorm (X, P, B, FALSE) ;
905 MAXERR (maxerr, r, 1) ;
906 CHOLMOD(free_dense) (&X, cm) ;
907
908 /* x=P'b and compare with sparse solve */
909 X = CHOLMOD(solve) (CHOLMOD_Pt, L, B, cm) ;
910 r = pnorm (X, P, B, TRUE) ;
911 MAXERR (maxerr, r, 1) ;
912 CHOLMOD(free_dense) (&X, cm) ;
913 }
914 cm->prefer_zomplex = save2 ;
915 }
916
917 /* turn memory tests back on, where we left off ] */
918 my_tries = save3 ;
919
920 CHOLMOD(free_dense) (&B, cm) ;
921 CHOLMOD(free_dense) (&B2, cm) ;
922
923 /* ---------------------------------------------------------------------- */
924 /* test the sparse solve */
925 /* ---------------------------------------------------------------------- */
926
927 /* turn off memory tests [ */
928 save3 = my_tries ;
929 my_tries = -1 ;
930
931 /* get ready for sparse solve */
932 Bset = CHOLMOD(allocate_sparse) (n, 1, 1, FALSE, TRUE, 0,
933 CHOLMOD_PATTERN, cm) ;
934 B = CHOLMOD(zeros) (n, 1, Lxtype, cm) ;
935
936 printf ("testing sparse solve, maxerr so far %g\n", maxerr) ;
937
938 if (Bset != NULL && B != NULL && n > 0)
939 {
940 cholmod_sparse *Xset ;
941 cholmod_dense *X2 ;
942 Int *Bseti, *Bsetp, *Xseti ;
943 double *X2x, *X1x, *Bx, *Bz ;
944
945 Xset = NULL ;
946 X2 = NULL ;
947
948 Bseti = Bset->i ;
949 Bsetp = Bset->p ;
950 Bsetp [0] = 0 ;
951 Bsetp [1] = 1 ;
952 Bseti [0] = 0 ; /* first entry in b is nonzero, all others zero */
953 Bx = B->x ;
954 Bz = B->z ;
955
956 if (B->xtype == CHOLMOD_REAL)
957 {
958 Bx [0] = 42 ;
959 }
960 else if (B->xtype == CHOLMOD_COMPLEX)
961 {
962 Bx [0] = -2 ;
963 Bx [1] = 1 ;
964 }
965 else
966 {
967 Bx [0] = 77 ;
968 Bz [0] = -2 ;
969 }
970
971 /* solve Ax=b and compare with the sparse solve */
972 X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
973 if (X != NULL)
974 {
975 r = resid (A, X, B) ;
976 MAXERR (maxerr, r, 1) ;
977 CHOLMOD(solve2) (CHOLMOD_A, L, B, Bset, &X2, &Xset,
978 &Ywork, &Ework, cm) ;
979 if (X2 != NULL)
980 {
981 X2x = X2->x ;
982 X1x = X->x ;
983 r = fabs (X2x [0] - X1x [0]) ;
984 printf ("solve2 Ax=b err %g\n", r) ;
985 MAXERR (maxerr, r, 1) ;
986 }
987 }
988 CHOLMOD(free_dense) (&X, cm) ;
989
990 /* solve LDL'x=b and compare with sparse solve */
991 X = CHOLMOD(solve) (CHOLMOD_LDLt, L, B, cm) ;
992 if (X != NULL)
993 {
994 /* printf ("LDL'x=b %p %p %p\n", Lo, X, B) ; */
995 r = resid3 (Lo, D, Up, X, B) ;
996 MAXERR (maxerr, r, 1) ;
997 CHOLMOD(solve2) (CHOLMOD_LDLt, L, B, Bset, &X2, &Xset,
998 &Ywork, &Ework, cm) ;
999 if (X2 != NULL)
1000 {
1001 X2x = X2->x ;
1002 X1x = X->x ;
1003 r = fabs (X2x [0] - X1x [0]) ;
1004 printf ("solve2 LDL'x=b err %g\n", r) ;
1005 MAXERR (maxerr, r, 1) ;
1006 }
1007 }
1008 CHOLMOD(free_dense) (&X, cm) ;
1009
1010 /* solve LDx=b and compare with sparse solve */
1011 X = CHOLMOD(solve) (CHOLMOD_LD, L, B, cm) ;
1012 if (X != NULL)
1013 {
1014 /* printf ("LDx=b %p %p %p\n", Lo, X, B) ; */
1015 r = resid3 (Lo, D, NULL, X, B) ;
1016 MAXERR (maxerr, r, 1) ;
1017 CHOLMOD(solve2) (CHOLMOD_LD, L, B, Bset, &X2, &Xset,
1018 &Ywork, &Ework, cm) ;
1019 if (X2 != NULL)
1020 {
1021 X2x = X2->x ;
1022 X1x = X->x ;
1023 r = fabs (X2x [0] - X1x [0]) ;
1024 printf ("solve2 LDx=b err %g\n", r) ;
1025 MAXERR (maxerr, r, 1) ;
1026 }
1027 }
1028 CHOLMOD(free_dense) (&X, cm) ;
1029
1030 /* solve DL'x=b and compare with sparse solve */
1031 X = CHOLMOD(solve) (CHOLMOD_DLt, L, B, cm) ;
1032 if (X != NULL)
1033 {
1034 /* printf ("DL'x=b %p %p %p\n", D, X, B) ; */
1035 r = resid3 (D, Up, NULL, X, B) ;
1036 MAXERR (maxerr, r, 1) ;
1037 CHOLMOD(solve2) (CHOLMOD_DLt, L, B, Bset, &X2, &Xset,
1038 &Ywork, &Ework, cm) ;
1039 if (X2 != NULL)
1040 {
1041 X2x = X2->x ;
1042 X1x = X->x ;
1043 r = fabs (X2x [0] - X1x [0]) ;
1044 printf ("solve2 DL'x=b err %g\n", r) ;
1045 MAXERR (maxerr, r, 1) ;
1046 }
1047 }
1048 CHOLMOD(free_dense) (&X, cm) ;
1049
1050 /* solve Lx=b and compare with sparse solve */
1051 X = CHOLMOD(solve) (CHOLMOD_L, L, B, cm) ;
1052 if (X != NULL)
1053 {
1054 /* printf ("Lx=b %p %p %p\n", Lo, X, B) ; */
1055 r = resid3 (Lo, NULL, NULL, X, B) ;
1056 MAXERR (maxerr, r, 1) ;
1057 CHOLMOD(solve2) (CHOLMOD_L, L, B, Bset, &X2, &Xset,
1058 &Ywork, &Ework, cm) ;
1059 if (X2 != NULL)
1060 {
1061 X2x = X2->x ;
1062 X1x = X->x ;
1063 r = fabs (X2x [0] - X1x [0]) ;
1064 printf ("solve2 Lx=b err %g\n", r) ;
1065 MAXERR (maxerr, r, 1) ;
1066 }
1067 }
1068 CHOLMOD(free_dense) (&X, cm) ;
1069
1070 /* solve L'x=b and compare with sparse solve */
1071 X = CHOLMOD(solve) (CHOLMOD_Lt, L, B, cm) ;
1072 if (X != NULL)
1073 {
1074 /* printf ("L'x=b %p %p %p\n", Up, X, B) ; */
1075 r = resid3 (Up, NULL, NULL, X, B) ;
1076 MAXERR (maxerr, r, 1) ;
1077 CHOLMOD(solve2) (CHOLMOD_Lt, L, B, Bset, &X2, &Xset,
1078 &Ywork, &Ework, cm) ;
1079 if (X2 != NULL)
1080 {
1081 X2x = X2->x ;
1082 X1x = X->x ;
1083 r = fabs (X2x [0] - X1x [0]) ;
1084 printf ("solve2 Ltx=b err %g\n", r) ;
1085 MAXERR (maxerr, r, 1) ;
1086 }
1087 }
1088 CHOLMOD(free_dense) (&X, cm) ;
1089
1090 /* solve Dx=b and compare with sparse solve */
1091 X = CHOLMOD(solve) (CHOLMOD_D, L, B, cm) ;
1092 if (X != NULL)
1093 {
1094 /* printf ("Dx=b %p %p %p\n", D, X, B) ; */
1095 r = resid3 (D, NULL, NULL, X, B) ;
1096 MAXERR (maxerr, r, 1) ;
1097 CHOLMOD(solve2) (CHOLMOD_D, L, B, Bset, &X2, &Xset,
1098 &Ywork, &Ework, cm) ;
1099 if (X2 != NULL)
1100 {
1101 X2x = X2->x ;
1102 X1x = X->x ;
1103 r = fabs (X2x [0] - X1x [0]) ;
1104 printf ("solve2 Dx=b err %g\n", r) ;
1105 MAXERR (maxerr, r, 1) ;
1106 }
1107 }
1108 CHOLMOD(free_dense) (&X, cm) ;
1109
1110 #if 0
1111 save2 = cm->prefer_zomplex ;
1112 for (k = 0 ; k <= 1 ; k++)
1113 {
1114 cm->prefer_zomplex = k ;
1115
1116 /* x=Pb and compare with sparse solve */
1117 X = CHOLMOD(solve) (CHOLMOD_P, L, B, cm) ;
1118 if (X != NULL)
1119 {
1120 r = pnorm (X, P, B, FALSE) ;
1121 MAXERR (maxerr, r, 1) ;
1122 CHOLMOD(solve2) (CHOLMOD_P, L, B, Bset, &X2, &Xset,
1123 &Ywork, &Ework, cm) ;
1124 if (X2 != NULL && Xset != NULL)
1125 {
1126 X2x = X2->x ;
1127 X1x = X->x ;
1128 Xseti = Xset->i ;
1129 i = Xseti [0] ;
1130 r = fabs (X2x [i] - X1x [i]) ;
1131 printf ("solve2 Px=b err %g\n", r) ;
1132 MAXERR (maxerr, r, 1) ;
1133 }
1134 }
1135 CHOLMOD(free_dense) (&X, cm) ;
1136
1137 /* x=P'b and compare with sparse solve */
1138 X = CHOLMOD(solve) (CHOLMOD_Pt, L, B, cm) ;
1139 if (X != NULL)
1140 {
1141 r = pnorm (X, P, B, TRUE) ;
1142 MAXERR (maxerr, r, 1) ;
1143 CHOLMOD(solve2) (CHOLMOD_Pt, L, B, Bset, &X2, &Xset,
1144 &Ywork, &Ework, cm) ;
1145 if (X2 != NULL && Xset != NULL)
1146 {
1147 X2x = X2->x ;
1148 X1x = X->x ;
1149 Xseti = Xset->i ;
1150 i = Xseti [0] ;
1151 r = fabs (X2x [i] - X1x [i]) ;
1152 printf ("solve2 P'x=b err %g\n", r) ;
1153 MAXERR (maxerr, r, 1) ;
1154 }
1155 }
1156 CHOLMOD(free_dense) (&X, cm) ;
1157 }
1158 cm->prefer_zomplex = save2 ;
1159 #endif
1160
1161 CHOLMOD(free_sparse) (&Xset, cm) ;
1162 CHOLMOD(free_dense) (&X2, cm) ;
1163 CHOLMOD(free_dense) (&Ywork, cm) ;
1164 CHOLMOD(free_dense) (&Ework, cm) ;
1165 }
1166
1167 CHOLMOD(free_sparse) (&Bset, cm) ;
1168 CHOLMOD(free_dense) (&B, cm) ;
1169
1170 /* turn memory tests back on, where we left off ] */
1171 my_tries = save3 ;
1172
1173 CHOLMOD(free_sparse) (&I, cm) ;
1174 CHOLMOD(free_sparse) (&D, cm) ;
1175 CHOLMOD(free_sparse) (&Lo, cm) ;
1176 CHOLMOD(free_sparse) (&Up, cm) ;
1177 CHOLMOD(free_sparse) (&Lmat, cm) ;
1178
1179 printf ("done testing sparse solve, maxerr so far %g\n", maxerr) ;
1180
1181 /* ---------------------------------------------------------------------- */
1182 /* update the factorization */
1183 /* ---------------------------------------------------------------------- */
1184
1185 /* turn off memory tests [ */
1186 save3 = my_tries ;
1187 my_tries = -1 ;
1188
1189 B = rhs (A, 1, n) ;
1190
1191 for (rank = 1 ; isreal && rank <= ((n < 100) ? 33 : 2) ; rank++)
1192 {
1193
1194 /* pick a random C */
1195 Cdense = CHOLMOD(zeros) (n, rank, CHOLMOD_REAL, cm) ;
1196
1197 if (Cdense != NULL)
1198 {
1199 Cx = Cdense->x ;
1200 for (k = 0 ; k < 10*rank ; k++)
1201 {
1202 i = nrand (n) ; /* RAND */
1203 j = nrand (rank) ; /* RAND */
1204 Cx [i+j*n] += xrand (1.) ; /* RAND */
1205 }
1206 }
1207
1208 C = CHOLMOD(dense_to_sparse) (Cdense, TRUE, cm) ;
1209 CHOLMOD(free_dense) (&Cdense, cm) ;
1210
1211 /* permute the rows according to L->Perm */
1212 Cperm = CHOLMOD(submatrix) (C, P, n, NULL, -1, TRUE, TRUE, cm) ;
1213
1214 /* update */
1215 CHOLMOD(updown) (TRUE, Cperm, L, cm) ;
1216 CHOLMOD(free_sparse) (&Cperm, cm) ;
1217
1218 /* solve (G+C*C')x=b */
1219 X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
1220
1221 /* an alternative method would be to compute (G*x + C*(C'*x) - b) */
1222
1223 /* compute S = G+C*C', with no sort */
1224 Ct = CHOLMOD(transpose) (C, 1, cm) ;
1225 CC = CHOLMOD(ssmult) (C, Ct, 0, TRUE, FALSE, cm) ;
1226 S = CHOLMOD(add) (G, CC, one, one, TRUE, TRUE, cm) ;
1227 Ssym = CHOLMOD(copy) (S, 1, 1, cm) ;
1228 CHOLMOD(free_sparse) (&CC, cm) ;
1229 CHOLMOD(free_sparse) (&Ct, cm) ;
1230
1231 /* compute norm (S*x-b) */
1232 r = resid (Ssym, X, B) ;
1233 MAXERR (maxerr, r, 1) ;
1234 CHOLMOD(free_dense) (&X, cm) ;
1235
1236 /* ------------------------------------------------------------------ */
1237 /* factor A+CC' from scratch, using same permutation */
1238 /* ------------------------------------------------------------------ */
1239
1240 if (rank == 1)
1241 {
1242 save = cm->nmethods ;
1243 save2 = cm->method [0].ordering ;
1244 cm->nmethods = 1 ;
1245 cm->method [0].ordering = CHOLMOD_GIVEN ;
1246 L2 = CHOLMOD(analyze_p) (Ssym, P, NULL, 0, cm) ;
1247 cm->nmethods = save ;
1248 cm->method [0].ordering = save2 ;
1249 CHOLMOD(factorize) (Ssym, L2, cm) ;
1250 X = CHOLMOD(solve) (CHOLMOD_A, L2, B, cm) ;
1251 r = resid (Ssym, X, B) ;
1252 MAXERR (maxerr, r, 1) ;
1253 CHOLMOD(free_dense) (&X, cm) ;
1254 CHOLMOD(free_factor) (&L2, cm) ;
1255 }
1256
1257 CHOLMOD(free_sparse) (&Ssym, cm) ;
1258
1259 /* ------------------------------------------------------------------ */
1260 /* downdate, with just the first half of C */
1261 /* ------------------------------------------------------------------ */
1262
1263 csize = MAX (1, rank / 2) ;
1264 cset = CHOLMOD(malloc) (csize, sizeof (Int), cm) ;
1265 if (cset != NULL)
1266 {
1267 for (i = 0 ; i < csize ; i++)
1268 {
1269 cset [i] = i ;
1270 }
1271 }
1272 C2 = CHOLMOD(submatrix) (C, NULL, -1, cset, csize, TRUE, TRUE, cm) ;
1273 Cperm = CHOLMOD(submatrix) (C2, P, n, NULL, -1, TRUE, TRUE, cm) ;
1274 CHOLMOD(free) (csize, sizeof (Int), cset, cm) ;
1275
1276 CHOLMOD(updown) (FALSE, Cperm, L, cm) ;
1277 CHOLMOD(free_sparse) (&Cperm, cm) ;
1278
1279 /* solve (G+C*C'-C2*C2')x=b */
1280 X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
1281
1282 /* This is an expensive way to compute the residual. A better
1283 * method would be to compute (G*x + C*(C'*x) - C2*(C2'*x) - b),
1284 * using sdmult. It is done just to test the ssmult
1285 * routine. */
1286
1287 /* compute S2 = G+C*C'-C2*C2', with no sort */
1288 Ct = CHOLMOD(transpose) (C2, 1, cm) ;
1289 CC = CHOLMOD(ssmult) (C2, Ct, 0, TRUE, FALSE, cm) ;
1290 S2 = CHOLMOD(add) (S, CC, one, minusone, TRUE, FALSE, cm) ;
1291 CHOLMOD(free_sparse) (&CC, cm) ;
1292 CHOLMOD(free_sparse) (&Ct, cm) ;
1293 CHOLMOD(free_sparse) (&C2, cm) ;
1294
1295 /* Ssym is a symmetric/upper copy of S2 */
1296 Ssym = CHOLMOD(copy) (S2, 1, 1, cm) ;
1297 CHOLMOD(free_sparse) (&S2, cm) ;
1298
1299 /* compute norm (S2*x-b) */
1300 r = resid (Ssym, X, B) ;
1301 MAXERR (maxerr, r, 1) ;
1302 CHOLMOD(free_dense) (&X, cm) ;
1303
1304 /* ------------------------------------------------------------------ */
1305 /* factor S2 scratch, using same permutation */
1306 /* ------------------------------------------------------------------ */
1307
1308 if (rank == 1)
1309 {
1310 save = cm->nmethods ;
1311 save2 = cm->method [0].ordering ;
1312 cm->nmethods = 1 ;
1313 cm->method [0].ordering = CHOLMOD_GIVEN ;
1314 L2 = CHOLMOD(analyze_p) (Ssym, P, NULL, 0, cm) ;
1315 cm->nmethods = save ;
1316 cm->method [0].ordering = save2 ;
1317 CHOLMOD(factorize) (Ssym, L2, cm) ;
1318 X = CHOLMOD(solve) (CHOLMOD_A, L2, B, cm) ;
1319 r = resid (Ssym, X, B) ;
1320 MAXERR (maxerr, r, 1) ;
1321 CHOLMOD(free_dense) (&X, cm) ;
1322 CHOLMOD(free_factor) (&L2, cm) ;
1323 }
1324
1325 /* ------------------------------------------------------------------ */
1326 /* re-do the symbolic factorization on L */
1327 /* ------------------------------------------------------------------ */
1328
1329 CHOLMOD(resymbol) (Ssym, NULL, 0, TRUE, L, cm) ;
1330
1331 /* solve (G+C*C'-C2*C2')x=b again */
1332 X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
1333
1334 /* compute norm (S2*x-b) */
1335 r = resid (Ssym, X, B) ;
1336 MAXERR (maxerr, r, 1) ;
1337 CHOLMOD(free_dense) (&X, cm) ;
1338
1339 CHOLMOD(free_sparse) (&Ssym, cm) ;
1340
1341 /* ------------------------------------------------------------------ */
1342 /* downdate, with the remaining part of C, to get original L */
1343 /* ------------------------------------------------------------------ */
1344
1345 if (rank > csize)
1346 {
1347 cset = CHOLMOD(malloc) (rank-csize, sizeof (Int), cm) ;
1348 if (cset != NULL)
1349 {
1350 for (i = csize ; i < rank ; i++)
1351 {
1352 cset [i-csize] = i ;
1353 }
1354 }
1355 if (rank - csize == 4)
1356 {
1357 /* test the subset print/check routine */
1358 CHOLMOD(print_subset) (cset, rank-csize, rank, "cset", cm) ;
1359 }
1360 C2 = CHOLMOD(submatrix) (C, NULL, -1, cset, rank-csize, TRUE, TRUE,
1361 cm) ;
1362 Cperm = CHOLMOD(submatrix) (C2, P, n, NULL, -1, TRUE, TRUE, cm) ;
1363
1364 CHOLMOD(free) (rank-csize, sizeof (Int), cset, cm) ;
1365 CHOLMOD(updown) (FALSE, Cperm, L, cm) ;
1366 CHOLMOD(free_sparse) (&Cperm, cm) ;
1367 CHOLMOD(free_sparse) (&C2, cm) ;
1368
1369 /* solve the original system */
1370 X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
1371
1372 /* compute the residual */
1373 r = resid (A, X, B) ;
1374 MAXERR (maxerr, r, 1) ;
1375 CHOLMOD(free_dense) (&X, cm) ;
1376 }
1377
1378 CHOLMOD(free_sparse) (&C, cm) ;
1379 CHOLMOD(free_sparse) (&S, cm) ;
1380 }
1381
1382 /* turn memory tests back on, where we left off ] */
1383 my_tries = save3 ;
1384
1385 CHOLMOD(free_dense) (&B, cm) ;
1386 CHOLMOD(free_sparse) (&G, cm) ;
1387 CHOLMOD(free_factor) (&L, cm) ;
1388
1389 /* ---------------------------------------------------------------------- */
1390 /* factor again, then change the factor type many times */
1391 /* ---------------------------------------------------------------------- */
1392
1393 C = CHOLMOD(copy_sparse) (A, cm) ;
1394 if (C != NULL)
1395 {
1396 C->sorted = FALSE ;
1397 }
1398 L = CHOLMOD(analyze) (C, cm) ;
1399 CHOLMOD(factorize) (C, L, cm) ;
1400
1401 if (L != NULL && !(L->is_super))
1402 {
1403 CHOLMOD(resymbol) (C, NULL, 0, TRUE, L, cm) ;
1404 }
1405
1406 B = rhs (C, 1, n) ;
1407 cm->prefer_zomplex = prefer_zomplex ;
1408 X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
1409 cm->prefer_zomplex = FALSE ;
1410 r = resid (C, X, B) ;
1411 MAXERR (maxerr, r, 1) ;
1412 CHOLMOD(free_dense) (&X, cm) ;
1413
1414 for (k = 0 ; k < NFTYPES ; k++)
1415 {
1416
1417 if (co_types [k] && n > 1)
1418 {
1419 /* reallocate column zero of L, to make it non-monotonic */
1420 CHOLMOD(reallocate_column) (0, n, L, cm) ;
1421 }
1422 Lxtype = (L == NULL) ? CHOLMOD_REAL : (L->xtype) ;
1423
1424 CHOLMOD(change_factor) (Lxtype, ll_types [k], FALSE, pk_types [k],
1425 mn_types [k], L, cm) ;
1426
1427 cm->prefer_zomplex = prefer_zomplex ;
1428 X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
1429 cm->prefer_zomplex = FALSE ;
1430 r = resid (C, X, B) ;
1431 MAXERR (maxerr, r, 1) ;
1432 CHOLMOD(free_dense) (&X, cm) ;
1433 }
1434
1435 /* reallocate a column and solve again */
1436 if (n > 3)
1437 {
1438 CHOLMOD(change_factor) (CHOLMOD_REAL, FALSE, FALSE, TRUE, TRUE, L, cm) ;
1439 CHOLMOD(reallocate_column) (0, n, L, cm) ;
1440 cm->prefer_zomplex = prefer_zomplex ;
1441 X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
1442 cm->prefer_zomplex = FALSE ;
1443 r = resid (C, X, B) ;
1444 MAXERR (maxerr, r, 1) ;
1445 CHOLMOD(free_dense) (&X, cm) ;
1446 }
1447
1448 CHOLMOD(free_sparse) (&C, cm) ;
1449 CHOLMOD(free_factor) (&L, cm) ;
1450 CHOLMOD(free_dense) (&B, cm) ;
1451 CHOLMOD(free_dense) (&Breal, cm) ;
1452 CHOLMOD(free_dense) (&Bcomplex, cm) ;
1453 CHOLMOD(free_dense) (&Bzomplex, cm) ;
1454
1455 /* ---------------------------------------------------------------------- */
1456 /* factorize and solve (A*A'+beta*I)x=b or A'x=b */
1457 /* ---------------------------------------------------------------------- */
1458
1459 if (A->stype == 0)
1460 {
1461 double *Rx, *Rz, *Xx, *Xz ;
1462 double beta [2] ;
1463 beta [0] = 3.14159 ;
1464 beta [1] = 0 ;
1465 L = CHOLMOD(analyze) (A, cm) ;
1466 CHOLMOD(factorize_p) (A, beta, NULL, 0, L, cm) ;
1467 B = rhs (A, 1, n) ;
1468 cm->prefer_zomplex = prefer_zomplex ;
1469 X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
1470 cm->prefer_zomplex = FALSE ;
1471
1472 /* compute the residual */
1473
1474 /* W = A'*X */
1475 W = CHOLMOD(zeros) (ncol, 1, xtype, cm) ;
1476 CHOLMOD(sdmult) (A, 2, one, zero, X, W, cm) ;
1477
1478 /* R = B */
1479 R = CHOLMOD(copy_dense) (B, cm) ;
1480
1481 /* R = A*W - R */
1482 CHOLMOD(sdmult) (A, 0, one, minusone, W, R, cm) ;
1483
1484 /* R = R + beta*X */
1485 if (R != NULL && X != NULL)
1486 {
1487 Rx = R->x ;
1488 Rz = R->z ;
1489 Xx = X->x ;
1490 Xz = X->z ;
1491
1492 for (i = 0 ; i < nrow ; i++)
1493 {
1494 switch (xtype)
1495 {
1496 case CHOLMOD_REAL:
1497 Rx [i] += beta [0] * Xx [i] ;
1498 break ;
1499
1500 case CHOLMOD_COMPLEX:
1501 Rx [2*i ] += beta [0] * Xx [2*i ] ;
1502 Rx [2*i+1] += beta [0] * Xx [2*i+1] ;
1503 break ;
1504
1505 case CHOLMOD_ZOMPLEX:
1506 Rx [i] += beta [0] * Xx [i] ;
1507 Rz [i] += beta [0] * Xz [i] ;
1508 break ;
1509 }
1510 }
1511 }
1512
1513 r = CHOLMOD(norm_dense) (R, 2, cm) ;
1514 bnorm = CHOLMOD(norm_dense) (B, 2, cm) ;
1515 xnorm = CHOLMOD(norm_dense) (X, 2, cm) ;
1516 norm = MAX (r, xnorm) ;
1517 if (norm > 0)
1518 {
1519 r /= norm ;
1520 }
1521 MAXERR (maxerr, r, 1) ;
1522
1523 CHOLMOD(free_dense) (&X, cm) ;
1524 CHOLMOD(free_dense) (&R, cm) ;
1525 CHOLMOD(free_dense) (&W, cm) ;
1526 CHOLMOD(free_factor) (&L, cm) ;
1527 CHOLMOD(free_dense) (&B, cm) ;
1528 }
1529
1530 /* ---------------------------------------------------------------------- */
1531 /* test rowdel and updown */
1532 /* ---------------------------------------------------------------------- */
1533
1534 if (isreal && A->stype == 1 && n > 0 && n < NLARGE)
1535 {
1536 Int save4, save5, save6 ;
1537 save4 = cm->nmethods ;
1538 save5 = cm->method [0].ordering ;
1539 save6 = cm->supernodal ;
1540
1541 cm->nmethods = 1 ;
1542 cm->method [0].ordering = CHOLMOD_NATURAL ;
1543 cm->supernodal = CHOLMOD_SUPERNODAL ;
1544
1545 B = rhs (A, 1, n) ;
1546 L = CHOLMOD(analyze) (A, cm) ;
1547 CHOLMOD(factorize) (A, L, cm) ;
1548
1549 /* solve Ax=b */
1550 X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
1551 r = resid (A, X, B) ;
1552 MAXERR (maxerr, r, 1) ;
1553 CHOLMOD(free_dense) (&X, cm) ;
1554
1555 /* determine the new system with row/column k missing */
1556 k = n/2 ;
1557 S = CHOLMOD(copy) (A, 0, 1, cm) ;
1558 RowK = CHOLMOD(submatrix) (S, NULL, -1, &k, 1, TRUE, TRUE, cm) ;
1559 CHOLMOD(print_sparse) (S, "S", cm) ;
1560 CHOLMOD(print_sparse) (RowK, "RowK of S", cm) ;
1561 CHOLMOD(free_sparse) (&RowK, cm) ;
1562 prune_row (S, k) ;
1563 if (S != NULL)
1564 {
1565 S->stype = 1 ;
1566 }
1567
1568 /* delete row k of L (converts to LDL') */
1569 /* printf ("rowdel here:\n") ; */
1570 CHOLMOD(rowdel) (k, NULL, L, cm) ;
1571 CHOLMOD(resymbol) (S, NULL, 0, TRUE, L, cm) ;
1572
1573 /* solve with row k missing */
1574 X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
1575 r = resid (S, X, B) ;
1576 MAXERR (maxerr, r, 1) ;
1577 CHOLMOD(free_dense) (&X, cm) ;
1578 CHOLMOD(free_sparse) (&S, cm) ;
1579
1580 /* factorize again */
1581 CHOLMOD(free_factor) (&L, cm) ;
1582 L = CHOLMOD(analyze) (A, cm) ;
1583 CHOLMOD(factorize) (A, L, cm) ;
1584
1585 /* rank-3 update (converts to LDL') and solve */
1586 C = CHOLMOD(speye) (n, 3, CHOLMOD_REAL, cm) ;
1587 CC = CHOLMOD(aat) (C, NULL, 0, 1, cm) ;
1588 S = CHOLMOD(add) (A, CC, one, one, TRUE, TRUE, cm) ;
1589 if (S != NULL)
1590 {
1591 S->stype = 1 ;
1592 }
1593 CHOLMOD(updown) (TRUE, C, L, cm) ;
1594 X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
1595 r = resid (S, X, B) ;
1596 MAXERR (maxerr, r, 1) ;
1597 CHOLMOD(free_dense) (&X, cm) ;
1598
1599 /* free everything */
1600 CHOLMOD(free_sparse) (&S, cm) ;
1601 CHOLMOD(free_sparse) (&CC, cm) ;
1602 CHOLMOD(free_sparse) (&C, cm) ;
1603 CHOLMOD(free_sparse) (&S, cm) ;
1604 CHOLMOD(free_factor) (&L, cm) ;
1605 CHOLMOD(free_dense) (&B, cm) ;
1606
1607 cm->nmethods = save4 ;
1608 cm->method [0].ordering = save5 ;
1609 cm->supernodal = save6 ;
1610 }
1611
1612 /* ---------------------------------------------------------------------- */
1613 /* free remaining workspace */
1614 /* ---------------------------------------------------------------------- */
1615
1616 OK (CHOLMOD(print_common) ("cm", cm)) ;
1617 progress (0, '.') ;
1618 return (maxerr) ;
1619 }
1620