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