1 /* ========================================================================== */
2 /* === Tcov/raw_factor ====================================================== */
3 /* ========================================================================== */
4 
5 /* -----------------------------------------------------------------------------
6  * CHOLMOD/Tcov Module.  Copyright (C) 2005-2006, Timothy A. Davis
7  * http://www.suitesparse.com
8  * -------------------------------------------------------------------------- */
9 
10 /* Factorize A using cholmod_rowfac for the simplicial case, and the
11  * cholmod_super_* routines for the supernodal case, and test the solution to
12  * linear systems. */
13 
14 #include "cm.h"
15 
16 
17 /* ========================================================================== */
18 /* === icomp ================================================================ */
19 /* ========================================================================== */
20 
21 /* for sorting by qsort */
icomp(Int * i,Int * j)22 static int icomp (Int *i, Int *j)
23 {
24     if (*i < *j)
25     {
26 	return (-1) ;
27     }
28     else
29     {
30 	return (1) ;
31     }
32 }
33 
34 
35 /* ========================================================================== */
36 /* === add_gunk ============================================================= */
37 /* ========================================================================== */
38 
add_gunk(cholmod_sparse * A)39 static cholmod_sparse *add_gunk (cholmod_sparse *A)
40 {
41     cholmod_sparse *S ;
42     double *Sx, *Sz ;
43     Int *Sp, *Si, nz, p, save3, j, n ;
44 
45     if (A == NULL) return (NULL) ;
46 
47     /* save3 = cm->print ; cm->print = 5 ; */
48 
49     A->nzmax++ ;
50     S = CHOLMOD(copy_sparse) (A, cm) ;
51     A->nzmax-- ;
52 
53     /* add a S(n,1)=1 entry to the matrix */
54     if (S != NULL)
55     {
56 	S->sorted = FALSE ;
57 	Sx = S->x ;
58 	Si = S->i ;
59 	Sp = S->p ;
60 	Sz = S->z ;
61 	n = S->ncol ;
62 	nz = Sp [n] ;
63 	for (j = 1 ; j <= n ; j++)
64 	{
65 	    Sp [j]++ ;
66 	}
67 	if (S->xtype == CHOLMOD_REAL)
68 	{
69 	    for (p = nz-1 ; p >= 0 ; p--)
70 	    {
71 		Si [p+1] = Si [p] ;
72 		Sx [p+1] = Sx [p] ;
73 	    }
74 	    Si [0] = n-1 ;
75 	    Sx [0] = 99999 ;
76 	}
77 	else if (S->xtype == CHOLMOD_COMPLEX)
78 	{
79 	    for (p = nz-1 ; p >= 0 ; p--)
80 	    {
81 		Si [p+1] = Si [p] ;
82 		Sx [2*p+2] = Sx [2*p] ;
83 		Sx [2*p+3] = Sx [2*p+1] ;
84 	    }
85 	    Si [0] = n-1 ;
86 	    Sx [0] = 99999 ;
87 	    Sx [1] = 0 ;
88 	}
89 	else if (S->xtype == CHOLMOD_ZOMPLEX)
90 	{
91 	    for (p = nz-1 ; p >= 0 ; p--)
92 	    {
93 		Si [p+1] = Si [p] ;
94 		Sx [p+1] = Sx [p] ;
95 		Sz [p+1] = Sz [p] ;
96 	    }
97 	    Si [0] = n-1 ;
98 	    Sx [0] = 99999 ;
99 	    Sz [0] = 0 ;
100 	}
101     }
102 
103     /* CHOLMOD(print_sparse) (A, "A for gunk", cm) ; */
104     /* CHOLMOD(print_sparse) (S, "S with gunk", cm) ; */
105     /* cm->print = save3 ; */
106 
107     return (S) ;
108 }
109 
110 
111 /* ========================================================================== */
112 /* === raw_factor =========================================================== */
113 /* ========================================================================== */
114 
115 /* Factor A, without using any fill-reducing permutation.  This may fail due
116  * to catastrophic fill-in (which is the desired test result for a large
117  * arrowhead matrix).
118  */
119 
raw_factor(cholmod_sparse * A,Int check_errors)120 double raw_factor (cholmod_sparse *A, Int check_errors)
121 {
122     double maxerr = 0, r, anorm ;
123     cholmod_sparse *AT, *C, *LT, *Lsparse, *S, *ST, *R, *A1 ;
124     cholmod_factor *L, *Lcopy ;
125     cholmod_dense *X, *W, *B, *X2 ;
126     Int i, k, n, ok, ok1, ok2, trial, rnz, lnz, Lxtype, Axtype, posdef,
127 	prefer_zomplex, Bxtype ;
128     Int *Parent, *Post, *First, *Level, *Ri, *Rp, *LTp = NULL, *LTi = NULL, *P,
129 	*mask, *RLinkUp ;
130     SuiteSparse_long lr ;
131     double beta [2] ;
132     unsigned SuiteSparse_long save ;
133 
134     /* ---------------------------------------------------------------------- */
135     /* create the problem */
136     /* ---------------------------------------------------------------------- */
137 
138     if (A == NULL || A->stype != 1)
139     {
140 	return (0) ;
141     }
142 
143     W = NULL ;
144     X2 = NULL ;
145     L = NULL ;
146     n = A->nrow ;
147     B = rhs (A, 1, n) ;
148     AT = CHOLMOD(transpose) (A, 2, cm) ;
149     Parent = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
150     Post = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
151     First = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
152     Level = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
153     beta [0] = 0 ;
154     beta [1] = 0 ;
155     anorm = CHOLMOD(norm_sparse) (A, 1, cm) ;
156 
157     prefer_zomplex = (A->xtype == CHOLMOD_ZOMPLEX) ;
158     Bxtype = A->xtype ;
159 
160     /* ---------------------------------------------------------------------- */
161     /* supernodal factorization */
162     /* ---------------------------------------------------------------------- */
163 
164     L = CHOLMOD(allocate_factor) (n, cm) ;
165     ok1 = CHOLMOD(etree) (A, Parent, cm) ;
166     lr = CHOLMOD(postorder) (Parent, n, NULL, Post, cm) ;
167     ok2 = CHOLMOD(rowcolcounts) (AT, NULL, 0, Parent, Post,
168 	NULL, (L != NULL) ? (L->ColCount) : NULL, First, Level, cm) ;
169 
170     if (ok2)
171     {
172 	printf ("raw_factor: cm->fl %g cm->lnz %g\n", cm->fl, cm->lnz) ;
173     }
174 
175     if (check_errors)
176     {
177 	OKP (AT) ;
178 	OKP (Parent) ;
179 	OKP (Post) ;
180 	OKP (First) ;
181 	OKP (Level) ;
182 	OK (AT->stype == -1) ;
183 	OKP (L) ;
184 	OK (ok1) ;
185 	OK (ok2) ;
186 	OK (lr >= 0) ;
187 
188 	/* rowcolcounts requires A in symmetric lower form */
189 	ok = CHOLMOD(rowcolcounts) (A, NULL, 0, Parent, Post,
190 	    NULL, L->ColCount, First, Level, cm) ;		    NOT (ok) ;
191     }
192 
193     /* super_symbolic needs A in upper form, so this will succeed
194      * unless the problem is huge */
195     ok = CHOLMOD(super_symbolic) (A, NULL, Parent, L, cm) ;
196 
197     /* super_symbolic should fail if lnz is too large */
198     if (cm->lnz > Size_max / 2)
199     {
200 	printf ("raw_factor: problem is huge\n") ;
201 	NOT (ok) ;
202 	OK (L->xtype == CHOLMOD_PATTERN && !(L->is_super)) ;
203 
204 	/* try changing to LDL packed, which should also fail */
205 	ok = CHOLMOD(change_factor) (CHOLMOD_REAL, FALSE, FALSE, TRUE, TRUE,
206 		L, cm) ;
207 	NOT (ok) ;
208 
209 	CHOLMOD(free_factor) (&L, cm) ;
210 	CHOLMOD(free_sparse) (&AT, cm) ;
211 	CHOLMOD(free_dense) (&B, cm) ;
212 	CHOLMOD(free) (n, sizeof (Int), First, cm) ;
213 	CHOLMOD(free) (n, sizeof (Int), Level, cm) ;
214 	CHOLMOD(free) (n, sizeof (Int), Parent, cm) ;
215 	CHOLMOD(free) (n, sizeof (Int), Post, cm) ;
216 	return (0) ;
217     }
218 
219     if (check_errors)
220     {
221 
222 	if (cm->status == CHOLMOD_OUT_OF_MEMORY
223 	 || cm->status == CHOLMOD_TOO_LARGE)
224 	{
225 	    /* no test case will reach here, but check just to be safe */
226 	    printf ("raw_factor: out of memory for symbolic case %d\n",
227 		cm->status) ;
228 	    CHOLMOD(free_factor) (&L, cm) ;
229 	    CHOLMOD(free_sparse) (&AT, cm) ;
230 	    CHOLMOD(free_dense) (&B, cm) ;
231 	    CHOLMOD(free) (n, sizeof (Int), First, cm) ;
232 	    CHOLMOD(free) (n, sizeof (Int), Level, cm) ;
233 	    CHOLMOD(free) (n, sizeof (Int), Parent, cm) ;
234 	    CHOLMOD(free) (n, sizeof (Int), Post, cm) ;
235 	    return (0) ;
236 	}
237 
238 	OK (ok) ;
239 	ok = CHOLMOD(super_symbolic)(A, NULL, Parent, L, cm) ;	    NOT (ok) ;
240 	ok = CHOLMOD(super_symbolic)(AT, NULL, Parent, L, cm) ;	    NOT (ok) ;
241 	ok = CHOLMOD(super_symbolic)(NULL, NULL, Parent, L, cm) ;   NOT (ok) ;
242 	ok = CHOLMOD(super_symbolic)(A, NULL, NULL, L, cm) ;	    NOT (ok) ;
243 	ok = CHOLMOD(super_symbolic)(A, NULL, Parent, NULL, cm) ;   NOT (ok) ;
244     }
245 
246     /* super_numeric needs A in lower form, so this will succeed unless
247      * the problem is huge */
248     ok = CHOLMOD(super_numeric) (AT, NULL, Zero, L, cm) ;
249 
250     if (check_errors)
251     {
252 
253 	if (cm->status == CHOLMOD_OUT_OF_MEMORY)
254 	{
255 	    /* For the 64-bit case, the Matrix/a1 problem will reach here */
256 	    printf ("raw_factor: out of memory for numeric case\n") ;
257 	    CHOLMOD(free_factor) (&L, cm) ;
258 	    CHOLMOD(free_sparse) (&AT, cm) ;
259 	    CHOLMOD(free_dense) (&B, cm) ;
260 	    CHOLMOD(free) (n, sizeof (Int), First, cm) ;
261 	    CHOLMOD(free) (n, sizeof (Int), Level, cm) ;
262 	    CHOLMOD(free) (n, sizeof (Int), Parent, cm) ;
263 	    CHOLMOD(free) (n, sizeof (Int), Post, cm) ;
264 	    return (0) ;
265 	}
266 
267 	OK (ok) ;
268 	ok = CHOLMOD(super_numeric)(A, NULL, Zero, L, cm) ;	    NOT (ok) ;
269 	ok = CHOLMOD(super_numeric)(NULL, NULL, Zero, L, cm) ;	    NOT (ok) ;
270 	ok = CHOLMOD(super_numeric)(AT, NULL, Zero, NULL, cm) ;	    NOT (ok) ;
271     }
272 
273     /* solve */
274     Lxtype = (L == NULL) ? CHOLMOD_REAL : L->xtype ;
275     W = CHOLMOD(zeros) (n, 1, Lxtype, cm) ;
276     X = CHOLMOD(copy_dense) (B, cm) ;
277     if (Bxtype == CHOLMOD_ZOMPLEX)
278     {
279 	CHOLMOD(dense_xtype) (CHOLMOD_COMPLEX, X, cm) ;
280     }
281 
282     CHOLMOD(print_factor) (L, "L for super l/ltsolve", cm) ;
283     CHOLMOD(print_dense) (W, "W", cm) ;
284     CHOLMOD(print_dense) (X, "X", cm) ;
285 
286     ok1 = CHOLMOD(super_lsolve) (L, X, W, cm) ;
287     CHOLMOD(print_dense) (X, "X", cm) ;
288 
289     ok2 = CHOLMOD(super_ltsolve) (L, X, W, cm) ;
290     CHOLMOD(print_dense) (X, "X", cm) ;
291 
292     if (Bxtype == CHOLMOD_ZOMPLEX)
293     {
294 	CHOLMOD(dense_xtype) (CHOLMOD_ZOMPLEX, X, cm) ;
295     }
296 
297     r = resid (A, X, B) ;
298     MAXERR (maxerr, r, 1) ;
299 
300     if (Bxtype == CHOLMOD_ZOMPLEX)
301     {
302 	CHOLMOD(dense_xtype) (CHOLMOD_COMPLEX, X, cm) ;
303     }
304 
305     if (check_errors)
306     {
307 	OKP (W) ;
308 	OKP (X) ;
309 	OK (ok1) ;
310 	OK (ok2) ;
311 	ok = CHOLMOD(super_lsolve) (NULL, X, W, cm) ;		    NOT (ok) ;
312 	ok = CHOLMOD(super_ltsolve) (NULL, X, W, cm) ;		    NOT (ok) ;
313 	ok = CHOLMOD(super_lsolve) (L, NULL, W, cm) ;		    NOT (ok) ;
314 	ok = CHOLMOD(super_ltsolve) (L, NULL, W, cm) ;		    NOT (ok) ;
315 	ok = CHOLMOD(super_lsolve) (L, X, NULL, cm) ;		    NOT (ok) ;
316 	ok = CHOLMOD(super_ltsolve) (L, X, NULL, cm) ;		    NOT (ok) ;
317 
318 	if (L != NULL && L->maxesize > 1)
319 	{
320 	    /* W is too small */
321 	    ok = CHOLMOD(free_dense) (&W, cm) ;			    OK (ok) ;
322 	    W = CHOLMOD(zeros) (1, 1, Lxtype, cm) ;		    OKP (W) ;
323 	    ok = CHOLMOD(super_lsolve) (L, X, W, cm) ;		    NOT (ok) ;
324 	    ok = CHOLMOD(super_ltsolve) (L, X, W, cm) ;		    NOT (ok) ;
325 	    ok = CHOLMOD(free_dense) (&W, cm) ;			    OK (ok) ;
326 	    W = CHOLMOD(zeros) (n, 1, Lxtype, cm) ;		    OKP (W) ;
327 	}
328 
329 	/* X2 has the wrong dimensions */
330 	X2 = CHOLMOD(zeros) (n+1, 1, Lxtype, cm) ;		    OKP (X2) ;
331 	ok = CHOLMOD(super_lsolve) (L, X2, W, cm) ;		    NOT (ok) ;
332 	ok = CHOLMOD(super_ltsolve) (L, X2, W, cm) ;		    NOT (ok) ;
333 	CHOLMOD(free_dense) (&X2, cm) ;
334     }
335 
336     CHOLMOD(free_dense) (&X, cm) ;
337 
338     /* X2 is n-by-0, which is OK */
339     X2 = CHOLMOD(zeros) (n, 0, Lxtype, cm) ;
340     ok1 = CHOLMOD(super_lsolve) (L, X2, W, cm) ;
341     ok2 = CHOLMOD(super_ltsolve) (L, X2, W, cm) ;
342     CHOLMOD(free_dense) (&W, cm) ;
343     CHOLMOD(free_dense) (&X2, cm) ;
344 
345     if (check_errors)
346     {
347 	OK (ok1) ;
348 	OK (ok2) ;
349 	test_memory_handler ( ) ;
350 	my_tries = 0 ;
351 	ok = CHOLMOD(super_symbolic) (A, NULL, Parent, L, cm) ;	    NOT (ok) ;
352 	ok = CHOLMOD(super_numeric) (AT, NULL, Zero, L, cm) ;	    NOT (ok) ;
353 	normal_memory_handler ( ) ;
354 	cm->error_handler = NULL ;
355     }
356 
357     /* R = space for result of row_subtree and row_lsubtree */
358     R = CHOLMOD(allocate_sparse)(n, 1, n, FALSE, TRUE, 0, CHOLMOD_PATTERN, cm) ;
359 
360     /* ---------------------------------------------------------------------- */
361     /* erroneous factorization */
362     /* ---------------------------------------------------------------------- */
363 
364     /* cannot use rowfac or row_lsubtree on a supernodal factorization */
365     if (check_errors && n > 0)
366     {
367 	ok = CHOLMOD(rowfac) (A, NULL, beta, 0, 0, L, cm) ;	    NOT (ok) ;
368 	ok = CHOLMOD(row_lsubtree) (A, &i, 0, n-1, L, R, cm) ;	    NOT (ok) ;
369     }
370 
371     /* ---------------------------------------------------------------------- */
372     /* convert to simplicial LDL' */
373     /* ---------------------------------------------------------------------- */
374 
375     CHOLMOD(change_factor) (Lxtype, FALSE, FALSE, TRUE, TRUE, L, cm) ;
376 
377     /* remove entries due to relaxed supernodal amalgamation */
378     CHOLMOD(resymbol) (A, NULL, 0, TRUE, L, cm) ;
379 
380     /* refactorize a numeric factor */
381     posdef = 0 ;   /* unknown */
382     if (A != NULL && A->stype >= 0)
383     {
384 	if (A->stype > 0 && A->packed)
385 	{
386 	    S = add_gunk (A) ;
387 	    CHOLMOD(rowfac) (S, NULL, beta, 0, n, L, cm) ;
388 	    if (S && S->xtype == CHOLMOD_COMPLEX)
389 	    {
390 		CHOLMOD(sparse_xtype) (CHOLMOD_ZOMPLEX, S, cm) ;
391 	    }
392 	    ok = CHOLMOD(free_sparse) (&S, cm) ;			    OK (ok) ;
393 	}
394 	else
395 	{
396 	    CHOLMOD(rowfac) (A, NULL, beta, 0, n, L, cm) ;
397 	}
398 	posdef = (cm->status == CHOLMOD_OK) ;
399     }
400 
401     /* convert to a sparse matrix, and transpose L */
402     Lcopy = CHOLMOD(copy_factor)(L, cm) ;
403     Lsparse = CHOLMOD(factor_to_sparse) (L, cm) ;
404 
405     LT = CHOLMOD(transpose) (Lsparse, 0, cm) ;
406     CHOLMOD(free_sparse) (&Lsparse, cm) ;
407 
408     if (LT != NULL)
409     {
410 	LTp = LT->p ;
411 	LTi = LT->i ;
412 	OK (LT->packed) ;
413     }
414 
415     /* remove the unit diagonal of LT */
416     CHOLMOD(band_inplace) (1, n, -1, LT, cm) ;
417 
418     /* ST = pattern of A(p,p)' */
419     P = (L == NULL) ? NULL : L->Perm ;
420     ST = CHOLMOD(ptranspose) (A, 0, P, NULL, 0, cm) ;
421 
422     /* S = pattern of A(p,p) */
423     S = CHOLMOD(transpose) (ST, 0, cm) ;
424     ok = CHOLMOD(free_sparse) (&ST, cm) ;
425 
426     if (R != NULL && LT != NULL && posdef && A != NULL && A->stype >= 0
427 	    && S != NULL && Lcopy != NULL)
428     {
429 	LTp = LT->p ;
430 	LTi = LT->i ;
431 
432 	Ri = R->i ;
433 	Rp = R->p ;
434 
435 	save = my_seed ( ) ;					/* RAND */
436 	for (trial = 0 ; trial < 30 ; trial++)
437 	{
438 	    /* pick a row at random */
439 	    i = nrand (n) ;					/* RAND */
440 
441 	    /* compute R = pattern of L(i,0:i-1), using row subtrees */
442 	    ok = CHOLMOD(row_subtree) (S, NULL, i, Parent, R, cm) ;
443 	    if (!ok)
444 	    {
445 		break ;
446 	    }
447 	    rnz = Rp [1] ;
448 
449 	    /* sort R */
450 	    qsort (Ri, rnz, sizeof (Int),
451 		    (int (*) (const void *, const void *)) icomp) ;
452 
453 	    /* compare with ith column of L transpose */
454 	    lnz = LTp [i+1] - LTp [i] ;
455 	    ok = TRUE ;
456 	    for (k = 0 ; k < MIN (rnz,lnz) ; k++)
457 	    {
458 		/* printf ("%d vs %d\n", Ri [k], LTi [LTp [i] + k]) ; */
459 		ok = ok && (Ri [k] == LTi [LTp [i] + k]) ;
460 	    }
461 	    OK (ok) ;
462 	    OK (rnz == lnz) ;
463 
464 	    /* compute R = pattern of L(i,0:i-1), using row lsubtrees */
465 	    ok = CHOLMOD(row_lsubtree) (S, NULL, 0, i, Lcopy, R, cm) ;
466 	    if (!ok)
467 	    {
468 		break ;
469 	    }
470 	    rnz = Rp [1] ;
471 
472 	    /* sort R */
473 	    qsort (Ri, rnz, sizeof (Int),
474 		    (int (*) (const void *, const void *)) icomp) ;
475 
476 	    /* compare with ith column of L transpose */
477 	    lnz = LTp [i+1] - LTp [i] ;
478 	    ok = TRUE ;
479 	    for (k = 0 ; k < MIN (rnz,lnz) ; k++)
480 	    {
481 		/* printf ("%d vs %d\n", Ri [k], LTi [LTp [i] + k]) ; */
482 		ok = ok && (Ri [k] == LTi [LTp [i] + k]) ;
483 	    }
484 	    OK (ok) ;
485 	    OK (rnz == lnz) ;
486 
487 	    /* L is symbolic, so cholmod_lsubtree will fail */
488 	    if (check_errors)
489 	    {
490 		ok = CHOLMOD(row_lsubtree) (S, NULL, 0, i, L, R, cm) ;
491 		NOT (ok) ;
492 	    }
493 
494 	}
495 	my_srand (save) ;					/* RAND */
496     }
497 
498     ok = CHOLMOD(free_factor) (&L, cm) ;			    OK (ok) ;
499     ok = CHOLMOD(free_factor) (&Lcopy, cm) ;			    OK (ok) ;
500     ok = CHOLMOD(free_sparse) (&LT, cm) ;			    OK (ok) ;
501     ok = CHOLMOD(free_sparse) (&R, cm) ;			    OK (ok) ;
502     ok = CHOLMOD(free_sparse) (&S, cm) ;			    OK (ok) ;
503 
504     /* ---------------------------------------------------------------------- */
505     /* simplicial LDL' or LL' factorization with no analysis */
506     /* ---------------------------------------------------------------------- */
507 
508     for (trial = 0 ; trial <= check_errors ; trial++)
509     {
510 
511 	/* create a simplicial symbolic factor */
512 	L = CHOLMOD(allocate_factor) (n, cm) ;
513 	ok = TRUE ;
514 	Axtype = (A == NULL) ? CHOLMOD_REAL : A->xtype ;
515 
516 	if (check_errors)
517 	{
518 	    OKP (L) ;
519 	    if (trial == 0)
520 	    {
521 		/* convert to packed LDL' first, then unpacked */
522 		ok = CHOLMOD(change_factor) (Axtype, FALSE, FALSE, TRUE,
523 			TRUE, L, cm) ;
524 		OK (ok);
525 		ok = CHOLMOD(change_factor) (Axtype, FALSE, FALSE, FALSE,
526 			TRUE, L, cm) ;
527 		OK (ok) ;
528 	    }
529 	    else if (trial == 1)
530 	    {
531 		ok = CHOLMOD(rowfac)(NULL, NULL, beta, 0, 0, L, cm) ; NOT (ok) ;
532 		ok = CHOLMOD(rowfac)(A, NULL, beta, 0, 0, NULL, cm) ; NOT (ok) ;
533 		ok = CHOLMOD(rowfac)(AT, NULL, beta, 0, 0, L, cm) ;   NOT (ok) ;
534 		if (n > 1)
535 		{
536 		    A1 = CHOLMOD(allocate_sparse)(1, 1, 1, TRUE, TRUE, 1,
537 			CHOLMOD_PATTERN, cm) ;			      OKP (A1) ;
538 		    ok = CHOLMOD(rowfac)(A1, NULL, beta, 0, 0, L, cm); NOT (ok);
539 		    ok = CHOLMOD(free_sparse)(&A1, cm) ;		OK (ok);
540 		}
541 	    }
542 	    else
543 	    {
544 		/* convert to symbolic LL' */
545 		ok = CHOLMOD(change_factor) (CHOLMOD_PATTERN, TRUE, FALSE, TRUE,
546 			TRUE, L, cm) ;
547 		OK (ok) ;
548 		OK (L->is_ll) ;
549 	    }
550 	}
551 
552 	/* factor */
553 	CHOLMOD(print_factor) (L, "L for rowfac", cm) ;
554 	CHOLMOD(print_sparse) (A, "A for rowfac", cm) ;
555 
556 	cm->dbound = 1e-15 ;
557 	for (k = 0 ; ok && k < n ; k++)
558 	{
559 	    if (!CHOLMOD(rowfac) (A, NULL, beta, k, k+1, L, cm))
560 	    {
561 		ok = FALSE ;
562 	    }
563 	    if (cm->status == CHOLMOD_NOT_POSDEF)
564 	    {
565 		/* LL' factorization failed; subsequent rowfac's should fail */
566 		k++ ;
567 		ok = CHOLMOD(rowfac) (A, NULL, beta, k, k+1, L, cm) ;
568 		NOT (ok) ;
569 		ok = TRUE ;
570 	    }
571 	}
572 	cm->dbound = 0 ;
573 
574 	if (check_errors)
575 	{
576 	    OK (ok) ;
577 	    ok = CHOLMOD(rowfac) (A, NULL, beta, n, n+1, L, cm) ;    NOT (ok) ;
578 	    ok = TRUE ;
579 	}
580 
581 	/* solve */
582 	if (ok)
583 	{
584 
585 /*
586 int saveit = cm->print ;
587 int saveit2 = cm->precise ;
588 cm->print = 5 ;
589 cm->precise = TRUE ;
590 
591 CHOLMOD (print_sparse) (A, "A here", cm) ;
592 CHOLMOD (print_factor) (L, "L here", cm) ;
593 CHOLMOD (print_dense) (B, "B here", cm) ;
594 */
595 
596 	    cm->prefer_zomplex = prefer_zomplex ;
597 	    X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
598 
599 /*
600 CHOLMOD (print_dense) (X, "X here", cm) ;
601 */
602 
603 	    cm->prefer_zomplex = FALSE ;
604 	    r = resid (A, X, B) ;
605 	    MAXERR (maxerr, r, 1) ;
606 	    CHOLMOD(free_dense) (&X, cm) ;
607 
608 /*
609 cm->print = saveit ;
610 cm->precise = saveit2 ;
611 fprintf (stderr, "solve %8.2e\n", r) ;
612 */
613 	}
614 
615 	CHOLMOD(free_factor) (&L, cm) ;
616     }
617 
618     /* ---------------------------------------------------------------------- */
619     /* factor again with entries in the (ignored) lower part A */
620     /* ---------------------------------------------------------------------- */
621 
622     if (A->packed)
623     {
624 	L = CHOLMOD(allocate_factor) (n, cm) ;
625 	C = add_gunk (A) ;
626 
627 /*
628 	C = CHOLMOD(copy) (A, 0, 1, cm) ;
629 	if (C != NULL)
630 	{
631 	    C->stype = 1 ;
632 	}
633 */
634 
635 	CHOLMOD(rowfac) (C, NULL, beta, 0, n, L, cm) ;
636 
637 	X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
638 
639 	r = resid (A, X, B) ;
640 	MAXERR (maxerr, r, 1) ;
641 	CHOLMOD(free_sparse) (&C, cm) ;
642 	CHOLMOD(free_factor) (&L, cm) ;
643 	CHOLMOD(free_dense) (&X, cm) ;
644     }
645 
646     /* ---------------------------------------------------------------------- */
647     /* factor again using rowfac_mask (for LPDASA only) */
648     /* ---------------------------------------------------------------------- */
649 
650     r = raw_factor2 (A, 0., 0) ;
651     MAXERR (maxerr, r, 1) ;
652 
653     r = raw_factor2 (A, 1e-16, 0) ;
654     MAXERR (maxerr, r, 1) ;
655 
656     /* ---------------------------------------------------------------------- */
657     /* free the problem */
658     /* ---------------------------------------------------------------------- */
659 
660     CHOLMOD(free_sparse) (&AT, cm) ;
661     CHOLMOD(free_dense) (&B, cm) ;
662     CHOLMOD(free) (n, sizeof (Int), First, cm) ;
663     CHOLMOD(free) (n, sizeof (Int), Level, cm) ;
664     CHOLMOD(free) (n, sizeof (Int), Parent, cm) ;
665     CHOLMOD(free) (n, sizeof (Int), Post, cm) ;
666     progress (0, '.') ;
667     return (maxerr) ;
668 }
669 
670 
671 /* ========================================================================== */
672 /* === raw_factor2 ========================================================== */
673 /* ========================================================================== */
674 
675 /* A->stype can be 0 (lower), 1 (upper) or 0 (unsymmetric).  In the first two
676  * cases, Ax=b is solved.  In the third, A*A'x=b is solved.  No analysis and no
677  * fill-reducing ordering is used.  Both simplicial LL' and LDL' factorizations
678  * are used (testing rowfac_mask, for LPDASA only). */
679 
raw_factor2(cholmod_sparse * A,double alpha,int domask)680 double raw_factor2 (cholmod_sparse *A, double alpha, int domask)
681 {
682     Int n, i, prefer_zomplex, is_ll, xtype, sorted, axtype, stype ;
683     Int *mask = NULL, *RLinkUp = NULL, nz = 0 ;
684     Int *Cp = NULL, added_gunk ;
685     double maxerr = 0, r = 0 ;
686     cholmod_sparse *AT = NULL, *C = NULL, *CT = NULL, *CC = NULL, *C2 = NULL ;
687     cholmod_factor *L = NULL ;
688     cholmod_dense *B = NULL, *X = NULL ;
689     double beta [2] ;
690 
691 /*
692 int saveit = cm->print ;
693 int saveit2 = cm->precise ;
694 cm->print = 5 ;
695 cm->precise = TRUE ;
696 */
697 
698     if (A == NULL)
699     {
700 	return (0) ;
701     }
702     n = A->nrow ;
703     if (n > 1000)
704     {
705 	printf ("\nSkipping rowfac, matrix too large\n") ;
706 	return (0) ;
707     }
708     axtype = A->xtype ;
709     beta [0] = alpha ;
710     beta [1] = 0 ;
711 
712     prefer_zomplex = (A->xtype == CHOLMOD_ZOMPLEX) ;
713     AT = CHOLMOD(transpose) (A, 2, cm) ;
714 
715     /* ensure C has stype of 0 or 1.  Do not prune any entries */
716     stype = A->stype ;
717     if (stype >= 0)
718     {
719 	A->stype = 0 ;
720 	C = CHOLMOD(copy_sparse) (A, cm) ;
721 	A->stype = stype ;
722 	if (C) C->stype = stype ;
723 	/*
724 	C = CHOLMOD(copy) (A, 0, 1, cm) ;
725 	if (C) C->stype = A->stype ;
726 	*/
727 	CT = AT ;
728 
729     }
730     else
731     {
732 	C = AT ;
733 	/*
734 	CT = CHOLMOD(copy_sparse) (A, cm) ;
735 	*/
736 	/*
737 	CT = CHOLMOD(copy) (A, 0, 1, cm) ;
738 	if (CT) CT->stype = A->stype ;
739 	*/
740 	A->stype = 0 ;
741 	CT = CHOLMOD(copy_sparse) (A, cm) ;
742 	A->stype = stype ;
743 	if (CT) CT->stype = stype ;
744 
745 	/* only domask if C is symmetric and upper part stored */
746 	domask = FALSE ;
747 
748     }
749 
750     mask = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
751     RLinkUp = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
752 
753     if (C && cm->status == CHOLMOD_OK)
754     {
755 	for (i = 0 ; i < n ; i++)
756 	{
757 	    mask [i] = -1 ;
758 	    RLinkUp [i] = i+1 ;
759 	}
760     }
761     else
762     {
763 	domask = FALSE ;
764     }
765 
766     if (C && !(C->packed) && !(C->sorted))
767     {
768 	/* do not do the unpacked or unsorted cases */
769 	domask = FALSE ;
770     }
771 
772     /* make a copy of C and add some gunk if stype > 0 */
773     added_gunk = (C && C->stype > 0) ;
774     if (added_gunk)
775     {
776 	C2 = add_gunk (C) ;
777     }
778     else
779     {
780 	C2 = CHOLMOD(copy_sparse) (C, cm) ;
781     }
782 
783     CC = CHOLMOD(copy_sparse) (C2, cm) ;
784 
785     if (CC && domask)
786     {
787 	Int *Cp, *Ci, p ;
788 	double *Cx, *Cz ;
789 
790 	/* this implicitly sets the first row/col of C to zero, except diag. */
791 	mask [0] = 1 ;
792 
793 	/* CC = C2, and then set the first row/col to zero, except diagonal */
794 	Cp = CC->p ;
795 	Ci = CC->i ;
796 	Cx = CC->x ;
797 	Cz = CC->z ;
798 	nz = Cp [n] ;
799 	switch (C->xtype)
800 	{
801 	    case CHOLMOD_REAL:
802 		for (p = 1 ; p < nz ; p++)
803 		{
804 		    if (Ci [p] == 0) Cx [p] = 0 ;
805 		}
806 		break ;
807 
808 	    case CHOLMOD_COMPLEX:
809 		for (p = 1 ; p < nz ; p++)
810 		{
811 		    if (Ci [p] == 0)
812 		    {
813 			Cx [2*p  ] = 0 ;
814 			Cx [2*p+1] = 0 ;
815 		    }
816 		}
817 		break ;
818 
819 	    case CHOLMOD_ZOMPLEX:
820 		for (p = 1 ; p < nz ; p++)
821 		{
822 		    if (Ci [p] == 0)
823 		    {
824 			Cx [p] = 0 ;
825 			Cz [p] = 0 ;
826 		    }
827 		}
828 		break ;
829 	}
830     }
831 
832     B = rhs (CC, 1, n) ;
833 
834     for (sorted = 1 ; sorted >= 0 ; sorted--)
835     {
836 
837 	if (!sorted)
838 	{
839 	    if (C2 && !added_gunk) C2->sorted = FALSE ;
840 	    if (C) C->sorted = FALSE ;
841 	    if (CT) CT->sorted = FALSE ;
842 	}
843 
844 	for (is_ll = 0 ; is_ll <= 1 ; is_ll++)
845 	{
846 	    for (xtype = 0 ; xtype <= 1 ; xtype++)
847 	    {
848 
849 		L = CHOLMOD(allocate_factor) (n, cm) ;
850 		if (L) L->is_ll = is_ll ;
851 
852 		if (xtype)
853 		{
854 		    CHOLMOD (change_factor) (axtype, is_ll, 0, 0, 1, L, cm) ;
855 		}
856 
857 		CHOLMOD(rowfac_mask) (sorted ? C : C2,
858 		    CT, beta, 0, n, mask, RLinkUp, L, cm) ;
859 
860 		cm->prefer_zomplex = prefer_zomplex ;
861 		X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
862 		cm->prefer_zomplex = FALSE ;
863 
864 		r = resid (CC, X, B) ;
865 		MAXERR (maxerr, r, 1) ;
866 
867 		printf ("rowfac mask: resid is %g\n", r) ;
868 
869 		CHOLMOD(free_factor) (&L, cm) ;
870 		CHOLMOD(free_dense) (&X, cm) ;
871 	    }
872 	}
873     }
874 
875     CHOLMOD(free) (n, sizeof (Int), mask, cm) ;
876     CHOLMOD(free) (n, sizeof (Int), RLinkUp, cm) ;
877 
878     CHOLMOD(free_sparse) (&C2, cm) ;
879     CHOLMOD(free_sparse) (&CC, cm) ;
880     CHOLMOD(free_sparse) (&CT, cm) ;
881     CHOLMOD(free_sparse) (&C, cm) ;
882     CHOLMOD(free_dense) (&B, cm) ;
883 
884     return (maxerr) ;
885 }
886