1 /* ========================================================================== */
2 /* === Cholesky/t_cholmod_lsolve ============================================ */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Cholesky Module. Copyright (C) 2005-2013, Timothy A. Davis
7 * -------------------------------------------------------------------------- */
8
9 /* Template routine to solve Lx=b with unit or non-unit diagonal, or
10 * solve LDx=b.
11 *
12 * The numeric xtype of L and Y must match. Y contains b on input and x on
13 * output, stored in row-form. Y is nrow-by-n, where nrow must equal 1 for the
14 * complex or zomplex cases, and nrow <= 4 for the real case.
15 *
16 * This file is not compiled separately. It is included in t_cholmod_solve.c
17 * instead. It contains no user-callable routines.
18 *
19 * workspace: none
20 *
21 * Supports real, complex, and zomplex factors.
22 */
23
24 /* undefine all prior definitions */
25 #undef FORM_NAME
26 #undef LSOLVE
27
28 /* -------------------------------------------------------------------------- */
29 /* define the method */
30 /* -------------------------------------------------------------------------- */
31
32 #ifdef LL
33 /* LL': solve Lx=b with non-unit diagonal */
34 #define FORM_NAME(prefix,rank) prefix ## ll_lsolve_ ## rank
35
36 #elif defined (LD)
37 /* LDL': solve LDx=b */
38 #define FORM_NAME(prefix,rank) prefix ## ldl_ldsolve_ ## rank
39
40 #else
41 /* LDL': solve Lx=b with unit diagonal */
42 #define FORM_NAME(prefix,rank) prefix ## ldl_lsolve_ ## rank
43
44 #endif
45
46 /* LSOLVE(k) defines the name of a routine for an n-by-k right-hand-side. */
47
48 #define LSOLVE(prefix,rank) FORM_NAME(prefix,rank)
49
50 #ifdef REAL
51
52 /* ========================================================================== */
53 /* === LSOLVE (1) =========================================================== */
54 /* ========================================================================== */
55
56 /* Solve Lx=b, where b has 1 column */
57
58 static void LSOLVE (PREFIX,1)
59 (
60 cholmod_factor *L,
61 double X [ ] /* n-by-1 in row form */
62 )
63 {
64 double *Lx = L->x ;
65 Int *Li = L->i ;
66 Int *Lp = L->p ;
67 Int *Lnz = L->nz ;
68 Int j, n = L->n ;
69
70 for (j = 0 ; j < n ; )
71 {
72 /* get the start, end, and length of column j */
73 Int p = Lp [j] ;
74 Int lnz = Lnz [j] ;
75 Int pend = p + lnz ;
76
77 /* find a chain of supernodes (up to j, j+1, and j+2) */
78 if (lnz < 4 || lnz != Lnz [j+1] + 1 || Li [p+1] != j+1)
79 {
80
81 /* -------------------------------------------------------------- */
82 /* solve with a single column of L */
83 /* -------------------------------------------------------------- */
84
85 double y = X [j] ;
86 #ifdef LL
87 y /= Lx [p] ;
88 X [j] = y ;
89 #elif defined (LD)
90 X [j] = y / Lx [p] ;
91 #endif
92 for (p++ ; p < pend ; p++)
93 {
94 X [Li [p]] -= Lx [p] * y ;
95 }
96 j++ ; /* advance to next column of L */
97
98 }
99 else if (lnz != Lnz [j+2] + 2 || Li [p+2] != j+2)
100 {
101
102 /* -------------------------------------------------------------- */
103 /* solve with a supernode of two columns of L */
104 /* -------------------------------------------------------------- */
105
106 double y [2] ;
107 Int q = Lp [j+1] ;
108 #ifdef LL
109 y [0] = X [j] / Lx [p] ;
110 y [1] = (X [j+1] - Lx [p+1] * y [0]) / Lx [q] ;
111 X [j ] = y [0] ;
112 X [j+1] = y [1] ;
113 #elif defined (LD)
114 y [0] = X [j] ;
115 y [1] = X [j+1] - Lx [p+1] * y [0] ;
116 X [j ] = y [0] / Lx [p] ;
117 X [j+1] = y [1] / Lx [q] ;
118 #else
119 y [0] = X [j] ;
120 y [1] = X [j+1] - Lx [p+1] * y [0] ;
121 X [j+1] = y [1] ;
122 #endif
123 for (p += 2, q++ ; p < pend ; p++, q++)
124 {
125 X [Li [p]] -= Lx [p] * y [0] + Lx [q] * y [1] ;
126 }
127 j += 2 ; /* advance to next column of L */
128
129 }
130 else
131 {
132
133 /* -------------------------------------------------------------- */
134 /* solve with a supernode of three columns of L */
135 /* -------------------------------------------------------------- */
136
137 double y [3] ;
138 Int q = Lp [j+1] ;
139 Int r = Lp [j+2] ;
140 #ifdef LL
141 y [0] = X [j] / Lx [p] ;
142 y [1] = (X [j+1] - Lx [p+1] * y [0]) / Lx [q] ;
143 y [2] = (X [j+2] - Lx [p+2] * y [0] - Lx [q+1] * y [1]) / Lx [r] ;
144 X [j ] = y [0] ;
145 X [j+1] = y [1] ;
146 X [j+2] = y [2] ;
147 #elif defined (LD)
148 y [0] = X [j] ;
149 y [1] = X [j+1] - Lx [p+1] * y [0] ;
150 y [2] = X [j+2] - Lx [p+2] * y [0] - Lx [q+1] * y [1] ;
151 X [j ] = y [0] / Lx [p] ;
152 X [j+1] = y [1] / Lx [q] ;
153 X [j+2] = y [2] / Lx [r] ;
154 #else
155 y [0] = X [j] ;
156 y [1] = X [j+1] - Lx [p+1] * y [0] ;
157 y [2] = X [j+2] - Lx [p+2] * y [0] - Lx [q+1] * y [1] ;
158 X [j+1] = y [1] ;
159 X [j+2] = y [2] ;
160 #endif
161 for (p += 3, q += 2, r++ ; p < pend ; p++, q++, r++)
162 {
163 X [Li [p]] -= Lx [p] * y [0] + Lx [q] * y [1] + Lx [r] * y [2] ;
164 }
165 j += 3 ; /* advance to next column of L */
166 }
167 }
168 }
169
170
171 /* ========================================================================== */
172 /* === LSOLVE (2) =========================================================== */
173 /* ========================================================================== */
174
175 /* Solve Lx=b, where b has 2 columns */
176
177 static void LSOLVE (PREFIX,2)
178 (
179 cholmod_factor *L,
180 double X [ ][2] /* n-by-2 in row form */
181 )
182 {
183 double *Lx = L->x ;
184 Int *Li = L->i ;
185 Int *Lp = L->p ;
186 Int *Lnz = L->nz ;
187 Int j, n = L->n ;
188
189 for (j = 0 ; j < n ; )
190 {
191 /* get the start, end, and length of column j */
192 Int p = Lp [j] ;
193 Int lnz = Lnz [j] ;
194 Int pend = p + lnz ;
195
196 /* find a chain of supernodes (up to j, j+1, and j+2) */
197 if (lnz < 4 || lnz != Lnz [j+1] + 1 || Li [p+1] != j+1)
198 {
199
200 /* -------------------------------------------------------------- */
201 /* solve with a single column of L */
202 /* -------------------------------------------------------------- */
203
204 double y [2] ;
205 y [0] = X [j][0] ;
206 y [1] = X [j][1] ;
207 #ifdef LL
208 y [0] /= Lx [p] ;
209 y [1] /= Lx [p] ;
210 X [j][0] = y [0] ;
211 X [j][1] = y [1] ;
212 #elif defined (LD)
213 X [j][0] = y [0] / Lx [p] ;
214 X [j][1] = y [1] / Lx [p] ;
215 #endif
216 for (p++ ; p < pend ; p++)
217 {
218 Int i = Li [p] ;
219 X [i][0] -= Lx [p] * y [0] ;
220 X [i][1] -= Lx [p] * y [1] ;
221 }
222 j++ ; /* advance to next column of L */
223
224 }
225 else if (lnz != Lnz [j+2] + 2 || Li [p+2] != j+2)
226 {
227
228 /* -------------------------------------------------------------- */
229 /* solve with a supernode of two columns of L */
230 /* -------------------------------------------------------------- */
231
232 double y [2][2] ;
233 Int q = Lp [j+1] ;
234 y [0][0] = X [j][0] ;
235 y [0][1] = X [j][1] ;
236 #ifdef LL
237 y [0][0] /= Lx [p] ;
238 y [0][1] /= Lx [p] ;
239 y [1][0] = (X [j+1][0] - Lx [p+1] * y [0][0]) / Lx [q] ;
240 y [1][1] = (X [j+1][1] - Lx [p+1] * y [0][1]) / Lx [q] ;
241 X [j ][0] = y [0][0] ;
242 X [j ][1] = y [0][1] ;
243 X [j+1][0] = y [1][0] ;
244 X [j+1][1] = y [1][1] ;
245 #elif defined (LD)
246 y [1][0] = X [j+1][0] - Lx [p+1] * y [0][0] ;
247 y [1][1] = X [j+1][1] - Lx [p+1] * y [0][1] ;
248 X [j ][0] = y [0][0] / Lx [p] ;
249 X [j ][1] = y [0][1] / Lx [p] ;
250 X [j+1][0] = y [1][0] / Lx [q] ;
251 X [j+1][1] = y [1][1] / Lx [q] ;
252 #else
253 y [1][0] = X [j+1][0] - Lx [p+1] * y [0][0] ;
254 y [1][1] = X [j+1][1] - Lx [p+1] * y [0][1] ;
255 X [j+1][0] = y [1][0] ;
256 X [j+1][1] = y [1][1] ;
257 #endif
258 for (p += 2, q++ ; p < pend ; p++, q++)
259 {
260 Int i = Li [p] ;
261 X [i][0] -= Lx [p] * y [0][0] + Lx [q] * y [1][0] ;
262 X [i][1] -= Lx [p] * y [0][1] + Lx [q] * y [1][1] ;
263 }
264 j += 2 ; /* advance to next column of L */
265
266 }
267 else
268 {
269
270 /* -------------------------------------------------------------- */
271 /* solve with a supernode of three columns of L */
272 /* -------------------------------------------------------------- */
273
274 double y [3][2] ;
275 Int q = Lp [j+1] ;
276 Int r = Lp [j+2] ;
277 y [0][0] = X [j][0] ;
278 y [0][1] = X [j][1] ;
279 #ifdef LL
280 y [0][0] /= Lx [p] ;
281 y [0][1] /= Lx [p] ;
282 y [1][0] = (X [j+1][0] - Lx[p+1] * y[0][0]) / Lx [q] ;
283 y [1][1] = (X [j+1][1] - Lx[p+1] * y[0][1]) / Lx [q] ;
284 y [2][0] = (X [j+2][0] - Lx[p+2] * y[0][0] - Lx[q+1]*y[1][0])/Lx[r];
285 y [2][1] = (X [j+2][1] - Lx[p+2] * y[0][1] - Lx[q+1]*y[1][1])/Lx[r];
286 X [j ][0] = y [0][0] ;
287 X [j ][1] = y [0][1] ;
288 X [j+1][0] = y [1][0] ;
289 X [j+1][1] = y [1][1] ;
290 X [j+2][0] = y [2][0] ;
291 X [j+2][1] = y [2][1] ;
292 #elif defined (LD)
293 y [1][0] = X [j+1][0] - Lx [p+1] * y [0][0] ;
294 y [1][1] = X [j+1][1] - Lx [p+1] * y [0][1] ;
295 y [2][0] = X [j+2][0] - Lx [p+2] * y [0][0] - Lx [q+1] * y [1][0] ;
296 y [2][1] = X [j+2][1] - Lx [p+2] * y [0][1] - Lx [q+1] * y [1][1] ;
297 X [j ][0] = y [0][0] / Lx [p] ;
298 X [j ][1] = y [0][1] / Lx [p] ;
299 X [j+1][0] = y [1][0] / Lx [q] ;
300 X [j+1][1] = y [1][1] / Lx [q] ;
301 X [j+2][0] = y [2][0] / Lx [r] ;
302 X [j+2][1] = y [2][1] / Lx [r] ;
303 #else
304 y [1][0] = X [j+1][0] - Lx [p+1] * y [0][0] ;
305 y [1][1] = X [j+1][1] - Lx [p+1] * y [0][1] ;
306 y [2][0] = X [j+2][0] - Lx [p+2] * y [0][0] - Lx [q+1] * y [1][0] ;
307 y [2][1] = X [j+2][1] - Lx [p+2] * y [0][1] - Lx [q+1] * y [1][1] ;
308 X [j+1][0] = y [1][0] ;
309 X [j+1][1] = y [1][1] ;
310 X [j+2][0] = y [2][0] ;
311 X [j+2][1] = y [2][1] ;
312 #endif
313 for (p += 3, q += 2, r++ ; p < pend ; p++, q++, r++)
314 {
315 Int i = Li [p] ;
316 X[i][0] -= Lx[p] * y[0][0] + Lx[q] * y[1][0] + Lx[r] * y[2][0] ;
317 X[i][1] -= Lx[p] * y[0][1] + Lx[q] * y[1][1] + Lx[r] * y[2][1] ;
318 }
319 j += 3 ; /* advance to next column of L */
320 }
321 }
322 }
323
324
325 /* ========================================================================== */
326 /* === LSOLVE (3) =========================================================== */
327 /* ========================================================================== */
328
329 /* Solve Lx=b, where b has 3 columns */
330
331 static void LSOLVE (PREFIX,3)
332 (
333 cholmod_factor *L,
334 double X [ ][3] /* n-by-3 in row form */
335 )
336 {
337 double *Lx = L->x ;
338 Int *Li = L->i ;
339 Int *Lp = L->p ;
340 Int *Lnz = L->nz ;
341 Int j, n = L->n ;
342
343 for (j = 0 ; j < n ; )
344 {
345 /* get the start, end, and length of column j */
346 Int p = Lp [j] ;
347 Int lnz = Lnz [j] ;
348 Int pend = p + lnz ;
349
350 /* find a chain of supernodes (up to j, j+1, and j+2) */
351 if (lnz < 4 || lnz != Lnz [j+1] + 1 || Li [p+1] != j+1)
352 {
353
354 /* -------------------------------------------------------------- */
355 /* solve with a single column of L */
356 /* -------------------------------------------------------------- */
357
358 double y [3] ;
359 y [0] = X [j][0] ;
360 y [1] = X [j][1] ;
361 y [2] = X [j][2] ;
362 #ifdef LL
363 y [0] /= Lx [p] ;
364 y [1] /= Lx [p] ;
365 y [2] /= Lx [p] ;
366 X [j][0] = y [0] ;
367 X [j][1] = y [1] ;
368 X [j][2] = y [2] ;
369 #elif defined (LD)
370 X [j][0] = y [0] / Lx [p] ;
371 X [j][1] = y [1] / Lx [p] ;
372 X [j][2] = y [2] / Lx [p] ;
373 #endif
374 for (p++ ; p < pend ; p++)
375 {
376 Int i = Li [p] ;
377 double lx = Lx [p] ;
378 X [i][0] -= lx * y [0] ;
379 X [i][1] -= lx * y [1] ;
380 X [i][2] -= lx * y [2] ;
381 }
382 j++ ; /* advance to next column of L */
383
384 }
385 else if (lnz != Lnz [j+2] + 2 || Li [p+2] != j+2)
386 {
387
388 /* -------------------------------------------------------------- */
389 /* solve with a supernode of two columns of L */
390 /* -------------------------------------------------------------- */
391
392 double y [2][3] ;
393 Int q = Lp [j+1] ;
394 y [0][0] = X [j][0] ;
395 y [0][1] = X [j][1] ;
396 y [0][2] = X [j][2] ;
397 #ifdef LL
398 y [0][0] /= Lx [p] ;
399 y [0][1] /= Lx [p] ;
400 y [0][2] /= Lx [p] ;
401 y [1][0] = (X [j+1][0] - Lx [p+1] * y [0][0]) / Lx [q] ;
402 y [1][1] = (X [j+1][1] - Lx [p+1] * y [0][1]) / Lx [q] ;
403 y [1][2] = (X [j+1][2] - Lx [p+1] * y [0][2]) / Lx [q] ;
404 X [j ][0] = y [0][0] ;
405 X [j ][1] = y [0][1] ;
406 X [j ][2] = y [0][2] ;
407 X [j+1][0] = y [1][0] ;
408 X [j+1][1] = y [1][1] ;
409 X [j+1][2] = y [1][2] ;
410 #elif defined (LD)
411 y [1][0] = X [j+1][0] - Lx [p+1] * y [0][0] ;
412 y [1][1] = X [j+1][1] - Lx [p+1] * y [0][1] ;
413 y [1][2] = X [j+1][2] - Lx [p+1] * y [0][2] ;
414 X [j ][0] = y [0][0] / Lx [p] ;
415 X [j ][1] = y [0][1] / Lx [p] ;
416 X [j ][2] = y [0][2] / Lx [p] ;
417 X [j+1][0] = y [1][0] / Lx [q] ;
418 X [j+1][1] = y [1][1] / Lx [q] ;
419 X [j+1][2] = y [1][2] / Lx [q] ;
420 #else
421 y [1][0] = X [j+1][0] - Lx [p+1] * y [0][0] ;
422 y [1][1] = X [j+1][1] - Lx [p+1] * y [0][1] ;
423 y [1][2] = X [j+1][2] - Lx [p+1] * y [0][2] ;
424 X [j+1][0] = y [1][0] ;
425 X [j+1][1] = y [1][1] ;
426 X [j+1][2] = y [1][2] ;
427 #endif
428 for (p += 2, q++ ; p < pend ; p++, q++)
429 {
430 Int i = Li [p] ;
431 double lx [2] ;
432 lx [0] = Lx [p] ;
433 lx [1] = Lx [q] ;
434 X [i][0] -= lx [0] * y [0][0] + lx [1] * y [1][0] ;
435 X [i][1] -= lx [0] * y [0][1] + lx [1] * y [1][1] ;
436 X [i][2] -= lx [0] * y [0][2] + lx [1] * y [1][2] ;
437 }
438 j += 2 ; /* advance to next column of L */
439
440 }
441 else
442 {
443
444 /* -------------------------------------------------------------- */
445 /* solve with a supernode of three columns of L */
446 /* -------------------------------------------------------------- */
447
448 double y [3][3] ;
449 Int q = Lp [j+1] ;
450 Int r = Lp [j+2] ;
451 y [0][0] = X [j][0] ;
452 y [0][1] = X [j][1] ;
453 y [0][2] = X [j][2] ;
454 #ifdef LL
455 y [0][0] /= Lx [p] ;
456 y [0][1] /= Lx [p] ;
457 y [0][2] /= Lx [p] ;
458 y [1][0] = (X [j+1][0] - Lx[p+1] * y[0][0]) / Lx [q] ;
459 y [1][1] = (X [j+1][1] - Lx[p+1] * y[0][1]) / Lx [q] ;
460 y [1][2] = (X [j+1][2] - Lx[p+1] * y[0][2]) / Lx [q] ;
461 y [2][0] = (X [j+2][0] - Lx[p+2] * y[0][0] - Lx[q+1]*y[1][0])/Lx[r];
462 y [2][1] = (X [j+2][1] - Lx[p+2] * y[0][1] - Lx[q+1]*y[1][1])/Lx[r];
463 y [2][2] = (X [j+2][2] - Lx[p+2] * y[0][2] - Lx[q+1]*y[1][2])/Lx[r];
464 X [j ][0] = y [0][0] ;
465 X [j ][1] = y [0][1] ;
466 X [j ][2] = y [0][2] ;
467 X [j+1][0] = y [1][0] ;
468 X [j+1][1] = y [1][1] ;
469 X [j+1][2] = y [1][2] ;
470 X [j+2][0] = y [2][0] ;
471 X [j+2][1] = y [2][1] ;
472 X [j+2][2] = y [2][2] ;
473 #elif defined (LD)
474 y [1][0] = X [j+1][0] - Lx [p+1] * y [0][0] ;
475 y [1][1] = X [j+1][1] - Lx [p+1] * y [0][1] ;
476 y [1][2] = X [j+1][2] - Lx [p+1] * y [0][2] ;
477 y [2][0] = X [j+2][0] - Lx [p+2] * y [0][0] - Lx [q+1] * y [1][0] ;
478 y [2][1] = X [j+2][1] - Lx [p+2] * y [0][1] - Lx [q+1] * y [1][1] ;
479 y [2][2] = X [j+2][2] - Lx [p+2] * y [0][2] - Lx [q+1] * y [1][2] ;
480 X [j ][0] = y [0][0] / Lx [p] ;
481 X [j ][1] = y [0][1] / Lx [p] ;
482 X [j ][2] = y [0][2] / Lx [p] ;
483 X [j+1][0] = y [1][0] / Lx [q] ;
484 X [j+1][1] = y [1][1] / Lx [q] ;
485 X [j+1][2] = y [1][2] / Lx [q] ;
486 X [j+2][0] = y [2][0] / Lx [r] ;
487 X [j+2][1] = y [2][1] / Lx [r] ;
488 X [j+2][2] = y [2][2] / Lx [r] ;
489 #else
490 y [1][0] = X [j+1][0] - Lx [p+1] * y [0][0] ;
491 y [1][1] = X [j+1][1] - Lx [p+1] * y [0][1] ;
492 y [1][2] = X [j+1][2] - Lx [p+1] * y [0][2] ;
493 y [2][0] = X [j+2][0] - Lx [p+2] * y [0][0] - Lx [q+1] * y [1][0] ;
494 y [2][1] = X [j+2][1] - Lx [p+2] * y [0][1] - Lx [q+1] * y [1][1] ;
495 y [2][2] = X [j+2][2] - Lx [p+2] * y [0][2] - Lx [q+1] * y [1][2] ;
496 X [j+1][0] = y [1][0] ;
497 X [j+1][1] = y [1][1] ;
498 X [j+1][2] = y [1][2] ;
499 X [j+2][0] = y [2][0] ;
500 X [j+2][1] = y [2][1] ;
501 X [j+2][2] = y [2][2] ;
502 #endif
503 for (p += 3, q += 2, r++ ; p < pend ; p++, q++, r++)
504 {
505 Int i = Li [p] ;
506 double lx [3] ;
507 lx [0] = Lx [p] ;
508 lx [1] = Lx [q] ;
509 lx [2] = Lx [r] ;
510 X [i][0] -= lx[0] * y[0][0] + lx[1] * y[1][0] + lx[2] * y[2][0];
511 X [i][1] -= lx[0] * y[0][1] + lx[1] * y[1][1] + lx[2] * y[2][1];
512 X [i][2] -= lx[0] * y[0][2] + lx[1] * y[1][2] + lx[2] * y[2][2];
513 }
514 j += 3 ; /* advance to next column of L */
515 }
516 }
517 }
518
519
520 /* ========================================================================== */
521 /* === LSOLVE (4) =========================================================== */
522 /* ========================================================================== */
523
524 /* Solve Lx=b, where b has 4 columns */
525
526 static void LSOLVE (PREFIX,4)
527 (
528 cholmod_factor *L,
529 double X [ ][4] /* n-by-4 in row form */
530 )
531 {
532 double *Lx = L->x ;
533 Int *Li = L->i ;
534 Int *Lp = L->p ;
535 Int *Lnz = L->nz ;
536 Int j, n = L->n ;
537
538 for (j = 0 ; j < n ; )
539 {
540 /* get the start, end, and length of column j */
541 Int p = Lp [j] ;
542 Int lnz = Lnz [j] ;
543 Int pend = p + lnz ;
544
545 /* find a chain of supernodes (up to j, j+1, and j+2) */
546 if (lnz < 4 || lnz != Lnz [j+1] + 1 || Li [p+1] != j+1)
547 {
548
549 /* -------------------------------------------------------------- */
550 /* solve with a single column of L */
551 /* -------------------------------------------------------------- */
552
553 double y [4] ;
554 y [0] = X [j][0] ;
555 y [1] = X [j][1] ;
556 y [2] = X [j][2] ;
557 y [3] = X [j][3] ;
558 #ifdef LL
559 y [0] /= Lx [p] ;
560 y [1] /= Lx [p] ;
561 y [2] /= Lx [p] ;
562 y [3] /= Lx [p] ;
563 X [j][0] = y [0] ;
564 X [j][1] = y [1] ;
565 X [j][2] = y [2] ;
566 X [j][3] = y [3] ;
567 #elif defined (LD)
568 X [j][0] = y [0] / Lx [p] ;
569 X [j][1] = y [1] / Lx [p] ;
570 X [j][2] = y [2] / Lx [p] ;
571 X [j][3] = y [3] / Lx [p] ;
572 #endif
573 for (p++ ; p < pend ; p++)
574 {
575 Int i = Li [p] ;
576 double lx = Lx [p] ;
577 X [i][0] -= lx * y [0] ;
578 X [i][1] -= lx * y [1] ;
579 X [i][2] -= lx * y [2] ;
580 X [i][3] -= lx * y [3] ;
581 }
582 j++ ; /* advance to next column of L */
583
584 }
585 else if (lnz != Lnz [j+2] + 2 || Li [p+2] != j+2)
586 {
587
588 /* -------------------------------------------------------------- */
589 /* solve with a supernode of two columns of L */
590 /* -------------------------------------------------------------- */
591
592 double y [2][4] ;
593 Int q = Lp [j+1] ;
594 y [0][0] = X [j][0] ;
595 y [0][1] = X [j][1] ;
596 y [0][2] = X [j][2] ;
597 y [0][3] = X [j][3] ;
598 #ifdef LL
599 y [0][0] /= Lx [p] ;
600 y [0][1] /= Lx [p] ;
601 y [0][2] /= Lx [p] ;
602 y [0][3] /= Lx [p] ;
603 y [1][0] = (X [j+1][0] - Lx [p+1] * y [0][0]) / Lx [q] ;
604 y [1][1] = (X [j+1][1] - Lx [p+1] * y [0][1]) / Lx [q] ;
605 y [1][2] = (X [j+1][2] - Lx [p+1] * y [0][2]) / Lx [q] ;
606 y [1][3] = (X [j+1][3] - Lx [p+1] * y [0][3]) / Lx [q] ;
607 X [j ][0] = y [0][0] ;
608 X [j ][1] = y [0][1] ;
609 X [j ][2] = y [0][2] ;
610 X [j ][3] = y [0][3] ;
611 X [j+1][0] = y [1][0] ;
612 X [j+1][1] = y [1][1] ;
613 X [j+1][2] = y [1][2] ;
614 X [j+1][3] = y [1][3] ;
615 #elif defined (LD)
616 y [1][0] = X [j+1][0] - Lx [p+1] * y [0][0] ;
617 y [1][1] = X [j+1][1] - Lx [p+1] * y [0][1] ;
618 y [1][2] = X [j+1][2] - Lx [p+1] * y [0][2] ;
619 y [1][3] = X [j+1][3] - Lx [p+1] * y [0][3] ;
620 X [j ][0] = y [0][0] / Lx [p] ;
621 X [j ][1] = y [0][1] / Lx [p] ;
622 X [j ][2] = y [0][2] / Lx [p] ;
623 X [j ][3] = y [0][3] / Lx [p] ;
624 X [j+1][0] = y [1][0] / Lx [q] ;
625 X [j+1][1] = y [1][1] / Lx [q] ;
626 X [j+1][2] = y [1][2] / Lx [q] ;
627 X [j+1][3] = y [1][3] / Lx [q] ;
628 #else
629 y [1][0] = X [j+1][0] - Lx [p+1] * y [0][0] ;
630 y [1][1] = X [j+1][1] - Lx [p+1] * y [0][1] ;
631 y [1][2] = X [j+1][2] - Lx [p+1] * y [0][2] ;
632 y [1][3] = X [j+1][3] - Lx [p+1] * y [0][3] ;
633 X [j+1][0] = y [1][0] ;
634 X [j+1][1] = y [1][1] ;
635 X [j+1][2] = y [1][2] ;
636 X [j+1][3] = y [1][3] ;
637 #endif
638 for (p += 2, q++ ; p < pend ; p++, q++)
639 {
640 Int i = Li [p] ;
641 double lx [2] ;
642 lx [0] = Lx [p] ;
643 lx [1] = Lx [q] ;
644 X [i][0] -= lx [0] * y [0][0] + lx [1] * y [1][0] ;
645 X [i][1] -= lx [0] * y [0][1] + lx [1] * y [1][1] ;
646 X [i][2] -= lx [0] * y [0][2] + lx [1] * y [1][2] ;
647 X [i][3] -= lx [0] * y [0][3] + lx [1] * y [1][3] ;
648 }
649 j += 2 ; /* advance to next column of L */
650
651 }
652 else
653 {
654
655 /* -------------------------------------------------------------- */
656 /* solve with a supernode of three columns of L */
657 /* -------------------------------------------------------------- */
658
659 double y [3][4] ;
660 Int q = Lp [j+1] ;
661 Int r = Lp [j+2] ;
662 y [0][0] = X [j][0] ;
663 y [0][1] = X [j][1] ;
664 y [0][2] = X [j][2] ;
665 y [0][3] = X [j][3] ;
666 #ifdef LL
667 y [0][0] /= Lx [p] ;
668 y [0][1] /= Lx [p] ;
669 y [0][2] /= Lx [p] ;
670 y [0][3] /= Lx [p] ;
671 y [1][0] = (X [j+1][0] - Lx[p+1] * y[0][0]) / Lx [q] ;
672 y [1][1] = (X [j+1][1] - Lx[p+1] * y[0][1]) / Lx [q] ;
673 y [1][2] = (X [j+1][2] - Lx[p+1] * y[0][2]) / Lx [q] ;
674 y [1][3] = (X [j+1][3] - Lx[p+1] * y[0][3]) / Lx [q] ;
675 y [2][0] = (X [j+2][0] - Lx[p+2] * y[0][0] - Lx[q+1]*y[1][0])/Lx[r];
676 y [2][1] = (X [j+2][1] - Lx[p+2] * y[0][1] - Lx[q+1]*y[1][1])/Lx[r];
677 y [2][2] = (X [j+2][2] - Lx[p+2] * y[0][2] - Lx[q+1]*y[1][2])/Lx[r];
678 y [2][3] = (X [j+2][3] - Lx[p+2] * y[0][3] - Lx[q+1]*y[1][3])/Lx[r];
679 X [j ][0] = y [0][0] ;
680 X [j ][1] = y [0][1] ;
681 X [j ][2] = y [0][2] ;
682 X [j ][3] = y [0][3] ;
683 X [j+1][0] = y [1][0] ;
684 X [j+1][1] = y [1][1] ;
685 X [j+1][2] = y [1][2] ;
686 X [j+1][3] = y [1][3] ;
687 X [j+2][0] = y [2][0] ;
688 X [j+2][1] = y [2][1] ;
689 X [j+2][2] = y [2][2] ;
690 X [j+2][3] = y [2][3] ;
691 #elif defined (LD)
692 y [1][0] = X [j+1][0] - Lx [p+1] * y [0][0] ;
693 y [1][1] = X [j+1][1] - Lx [p+1] * y [0][1] ;
694 y [1][2] = X [j+1][2] - Lx [p+1] * y [0][2] ;
695 y [1][3] = X [j+1][3] - Lx [p+1] * y [0][3] ;
696 y [2][0] = X [j+2][0] - Lx [p+2] * y [0][0] - Lx [q+1] * y [1][0] ;
697 y [2][1] = X [j+2][1] - Lx [p+2] * y [0][1] - Lx [q+1] * y [1][1] ;
698 y [2][2] = X [j+2][2] - Lx [p+2] * y [0][2] - Lx [q+1] * y [1][2] ;
699 y [2][3] = X [j+2][3] - Lx [p+2] * y [0][3] - Lx [q+1] * y [1][3] ;
700 X [j ][0] = y [0][0] / Lx [p] ;
701 X [j ][1] = y [0][1] / Lx [p] ;
702 X [j ][2] = y [0][2] / Lx [p] ;
703 X [j ][3] = y [0][3] / Lx [p] ;
704 X [j+1][0] = y [1][0] / Lx [q] ;
705 X [j+1][1] = y [1][1] / Lx [q] ;
706 X [j+1][2] = y [1][2] / Lx [q] ;
707 X [j+1][3] = y [1][3] / Lx [q] ;
708 X [j+2][0] = y [2][0] / Lx [r] ;
709 X [j+2][1] = y [2][1] / Lx [r] ;
710 X [j+2][2] = y [2][2] / Lx [r] ;
711 X [j+2][3] = y [2][3] / Lx [r] ;
712 #else
713 y [1][0] = X [j+1][0] - Lx [p+1] * y [0][0] ;
714 y [1][1] = X [j+1][1] - Lx [p+1] * y [0][1] ;
715 y [1][2] = X [j+1][2] - Lx [p+1] * y [0][2] ;
716 y [1][3] = X [j+1][3] - Lx [p+1] * y [0][3] ;
717 y [2][0] = X [j+2][0] - Lx [p+2] * y [0][0] - Lx [q+1] * y [1][0] ;
718 y [2][1] = X [j+2][1] - Lx [p+2] * y [0][1] - Lx [q+1] * y [1][1] ;
719 y [2][2] = X [j+2][2] - Lx [p+2] * y [0][2] - Lx [q+1] * y [1][2] ;
720 y [2][3] = X [j+2][3] - Lx [p+2] * y [0][3] - Lx [q+1] * y [1][3] ;
721 X [j+1][0] = y [1][0] ;
722 X [j+1][1] = y [1][1] ;
723 X [j+1][2] = y [1][2] ;
724 X [j+1][3] = y [1][3] ;
725 X [j+2][0] = y [2][0] ;
726 X [j+2][1] = y [2][1] ;
727 X [j+2][2] = y [2][2] ;
728 X [j+2][3] = y [2][3] ;
729 #endif
730 for (p += 3, q += 2, r++ ; p < pend ; p++, q++, r++)
731 {
732 Int i = Li [p] ;
733 double lx [3] ;
734 lx [0] = Lx [p] ;
735 lx [1] = Lx [q] ;
736 lx [2] = Lx [r] ;
737 X [i][0] -= lx[0] * y[0][0] + lx[1] * y[1][0] + lx[2] * y[2][0];
738 X [i][1] -= lx[0] * y[0][1] + lx[1] * y[1][1] + lx[2] * y[2][1];
739 X [i][2] -= lx[0] * y[0][2] + lx[1] * y[1][2] + lx[2] * y[2][2];
740 X [i][3] -= lx[0] * y[0][3] + lx[1] * y[1][3] + lx[2] * y[2][3];
741 }
742 j += 3 ; /* advance to next column of L */
743 }
744 }
745 }
746
747 #endif
748
749
750 /* ========================================================================== */
751 /* === LSOLVE (k) =========================================================== */
752 /* ========================================================================== */
753
LSOLVE(PREFIX,k)754 static void LSOLVE (PREFIX,k)
755 (
756 cholmod_factor *L,
757 cholmod_dense *Y, /* nr-by-n where nr is 1 to 4 */
758 Int *Yseti, Int ysetlen
759 )
760 {
761
762 double yx [2] ;
763 #ifdef ZOMPLEX
764 double yz [1] ;
765 double *Lz = L->z ;
766 double *Xz = Y->z ;
767 #endif
768 double *Lx = L->x ;
769 double *Xx = Y->x ;
770 Int *Li = L->i ;
771 Int *Lp = L->p ;
772 Int *Lnz = L->nz ;
773 Int n = L->n, jj, jjiters ;
774
775 ASSERT (L->xtype == Y->xtype) ; /* L and Y must have the same xtype */
776 ASSERT (L->n == Y->ncol) ; /* dimensions must match */
777 ASSERT (Y->nrow == Y->d) ; /* leading dimension of Y = # rows of Y */
778 ASSERT (L->xtype != CHOLMOD_PATTERN) ; /* L is not symbolic */
779 ASSERT (!(L->is_super)) ; /* L is simplicial LL' or LDL' */
780
781 #ifdef REAL
782
783 if (Yseti == NULL)
784 {
785
786 /* ------------------------------------------------------------------ */
787 /* real case, no Yseti, with 1 to 4 RHS's and dynamic supernodes */
788 /* ------------------------------------------------------------------ */
789
790 ASSERT (Y->nrow <= 4) ;
791
792 switch (Y->nrow)
793 {
794 case 1: LSOLVE (PREFIX,1) (L, Y->x) ; break ;
795 case 2: LSOLVE (PREFIX,2) (L, Y->x) ; break ;
796 case 3: LSOLVE (PREFIX,3) (L, Y->x) ; break ;
797 case 4: LSOLVE (PREFIX,4) (L, Y->x) ; break ;
798 }
799
800 }
801 else
802 #endif
803 {
804
805 /* ------------------------------------------------------------------ */
806 /* solve a complex linear system or solve with Yseti */
807 /* ------------------------------------------------------------------ */
808
809 ASSERT (Y->nrow == 1) ;
810
811 jjiters = Yseti ? ysetlen : n ;
812
813 for (jj = 0 ; jj < jjiters ; jj++)
814 {
815 Int j = Yseti ? Yseti [jj] : jj ;
816
817 /* get the start, end, and length of column j */
818 Int p = Lp [j] ;
819 Int lnz = Lnz [j] ;
820 Int pend = p + lnz ;
821
822 /* y = X [j] ; */
823 ASSIGN (yx,yz,0, Xx,Xz,j) ;
824
825 #ifdef LL
826 /* y /= Lx [p] ; */
827 /* X [j] = y ; */
828 DIV_REAL (yx,yz,0, yx,yz,0, Lx,p) ;
829 ASSIGN (Xx,Xz,j, yx,yz,0) ;
830 #elif defined (LD)
831 /* X [j] = y / Lx [p] ; */
832 DIV_REAL (Xx,Xz,j, yx,yz,0, Lx,p) ;
833 #endif
834
835 for (p++ ; p < pend ; p++)
836 {
837 /* X [Li [p]] -= Lx [p] * y ; */
838 Int i = Li [p] ;
839 MULTSUB (Xx,Xz,i, Lx,Lz,p, yx,yz,0) ;
840 }
841 }
842 }
843 }
844
845 /* prepare for the next inclusion of this file in cholmod_solve.c */
846 #undef LL
847 #undef LD
848