1 /* ========================================================================== */
2 /* === Cholesky/cholmod_solve =============================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Cholesky Module. Copyright (C) 2005-2013, Timothy A. Davis
7 * -------------------------------------------------------------------------- */
8
9 /* Solve one of the following systems. D is identity for an LL' factorization,
10 * in which the D operation is skipped:
11 *
12 * Ax=b 0: CHOLMOD_A x = P' * (L' \ (D \ (L \ (P * b))))
13 * LDL'x=b 1: CHOLMOD_LDLt x = (L' \ (D \ (L \ ( b))))
14 * LDx=b 2: CHOLMOD_LD x = ( (D \ (L \ ( b))))
15 * DL'x=b 3: CHOLMOD_DLt x = (L' \ (D \ ( ( b))))
16 * Lx=b 4: CHOLMOD_L x = ( ( (L \ ( b))))
17 * L'x=b 5: CHOLMOD_Lt x = (L' \ ( ( ( b))))
18 * Dx=b 6: CHOLMOD_D x = ( (D \ ( ( b))))
19 * x=Pb 7: CHOLMOD_P x = ( ( ( (P * b))))
20 * x=P'b 8: CHOLMOD_Pt x = P' * ( ( ( ( b))))
21 *
22 * The factorization can be simplicial LDL', simplicial LL', or supernodal LL'.
23 * For an LL' factorization, D is the identity matrix. Thus CHOLMOD_LD and
24 * CHOLMOD_L solve the same system if an LL' factorization was performed,
25 * for example.
26 *
27 * The supernodal solver uses BLAS routines dtrsv, dgemv, dtrsm, and dgemm,
28 * or their complex counterparts ztrsv, zgemv, ztrsm, and zgemm.
29 *
30 * If both L and B are real, then X is returned real. If either is complex
31 * or zomplex, X is returned as either complex or zomplex, depending on the
32 * Common->prefer_zomplex parameter.
33 *
34 * Supports any numeric xtype (pattern-only matrices not supported).
35 *
36 * This routine does not check to see if the diagonal of L or D is zero,
37 * because sometimes a partial solve can be done with indefinite or singular
38 * matrix. If you wish to check in your own code, test L->minor. If
39 * L->minor == L->n, then the matrix has no zero diagonal entries.
40 * If k = L->minor < L->n, then L(k,k) is zero for an LL' factorization, or
41 * D(k,k) is zero for an LDL' factorization.
42 *
43 * This routine returns X as NULL only if it runs out of memory. If L is
44 * indefinite or singular, then X may contain Inf's or NaN's, but it will
45 * exist on output.
46 */
47
48 #ifndef NCHOLESKY
49
50 #include "cholmod_internal.h"
51 #include "cholmod_cholesky.h"
52
53 #ifndef NSUPERNODAL
54 #include "cholmod_supernodal.h"
55 #endif
56
57
58 /* ========================================================================== */
59 /* === TEMPLATE ============================================================= */
60 /* ========================================================================== */
61
62 #define REAL
63 #include "t_cholmod_solve.c"
64
65 #define COMPLEX
66 #include "t_cholmod_solve.c"
67
68 #define ZOMPLEX
69 #include "t_cholmod_solve.c"
70
71 /* ========================================================================== */
72 /* === Permutation macro ==================================================== */
73 /* ========================================================================== */
74
75 /* If Perm is NULL, it is interpretted as the identity permutation */
76
77 #define P(k) ((Perm == NULL) ? (k) : Perm [k])
78
79
80 /* ========================================================================== */
81 /* === perm ================================================================= */
82 /* ========================================================================== */
83
84 /* Y = B (P (1:nrow), k1 : min (k1+ncols,ncol)-1) where B is nrow-by-ncol.
85 *
86 * Creates a permuted copy of a contiguous set of columns of B.
87 * Y is already allocated on input. Y must be of sufficient size. Let nk be
88 * the number of columns accessed in B. Y->xtype determines the complexity of
89 * the result.
90 *
91 * If B is real and Y is complex (or zomplex), only the real part of B is
92 * copied into Y. The imaginary part of Y is set to zero.
93 *
94 * If B is complex (or zomplex) and Y is real, both the real and imaginary and
95 * parts of B are returned in Y. Y is returned as nrow-by-2*nk. The even
96 * columns of Y contain the real part of B and the odd columns contain the
97 * imaginary part of B. Y->nzmax must be >= 2*nrow*nk. Otherise, Y is
98 * returned as nrow-by-nk with leading dimension nrow. Y->nzmax must be >=
99 * nrow*nk.
100 *
101 * The case where the input (B) is real and the output (Y) is zomplex is
102 * not used.
103 */
104
perm(cholmod_dense * B,Int * Perm,Int k1,Int ncols,cholmod_dense * Y)105 static void perm
106 (
107 /* ---- input ---- */
108 cholmod_dense *B, /* input matrix B */
109 Int *Perm, /* optional input permutation (can be NULL) */
110 Int k1, /* first column of B to copy */
111 Int ncols, /* last column to copy is min(k1+ncols,B->ncol)-1 */
112 /* ---- in/out --- */
113 cholmod_dense *Y /* output matrix Y, already allocated */
114 )
115 {
116 double *Yx, *Yz, *Bx, *Bz ;
117 Int k2, nk, p, k, j, nrow, ncol, d, dual, dj, j2 ;
118
119 /* ---------------------------------------------------------------------- */
120 /* get inputs */
121 /* ---------------------------------------------------------------------- */
122
123 ncol = B->ncol ;
124 nrow = B->nrow ;
125 k2 = MIN (k1+ncols, ncol) ;
126 nk = MAX (k2 - k1, 0) ;
127 dual = (Y->xtype == CHOLMOD_REAL && B->xtype != CHOLMOD_REAL) ? 2 : 1 ;
128 d = B->d ;
129 Bx = B->x ;
130 Bz = B->z ;
131 Yx = Y->x ;
132 Yz = Y->z ;
133 Y->nrow = nrow ;
134 Y->ncol = dual*nk ;
135 Y->d = nrow ;
136 ASSERT (((Int) Y->nzmax) >= nrow*nk*dual) ;
137
138 /* ---------------------------------------------------------------------- */
139 /* Y = B (P (1:nrow), k1:k2-1) */
140 /* ---------------------------------------------------------------------- */
141
142 switch (Y->xtype)
143 {
144
145 case CHOLMOD_REAL:
146
147 switch (B->xtype)
148 {
149
150 case CHOLMOD_REAL:
151 /* Y real, B real */
152 for (j = k1 ; j < k2 ; j++)
153 {
154 dj = d*j ;
155 j2 = nrow * (j-k1) ;
156 for (k = 0 ; k < nrow ; k++)
157 {
158 p = P(k) + dj ;
159 Yx [k + j2] = Bx [p] ; /* real */
160 }
161 }
162 break ;
163
164 case CHOLMOD_COMPLEX:
165 /* Y real, B complex. Y is nrow-by-2*nk */
166 for (j = k1 ; j < k2 ; j++)
167 {
168 dj = d*j ;
169 j2 = nrow * 2 * (j-k1) ;
170 for (k = 0 ; k < nrow ; k++)
171 {
172 p = P(k) + dj ;
173 Yx [k + j2 ] = Bx [2*p ] ; /* real */
174 Yx [k + j2 + nrow] = Bx [2*p+1] ; /* imag */
175 }
176 }
177 break ;
178
179 case CHOLMOD_ZOMPLEX:
180 /* Y real, B zomplex. Y is nrow-by-2*nk */
181 for (j = k1 ; j < k2 ; j++)
182 {
183 dj = d*j ;
184 j2 = nrow * 2 * (j-k1) ;
185 for (k = 0 ; k < nrow ; k++)
186 {
187 p = P(k) + dj ;
188 Yx [k + j2 ] = Bx [p] ; /* real */
189 Yx [k + j2 + nrow] = Bz [p] ; /* imag */
190 }
191 }
192 break ;
193
194 }
195 break ;
196
197 case CHOLMOD_COMPLEX:
198
199 switch (B->xtype)
200 {
201
202 case CHOLMOD_REAL:
203 /* Y complex, B real */
204 for (j = k1 ; j < k2 ; j++)
205 {
206 dj = d*j ;
207 j2 = nrow * 2 * (j-k1) ;
208 for (k = 0 ; k < nrow ; k++)
209 {
210 p = P(k) + dj ;
211 Yx [2*k + j2] = Bx [p] ; /* real */
212 Yx [2*k+1 + j2] = 0 ; /* imag */
213 }
214 }
215 break ;
216
217 case CHOLMOD_COMPLEX:
218 /* Y complex, B complex */
219 for (j = k1 ; j < k2 ; j++)
220 {
221 dj = d*j ;
222 j2 = nrow * 2 * (j-k1) ;
223 for (k = 0 ; k < nrow ; k++)
224 {
225 p = P(k) + dj ;
226 Yx [2*k + j2] = Bx [2*p ] ; /* real */
227 Yx [2*k+1 + j2] = Bx [2*p+1] ; /* imag */
228 }
229 }
230 break ;
231
232 case CHOLMOD_ZOMPLEX:
233 /* Y complex, B zomplex */
234 for (j = k1 ; j < k2 ; j++)
235 {
236 dj = d*j ;
237 j2 = nrow * 2 * (j-k1) ;
238 for (k = 0 ; k < nrow ; k++)
239 {
240 p = P(k) + dj ;
241 Yx [2*k + j2] = Bx [p] ; /* real */
242 Yx [2*k+1 + j2] = Bz [p] ; /* imag */
243 }
244 }
245 break ;
246
247 }
248 break ;
249
250 case CHOLMOD_ZOMPLEX:
251
252 switch (B->xtype)
253 {
254
255 #if 0
256 case CHOLMOD_REAL:
257 /* this case is not used */
258 break ;
259 #endif
260
261 case CHOLMOD_COMPLEX:
262 /* Y zomplex, B complex */
263 for (j = k1 ; j < k2 ; j++)
264 {
265 dj = d*j ;
266 j2 = nrow * (j-k1) ;
267 for (k = 0 ; k < nrow ; k++)
268 {
269 p = P(k) + dj ;
270 Yx [k + j2] = Bx [2*p ] ; /* real */
271 Yz [k + j2] = Bx [2*p+1] ; /* imag */
272 }
273 }
274 break ;
275
276 case CHOLMOD_ZOMPLEX:
277 /* Y zomplex, B zomplex */
278 for (j = k1 ; j < k2 ; j++)
279 {
280 dj = d*j ;
281 j2 = nrow * (j-k1) ;
282 for (k = 0 ; k < nrow ; k++)
283 {
284 p = P(k) + dj ;
285 Yx [k + j2] = Bx [p] ; /* real */
286 Yz [k + j2] = Bz [p] ; /* imag */
287 }
288 }
289 break ;
290
291 }
292 break ;
293
294 }
295 }
296
297
298 /* ========================================================================== */
299 /* === iperm ================================================================ */
300 /* ========================================================================== */
301
302 /* X (P (1:nrow), k1 : min (k1+ncols,ncol)-1) = Y where X is nrow-by-ncol.
303 *
304 * Copies and permutes Y into a contiguous set of columns of X. X is already
305 * allocated on input. Y must be of sufficient size. Let nk be the number
306 * of columns accessed in X. X->xtype determines the complexity of the result.
307 *
308 * If X is real and Y is complex (or zomplex), only the real part of B is
309 * copied into X. The imaginary part of Y is ignored.
310 *
311 * If X is complex (or zomplex) and Y is real, both the real and imaginary and
312 * parts of Y are returned in X. Y is nrow-by-2*nk. The even
313 * columns of Y contain the real part of B and the odd columns contain the
314 * imaginary part of B. Y->nzmax must be >= 2*nrow*nk. Otherise, Y is
315 * nrow-by-nk with leading dimension nrow. Y->nzmax must be >= nrow*nk.
316 *
317 * The case where the input (Y) is complex and the output (X) is real,
318 * and the case where the input (Y) is zomplex and the output (X) is real,
319 * are not used.
320 */
321
iperm(cholmod_dense * Y,Int * Perm,Int k1,Int ncols,cholmod_dense * X)322 static void iperm
323 (
324 /* ---- input ---- */
325 cholmod_dense *Y, /* input matrix Y */
326 Int *Perm, /* optional input permutation (can be NULL) */
327 Int k1, /* first column of B to copy */
328 Int ncols, /* last column to copy is min(k1+ncols,B->ncol)-1 */
329 /* ---- in/out --- */
330 cholmod_dense *X /* output matrix X, already allocated */
331 )
332 {
333 double *Yx, *Yz, *Xx, *Xz ;
334 Int k2, nk, p, k, j, nrow, ncol, d, dj, j2 ;
335
336 /* ---------------------------------------------------------------------- */
337 /* get inputs */
338 /* ---------------------------------------------------------------------- */
339
340 ncol = X->ncol ;
341 nrow = X->nrow ;
342 k2 = MIN (k1+ncols, ncol) ;
343 nk = MAX (k2 - k1, 0) ;
344 d = X->d ;
345 Xx = X->x ;
346 Xz = X->z ;
347 Yx = Y->x ;
348 Yz = Y->z ;
349 ASSERT (((Int) Y->nzmax) >= nrow*nk*
350 ((X->xtype != CHOLMOD_REAL && Y->xtype == CHOLMOD_REAL) ? 2:1)) ;
351
352 /* ---------------------------------------------------------------------- */
353 /* X (P (1:nrow), k1:k2-1) = Y */
354 /* ---------------------------------------------------------------------- */
355
356 switch (Y->xtype)
357 {
358
359 case CHOLMOD_REAL:
360
361 switch (X->xtype)
362 {
363
364 case CHOLMOD_REAL:
365 /* Y real, X real */
366 for (j = k1 ; j < k2 ; j++)
367 {
368 dj = d*j ;
369 j2 = nrow * (j-k1) ;
370 for (k = 0 ; k < nrow ; k++)
371 {
372 p = P(k) + dj ;
373 Xx [p] = Yx [k + j2] ; /* real */
374 }
375 }
376 break ;
377
378 case CHOLMOD_COMPLEX:
379 /* Y real, X complex. Y is nrow-by-2*nk */
380 for (j = k1 ; j < k2 ; j++)
381 {
382 dj = d*j ;
383 j2 = nrow * 2 * (j-k1) ;
384 for (k = 0 ; k < nrow ; k++)
385 {
386 p = P(k) + dj ;
387 Xx [2*p ] = Yx [k + j2 ] ; /* real */
388 Xx [2*p+1] = Yx [k + j2 + nrow] ; /* imag */
389 }
390 }
391 break ;
392
393 case CHOLMOD_ZOMPLEX:
394 /* Y real, X zomplex. Y is nrow-by-2*nk */
395 for (j = k1 ; j < k2 ; j++)
396 {
397 dj = d*j ;
398 j2 = nrow * 2 * (j-k1) ;
399 for (k = 0 ; k < nrow ; k++)
400 {
401 p = P(k) + dj ;
402 Xx [p] = Yx [k + j2 ] ; /* real */
403 Xz [p] = Yx [k + j2 + nrow] ; /* imag */
404 }
405 }
406 break ;
407
408 }
409 break ;
410
411 case CHOLMOD_COMPLEX:
412
413 switch (X->xtype)
414 {
415
416 #if 0
417 case CHOLMOD_REAL:
418 /* this case is not used */
419 break ;
420 #endif
421
422 case CHOLMOD_COMPLEX:
423 /* Y complex, X complex */
424 for (j = k1 ; j < k2 ; j++)
425 {
426 dj = d*j ;
427 j2 = nrow * 2 * (j-k1) ;
428 for (k = 0 ; k < nrow ; k++)
429 {
430 p = P(k) + dj ;
431 Xx [2*p ] = Yx [2*k + j2] ; /* real */
432 Xx [2*p+1] = Yx [2*k+1 + j2] ; /* imag */
433 }
434 }
435 break ;
436
437 case CHOLMOD_ZOMPLEX:
438 /* Y complex, X zomplex */
439 for (j = k1 ; j < k2 ; j++)
440 {
441 dj = d*j ;
442 j2 = nrow * 2 * (j-k1) ;
443 for (k = 0 ; k < nrow ; k++)
444 {
445 p = P(k) + dj ;
446 Xx [p] = Yx [2*k + j2] ; /* real */
447 Xz [p] = Yx [2*k+1 + j2] ; /* imag */
448 }
449 }
450 break ;
451
452 }
453 break ;
454
455 case CHOLMOD_ZOMPLEX:
456
457 switch (X->xtype)
458 {
459
460 #if 0
461 case CHOLMOD_REAL:
462 /* this case is not used */
463 break ;
464 #endif
465
466 case CHOLMOD_COMPLEX:
467 /* Y zomplex, X complex */
468 for (j = k1 ; j < k2 ; j++)
469 {
470 dj = d*j ;
471 j2 = nrow * (j-k1) ;
472 for (k = 0 ; k < nrow ; k++)
473 {
474 p = P(k) + dj ;
475 Xx [2*p ] = Yx [k + j2] ; /* real */
476 Xx [2*p+1] = Yz [k + j2] ; /* imag */
477 }
478 }
479 break ;
480
481 case CHOLMOD_ZOMPLEX:
482 /* Y zomplex, X zomplex */
483 for (j = k1 ; j < k2 ; j++)
484 {
485 dj = d*j ;
486 j2 = nrow * (j-k1) ;
487 for (k = 0 ; k < nrow ; k++)
488 {
489 p = P(k) + dj ;
490 Xx [p] = Yx [k + j2] ; /* real */
491 Xz [p] = Yz [k + j2] ; /* imag */
492 }
493 }
494 break ;
495
496 }
497 break ;
498
499 }
500 }
501
502
503 /* ========================================================================== */
504 /* === ptrans =============================================================== */
505 /* ========================================================================== */
506
507 /* Y = B (P (1:nrow), k1 : min (k1+ncols,ncol)-1)' where B is nrow-by-ncol.
508 *
509 * Creates a permuted and transposed copy of a contiguous set of columns of B.
510 * Y is already allocated on input. Y must be of sufficient size. Let nk be
511 * the number of columns accessed in B. Y->xtype determines the complexity of
512 * the result.
513 *
514 * If B is real and Y is complex (or zomplex), only the real part of B is
515 * copied into Y. The imaginary part of Y is set to zero.
516 *
517 * If B is complex (or zomplex) and Y is real, both the real and imaginary and
518 * parts of B are returned in Y. Y is returned as 2*nk-by-nrow. The even
519 * rows of Y contain the real part of B and the odd rows contain the
520 * imaginary part of B. Y->nzmax must be >= 2*nrow*nk. Otherise, Y is
521 * returned as nk-by-nrow with leading dimension nk. Y->nzmax must be >=
522 * nrow*nk.
523 *
524 * The array transpose is performed, not the complex conjugate transpose.
525 */
526
ptrans(cholmod_dense * B,Int * Perm,Int k1,Int ncols,cholmod_dense * Y)527 static void ptrans
528 (
529 /* ---- input ---- */
530 cholmod_dense *B, /* input matrix B */
531 Int *Perm, /* optional input permutation (can be NULL) */
532 Int k1, /* first column of B to copy */
533 Int ncols, /* last column to copy is min(k1+ncols,B->ncol)-1 */
534 /* ---- in/out --- */
535 cholmod_dense *Y /* output matrix Y, already allocated */
536 )
537 {
538 double *Yx, *Yz, *Bx, *Bz ;
539 Int k2, nk, p, k, j, nrow, ncol, d, dual, dj, j2 ;
540
541 /* ---------------------------------------------------------------------- */
542 /* get inputs */
543 /* ---------------------------------------------------------------------- */
544
545 ncol = B->ncol ;
546 nrow = B->nrow ;
547 k2 = MIN (k1+ncols, ncol) ;
548 nk = MAX (k2 - k1, 0) ;
549 dual = (Y->xtype == CHOLMOD_REAL && B->xtype != CHOLMOD_REAL) ? 2 : 1 ;
550 d = B->d ;
551 Bx = B->x ;
552 Bz = B->z ;
553 Yx = Y->x ;
554 Yz = Y->z ;
555 Y->nrow = dual*nk ;
556 Y->ncol = nrow ;
557 Y->d = dual*nk ;
558 ASSERT (((Int) Y->nzmax) >= nrow*nk*dual) ;
559
560 /* ---------------------------------------------------------------------- */
561 /* Y = B (P (1:nrow), k1:k2-1)' */
562 /* ---------------------------------------------------------------------- */
563
564 switch (Y->xtype)
565 {
566
567 case CHOLMOD_REAL:
568
569 switch (B->xtype)
570 {
571
572 case CHOLMOD_REAL:
573 /* Y real, B real */
574 for (j = k1 ; j < k2 ; j++)
575 {
576 dj = d*j ;
577 j2 = j-k1 ;
578 for (k = 0 ; k < nrow ; k++)
579 {
580 p = P(k) + dj ;
581 Yx [j2 + k*nk] = Bx [p] ; /* real */
582 }
583 }
584 break ;
585
586 case CHOLMOD_COMPLEX:
587 /* Y real, B complex. Y is 2*nk-by-nrow */
588 for (j = k1 ; j < k2 ; j++)
589 {
590 dj = d*j ;
591 j2 = 2*(j-k1) ;
592 for (k = 0 ; k < nrow ; k++)
593 {
594 p = P(k) + dj ;
595 Yx [j2 + k*2*nk] = Bx [2*p ] ; /* real */
596 Yx [j2+1 + k*2*nk] = Bx [2*p+1] ; /* imag */
597 }
598 }
599 break ;
600
601 case CHOLMOD_ZOMPLEX:
602 /* Y real, B zomplex. Y is 2*nk-by-nrow */
603 for (j = k1 ; j < k2 ; j++)
604 {
605 dj = d*j ;
606 j2 = 2*(j-k1) ;
607 for (k = 0 ; k < nrow ; k++)
608 {
609 p = P(k) + dj ;
610 Yx [j2 + k*2*nk] = Bx [p] ; /* real */
611 Yx [j2+1 + k*2*nk] = Bz [p] ; /* imag */
612 }
613 }
614 break ;
615
616 }
617 break ;
618
619 case CHOLMOD_COMPLEX:
620
621 switch (B->xtype)
622 {
623
624 case CHOLMOD_REAL:
625 /* Y complex, B real */
626 for (j = k1 ; j < k2 ; j++)
627 {
628 dj = d*j ;
629 j2 = 2*(j-k1) ;
630 for (k = 0 ; k < nrow ; k++)
631 {
632 p = P(k) + dj ;
633 Yx [j2 + k*2*nk] = Bx [p] ; /* real */
634 Yx [j2+1 + k*2*nk] = 0 ; /* imag */
635 }
636 }
637 break ;
638
639 case CHOLMOD_COMPLEX:
640 /* Y complex, B complex */
641 for (j = k1 ; j < k2 ; j++)
642 {
643 dj = d*j ;
644 j2 = 2*(j-k1) ;
645 for (k = 0 ; k < nrow ; k++)
646 {
647 p = P(k) + dj ;
648 Yx [j2 + k*2*nk] = Bx [2*p ] ; /* real */
649 Yx [j2+1 + k*2*nk] = Bx [2*p+1] ; /* imag */
650 }
651 }
652 break ;
653
654 case CHOLMOD_ZOMPLEX:
655 /* Y complex, B zomplex */
656 for (j = k1 ; j < k2 ; j++)
657 {
658 dj = d*j ;
659 j2 = 2*(j-k1) ;
660 for (k = 0 ; k < nrow ; k++)
661 {
662 p = P(k) + dj ;
663 Yx [j2 + k*2*nk] = Bx [p] ; /* real */
664 Yx [j2+1 + k*2*nk] = Bz [p] ; /* imag */
665 }
666 }
667 break ;
668
669 }
670 break ;
671
672 case CHOLMOD_ZOMPLEX:
673
674 switch (B->xtype)
675 {
676
677 case CHOLMOD_REAL:
678 /* Y zomplex, B real */
679 for (j = k1 ; j < k2 ; j++)
680 {
681 dj = d*j ;
682 j2 = j-k1 ;
683 for (k = 0 ; k < nrow ; k++)
684 {
685 p = P(k) + dj ;
686 Yx [j2 + k*nk] = Bx [p] ; /* real */
687 Yz [j2 + k*nk] = 0 ; /* imag */
688 }
689 }
690 break ;
691
692 case CHOLMOD_COMPLEX:
693 /* Y zomplex, B complex */
694 for (j = k1 ; j < k2 ; j++)
695 {
696 dj = d*j ;
697 j2 = j-k1 ;
698 for (k = 0 ; k < nrow ; k++)
699 {
700 p = P(k) + dj ;
701 Yx [j2 + k*nk] = Bx [2*p ] ; /* real */
702 Yz [j2 + k*nk] = Bx [2*p+1] ; /* imag */
703 }
704 }
705 break ;
706
707 case CHOLMOD_ZOMPLEX:
708 /* Y zomplex, B zomplex */
709 for (j = k1 ; j < k2 ; j++)
710 {
711 dj = d*j ;
712 j2 = j-k1 ;
713 for (k = 0 ; k < nrow ; k++)
714 {
715 p = P(k) + dj ;
716 Yx [j2 + k*nk] = Bx [p] ; /* real */
717 Yz [j2 + k*nk] = Bz [p] ; /* imag */
718 }
719 }
720 break ;
721
722 }
723 break ;
724
725 }
726 }
727
728
729 /* ========================================================================== */
730 /* === iptrans ============================================================== */
731 /* ========================================================================== */
732
733 /* X (P (1:nrow), k1 : min (k1+ncols,ncol)-1) = Y' where X is nrow-by-ncol.
734 *
735 * Copies into a permuted and transposed contiguous set of columns of X.
736 * X is already allocated on input. Y must be of sufficient size. Let nk be
737 * the number of columns accessed in X. X->xtype determines the complexity of
738 * the result.
739 *
740 * If X is real and Y is complex (or zomplex), only the real part of Y is
741 * copied into X. The imaginary part of Y is ignored.
742 *
743 * If X is complex (or zomplex) and Y is real, both the real and imaginary and
744 * parts of X are returned in Y. Y is 2*nk-by-nrow. The even
745 * rows of Y contain the real part of X and the odd rows contain the
746 * imaginary part of X. Y->nzmax must be >= 2*nrow*nk. Otherise, Y is
747 * nk-by-nrow with leading dimension nk. Y->nzmax must be >= nrow*nk.
748 *
749 * The case where Y is complex or zomplex, and X is real, is not used.
750 *
751 * The array transpose is performed, not the complex conjugate transpose.
752 */
753
iptrans(cholmod_dense * Y,Int * Perm,Int k1,Int ncols,cholmod_dense * X)754 static void iptrans
755 (
756 /* ---- input ---- */
757 cholmod_dense *Y, /* input matrix Y */
758 Int *Perm, /* optional input permutation (can be NULL) */
759 Int k1, /* first column of X to copy into */
760 Int ncols, /* last column to copy is min(k1+ncols,X->ncol)-1 */
761 /* ---- in/out --- */
762 cholmod_dense *X /* output matrix X, already allocated */
763 )
764 {
765 double *Yx, *Yz, *Xx, *Xz ;
766 Int k2, nk, p, k, j, nrow, ncol, d, dj, j2 ;
767
768 /* ---------------------------------------------------------------------- */
769 /* get inputs */
770 /* ---------------------------------------------------------------------- */
771
772 ncol = X->ncol ;
773 nrow = X->nrow ;
774 k2 = MIN (k1+ncols, ncol) ;
775 nk = MAX (k2 - k1, 0) ;
776 d = X->d ;
777 Xx = X->x ;
778 Xz = X->z ;
779 Yx = Y->x ;
780 Yz = Y->z ;
781 ASSERT (((Int) Y->nzmax) >= nrow*nk*
782 ((X->xtype != CHOLMOD_REAL && Y->xtype == CHOLMOD_REAL) ? 2:1)) ;
783
784 /* ---------------------------------------------------------------------- */
785 /* X (P (1:nrow), k1:k2-1) = Y' */
786 /* ---------------------------------------------------------------------- */
787
788 switch (Y->xtype)
789 {
790
791 case CHOLMOD_REAL:
792
793 switch (X->xtype)
794 {
795
796 case CHOLMOD_REAL:
797 /* Y real, X real */
798 for (j = k1 ; j < k2 ; j++)
799 {
800 dj = d*j ;
801 j2 = j-k1 ;
802 for (k = 0 ; k < nrow ; k++)
803 {
804 p = P(k) + dj ;
805 Xx [p] = Yx [j2 + k*nk] ; /* real */
806 }
807 }
808 break ;
809
810 case CHOLMOD_COMPLEX:
811 /* Y real, X complex. Y is 2*nk-by-nrow */
812 for (j = k1 ; j < k2 ; j++)
813 {
814 dj = d*j ;
815 j2 = 2*(j-k1) ;
816 for (k = 0 ; k < nrow ; k++)
817 {
818 p = P(k) + dj ;
819 Xx [2*p ] = Yx [j2 + k*2*nk] ; /* real */
820 Xx [2*p+1] = Yx [j2+1 + k*2*nk] ; /* imag */
821 }
822 }
823 break ;
824
825 case CHOLMOD_ZOMPLEX:
826 /* Y real, X zomplex. Y is 2*nk-by-nrow */
827 for (j = k1 ; j < k2 ; j++)
828 {
829 dj = d*j ;
830 j2 = 2*(j-k1) ;
831 for (k = 0 ; k < nrow ; k++)
832 {
833 p = P(k) + dj ;
834 Xx [p] = Yx [j2 + k*2*nk] ; /* real */
835 Xz [p] = Yx [j2+1 + k*2*nk] ; /* imag */
836 }
837 }
838 break ;
839
840 }
841 break ;
842
843 case CHOLMOD_COMPLEX:
844
845 switch (X->xtype)
846 {
847
848 #if 0
849 case CHOLMOD_REAL:
850 /* this case is not used */
851 break ;
852 #endif
853
854 case CHOLMOD_COMPLEX:
855 /* Y complex, X complex */
856 for (j = k1 ; j < k2 ; j++)
857 {
858 dj = d*j ;
859 j2 = 2*(j-k1) ;
860 for (k = 0 ; k < nrow ; k++)
861 {
862 p = P(k) + dj ;
863 Xx [2*p ] = Yx [j2 + k*2*nk] ; /* real */
864 Xx [2*p+1] = Yx [j2+1 + k*2*nk] ; /* imag */
865 }
866 }
867 break ;
868
869 case CHOLMOD_ZOMPLEX:
870 /* Y complex, X zomplex */
871 for (j = k1 ; j < k2 ; j++)
872 {
873 dj = d*j ;
874 j2 = 2*(j-k1) ;
875 for (k = 0 ; k < nrow ; k++)
876 {
877 p = P(k) + dj ;
878 Xx [p] = Yx [j2 + k*2*nk] ; /* real */
879 Xz [p] = Yx [j2+1 + k*2*nk] ; /* imag */
880 }
881 }
882 break ;
883
884 }
885 break ;
886
887 case CHOLMOD_ZOMPLEX:
888
889 switch (X->xtype)
890 {
891
892 #if 0
893 case CHOLMOD_REAL:
894 /* this case is not used */
895 break ;
896 #endif
897
898 case CHOLMOD_COMPLEX:
899 /* Y zomplex, X complex */
900 for (j = k1 ; j < k2 ; j++)
901 {
902 dj = d*j ;
903 j2 = j-k1 ;
904 for (k = 0 ; k < nrow ; k++)
905 {
906 p = P(k) + dj ;
907 Xx [2*p ] = Yx [j2 + k*nk] ; /* real */
908 Xx [2*p+1] = Yz [j2 + k*nk] ; /* imag */
909 }
910 }
911 break ;
912
913 case CHOLMOD_ZOMPLEX:
914 /* Y zomplex, X zomplex */
915 for (j = k1 ; j < k2 ; j++)
916 {
917 dj = d*j ;
918 j2 = j-k1 ;
919 for (k = 0 ; k < nrow ; k++)
920 {
921 p = P(k) + dj ;
922 Xx [p] = Yx [j2 + k*nk] ; /* real */
923 Xz [p] = Yz [j2 + k*nk] ; /* imag */
924 }
925 }
926 break ;
927
928 }
929 break ;
930
931 }
932 }
933
934
935 /* ========================================================================== */
936 /* === cholmod_solve ======================================================== */
937 /* ========================================================================== */
938
939 /* Solve a linear system.
940 *
941 * The factorization can be simplicial LDL', simplicial LL', or supernodal LL'.
942 * The Dx=b solve returns silently for the LL' factorizations (it is implicitly
943 * identity).
944 */
945
CHOLMOD(solve)946 cholmod_dense *CHOLMOD(solve)
947 (
948 /* ---- input ---- */
949 int sys, /* system to solve */
950 cholmod_factor *L, /* factorization to use */
951 cholmod_dense *B, /* right-hand-side */
952 /* --------------- */
953 cholmod_common *Common
954 )
955 {
956 cholmod_dense *Y = NULL, *X = NULL ;
957 cholmod_dense *E = NULL ;
958 int ok ;
959
960 /* do the solve, allocating workspaces as needed */
961 ok = CHOLMOD (solve2) (sys, L, B, NULL, &X, NULL, &Y, &E, Common) ;
962
963 /* free workspaces if allocated, and free result if an error occured */
964 CHOLMOD(free_dense) (&Y, Common) ;
965 CHOLMOD(free_dense) (&E, Common) ;
966 if (!ok)
967 {
968 CHOLMOD(free_dense) (&X, Common) ;
969 }
970 return (X) ;
971 }
972
973
974 /* ========================================================================== */
975 /* === cholmod_solve2 ======================================================= */
976 /* ========================================================================== */
977
978 /* This function acts just like cholmod_solve, except that the solution X and
979 * the internal workspace (Y and E) can be passed in preallocated. If the
980 * solution X or any required workspaces are not allocated on input, or if they
981 * are the wrong size or type, then this function frees them and reallocates
982 * them as the proper size and type. Thus, if you have a sequence of solves to
983 * do, you can let this function allocate X, Y, and E on the first call.
984 * Subsequent calls to cholmod_solve2 can then reuse this space. You must
985 * then free the workspaces Y and E (and X if desired) when you are finished.
986 * For example, the first call to cholmod_l_solve2, below, will solve the
987 * requested system. The next 2 calls (with different right-hand-sides but
988 * the same value of "sys") will resuse the workspace and solution X from the
989 * first call. Finally, when all solves are done, you must free the workspaces
990 * Y and E (otherwise you will have a memory leak), and you should also free X
991 * when you are done with it. Note that on input, X, Y, and E must be either
992 * valid cholmod_dense matrices, or initialized to NULL. You cannot pass in an
993 * uninitialized X, Y, or E.
994 *
995 * cholmod_dense *X = NULL, *Y = NULL, *E = NULL ;
996 * ...
997 * cholmod_l_solve2 (sys, L, B1, NULL, &X, NULL, &Y, &E, Common) ;
998 * cholmod_l_solve2 (sys, L, B2, NULL, &X, NULL, &Y, &E, Common) ;
999 * cholmod_l_solve2 (sys, L, B3, NULL, &X, NULL, &Y, &E, Common) ;
1000 * cholmod_l_free_dense (&X, Common) ;
1001 * cholmod_l_free_dense (&Y, Common) ;
1002 * cholmod_l_free_dense (&E, Common) ;
1003 *
1004 * The equivalent when using cholmod_l_solve is:
1005 *
1006 * cholmod_dense *X = NULL, *Y = NULL, *E = NULL ;
1007 * ...
1008 * X = cholmod_l_solve (sys, L, B1, Common) ;
1009 * cholmod_l_free_dense (&X, Common) ;
1010 * X = cholmod_l_solve (sys, L, B2, Common) ;
1011 * cholmod_l_free_dense (&X, Common) ;
1012 * X = cholmod_l_solve (sys, L, B3, Common) ;
1013 * cholmod_l_free_dense (&X, Common) ;
1014 *
1015 * Both methods work fine, but in the 2nd method with cholmod_solve, the
1016 * internal workspaces (Y and E) are allocated and freed on each call.
1017 *
1018 * Bset is an optional sparse column (pattern only) that specifies a set
1019 * of row indices. It is ignored if NULL, or if sys is CHOLMOD_P or
1020 * CHOLMOD_Pt. If it is present and not ignored, B must be a dense column
1021 * vector, and only entries B(i) where i is in the pattern of Bset are
1022 * considered. All others are treated as if they were zero (they are not
1023 * accessed). L must be a simplicial factorization, not supernodal. L is
1024 * converted from supernodal to simplicial if necessary. The solution X is
1025 * defined only for entries in the output sparse pattern of Xset.
1026 * The xtype (real/complex/zomplex) of L and B must match.
1027 *
1028 * NOTE: If Bset is present and L is supernodal, it is converted to simplicial
1029 * on output.
1030 */
1031
CHOLMOD(solve2)1032 int CHOLMOD(solve2) /* returns TRUE on success, FALSE on failure */
1033 (
1034 /* ---- input ---- */
1035 int sys, /* system to solve */
1036 cholmod_factor *L, /* factorization to use */
1037 cholmod_dense *B, /* right-hand-side */
1038 cholmod_sparse *Bset,
1039 /* ---- output --- */
1040 cholmod_dense **X_Handle, /* solution, allocated if need be */
1041 cholmod_sparse **Xset_Handle,
1042 /* ---- workspace */
1043 cholmod_dense **Y_Handle, /* workspace, or NULL */
1044 cholmod_dense **E_Handle, /* workspace, or NULL */
1045 /* --------------- */
1046 cholmod_common *Common
1047 )
1048 {
1049 double *Yx, *Yz, *Bx, *Bz, *Xx, *Xz ;
1050 cholmod_dense *Y = NULL, *X = NULL ;
1051 cholmod_sparse *C, *Yset, C_header, Yset_header, *Xset ;
1052 Int *Perm = NULL, *IPerm = NULL ;
1053 Int n, nrhs, ncols, ctype, xtype, k1, nr, ytype, k, blen, p, i, d, nrow ;
1054 Int Cp [2], Ysetp [2], *Ci, *Yseti, ysetlen ;
1055 Int *Bsetp, *Bseti, *Bsetnz, *Xseti, *Xsetp, *Iwork ;
1056
1057 /* ---------------------------------------------------------------------- */
1058 /* check inputs */
1059 /* ---------------------------------------------------------------------- */
1060
1061 RETURN_IF_NULL_COMMON (FALSE) ;
1062 RETURN_IF_NULL (L, FALSE) ;
1063 RETURN_IF_NULL (B, FALSE) ;
1064 RETURN_IF_XTYPE_INVALID (L, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
1065 RETURN_IF_XTYPE_INVALID (B, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
1066 if (sys < CHOLMOD_A || sys > CHOLMOD_Pt)
1067 {
1068 ERROR (CHOLMOD_INVALID, "invalid system") ;
1069 return (FALSE) ;
1070 }
1071 DEBUG (CHOLMOD(dump_factor) (L, "L", Common)) ;
1072 DEBUG (CHOLMOD(dump_dense) (B, "B", Common)) ;
1073 nrhs = B->ncol ;
1074 n = (Int) L->n ;
1075 d = (Int) B->d ;
1076 nrow = (Int) B->nrow ;
1077 if (d < n || nrow != n)
1078 {
1079 ERROR (CHOLMOD_INVALID, "dimensions of L and B do not match") ;
1080 return (FALSE) ;
1081 }
1082 if (Bset)
1083 {
1084 if (nrhs != 1)
1085 {
1086 ERROR (CHOLMOD_INVALID, "Bset requires a single right-hand side") ;
1087 return (FALSE) ;
1088 }
1089 if (L->xtype != B->xtype)
1090 {
1091 ERROR (CHOLMOD_INVALID, "Bset requires xtype of L and B to match") ;
1092 return (FALSE) ;
1093 }
1094 DEBUG (CHOLMOD(dump_sparse) (Bset, "Bset", Common)) ;
1095 }
1096 Common->status = CHOLMOD_OK ;
1097
1098 /* ---------------------------------------------------------------------- */
1099 /* get inputs */
1100 /* ---------------------------------------------------------------------- */
1101
1102 if ((sys == CHOLMOD_P || sys == CHOLMOD_Pt || sys == CHOLMOD_A)
1103 && L->ordering != CHOLMOD_NATURAL)
1104 {
1105 /* otherwise, Perm is NULL, and the identity permutation is used */
1106 Perm = L->Perm ;
1107 }
1108
1109 /* ---------------------------------------------------------------------- */
1110 /* allocate the result X (or resuse the space from a prior call) */
1111 /* ---------------------------------------------------------------------- */
1112
1113 ctype = (Common->prefer_zomplex) ? CHOLMOD_ZOMPLEX : CHOLMOD_COMPLEX ;
1114
1115 if (Bset)
1116 {
1117 xtype = L->xtype ;
1118 }
1119 else if (sys == CHOLMOD_P || sys == CHOLMOD_Pt)
1120 {
1121 /* x=Pb and x=P'b return X real if B is real; X is the preferred
1122 * complex/zcomplex type if B is complex or zomplex */
1123 xtype = (B->xtype == CHOLMOD_REAL) ? CHOLMOD_REAL : ctype ;
1124 }
1125 else if (L->xtype == CHOLMOD_REAL && B->xtype == CHOLMOD_REAL)
1126 {
1127 /* X is real if both L and B are real */
1128 xtype = CHOLMOD_REAL ;
1129 }
1130 else
1131 {
1132 /* X is complex, use the preferred complex/zomplex type */
1133 xtype = ctype ;
1134 }
1135
1136 /* ensure X has the right size and type */
1137 X = CHOLMOD(ensure_dense) (X_Handle, n, nrhs, n, xtype, Common) ;
1138 if (Common->status < CHOLMOD_OK)
1139 {
1140 return (FALSE) ;
1141 }
1142
1143 /* ---------------------------------------------------------------------- */
1144 /* solve using L, D, L', P, or some combination */
1145 /* ---------------------------------------------------------------------- */
1146
1147 if (Bset)
1148 {
1149
1150 /* ------------------------------------------------------------------ */
1151 /* solve for a subset of x, with a sparse b */
1152 /* ------------------------------------------------------------------ */
1153
1154 Int save_realloc_state ;
1155
1156 #ifndef NSUPERNODAL
1157 /* convert a supernodal L to simplicial when using Bset */
1158 if (L->is_super)
1159 {
1160 /* Can only use Bset on a simplicial factorization. The supernodal
1161 * factor L is converted to simplicial, leaving the xtype unchanged
1162 * (real, complex, or zomplex). Since the supernodal factorization
1163 * is already LL', it is left in that form. This conversion uses
1164 * the ll_super_to_simplicial_numeric function in
1165 * cholmod_change_factor.
1166 */
1167 CHOLMOD(change_factor) (
1168 CHOLMOD_REAL, /* ignored, since L is already numeric */
1169 TRUE, /* convert to LL' (no change to num. values) */
1170 FALSE, /* convert to simplicial */
1171 FALSE, /* do not pack the columns of L */
1172 FALSE, /* (ignored) */
1173 L, Common) ;
1174 if (Common->status < CHOLMOD_OK)
1175 {
1176 /* out of memory, L is returned unchanged */
1177 return (FALSE) ;
1178 }
1179 }
1180 #endif
1181
1182 /* L, X, and B are all the same xtype */
1183 /* ensure Y is the the right size */
1184 Y = CHOLMOD(ensure_dense) (Y_Handle, 1, n, 1, L->xtype, Common) ;
1185 if (Common->status < CHOLMOD_OK)
1186 {
1187 /* out of memory */
1188 return (FALSE) ;
1189 }
1190
1191 /* ------------------------------------------------------------------ */
1192 /* get the inverse permutation, constructing it if needed */
1193 /* ------------------------------------------------------------------ */
1194
1195 DEBUG (CHOLMOD (dump_perm) (Perm, n,n, "Perm", Common)) ;
1196
1197 if ((sys == CHOLMOD_A || sys == CHOLMOD_P) && Perm != NULL)
1198 {
1199 /* The inverse permutation IPerm is used for the c=Pb step,
1200 which is needed only for solving Ax=b or x=Pb. No other
1201 steps should use IPerm */
1202 if (L->IPerm == NULL)
1203 {
1204 /* construct the inverse permutation. This is done only once
1205 * and then stored in L permanently. */
1206 L->IPerm = CHOLMOD(malloc) (n, sizeof (Int), Common) ;
1207 if (Common->status < CHOLMOD_OK)
1208 {
1209 /* out of memory */
1210 return (FALSE) ;
1211 }
1212 IPerm = L->IPerm ;
1213 for (k = 0 ; k < n ; k++)
1214 {
1215 IPerm [Perm [k]] = k ;
1216 }
1217 }
1218 /* x=A\b and x=Pb both need IPerm */
1219 IPerm = L->IPerm ;
1220 }
1221
1222 if (sys == CHOLMOD_P)
1223 {
1224 /* x=Pb needs to turn off the subsequent x=P'b permutation */
1225 Perm = NULL ;
1226 }
1227
1228 DEBUG (CHOLMOD (dump_perm) (Perm, n,n, "Perm", Common)) ;
1229 DEBUG (CHOLMOD (dump_perm) (IPerm, n,n, "IPerm", Common)) ;
1230
1231 /* ------------------------------------------------------------------ */
1232 /* ensure Xset is the right size and type */
1233 /* ------------------------------------------------------------------ */
1234
1235 /* Xset is n-by-1, nzmax >= n, pattern-only, packed, unsorted */
1236 Xset = *Xset_Handle ;
1237 if (Xset == NULL || (Int) Xset->nrow != n || (Int) Xset->ncol != 1 ||
1238 (Int) Xset->nzmax < n || Xset->itype != CHOLMOD_PATTERN)
1239 {
1240 /* this is done only once, for the 1st call to cholmod_solve */
1241 CHOLMOD(free_sparse) (Xset_Handle, Common) ;
1242 Xset = CHOLMOD(allocate_sparse) (n, 1, n, FALSE, TRUE, 0,
1243 CHOLMOD_PATTERN, Common) ;
1244 *Xset_Handle = Xset ;
1245 }
1246 Xset->sorted = FALSE ;
1247 Xset->stype = 0 ;
1248 if (Common->status < CHOLMOD_OK)
1249 {
1250 /* out of memory */
1251 return (FALSE) ;
1252 }
1253
1254 /* -------------------------------------------------------------- */
1255 /* ensure Flag of size n, and 3*n Int workspace is available */
1256 /* -------------------------------------------------------------- */
1257
1258 /* does no work if prior calls already allocated enough space */
1259 CHOLMOD(allocate_work) (n, 3*n, 0, Common) ;
1260 if (Common->status < CHOLMOD_OK)
1261 {
1262 /* out of memory */
1263 return (FALSE) ;
1264 }
1265
1266 /* [ use Iwork (n:3n-1) for Ci and Yseti */
1267 Iwork = Common->Iwork ;
1268 /* Iwork (0:n-1) is not used because it is used by check_perm,
1269 print_perm, check_sparse, and print_sparse */
1270 Ci = Iwork + n ;
1271 Yseti = Ci + n ;
1272
1273 /* reallocating workspace would break Ci and Yseti */
1274 save_realloc_state = Common->no_workspace_reallocate ;
1275 Common->no_workspace_reallocate = TRUE ;
1276
1277 /* -------------------------------------------------------------- */
1278 /* C = permuted Bset, to correspond to the permutation of L */
1279 /* -------------------------------------------------------------- */
1280
1281 /* C = IPerm (Bset) */
1282 DEBUG (CHOLMOD(dump_sparse) (Bset, "Bset", Common)) ;
1283
1284 Bsetp = Bset->p ;
1285 Bseti = Bset->i ;
1286 Bsetnz = Bset->nz ;
1287 blen = (Bset->packed) ? Bsetp [1] : Bsetnz [0] ;
1288
1289 /* C = spones (P*B) or C = spones (B) if IPerm is NULL */
1290 C = &C_header ;
1291 C->nrow = n ;
1292 C->ncol = 1 ;
1293 C->nzmax = n ;
1294 C->packed = TRUE ;
1295 C->stype = 0 ;
1296 C->itype = ITYPE ;
1297 C->xtype = CHOLMOD_PATTERN ;
1298 C->dtype = CHOLMOD_DOUBLE ;
1299 C->nz = NULL ;
1300 C->p = Cp ;
1301 C->i = Ci ;
1302 C->x = NULL ;
1303 C->z = NULL ;
1304 C->sorted = FALSE ;
1305 Cp [0] = 0 ;
1306 Cp [1] = blen ;
1307 for (p = 0 ; p < blen ; p++)
1308 {
1309 Int iold = Bseti [p] ;
1310 Ci [p] = IPerm ? IPerm [iold] : iold ;
1311 }
1312 DEBUG (CHOLMOD (dump_sparse) (C, "C", Common)) ;
1313
1314 /* create a sparse column Yset from Iwork (n:2n-1) */
1315 Yset = &Yset_header ;
1316 Yset->nrow = n ;
1317 Yset->ncol = 1 ;
1318 Yset->nzmax = n ;
1319 Yset->packed = TRUE ;
1320 Yset->stype = 0 ;
1321 Yset->itype = ITYPE ;
1322 Yset->xtype = CHOLMOD_PATTERN ;
1323 Yset->dtype = CHOLMOD_DOUBLE ;
1324 Yset->nz = NULL ;
1325 Yset->p = Ysetp ;
1326 Yset->i = Yseti ;
1327 Yset->x = NULL ;
1328 Yset->z = NULL ;
1329 Yset->sorted = FALSE ;
1330 Ysetp [0] = 0 ;
1331 Ysetp [1] = 0 ;
1332 DEBUG (CHOLMOD (dump_sparse) (Yset, "Yset empty", Common)) ;
1333
1334 /* -------------------------------------------------------------- */
1335 /* Yset = nonzero pattern of L\C, or just C itself */
1336 /* -------------------------------------------------------------- */
1337
1338 /* this takes O(ysetlen) time */
1339 if (sys == CHOLMOD_P || sys == CHOLMOD_Pt || sys == CHOLMOD_D)
1340 {
1341 Ysetp [1] = blen ;
1342 for (p = 0 ; p < blen ; p++)
1343 {
1344 Yseti [p] = Ci [p] ;
1345 }
1346 }
1347 else
1348 {
1349 if (!CHOLMOD(lsolve_pattern) (C, L, Yset, Common))
1350 {
1351 Common->no_workspace_reallocate = save_realloc_state ;
1352 return (FALSE) ;
1353 }
1354 }
1355 DEBUG (CHOLMOD (dump_sparse) (Yset, "Yset", Common)) ;
1356
1357 /* -------------------------------------------------------------- */
1358 /* clear the parts of Y that we will use in the solve */
1359 /* -------------------------------------------------------------- */
1360
1361 Yx = Y->x ;
1362 Yz = Y->z ;
1363 ysetlen = Ysetp [1] ;
1364
1365 switch (L->xtype)
1366 {
1367
1368 case CHOLMOD_REAL:
1369 for (p = 0 ; p < ysetlen ; p++)
1370 {
1371 i = Yseti [p] ;
1372 Yx [i] = 0 ;
1373 }
1374 break ;
1375
1376 case CHOLMOD_COMPLEX:
1377 for (p = 0 ; p < ysetlen ; p++)
1378 {
1379 i = Yseti [p] ;
1380 Yx [2*i ] = 0 ;
1381 Yx [2*i+1] = 0 ;
1382 }
1383 break ;
1384
1385 case CHOLMOD_ZOMPLEX:
1386 for (p = 0 ; p < ysetlen ; p++)
1387 {
1388 i = Yseti [p] ;
1389 Yx [i] = 0 ;
1390 Yz [i] = 0 ;
1391 }
1392 break ;
1393 }
1394
1395 DEBUG (CHOLMOD (dump_dense) (Y, "Y (Yset) = 0", Common)) ;
1396
1397 /* -------------------------------------------------------------- */
1398 /* scatter and permute B into Y */
1399 /* -------------------------------------------------------------- */
1400
1401 /* Y (C) = B (Bset) */
1402 Bx = B->x ;
1403 Bz = B->z ;
1404
1405 switch (L->xtype)
1406 {
1407
1408 case CHOLMOD_REAL:
1409 for (p = 0 ; p < blen ; p++)
1410 {
1411 Int iold = Bseti [p] ;
1412 Int inew = Ci [p] ;
1413 Yx [inew] = Bx [iold] ;
1414 }
1415 break ;
1416
1417 case CHOLMOD_COMPLEX:
1418 for (p = 0 ; p < blen ; p++)
1419 {
1420 Int iold = Bseti [p] ;
1421 Int inew = Ci [p] ;
1422 Yx [2*inew ] = Bx [2*iold ] ;
1423 Yx [2*inew+1] = Bx [2*iold+1] ;
1424 }
1425 break ;
1426
1427 case CHOLMOD_ZOMPLEX:
1428 for (p = 0 ; p < blen ; p++)
1429 {
1430 Int iold = Bseti [p] ;
1431 Int inew = Ci [p] ;
1432 Yx [inew] = Bx [iold] ;
1433 Yz [inew] = Bz [iold] ;
1434 }
1435 break ;
1436 }
1437
1438 DEBUG (CHOLMOD (dump_dense) (Y, "Y (C) = B (Bset)", Common)) ;
1439
1440 /* -------------------------------------------------------------- */
1441 /* solve Y = (L' \ (L \ Y'))', or other system, with template */
1442 /* -------------------------------------------------------------- */
1443
1444 /* the solve only iterates over columns in Yseti [0...ysetlen-1] */
1445
1446 if (! (sys == CHOLMOD_P || sys == CHOLMOD_Pt))
1447 {
1448 switch (L->xtype)
1449 {
1450 case CHOLMOD_REAL:
1451 r_simplicial_solver (sys, L, Y, Yseti, ysetlen) ;
1452 break ;
1453
1454 case CHOLMOD_COMPLEX:
1455 c_simplicial_solver (sys, L, Y, Yseti, ysetlen) ;
1456 break ;
1457
1458 case CHOLMOD_ZOMPLEX:
1459 z_simplicial_solver (sys, L, Y, Yseti, ysetlen) ;
1460 break ;
1461 }
1462 }
1463
1464 DEBUG (CHOLMOD (dump_dense) (Y, "Y after solve", Common)) ;
1465
1466 /* -------------------------------------------------------------- */
1467 /* X = P'*Y, but only for rows in Yset, and create Xset */
1468 /* -------------------------------------------------------------- */
1469
1470 /* X (Perm (Yset)) = Y (Yset) */
1471 Xx = X->x ;
1472 Xz = X->z ;
1473 Xseti = Xset->i ;
1474 Xsetp = Xset->p ;
1475
1476 switch (L->xtype)
1477 {
1478
1479 case CHOLMOD_REAL:
1480 for (p = 0 ; p < ysetlen ; p++)
1481 {
1482 Int inew = Yseti [p] ;
1483 Int iold = Perm ? Perm [inew] : inew ;
1484 Xx [iold] = Yx [inew] ;
1485 Xseti [p] = iold ;
1486 }
1487 break ;
1488
1489 case CHOLMOD_COMPLEX:
1490 for (p = 0 ; p < ysetlen ; p++)
1491 {
1492 Int inew = Yseti [p] ;
1493 Int iold = Perm ? Perm [inew] : inew ;
1494 Xx [2*iold ] = Yx [2*inew] ;
1495 Xx [2*iold+1] = Yx [2*inew+1] ;
1496 Xseti [p] = iold ;
1497 }
1498 break ;
1499
1500 case CHOLMOD_ZOMPLEX:
1501 for (p = 0 ; p < ysetlen ; p++)
1502 {
1503 Int inew = Yseti [p] ;
1504 Int iold = Perm ? Perm [inew] : inew ;
1505 Xx [iold] = Yx [inew] ;
1506 Xz [iold] = Yz [inew] ;
1507 Xseti [p] = iold ;
1508 }
1509 break ;
1510 }
1511
1512 Xsetp [0] = 0 ;
1513 Xsetp [1] = ysetlen ;
1514
1515 DEBUG (CHOLMOD(dump_sparse) (Xset, "Xset", Common)) ;
1516 DEBUG (CHOLMOD(dump_dense) (X, "X", Common)) ;
1517 Common->no_workspace_reallocate = save_realloc_state ;
1518 /* done using Iwork (n:3n-1) for Ci and Yseti ] */
1519
1520 }
1521 else if (sys == CHOLMOD_P)
1522 {
1523
1524 /* ------------------------------------------------------------------ */
1525 /* x = P*b */
1526 /* ------------------------------------------------------------------ */
1527
1528 perm (B, Perm, 0, nrhs, X) ;
1529
1530 }
1531 else if (sys == CHOLMOD_Pt)
1532 {
1533
1534 /* ------------------------------------------------------------------ */
1535 /* x = P'*b */
1536 /* ------------------------------------------------------------------ */
1537
1538 iperm (B, Perm, 0, nrhs, X) ;
1539
1540 }
1541 else if (L->is_super)
1542 {
1543
1544 /* ------------------------------------------------------------------ */
1545 /* solve using a supernodal LL' factorization */
1546 /* ------------------------------------------------------------------ */
1547
1548 #ifndef NSUPERNODAL
1549 /* allocate workspace */
1550 cholmod_dense *E ;
1551 Int dual ;
1552 Common->blas_ok = TRUE ;
1553 dual = (L->xtype == CHOLMOD_REAL && B->xtype != CHOLMOD_REAL) ? 2 : 1 ;
1554 Y = CHOLMOD(ensure_dense) (Y_Handle, n, dual*nrhs, n, L->xtype, Common);
1555 E = CHOLMOD(ensure_dense) (E_Handle, dual*nrhs, L->maxesize, dual*nrhs,
1556 L->xtype, Common) ;
1557
1558 if (Common->status < CHOLMOD_OK)
1559 {
1560 /* out of memory */
1561 return (FALSE) ;
1562 }
1563
1564 perm (B, Perm, 0, nrhs, Y) ; /* Y = P*B */
1565
1566 if (sys == CHOLMOD_A || sys == CHOLMOD_LDLt)
1567 {
1568 CHOLMOD(super_lsolve) (L, Y, E, Common) ; /* Y = L\Y */
1569 CHOLMOD(super_ltsolve) (L, Y, E, Common) ; /* Y = L'\Y*/
1570 }
1571 else if (sys == CHOLMOD_L || sys == CHOLMOD_LD)
1572 {
1573 CHOLMOD(super_lsolve) (L, Y, E, Common) ; /* Y = L\Y */
1574 }
1575 else if (sys == CHOLMOD_Lt || sys == CHOLMOD_DLt)
1576 {
1577 CHOLMOD(super_ltsolve) (L, Y, E, Common) ; /* Y = L'\Y*/
1578 }
1579
1580 iperm (Y, Perm, 0, nrhs, X) ; /* X = P'*Y */
1581
1582 if (CHECK_BLAS_INT && !Common->blas_ok)
1583 {
1584 /* Integer overflow in the BLAS. This is probably impossible,
1585 * since the BLAS were used to create the supernodal factorization.
1586 * It might be possible for the calls to the BLAS to differ between
1587 * factorization and forward/backsolves, however. This statement
1588 * is untested; it does not appear in the compiled code if
1589 * CHECK_BLAS_INT is true (when the same integer is used in
1590 * CHOLMOD and the BLAS. */
1591 return (FALSE) ;
1592 }
1593
1594 #else
1595 /* CHOLMOD Supernodal module not installed */
1596 ERROR (CHOLMOD_NOT_INSTALLED,"Supernodal module not installed") ;
1597 #endif
1598
1599 }
1600 else
1601 {
1602
1603 /* ------------------------------------------------------------------ */
1604 /* solve using a simplicial LL' or LDL' factorization */
1605 /* ------------------------------------------------------------------ */
1606
1607 if (L->xtype == CHOLMOD_REAL && B->xtype == CHOLMOD_REAL)
1608 {
1609 /* L, B, and Y are all real */
1610 /* solve with up to 4 columns of B at a time */
1611 ncols = 4 ;
1612 nr = MAX (4, nrhs) ;
1613 ytype = CHOLMOD_REAL ;
1614 }
1615 else if (L->xtype == CHOLMOD_REAL)
1616 {
1617 /* L is real and B is complex or zomplex */
1618 /* solve with one column of B (real/imag), at a time */
1619 ncols = 1 ;
1620 nr = 2 ;
1621 ytype = CHOLMOD_REAL ;
1622 }
1623 else
1624 {
1625 /* L is complex or zomplex, B is real/complex/zomplex, Y has the
1626 * same complexity as L. Solve with one column of B at a time. */
1627 ncols = 1 ;
1628 nr = 1 ;
1629 ytype = L->xtype ;
1630 }
1631
1632 Y = CHOLMOD(ensure_dense) (Y_Handle, nr, n, nr, ytype, Common) ;
1633 if (Common->status < CHOLMOD_OK)
1634 {
1635 /* out of memory */
1636 return (FALSE) ;
1637 }
1638
1639 for (k1 = 0 ; k1 < nrhs ; k1 += ncols)
1640 {
1641
1642 /* -------------------------------------------------------------- */
1643 /* Y = B (P, k1:k1+ncols-1)' = (P * B (:,...))' */
1644 /* -------------------------------------------------------------- */
1645
1646 ptrans (B, Perm, k1, ncols, Y) ;
1647
1648 /* -------------------------------------------------------------- */
1649 /* solve Y = (L' \ (L \ Y'))', or other system, with template */
1650 /* -------------------------------------------------------------- */
1651
1652 switch (L->xtype)
1653 {
1654 case CHOLMOD_REAL:
1655 r_simplicial_solver (sys, L, Y, NULL, 0) ;
1656 break ;
1657
1658 case CHOLMOD_COMPLEX:
1659 c_simplicial_solver (sys, L, Y, NULL, 0) ;
1660 break ;
1661
1662 case CHOLMOD_ZOMPLEX:
1663 z_simplicial_solver (sys, L, Y, NULL, 0) ;
1664 break ;
1665 }
1666
1667 /* -------------------------------------------------------------- */
1668 /* X (P, k1:k2+ncols-1) = Y' */
1669 /* -------------------------------------------------------------- */
1670
1671 iptrans (Y, Perm, k1, ncols, X) ;
1672 }
1673 }
1674
1675 DEBUG (CHOLMOD(dump_dense) (X, "X result", Common)) ;
1676 return (TRUE) ;
1677 }
1678 #endif
1679