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