1 /* ========================================================================== */
2 /* === Core/cholmod_change_factor =========================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Core Module. Copyright (C) 2005-2006,
7 * Univ. of Florida. Author: Timothy A. Davis
8 * -------------------------------------------------------------------------- */
9
10 /* Change the numeric/symbolic, LL/LDL, simplicial/super, packed/unpacked,
11 * monotonic/non-monotonic status of a cholmod_factor object.
12 *
13 * There are four basic classes of factor types:
14 *
15 * (1) simplicial symbolic: Consists of two size-n arrays: the fill-reducing
16 * permutation (L->Perm) and the nonzero count for each column of L
17 * (L->ColCount). All other factor types also include this information.
18 * L->ColCount may be exact (obtained from the analysis routines), or
19 * it may be a guess. During factorization, and certainly after update/
20 * downdate, the columns of L can have a different number of nonzeros.
21 * L->ColCount is used to allocate space. L->ColCount is exact for the
22 * supernodal factorizations. The nonzero pattern of L is not kept.
23 *
24 * (2) simplicial numeric: These represent L in a compressed column form. The
25 * variants of this type are:
26 *
27 * LDL': L is unit diagonal. Row indices in column j are located in
28 * L->i [L->p [j] ... L->p [j] + L->nz [j]], and corresponding numeric
29 * values are in the same locations in L->x. The total number of
30 * entries is the sum of L->nz [j]. The unit diagonal is not stored;
31 * D is stored on the diagonal of L instead. L->p may or may not be
32 * monotonic. The order of storage of the columns in L->i and L->x is
33 * given by a doubly-linked list (L->prev and L->next). L->p is of
34 * size n+1, but only the first n entries are used (it is used if L
35 * is converted to a sparse matrix via cholmod_factor_to_sparse).
36 *
37 * For the complex case, L->x is stored interleaved with real/imag
38 * parts, and is of size 2*lnz*sizeof(double). For the zomplex case,
39 * L->x is of size lnz*sizeof(double) and holds the real part; L->z
40 * is the same size and holds the imaginary part.
41 *
42 * LL': This is identical to the LDL' form, except that the non-unit
43 * diagonal of L is stored as the first entry in each column of L.
44 *
45 * (3) supernodal symbolic: A representation of the nonzero pattern of the
46 * supernodes for a supernodal factorization. There are L->nsuper
47 * supernodes. Columns L->super [k] to L->super [k+1]-1 are in the kth
48 * supernode. The row indices for the kth supernode are in
49 * L->s [L->pi [k] ... L->pi [k+1]-1]. The numerical values are not
50 * allocated (L->x), but when they are they will be located in
51 * L->x [L->px [k] ... L->px [k+1]-1], and the L->px array is defined
52 * in this factor type.
53 *
54 * For the complex case, L->x is stored interleaved with real/imag parts,
55 * and is of size 2*L->xsize*sizeof(double). The zomplex supernodal case
56 * is not supported, since it is not compatible with LAPACK and the BLAS.
57 *
58 * (4) supernodal numeric: Always an LL' factorization. L is non-unit
59 * diagonal. L->x contains the numerical values of the supernodes, as
60 * described above for the supernodal symbolic factor.
61 * For the complex case, L->x is stored interleaved, and is of size
62 * 2*L->xsize*sizeof(double). The zomplex supernodal case is not
63 * supported, since it is not compatible with LAPACK and the BLAS.
64 *
65 * FUTURE WORK: support a supernodal LDL' factor.
66 *
67 *
68 * In all cases, the row indices in each column (L->i for simplicial L and
69 * L->s for supernodal L) are kept sorted from low indices to high indices.
70 * This means the diagonal of L (or D for LDL' factors) is always kept as the
71 * first entry in each column.
72 *
73 * The cholmod_change_factor routine can do almost all possible conversions.
74 * It cannot do the following conversions:
75 *
76 * (1) Simplicial numeric types cannot be converted to a supernodal
77 * symbolic type. This would simultaneously deallocate the
78 * simplicial pattern and numeric values and reallocate uninitialized
79 * space for the supernodal pattern. This isn't useful for the user,
80 * and not needed by CHOLMOD's own routines either.
81 *
82 * (2) Only a symbolic factor (simplicial to supernodal) can be converted
83 * to a supernodal numeric factor.
84 *
85 * Some conversions are meant only to be used internally by other CHOLMOD
86 * routines, and should not be performed by the end user. They allocate space
87 * whose contents are undefined:
88 *
89 * (1) converting from simplicial symbolic to supernodal symbolic.
90 * (2) converting any factor to supernodal numeric.
91 *
92 * workspace: no conversion routine uses workspace in Common. No temporary
93 * workspace is allocated.
94 *
95 * Supports all xtypes, except that there is no supernodal zomplex L.
96 *
97 * The to_xtype parameter is used only when converting from symbolic to numeric
98 * or numeric to symbolic. It cannot be used to convert a numeric xtype (real,
99 * complex, or zomplex) to a different numeric xtype. For that conversion,
100 * use cholmod_factor_xtype instead.
101 */
102
103 #include "cholmod_internal.h"
104 #include "cholmod_core.h"
105
106 static void natural_list (cholmod_factor *L) ;
107
108 /* ========================================================================== */
109 /* === TEMPLATE ============================================================= */
110 /* ========================================================================== */
111
112 #define REAL
113 #include "t_cholmod_change_factor.c"
114 #define COMPLEX
115 #include "t_cholmod_change_factor.c"
116 #define ZOMPLEX
117 #include "t_cholmod_change_factor.c"
118
119
120 /* ========================================================================== */
121 /* === L_is_packed ========================================================== */
122 /* ========================================================================== */
123
124 /* Return TRUE if the columns of L are packed, FALSE otherwise. For debugging
125 * only. */
126
127 #ifndef NDEBUG
L_is_packed(cholmod_factor * L,cholmod_common * Common)128 static int L_is_packed (cholmod_factor *L, cholmod_common *Common)
129 {
130 Int j ;
131 Int *Lnz = L->nz ;
132 Int *Lp = L->p ;
133 Int n = L->n ;
134
135 if (L->xtype == CHOLMOD_PATTERN || L->is_super)
136 {
137 return (TRUE) ;
138 }
139
140 if (Lnz == NULL || Lp == NULL)
141 {
142 return (TRUE) ;
143 }
144
145 for (j = 0 ; j < n ; j++)
146 {
147 PRINT3 (("j: "ID" Lnz "ID" Lp[j+1] "ID" Lp[j] "ID"\n", j, Lnz [j],
148 Lp [j+1], Lp [j])) ;
149 if (Lnz [j] != (Lp [j+1] - Lp [j]))
150 {
151 PRINT2 (("L is not packed\n")) ;
152 return (FALSE) ;
153 }
154 }
155 return (TRUE) ;
156 }
157 #endif
158
159
160 /* ========================================================================== */
161 /* === natural_list ========================================================= */
162 /* ========================================================================== */
163
164 /* Create a naturally-ordered doubly-linked list of columns. */
165
natural_list(cholmod_factor * L)166 static void natural_list (cholmod_factor *L)
167 {
168 Int head, tail, n, j ;
169 Int *Lnext, *Lprev ;
170 Lnext = L->next ;
171 Lprev = L->prev ;
172 ASSERT (Lprev != NULL && Lnext != NULL) ;
173 n = L->n ;
174 head = n+1 ;
175 tail = n ;
176 Lnext [head] = 0 ;
177 Lprev [head] = EMPTY ;
178 Lnext [tail] = EMPTY ;
179 Lprev [tail] = n-1 ;
180 for (j = 0 ; j < n ; j++)
181 {
182 Lnext [j] = j+1 ;
183 Lprev [j] = j-1 ;
184 }
185 Lprev [0] = head ;
186 L->is_monotonic = TRUE ;
187 }
188
189
190 /* ========================================================================== */
191 /* === allocate_simplicial_numeric ========================================== */
192 /* ========================================================================== */
193
194 /* Allocate O(n) arrays for simplicial numeric factorization. Initializes
195 * the link lists only. Does not allocate the L->i, L->x, or L->z arrays. */
196
allocate_simplicial_numeric(cholmod_factor * L,cholmod_common * Common)197 static int allocate_simplicial_numeric
198 (
199 cholmod_factor *L,
200 cholmod_common *Common
201 )
202 {
203 Int n ;
204 Int *Lp, *Lnz, *Lprev, *Lnext ;
205 size_t n1, n2 ;
206
207 PRINT1 (("Allocate simplicial\n")) ;
208
209 ASSERT (L->xtype == CHOLMOD_PATTERN || L->is_super) ;
210 ASSERT (L->p == NULL) ;
211 ASSERT (L->nz == NULL) ;
212 ASSERT (L->prev == NULL) ;
213 ASSERT (L->next == NULL) ;
214
215 n = L->n ;
216
217 /* this cannot cause size_t overflow */
218 n1 = ((size_t) n) + 1 ;
219 n2 = ((size_t) n) + 2 ;
220
221 Lp = CHOLMOD(malloc) (n1, sizeof (Int), Common) ;
222 Lnz = CHOLMOD(malloc) (n, sizeof (Int), Common) ;
223 Lprev = CHOLMOD(malloc) (n2, sizeof (Int), Common) ;
224 Lnext = CHOLMOD(malloc) (n2, sizeof (Int), Common) ;
225
226 if (Common->status < CHOLMOD_OK)
227 {
228 CHOLMOD(free) (n1, sizeof (Int), Lp, Common) ;
229 CHOLMOD(free) (n, sizeof (Int), Lnz, Common) ;
230 CHOLMOD(free) (n2, sizeof (Int), Lprev, Common) ;
231 CHOLMOD(free) (n2, sizeof (Int), Lnext, Common) ;
232 PRINT1 (("Allocate simplicial failed\n")) ;
233 return (FALSE) ; /* out of memory */
234 }
235
236 /* ============================================== commit the changes to L */
237
238 L->p = Lp ;
239 L->nz = Lnz ;
240 L->prev = Lprev ;
241 L->next = Lnext ;
242 /* initialize a doubly linked list for columns in natural order */
243 natural_list (L) ;
244 PRINT1 (("Allocate simplicial done\n")) ;
245 return (TRUE) ;
246 }
247
248
249 /* ========================================================================== */
250 /* === simplicial_symbolic_to_super_symbolic ================================ */
251 /* ========================================================================== */
252
253 /* Convert a simplicial symbolic factor supernodal symbolic factor. Does not
254 * initialize the new space. */
255
simplicial_symbolic_to_super_symbolic(cholmod_factor * L,cholmod_common * Common)256 static int simplicial_symbolic_to_super_symbolic
257 (
258 cholmod_factor *L,
259 cholmod_common *Common
260 )
261 {
262 Int nsuper, xsize, ssize ;
263 Int *Lsuper, *Lpi, *Lpx, *Ls ;
264 size_t nsuper1 ;
265
266 ASSERT (L->xtype == CHOLMOD_PATTERN && !(L->is_super)) ;
267
268 xsize = L->xsize ;
269 ssize = L->ssize ;
270 nsuper = L->nsuper ;
271 nsuper1 = ((size_t) nsuper) + 1 ;
272
273 PRINT1 (("simple sym to super sym: ssize "ID" xsize "ID" nsuper "ID""
274 " status %d\n", ssize, xsize, nsuper, Common->status)) ;
275
276 /* O(nsuper) arrays, where nsuper <= n */
277 Lsuper = CHOLMOD(malloc) (nsuper1, sizeof (Int), Common) ;
278 Lpi = CHOLMOD(malloc) (nsuper1, sizeof (Int), Common) ;
279 Lpx = CHOLMOD(malloc) (nsuper1, sizeof (Int), Common) ;
280
281 /* O(ssize) array, where ssize <= nnz(L), and usually much smaller */
282 Ls = CHOLMOD(malloc) (ssize, sizeof (Int), Common) ;
283
284 if (Common->status < CHOLMOD_OK)
285 {
286 CHOLMOD(free) (nsuper1, sizeof (Int), Lsuper, Common) ;
287 CHOLMOD(free) (nsuper1, sizeof (Int), Lpi, Common) ;
288 CHOLMOD(free) (nsuper1, sizeof (Int), Lpx, Common) ;
289 CHOLMOD(free) (ssize, sizeof (Int), Ls, Common) ;
290 return (FALSE) ; /* out of memory */
291 }
292
293 /* ============================================== commit the changes to L */
294
295 ASSERT (Lsuper != NULL && Lpi != NULL && Lpx != NULL && Ls != NULL) ;
296
297 L->maxcsize = 0 ;
298 L->maxesize = 0 ;
299
300 L->super = Lsuper ;
301 L->pi = Lpi ;
302 L->px = Lpx ;
303 L->s = Ls ;
304 Ls [0] = EMPTY ; /* supernodal pattern undefined */
305
306 L->is_super = TRUE ;
307 L->is_ll = TRUE ; /* supernodal LDL' not supported */
308 L->xtype = CHOLMOD_PATTERN ;
309 L->dtype = DTYPE ;
310 L->minor = L->n ;
311 return (TRUE) ;
312 }
313
314
315 /* ========================================================================== */
316 /* === any_to_simplicial_symbolic =========================================== */
317 /* ========================================================================== */
318
319 /* Convert any factor L to a simplicial symbolic factor, leaving only L->Perm
320 * and L->ColCount. Cannot fail. Any of the components of L (except Perm and
321 * ColCount) may already be free'd. */
322
any_to_simplicial_symbolic(cholmod_factor * L,int to_ll,cholmod_common * Common)323 static void any_to_simplicial_symbolic
324 (
325 cholmod_factor *L,
326 int to_ll,
327 cholmod_common *Common
328 )
329 {
330 Int n, lnz, xs, ss, s, e ;
331 size_t n1, n2 ;
332
333 /* ============================================== commit the changes to L */
334
335 n = L->n ;
336 lnz = L->nzmax ;
337 s = L->nsuper + 1 ;
338 xs = (L->is_super) ? ((Int) (L->xsize)) : (lnz) ;
339 e = (L->xtype == CHOLMOD_COMPLEX ? 2 : 1) ;
340 ss = L->ssize ;
341
342 /* this cannot cause size_t overflow */
343 n1 = ((size_t) n) + 1 ;
344 n2 = ((size_t) n) + 2 ;
345
346 /* free all but the symbolic analysis (Perm and ColCount) */
347 L->p = CHOLMOD(free) (n1, sizeof (Int), L->p, Common) ;
348 L->i = CHOLMOD(free) (lnz, sizeof (Int), L->i, Common) ;
349 L->x = CHOLMOD(free) (xs, e*sizeof (double), L->x, Common) ;
350 L->z = CHOLMOD(free) (lnz, sizeof (double), L->z, Common) ;
351 L->nz = CHOLMOD(free) (n, sizeof (Int), L->nz, Common) ;
352 L->next = CHOLMOD(free) (n2, sizeof (Int), L->next, Common) ;
353 L->prev = CHOLMOD(free) (n2, sizeof (Int), L->prev, Common) ;
354 L->super = CHOLMOD(free) (s, sizeof (Int), L->super, Common) ;
355 L->pi = CHOLMOD(free) (s, sizeof (Int), L->pi, Common) ;
356 L->px = CHOLMOD(free) (s, sizeof (Int), L->px, Common) ;
357 L->s = CHOLMOD(free) (ss, sizeof (Int), L->s, Common) ;
358 L->nzmax = 0 ;
359 L->is_super = FALSE ;
360 L->xtype = CHOLMOD_PATTERN ;
361 L->dtype = DTYPE ;
362 L->minor = n ;
363 L->is_ll = to_ll ;
364 }
365
366
367 /* ========================================================================== */
368 /* === ll_super_to_super_symbolic =========================================== */
369 /* ========================================================================== */
370
371 /* Convert a numerical supernodal L to symbolic supernodal. Cannot fail. */
372
ll_super_to_super_symbolic(cholmod_factor * L,cholmod_common * Common)373 static void ll_super_to_super_symbolic
374 (
375 cholmod_factor *L,
376 cholmod_common *Common
377 )
378 {
379
380 /* ============================================== commit the changes to L */
381
382 /* free all but the supernodal numerical factor */
383 ASSERT (L->xtype != CHOLMOD_PATTERN && L->is_super && L->is_ll) ;
384 DEBUG (CHOLMOD(dump_factor) (L, "start to super symbolic", Common)) ;
385 L->x = CHOLMOD(free) (L->xsize,
386 (L->xtype == CHOLMOD_COMPLEX ? 2 : 1) * sizeof (double), L->x,
387 Common) ;
388 L->xtype = CHOLMOD_PATTERN ;
389 L->dtype = DTYPE ;
390 L->minor = L->n ;
391 L->is_ll = TRUE ; /* supernodal LDL' not supported */
392 DEBUG (CHOLMOD(dump_factor) (L, "done to super symbolic", Common)) ;
393 }
394
395
396 /* ========================================================================== */
397 /* === simplicial_symbolic_to_simplicial_numeric ============================ */
398 /* ========================================================================== */
399
400 /* Convert a simplicial symbolic L to a simplicial numeric L; allocate space
401 * for L using L->ColCount from symbolic analysis, and set L to identity.
402 *
403 * If packed < 0, then this routine is creating a copy of another factor
404 * (via cholmod_copy_factor). In this case, the space is not initialized. */
405
simplicial_symbolic_to_simplicial_numeric(cholmod_factor * L,int to_ll,int packed,int to_xtype,cholmod_common * Common)406 static void simplicial_symbolic_to_simplicial_numeric
407 (
408 cholmod_factor *L,
409 int to_ll,
410 int packed,
411 int to_xtype,
412 cholmod_common *Common
413 )
414 {
415 double grow0, grow1, xlen, xlnz ;
416 double *Lx, *Lz ;
417 Int *Li, *Lp, *Lnz, *ColCount ;
418 Int n, grow, grow2, p, j, lnz, len, ok, e ;
419
420 ASSERT (L->xtype == CHOLMOD_PATTERN && !(L->is_super)) ;
421 if (!allocate_simplicial_numeric (L, Common))
422 {
423 PRINT1 (("out of memory, allocate simplicial numeric\n")) ;
424 return ; /* out of memory */
425 }
426 ASSERT (L->ColCount != NULL && L->nz != NULL && L->p != NULL) ;
427 ASSERT (L->x == NULL && L->z == NULL && L->i == NULL) ;
428
429 ColCount = L->ColCount ;
430 Lnz = L->nz ;
431 Lp = L->p ;
432 ok = TRUE ;
433 n = L->n ;
434
435 if (packed < 0)
436 {
437
438 /* ------------------------------------------------------------------ */
439 /* used by cholmod_copy_factor to allocate a copy of a factor object */
440 /* ------------------------------------------------------------------ */
441
442 lnz = L->nzmax ;
443 L->nzmax = 0 ;
444
445 }
446 else if (packed)
447 {
448
449 /* ------------------------------------------------------------------ */
450 /* LDL' or LL' packed */
451 /* ------------------------------------------------------------------ */
452
453 PRINT1 (("convert to packed LL' or LDL'\n")) ;
454 lnz = 0 ;
455 for (j = 0 ; ok && j < n ; j++)
456 {
457 /* ensure len is in the range 1 to n-j */
458 len = ColCount [j] ;
459 len = MAX (1, len) ;
460 len = MIN (len, n-j) ;
461 lnz += len ;
462 ok = (lnz >= 0) ;
463 }
464 for (j = 0 ; j <= n ; j++)
465 {
466 Lp [j] = j ;
467 }
468 for (j = 0 ; j < n ; j++)
469 {
470 Lnz [j] = 1 ;
471 }
472
473 }
474 else
475 {
476
477 /* ------------------------------------------------------------------ */
478 /* LDL' unpacked */
479 /* ------------------------------------------------------------------ */
480
481 PRINT1 (("convert to unpacked\n")) ;
482 /* compute new lnzmax */
483 /* if any parameter is NaN, grow is false */
484 grow0 = Common->grow0 ;
485 grow1 = Common->grow1 ;
486 grow2 = Common->grow2 ;
487 grow0 = IS_NAN (grow0) ? 1 : grow0 ;
488 grow1 = IS_NAN (grow1) ? 1 : grow1 ;
489 /* fl.pt. compare, but no NaN's: */
490 grow = (grow0 >= 1.0) && (grow1 >= 1.0) && (grow2 > 0) ;
491 PRINT1 (("init, grow1 %g grow2 "ID"\n", grow1, grow2)) ;
492 /* initialize Lp and Lnz for each column */
493 lnz = 0 ;
494 for (j = 0 ; ok && j < n ; j++)
495 {
496 Lp [j] = lnz ;
497 Lnz [j] = 1 ;
498
499 /* ensure len is in the range 1 to n-j */
500 len = ColCount [j] ;
501 len = MAX (1, len) ;
502 len = MIN (len, n-j) ;
503
504 /* compute len in double to avoid integer overflow */
505 PRINT1 (("ColCount ["ID"] = "ID"\n", j, len)) ;
506 if (grow)
507 {
508 xlen = (double) len ;
509 xlen = grow1 * xlen + grow2 ;
510 xlen = MIN (xlen, n-j) ;
511 len = (Int) xlen ;
512 }
513 ASSERT (len >= 1 && len <= n-j) ;
514 lnz += len ;
515 ok = (lnz >= 0) ;
516 }
517 if (ok)
518 {
519 Lp [n] = lnz ;
520 if (grow)
521 {
522 /* add extra space */
523 xlnz = (double) lnz ;
524 xlnz *= grow0 ;
525 xlnz = MIN (xlnz, Size_max) ;
526 xlnz = MIN (xlnz, ((double) n * (double) n + (double) n) / 2) ;
527 lnz = (Int) xlnz ;
528 }
529 }
530 }
531
532 lnz = MAX (1, lnz) ;
533
534 if (!ok)
535 {
536 ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
537 }
538
539 /* allocate L->i, L->x, and L->z */
540 PRINT1 (("resizing from zero size to lnz "ID"\n", lnz)) ;
541 ASSERT (L->nzmax == 0) ;
542 e = (to_xtype == CHOLMOD_COMPLEX ? 2 : 1) ;
543 if (!ok || !CHOLMOD(realloc_multiple) (lnz, 1, to_xtype, &(L->i), NULL,
544 &(L->x), &(L->z), &(L->nzmax), Common))
545 {
546 L->p = CHOLMOD(free) (n+1, sizeof (Int), L->p, Common) ;
547 L->nz = CHOLMOD(free) (n, sizeof (Int), L->nz, Common) ;
548 L->prev = CHOLMOD(free) (n+2, sizeof (Int), L->prev, Common) ;
549 L->next = CHOLMOD(free) (n+2, sizeof (Int), L->next, Common) ;
550 L->i = CHOLMOD(free) (lnz, sizeof (Int), L->i, Common) ;
551 L->x = CHOLMOD(free) (lnz, e*sizeof (double), L->x, Common) ;
552 L->z = CHOLMOD(free) (lnz, sizeof (double), L->z, Common) ;
553 PRINT1 (("cannot realloc simplicial numeric\n")) ;
554 return ; /* out of memory */
555 }
556
557 /* ============================================== commit the changes to L */
558
559 /* initialize L to be the identity matrix */
560 L->xtype = to_xtype ;
561 L->dtype = DTYPE ;
562 L->minor = n ;
563
564 Li = L->i ;
565 Lx = L->x ;
566 Lz = L->z ;
567
568 #if 0
569 if (lnz == 1)
570 {
571 /* the user won't expect to access this entry, but some CHOLMOD
572 * routines may. Set it to zero so that valgrind doesn't complain. */
573 switch (to_xtype)
574 {
575 case CHOLMOD_REAL:
576 Lx [0] = 0 ;
577 break ;
578
579 case CHOLMOD_COMPLEX:
580 Lx [0] = 0 ;
581 Lx [1] = 0 ;
582 break ;
583
584 case CHOLMOD_ZOMPLEX:
585 Lx [0] = 0 ;
586 Lz [0] = 0 ;
587 break ;
588 }
589 }
590 #endif
591
592 if (packed >= 0)
593 {
594 /* create the unit diagonal for either the LL' or LDL' case */
595
596 switch (L->xtype)
597 {
598 case CHOLMOD_REAL:
599 for (j = 0 ; j < n ; j++)
600 {
601 ASSERT (Lp [j] < Lp [j+1]) ;
602 p = Lp [j] ;
603 Li [p] = j ;
604 Lx [p] = 1 ;
605 }
606 break ;
607
608 case CHOLMOD_COMPLEX:
609 for (j = 0 ; j < n ; j++)
610 {
611 ASSERT (Lp [j] < Lp [j+1]) ;
612 p = Lp [j] ;
613 Li [p] = j ;
614 Lx [2*p ] = 1 ;
615 Lx [2*p+1] = 0 ;
616 }
617 break ;
618
619 case CHOLMOD_ZOMPLEX:
620 for (j = 0 ; j < n ; j++)
621 {
622 ASSERT (Lp [j] < Lp [j+1]) ;
623 p = Lp [j] ;
624 Li [p] = j ;
625 Lx [p] = 1 ;
626 Lz [p] = 0 ;
627 }
628 break ;
629 }
630 }
631
632 L->is_ll = to_ll ;
633
634 PRINT1 (("done convert simplicial symbolic to numeric\n")) ;
635 }
636
637
638 /* ========================================================================== */
639 /* === change_simplicial_numeric ============================================ */
640 /* ========================================================================== */
641
642 /* Change LL' to LDL', LDL' to LL', or leave as-is.
643 *
644 * If to_packed is TRUE, then the columns of L are packed and made monotonic
645 * (to_monotonic is ignored; it is implicitly TRUE).
646 *
647 * If to_monotonic is TRUE but to_packed is FALSE, the columns of L are made
648 * monotonic but not packed.
649 *
650 * If both to_packed and to_monotonic are FALSE, then the columns of L are
651 * left as-is, and the conversion is done in place.
652 *
653 * If L is already monotonic, or if it is to be left non-monotonic, then this
654 * conversion always succeeds.
655 *
656 * When converting an LDL' to LL' factorization, any column with a negative
657 * or zero diagonal entry is not modified so that conversion back to LDL' will
658 * succeed. This can result in a matrix L with a negative entry on the diagonal
659 * If the kth entry on the diagonal of D is negative, it and the kth column of
660 * L are left unchanged. A subsequent conversion back to an LDL' form will also
661 * leave the column unchanged, so the correct LDL' factorization will be
662 * restored. L->minor is set to the smallest k for which D (k,k) is negative.
663 */
664
change_simplicial_numeric(cholmod_factor * L,int to_ll,int to_packed,int to_monotonic,cholmod_common * Common)665 static void change_simplicial_numeric
666 (
667 cholmod_factor *L,
668 int to_ll,
669 int to_packed,
670 int to_monotonic,
671 cholmod_common *Common
672 )
673 {
674 double grow0, grow1, xlen, xlnz ;
675 void *newLi, *newLx, *newLz ;
676 double *Lx, *Lz ;
677 Int *Lp, *Li, *Lnz ;
678 Int make_monotonic, grow2, n, j, lnz, len, grow, ok, make_ll, make_ldl ;
679 size_t nzmax0 ;
680
681 PRINT1 (("\n===Change simplicial numeric: %d %d %d\n",
682 to_ll, to_packed, to_monotonic)) ;
683 DEBUG (CHOLMOD(dump_factor) (L, "change simplicial numeric", Common)) ;
684 ASSERT (L->xtype != CHOLMOD_PATTERN && !(L->is_super)) ;
685
686 make_monotonic = ((to_packed || to_monotonic) && !(L->is_monotonic)) ;
687 make_ll = (to_ll && !(L->is_ll)) ;
688 make_ldl = (!to_ll && L->is_ll) ;
689
690 n = L->n ;
691 Lp = L->p ;
692 Li = L->i ;
693 Lx = L->x ;
694 Lz = L->z ;
695 Lnz = L->nz ;
696
697 grow = FALSE ;
698 grow0 = Common->grow0 ;
699 grow1 = Common->grow1 ;
700 grow2 = Common->grow2 ;
701 grow0 = IS_NAN (grow0) ? 1 : grow0 ;
702 grow1 = IS_NAN (grow1) ? 1 : grow1 ;
703 ok = TRUE ;
704 newLi = NULL ;
705 newLx = NULL ;
706 newLz = NULL ;
707 lnz = 0 ;
708
709 if (make_monotonic)
710 {
711
712 /* ------------------------------------------------------------------ */
713 /* Columns out of order, but will be reordered and optionally packed. */
714 /* ------------------------------------------------------------------ */
715
716 PRINT1 (("L is non-monotonic\n")) ;
717
718 /* compute new L->nzmax */
719 if (!to_packed)
720 {
721 /* if any parameter is NaN, grow is false */
722 /* fl.pt. comparisons below are false if any parameter is NaN */
723 grow = (grow0 >= 1.0) && (grow1 >= 1.0) && (grow2 > 0) ;
724 }
725 for (j = 0 ; ok && j < n ; j++)
726 {
727 len = Lnz [j] ;
728 ASSERT (len >= 1 && len <= n-j) ;
729
730 /* compute len in double to avoid integer overflow */
731 if (grow)
732 {
733 xlen = (double) len ;
734 xlen = grow1 * xlen + grow2 ;
735 xlen = MIN (xlen, n-j) ;
736 len = (Int) xlen ;
737 }
738 ASSERT (len >= Lnz [j] && len <= n-j) ;
739
740 PRINT2 (("j: "ID" Lnz[j] "ID" len "ID" p "ID"\n",
741 j, Lnz [j], len, lnz)) ;
742
743 lnz += len ;
744 ok = (lnz >= 0) ;
745 }
746
747 if (!ok)
748 {
749 ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
750 return ;
751 }
752
753 if (grow)
754 {
755 xlnz = (double) lnz ;
756 xlnz *= grow0 ;
757 xlnz = MIN (xlnz, Size_max) ;
758 xlnz = MIN (xlnz, ((double) n * (double) n + (double) n) / 2) ;
759 lnz = (Int) xlnz ;
760 }
761
762 lnz = MAX (1, lnz) ;
763 PRINT1 (("final lnz "ID"\n", lnz)) ;
764 nzmax0 = 0 ;
765
766 CHOLMOD(realloc_multiple) (lnz, 1, L->xtype, &newLi, NULL,
767 &newLx, &newLz, &nzmax0, Common) ;
768
769 if (Common->status < CHOLMOD_OK)
770 {
771 return ; /* out of memory */
772 }
773 }
774
775 /* ============================================== commit the changes to L */
776
777 /* ---------------------------------------------------------------------- */
778 /* convert the simplicial L, using template routine */
779 /* ---------------------------------------------------------------------- */
780
781 switch (L->xtype)
782 {
783
784 case CHOLMOD_REAL:
785 r_change_simplicial_numeric (L, to_ll, to_packed,
786 newLi, newLx, newLz, lnz, grow, grow1, grow2,
787 make_ll, make_monotonic, make_ldl, Common) ;
788 break ;
789
790 case CHOLMOD_COMPLEX:
791 c_change_simplicial_numeric (L, to_ll, to_packed,
792 newLi, newLx, newLz, lnz, grow, grow1, grow2,
793 make_ll, make_monotonic, make_ldl, Common) ;
794 break ;
795
796 case CHOLMOD_ZOMPLEX:
797 z_change_simplicial_numeric (L, to_ll, to_packed,
798 newLi, newLx, newLz, lnz, grow, grow1, grow2,
799 make_ll, make_monotonic, make_ldl, Common) ;
800 break ;
801 }
802
803 DEBUG (CHOLMOD(dump_factor) (L, "L simplicial changed", Common)) ;
804 }
805
806
807 /* ========================================================================== */
808 /* === ll_super_to_simplicial_numeric ======================================= */
809 /* ========================================================================== */
810
811 /* Convert a supernodal numeric factorization to any simplicial numeric one.
812 * Leaves L->xtype unchanged (real or complex, not zomplex since there is
813 * no supernodal zomplex L). */
814
ll_super_to_simplicial_numeric(cholmod_factor * L,int to_packed,int to_ll,cholmod_common * Common)815 static void ll_super_to_simplicial_numeric
816 (
817 cholmod_factor *L,
818 int to_packed,
819 int to_ll,
820 cholmod_common *Common
821 )
822 {
823 Int *Ls, *Lpi, *Lpx, *Super, *Li ;
824 Int n, lnz, s, nsuper, psi, psend, nsrow, nscol, k1, k2, erows ;
825
826 DEBUG (CHOLMOD(dump_factor) (L, "start LL super to simplicial", Common)) ;
827 PRINT1 (("super -> simplicial (%d %d)\n", to_packed, to_ll)) ;
828 ASSERT (L->xtype != CHOLMOD_PATTERN && L->is_ll && L->is_super) ;
829 ASSERT (L->x != NULL && L->i == NULL) ;
830
831 n = L->n ;
832 nsuper = L->nsuper ;
833 Lpi = L->pi ;
834 Lpx = L->px ;
835 Ls = L->s ;
836 Super = L->super ;
837
838 /* Int overflow cannot occur since supernodal L already exists */
839
840 if (to_packed)
841 {
842 /* count the number of nonzeros in L. Each supernode is of the form
843 *
844 * l . . . For this example, nscol = 4 (# columns). nsrow = 9.
845 * l l . . The "." entries are allocated in the supernodal
846 * l l l . factor, but not used. They are not copied to the
847 * l l l l simplicial factor. Some "l" and "e" entries may be
848 * e e e e numerically zero and even symbolically zero if a
849 * e e e e tight simplicial factorization or resymbol were
850 * e e e e done, because of numerical cancellation and relaxed
851 * e e e e supernode amalgamation, respectively.
852 * e e e e
853 */
854 lnz = 0 ;
855 for (s = 0 ; s < nsuper ; s++)
856 {
857 k1 = Super [s] ;
858 k2 = Super [s+1] ;
859 psi = Lpi [s] ;
860 psend = Lpi [s+1] ;
861 nsrow = psend - psi ;
862 nscol = k2 - k1 ;
863 ASSERT (nsrow >= nscol) ;
864 erows = nsrow - nscol ;
865
866 /* lower triangular part, including the diagonal,
867 * counting the "l" terms in the figure above. */
868 lnz += nscol * (nscol+1) / 2 ;
869
870 /* rectangular part, below the diagonal block (the "e" terms) */
871 lnz += nscol * erows ;
872 }
873 ASSERT (lnz <= (Int) (L->xsize)) ;
874 }
875 else
876 {
877 /* Li will be the same size as Lx */
878 lnz = L->xsize ;
879 }
880 ASSERT (lnz >= 0) ;
881 PRINT1 (("simplicial lnz = "ID" to_packed: %d to_ll: %d L->xsize %g\n",
882 lnz, to_ll, to_packed, (double) L->xsize)) ;
883
884 Li = CHOLMOD(malloc) (lnz, sizeof (Int), Common) ;
885 if (Common->status < CHOLMOD_OK)
886 {
887 return ; /* out of memory */
888 }
889
890 if (!allocate_simplicial_numeric (L, Common))
891 {
892 CHOLMOD(free) (lnz, sizeof (Int), Li, Common) ;
893 return ; /* out of memory */
894 }
895
896 /* ============================================== commit the changes to L */
897
898 L->i = Li ;
899 L->nzmax = lnz ;
900
901 /* ---------------------------------------------------------------------- */
902 /* convert the supernodal L, using template routine */
903 /* ---------------------------------------------------------------------- */
904
905 switch (L->xtype)
906 {
907
908 case CHOLMOD_REAL:
909 r_ll_super_to_simplicial_numeric (L, to_packed, to_ll, Common) ;
910 break ;
911
912 case CHOLMOD_COMPLEX:
913 c_ll_super_to_simplicial_numeric (L, to_packed, to_ll, Common) ;
914 break ;
915 }
916
917 /* ---------------------------------------------------------------------- */
918 /* free unused parts of L */
919 /* ---------------------------------------------------------------------- */
920
921 L->super = CHOLMOD(free) (nsuper+1, sizeof (Int), L->super, Common) ;
922 L->pi = CHOLMOD(free) (nsuper+1, sizeof (Int), L->pi, Common) ;
923 L->px = CHOLMOD(free) (nsuper+1, sizeof (Int), L->px, Common) ;
924 L->s = CHOLMOD(free) (L->ssize, sizeof (Int), L->s, Common) ;
925
926 L->ssize = 0 ;
927 L->xsize = 0 ;
928 L->nsuper = 0 ;
929 L->maxesize = 0 ;
930 L->maxcsize = 0 ;
931
932 L->is_super = FALSE ;
933
934 DEBUG (CHOLMOD(dump_factor) (L, "done LL super to simplicial", Common)) ;
935 }
936
937
938 /* ========================================================================== */
939 /* === super_symbolic_to_ll_super =========================================== */
940 /* ========================================================================== */
941
942 /* Convert a supernodal symbolic factorization to a supernodal numeric
943 * factorization by allocating L->x. Contents of L->x are undefined.
944 */
945
super_symbolic_to_ll_super(int to_xtype,cholmod_factor * L,cholmod_common * Common)946 static int super_symbolic_to_ll_super
947 (
948 int to_xtype,
949 cholmod_factor *L,
950 cholmod_common *Common
951 )
952 {
953 double *Lx ;
954 Int wentry = (to_xtype == CHOLMOD_REAL) ? 1 : 2 ;
955 PRINT1 (("convert super sym to num\n")) ;
956 ASSERT (L->xtype == CHOLMOD_PATTERN && L->is_super) ;
957 Lx = CHOLMOD(malloc) (L->xsize, wentry * sizeof (double), Common) ;
958 PRINT1 (("xsize %g\n", (double) L->xsize)) ;
959 if (Common->status < CHOLMOD_OK)
960 {
961 return (FALSE) ; /* out of memory */
962 }
963
964 /* ============================================== commit the changes to L */
965
966 if (L->xsize == 1)
967 {
968 /* the caller won't expect to access this entry, but some CHOLMOD
969 * routines may. Set it to zero so that valgrind doesn't complain. */
970 switch (to_xtype)
971 {
972 case CHOLMOD_REAL:
973 Lx [0] = 0 ;
974 break ;
975
976 case CHOLMOD_COMPLEX:
977 Lx [0] = 0 ;
978 Lx [1] = 0 ;
979 break ;
980 }
981 }
982
983 L->x = Lx ;
984 L->xtype = to_xtype ;
985 L->dtype = DTYPE ;
986 L->minor = L->n ;
987 return (TRUE) ;
988 }
989
990
991 /* ========================================================================== */
992 /* === cholmod_change_factor ================================================ */
993 /* ========================================================================== */
994
995 /* Convert a factor L. Some conversions simply allocate uninitialized space
996 * that meant to be filled later.
997 *
998 * If the conversion fails, the factor is left in its original form, with one
999 * exception. Converting a supernodal symbolic factor to a simplicial numeric
1000 * one (with L=D=I) may leave the factor in simplicial symbolic form.
1001 *
1002 * Memory allocated for each conversion is listed below.
1003 */
1004
CHOLMOD(change_factor)1005 int CHOLMOD(change_factor)
1006 (
1007 /* ---- input ---- */
1008 int to_xtype, /* convert to CHOLMOD_PATTERN, _REAL, _COMPLEX, or
1009 * _ZOMPLEX */
1010 int to_ll, /* TRUE: convert to LL', FALSE: LDL' */
1011 int to_super, /* TRUE: convert to supernodal, FALSE: simplicial */
1012 int to_packed, /* TRUE: pack simplicial columns, FALSE: do not pack */
1013 int to_monotonic, /* TRUE: put simplicial columns in order, FALSE: not */
1014 /* ---- in/out --- */
1015 cholmod_factor *L, /* factor to modify */
1016 /* --------------- */
1017 cholmod_common *Common
1018 )
1019 {
1020
1021 /* ---------------------------------------------------------------------- */
1022 /* get inputs */
1023 /* ---------------------------------------------------------------------- */
1024
1025 RETURN_IF_NULL_COMMON (FALSE) ;
1026 RETURN_IF_NULL (L, FALSE) ;
1027 RETURN_IF_XTYPE_INVALID (L, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
1028 if (to_xtype < CHOLMOD_PATTERN || to_xtype > CHOLMOD_ZOMPLEX)
1029 {
1030 ERROR (CHOLMOD_INVALID, "xtype invalid") ;
1031 return (FALSE) ;
1032 }
1033 Common->status = CHOLMOD_OK ;
1034
1035 PRINT1 (("-----convert from (%d,%d,%d,%d,%d) to (%d,%d,%d,%d,%d)\n",
1036 L->xtype, L->is_ll, L->is_super, L_is_packed (L, Common), L->is_monotonic,
1037 to_xtype, to_ll, to_super, to_packed, to_monotonic)) ;
1038
1039 /* ensure all parameters are TRUE/FALSE */
1040 to_ll = BOOLEAN (to_ll) ;
1041 to_super = BOOLEAN (to_super) ;
1042
1043 ASSERT (BOOLEAN (L->is_ll) == L->is_ll) ;
1044 ASSERT (BOOLEAN (L->is_super) == L->is_super) ;
1045
1046 if (to_super && to_xtype == CHOLMOD_ZOMPLEX)
1047 {
1048 ERROR (CHOLMOD_INVALID, "supernodal zomplex L not supported") ;
1049 return (FALSE) ;
1050 }
1051
1052 /* ---------------------------------------------------------------------- */
1053 /* convert */
1054 /* ---------------------------------------------------------------------- */
1055
1056 if (to_xtype == CHOLMOD_PATTERN)
1057 {
1058
1059 /* ------------------------------------------------------------------ */
1060 /* convert to symbolic */
1061 /* ------------------------------------------------------------------ */
1062
1063 if (!to_super)
1064 {
1065
1066 /* -------------------------------------------------------------- */
1067 /* convert any factor into a simplicial symbolic factor */
1068 /* -------------------------------------------------------------- */
1069
1070 any_to_simplicial_symbolic (L, to_ll, Common) ; /* cannot fail */
1071
1072 }
1073 else
1074 {
1075
1076 /* -------------------------------------------------------------- */
1077 /* convert to a supernodal symbolic factor */
1078 /* -------------------------------------------------------------- */
1079
1080 if (L->xtype != CHOLMOD_PATTERN && L->is_super)
1081 {
1082 /* convert from supernodal numeric to supernodal symbolic.
1083 * this preserves symbolic pattern of L, discards numeric
1084 * values */
1085 ll_super_to_super_symbolic (L, Common) ; /* cannot fail */
1086 }
1087 else if (L->xtype == CHOLMOD_PATTERN && !(L->is_super))
1088 {
1089 /* convert from simplicial symbolic to supernodal symbolic.
1090 * contents of supernodal pattern are uninitialized. Not meant
1091 * for the end user. */
1092 simplicial_symbolic_to_super_symbolic (L, Common) ;
1093 }
1094 else
1095 {
1096 /* cannot convert from simplicial numeric to supernodal
1097 * symbolic */
1098 ERROR (CHOLMOD_INVALID,
1099 "cannot convert L to supernodal symbolic") ;
1100 }
1101 }
1102
1103 }
1104 else
1105 {
1106
1107 /* ------------------------------------------------------------------ */
1108 /* convert to numeric */
1109 /* ------------------------------------------------------------------ */
1110
1111 if (to_super)
1112 {
1113
1114 /* -------------------------------------------------------------- */
1115 /* convert to supernodal numeric factor */
1116 /* -------------------------------------------------------------- */
1117
1118 if (L->xtype == CHOLMOD_PATTERN)
1119 {
1120 if (L->is_super)
1121 {
1122 /* Convert supernodal symbolic to supernodal numeric.
1123 * Contents of supernodal numeric values are uninitialized.
1124 * This is used by cholmod_super_numeric. Not meant for
1125 * the end user. */
1126 super_symbolic_to_ll_super (to_xtype, L, Common) ;
1127 }
1128 else
1129 {
1130 /* Convert simplicial symbolic to supernodal numeric.
1131 * Contents not defined. This is used by
1132 * Core/cholmod_copy_factor only. Not meant for the end
1133 * user. */
1134 if (!simplicial_symbolic_to_super_symbolic (L, Common))
1135 {
1136 /* failure, convert back to simplicial symbolic */
1137 any_to_simplicial_symbolic (L, to_ll, Common) ;
1138 }
1139 else
1140 {
1141 /* conversion to super symbolic OK, allocate numeric
1142 * part */
1143 super_symbolic_to_ll_super (to_xtype, L, Common) ;
1144 }
1145 }
1146 }
1147 else
1148 {
1149 /* nothing to do if L is already in supernodal numeric form */
1150 if (!(L->is_super))
1151 {
1152 ERROR (CHOLMOD_INVALID,
1153 "cannot convert simplicial L to supernodal") ;
1154 }
1155 /* FUTURE WORK: convert to/from supernodal LL' and LDL' */
1156 }
1157
1158 }
1159 else
1160 {
1161
1162 /* -------------------------------------------------------------- */
1163 /* convert any factor to simplicial numeric */
1164 /* -------------------------------------------------------------- */
1165
1166 if (L->xtype == CHOLMOD_PATTERN && !(L->is_super))
1167 {
1168
1169 /* ---------------------------------------------------------- */
1170 /* convert simplicial symbolic to simplicial numeric (L=I,D=I)*/
1171 /* ---------------------------------------------------------- */
1172
1173 simplicial_symbolic_to_simplicial_numeric (L, to_ll, to_packed,
1174 to_xtype, Common) ;
1175
1176 }
1177 else if (L->xtype != CHOLMOD_PATTERN && L->is_super)
1178 {
1179
1180 /* ---------------------------------------------------------- */
1181 /* convert a supernodal LL' to simplicial numeric */
1182 /* ---------------------------------------------------------- */
1183
1184 ll_super_to_simplicial_numeric (L, to_packed, to_ll, Common) ;
1185
1186 }
1187 else if (L->xtype == CHOLMOD_PATTERN && L->is_super)
1188 {
1189
1190 /* ---------------------------------------------------------- */
1191 /* convert a supernodal symbolic to simplicial numeric (L=D=I)*/
1192 /* ---------------------------------------------------------- */
1193
1194 any_to_simplicial_symbolic (L, to_ll, Common) ;
1195 /* if the following fails, it leaves the factor in simplicial
1196 * symbolic form */
1197 simplicial_symbolic_to_simplicial_numeric (L, to_ll, to_packed,
1198 to_xtype, Common) ;
1199
1200 }
1201 else
1202 {
1203
1204 /* ---------------------------------------------------------- */
1205 /* change a simplicial numeric factor */
1206 /* ---------------------------------------------------------- */
1207
1208 /* change LL' to LDL', LDL' to LL', or leave as-is. pack the
1209 * columns of L, or leave as-is. Ensure the columns are
1210 * monotonic, or leave as-is. */
1211
1212 change_simplicial_numeric (L, to_ll, to_packed, to_monotonic,
1213 Common) ;
1214 }
1215 }
1216 }
1217
1218 /* ---------------------------------------------------------------------- */
1219 /* return result */
1220 /* ---------------------------------------------------------------------- */
1221
1222 return (Common->status >= CHOLMOD_OK) ;
1223 }
1224