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