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