Amesos Package Browser (Single Doxygen Collection) Development
amesos_t_cholmod_ltsolve.c
Go to the documentation of this file.
```00001 /* ========================================================================== */
00002 /* === Cholesky/t_cholmod_ltsolve =========================================== */
00003 /* ========================================================================== */
00004
00005 /* -----------------------------------------------------------------------------
00006  * CHOLMOD/Cholesky Module.  Copyright (C) 2005-2006, Timothy A. Davis
00007  * The CHOLMOD/Cholesky Module is licensed under Version 2.1 of the GNU
00008  * Lesser General Public License.  See lesser.txt for a text of the license.
00009  * CHOLMOD is also available under other licenses; contact authors for details.
00010  * http://www.cise.ufl.edu/research/sparse
00011  * -------------------------------------------------------------------------- */
00012
00013 /* Template routine to solve L'x=b with unit or non-unit diagonal, or
00014  * solve DL'x=b.
00015  *
00016  * The numeric xtype of L and Y must match.  Y contains b on input and x on
00017  * output, stored in row-form.  Y is nrow-by-n, where nrow must equal 1 for the
00018  * complex or zomplex cases, and nrow <= 4 for the real case.
00019  *
00020  * This file is not compiled separately.  It is included in t_cholmod_solve.c
00021  * instead.  It contains no user-callable routines.
00022  *
00023  * workspace: none
00024  *
00025  * Supports real, complex, and zomplex factors.
00026  */
00027
00028 /* undefine all prior definitions */
00029 #undef FORM_NAME
00030 #undef LSOLVE
00031 #undef DIAG
00032
00033 /* -------------------------------------------------------------------------- */
00034 /* define the method */
00035 /* -------------------------------------------------------------------------- */
00036
00037 #ifdef LL
00038 /* LL': solve Lx=b with non-unit diagonal */
00039 #define FORM_NAME(prefix,rank) prefix ## amesos_ll_ltsolve_ ## rank
00040 #define DIAG
00041
00042 #elif defined (LD)
00043 /* LDL': solve LDx=b */
00044 #define FORM_NAME(prefix,rank) prefix ## amesos_ldl_dltsolve_ ## rank
00045 #define DIAG
00046
00047 #else
00048 /* LDL': solve Lx=b with unit diagonal */
00049 #define FORM_NAME(prefix,rank) prefix ## amesos_ldl_ltsolve_ ## rank
00050
00051 #endif
00052
00053 /* LSOLVE(k) defines the name of a routine for an n-by-k right-hand-side. */
00054 #define LSOLVE(prefix,rank) FORM_NAME(prefix,rank)
00055
00056 #ifdef REAL
00057
00058 /* ========================================================================== */
00059 /* === LSOLVE (1) =========================================================== */
00060 /* ========================================================================== */
00061
00062 /* Solve L'x=b, where b has 1 column  */
00063
00064 static void LSOLVE (PREFIX,1)
00065 (
00066     cholmod_factor *L,
00067     double X [ ]        /* n-by-1 in row form */
00068 )
00069 {
00070     double *Lx = L->x ;
00071     Int *Li = L->i ;
00072     Int *Lp = L->p ;
00073     Int *Lnz = L->nz ;
00074     Int j, n = L->n ;
00075
00076     for (j = n-1 ; j >= 0 ; )
00077     {
00078   /* get the start, end, and length of column j */
00079   Int p = Lp [j] ;
00080   Int lnz = Lnz [j] ;
00081   Int pend = p + lnz ;
00082
00083   /* find a chain of supernodes (up to j, j-1, and j-2) */
00084   if (j < 4 || lnz != Lnz [j-1] - 1 || Li [Lp [j-1]+1] != j)
00085   {
00086
00087       /* -------------------------------------------------------------- */
00088       /* solve with a single column of L */
00089       /* -------------------------------------------------------------- */
00090
00091       double y = X [j] ;
00092 #ifdef DIAG
00093       double d = Lx [p] ;
00094 #endif
00095 #ifdef LD
00096       y /= d ;
00097 #endif
00098       for (p++ ; p < pend ; p++)
00099       {
00100     y -= Lx [p] * X [Li [p]] ;
00101       }
00102 #ifdef LL
00103       X [j] = y / d ;
00104 #else
00105       X [j] = y ;
00106 #endif
00107       j-- ; /* advance to the next column of L */
00108
00109   }
00110   else if (lnz != Lnz [j-2]-2 || Li [Lp [j-2]+2] != j)
00111   {
00112
00113       /* -------------------------------------------------------------- */
00114       /* solve with a supernode of two columns of L */
00115       /* -------------------------------------------------------------- */
00116
00117       double y [2], t ;
00118       Int q = Lp [j-1] ;
00119 #ifdef DIAG
00120       double d [2] ;
00121       d [0] = Lx [p] ;
00122       d [1] = Lx [q] ;
00123 #endif
00124       t = Lx [q+1] ;
00125 #ifdef LD
00126       y [0] = X [j  ] / d [0] ;
00127       y [1] = X [j-1] / d [1] ;
00128 #else
00129       y [0] = X [j  ] ;
00130       y [1] = X [j-1] ;
00131 #endif
00132       for (p++, q += 2 ; p < pend ; p++, q++)
00133       {
00134     Int i = Li [p] ;
00135     y [0] -= Lx [p] * X [i] ;
00136     y [1] -= Lx [q] * X [i] ;
00137       }
00138 #ifdef LL
00139       y [0] /= d [0] ;
00140       y [1] = (y [1] - t * y [0]) / d [1] ;
00141 #else
00142       y [1] -= t * y [0] ;
00143 #endif
00144       X [j  ] = y [0] ;
00145       X [j-1] = y [1] ;
00146       j -= 2 ;      /* advance to the next column of L */
00147
00148   }
00149   else
00150   {
00151
00152       /* -------------------------------------------------------------- */
00153       /* solve with a supernode of three columns of L */
00154       /* -------------------------------------------------------------- */
00155
00156       double y [3], t [3] ;
00157       Int q = Lp [j-1] ;
00158       Int r = Lp [j-2] ;
00159 #ifdef DIAG
00160       double d [3] ;
00161       d [0] = Lx [p] ;
00162       d [1] = Lx [q] ;
00163       d [2] = Lx [r] ;
00164 #endif
00165       t [0] = Lx [q+1] ;
00166       t [1] = Lx [r+1] ;
00167       t [2] = Lx [r+2] ;
00168 #ifdef LD
00169       y [0] = X [j]   / d [0] ;
00170       y [1] = X [j-1] / d [1] ;
00171       y [2] = X [j-2] / d [2] ;
00172 #else
00173       y [0] = X [j] ;
00174       y [1] = X [j-1] ;
00175       y [2] = X [j-2] ;
00176 #endif
00177       for (p++, q += 2, r += 3 ; p < pend ; p++, q++, r++)
00178       {
00179     Int i = Li [p] ;
00180     y [0] -= Lx [p] * X [i] ;
00181     y [1] -= Lx [q] * X [i] ;
00182     y [2] -= Lx [r] * X [i] ;
00183       }
00184 #ifdef LL
00185       y [0] /= d [0] ;
00186       y [1] = (y [1] - t [0] * y [0]) / d [1] ;
00187       y [2] = (y [2] - t [2] * y [0] - t [1] * y [1]) / d [2] ;
00188 #else
00189       y [1] -= t [0] * y [0] ;
00190       y [2] -= t [2] * y [0] + t [1] * y [1] ;
00191 #endif
00192       X [j-2] = y [2] ;
00193       X [j-1] = y [1] ;
00194       X [j  ] = y [0] ;
00195       j -= 3 ;      /* advance to the next column of L */
00196   }
00197     }
00198 }
00199
00200
00201 /* ========================================================================== */
00202 /* === LSOLVE (2) =========================================================== */
00203 /* ========================================================================== */
00204
00205 /* Solve L'x=b, where b has 2 columns */
00206
00207 static void LSOLVE (PREFIX,2)
00208 (
00209     cholmod_factor *L,
00210     double X [ ][2]       /* n-by-2 in row form */
00211 )
00212 {
00213     double *Lx = L->x ;
00214     Int *Li = L->i ;
00215     Int *Lp = L->p ;
00216     Int *Lnz = L->nz ;
00217     Int j, n = L->n ;
00218
00219     for (j = n-1 ; j >= 0 ; )
00220     {
00221   /* get the start, end, and length of column j */
00222   Int p = Lp [j] ;
00223   Int lnz = Lnz [j] ;
00224   Int pend = p + lnz ;
00225
00226   /* find a chain of supernodes (up to j, j-1, and j-2) */
00227   if (j < 4 || lnz != Lnz [j-1] - 1 || Li [Lp [j-1]+1] != j)
00228   {
00229
00230       /* -------------------------------------------------------------- */
00231       /* solve with a single column of L */
00232       /* -------------------------------------------------------------- */
00233
00234       double y [2] ;
00235 #ifdef DIAG
00236       double d = Lx [p] ;
00237 #endif
00238 #ifdef LD
00239       y [0] = X [j][0] / d ;
00240       y [1] = X [j][1] / d ;
00241 #else
00242       y [0] = X [j][0] ;
00243       y [1] = X [j][1] ;
00244 #endif
00245       for (p++ ; p < pend ; p++)
00246       {
00247     Int i = Li [p] ;
00248     y [0] -= Lx [p] * X [i][0] ;
00249     y [1] -= Lx [p] * X [i][1] ;
00250       }
00251 #ifdef LL
00252       X [j][0] = y [0] / d ;
00253       X [j][1] = y [1] / d ;
00254 #else
00255       X [j][0] = y [0] ;
00256       X [j][1] = y [1] ;
00257 #endif
00258       j-- ; /* advance to the next column of L */
00259
00260   }
00261   else if (lnz != Lnz [j-2]-2 || Li [Lp [j-2]+2] != j)
00262   {
00263
00264       /* -------------------------------------------------------------- */
00265       /* solve with a supernode of two columns of L */
00266       /* -------------------------------------------------------------- */
00267
00268       double y [2][2], t ;
00269       Int q = Lp [j-1] ;
00270 #ifdef DIAG
00271       double d [2] ;
00272       d [0] = Lx [p] ;
00273       d [1] = Lx [q] ;
00274 #endif
00275       t = Lx [q+1] ;
00276 #ifdef LD
00277       y [0][0] = X [j  ][0] / d [0] ;
00278       y [0][1] = X [j  ][1] / d [0] ;
00279       y [1][0] = X [j-1][0] / d [1] ;
00280       y [1][1] = X [j-1][1] / d [1] ;
00281 #else
00282       y [0][0] = X [j  ][0] ;
00283       y [0][1] = X [j  ][1] ;
00284       y [1][0] = X [j-1][0] ;
00285       y [1][1] = X [j-1][1] ;
00286 #endif
00287       for (p++, q += 2 ; p < pend ; p++, q++)
00288       {
00289     Int i = Li [p] ;
00290     y [0][0] -= Lx [p] * X [i][0] ;
00291     y [0][1] -= Lx [p] * X [i][1] ;
00292     y [1][0] -= Lx [q] * X [i][0] ;
00293     y [1][1] -= Lx [q] * X [i][1] ;
00294       }
00295 #ifdef LL
00296       y [0][0] /= d [0] ;
00297       y [0][1] /= d [0] ;
00298       y [1][0] = (y [1][0] - t * y [0][0]) / d [1] ;
00299       y [1][1] = (y [1][1] - t * y [0][1]) / d [1] ;
00300 #else
00301       y [1][0] -= t * y [0][0] ;
00302       y [1][1] -= t * y [0][1] ;
00303 #endif
00304       X [j  ][0] = y [0][0] ;
00305       X [j  ][1] = y [0][1] ;
00306       X [j-1][0] = y [1][0] ;
00307       X [j-1][1] = y [1][1] ;
00308       j -= 2 ;      /* advance to the next column of L */
00309
00310   }
00311   else
00312   {
00313
00314       /* -------------------------------------------------------------- */
00315       /* solve with a supernode of three columns of L */
00316       /* -------------------------------------------------------------- */
00317
00318       double y [3][2], t [3] ;
00319       Int q = Lp [j-1] ;
00320       Int r = Lp [j-2] ;
00321 #ifdef DIAG
00322       double d [3] ;
00323       d [0] = Lx [p] ;
00324       d [1] = Lx [q] ;
00325       d [2] = Lx [r] ;
00326 #endif
00327       t [0] = Lx [q+1] ;
00328       t [1] = Lx [r+1] ;
00329       t [2] = Lx [r+2] ;
00330 #ifdef LD
00331       y [0][0] = X [j  ][0] / d [0] ;
00332       y [0][1] = X [j  ][1] / d [0] ;
00333       y [1][0] = X [j-1][0] / d [1] ;
00334       y [1][1] = X [j-1][1] / d [1] ;
00335       y [2][0] = X [j-2][0] / d [2] ;
00336       y [2][1] = X [j-2][1] / d [2] ;
00337 #else
00338       y [0][0] = X [j  ][0] ;
00339       y [0][1] = X [j  ][1] ;
00340       y [1][0] = X [j-1][0] ;
00341       y [1][1] = X [j-1][1] ;
00342       y [2][0] = X [j-2][0] ;
00343       y [2][1] = X [j-2][1] ;
00344 #endif
00345       for (p++, q += 2, r += 3 ; p < pend ; p++, q++, r++)
00346       {
00347     Int i = Li [p] ;
00348     y [0][0] -= Lx [p] * X [i][0] ;
00349     y [0][1] -= Lx [p] * X [i][1] ;
00350     y [1][0] -= Lx [q] * X [i][0] ;
00351     y [1][1] -= Lx [q] * X [i][1] ;
00352     y [2][0] -= Lx [r] * X [i][0] ;
00353     y [2][1] -= Lx [r] * X [i][1] ;
00354       }
00355 #ifdef LL
00356       y [0][0] /= d [0] ;
00357       y [0][1] /= d [0] ;
00358       y [1][0] = (y [1][0] - t [0] * y [0][0]) / d [1] ;
00359       y [1][1] = (y [1][1] - t [0] * y [0][1]) / d [1] ;
00360       y [2][0] = (y [2][0] - t [2] * y [0][0] - t [1] * y [1][0]) / d [2];
00361       y [2][1] = (y [2][1] - t [2] * y [0][1] - t [1] * y [1][1]) / d [2];
00362 #else
00363       y [1][0] -= t [0] * y [0][0] ;
00364       y [1][1] -= t [0] * y [0][1] ;
00365       y [2][0] -= t [2] * y [0][0] + t [1] * y [1][0] ;
00366       y [2][1] -= t [2] * y [0][1] + t [1] * y [1][1] ;
00367 #endif
00368       X [j  ][0] = y [0][0] ;
00369       X [j  ][1] = y [0][1] ;
00370       X [j-1][0] = y [1][0] ;
00371       X [j-1][1] = y [1][1] ;
00372       X [j-2][0] = y [2][0] ;
00373       X [j-2][1] = y [2][1] ;
00374       j -= 3 ;      /* advance to the next column of L */
00375   }
00376     }
00377 }
00378
00379
00380 /* ========================================================================== */
00381 /* === LSOLVE (3) =========================================================== */
00382 /* ========================================================================== */
00383
00384 /* Solve L'x=b, where b has 3 columns */
00385
00386 static void LSOLVE (PREFIX,3)
00387 (
00388     cholmod_factor *L,
00389     double X [ ][3]       /* n-by-3 in row form */
00390 )
00391 {
00392     double *Lx = L->x ;
00393     Int *Li = L->i ;
00394     Int *Lp = L->p ;
00395     Int *Lnz = L->nz ;
00396     Int j, n = L->n ;
00397
00398     for (j = n-1 ; j >= 0 ; )
00399     {
00400   /* get the start, end, and length of column j */
00401   Int p = Lp [j] ;
00402   Int lnz = Lnz [j] ;
00403   Int pend = p + lnz ;
00404
00405   /* find a chain of supernodes (up to j, j-1, and j-2) */
00406   if (j < 4 || lnz != Lnz [j-1] - 1 || Li [Lp [j-1]+1] != j)
00407   {
00408
00409       /* -------------------------------------------------------------- */
00410       /* solve with a single column of L */
00411       /* -------------------------------------------------------------- */
00412
00413       double y [3] ;
00414 #ifdef DIAG
00415       double d = Lx [p] ;
00416 #endif
00417 #ifdef LD
00418       y [0] = X [j][0] / d ;
00419       y [1] = X [j][1] / d ;
00420       y [2] = X [j][2] / d ;
00421 #else
00422       y [0] = X [j][0] ;
00423       y [1] = X [j][1] ;
00424       y [2] = X [j][2] ;
00425 #endif
00426       for (p++ ; p < pend ; p++)
00427       {
00428     Int i = Li [p] ;
00429     y [0] -= Lx [p] * X [i][0] ;
00430     y [1] -= Lx [p] * X [i][1] ;
00431     y [2] -= Lx [p] * X [i][2] ;
00432       }
00433 #ifdef LL
00434       X [j][0] = y [0] / d ;
00435       X [j][1] = y [1] / d ;
00436       X [j][2] = y [2] / d ;
00437 #else
00438       X [j][0] = y [0] ;
00439       X [j][1] = y [1] ;
00440       X [j][2] = y [2] ;
00441 #endif
00442       j-- ; /* advance to the next column of L */
00443
00444   }
00445   else if (lnz != Lnz [j-2]-2 || Li [Lp [j-2]+2] != j)
00446   {
00447
00448       /* -------------------------------------------------------------- */
00449       /* solve with a supernode of two columns of L */
00450       /* -------------------------------------------------------------- */
00451
00452       double y [2][3], t ;
00453       Int q = Lp [j-1] ;
00454 #ifdef DIAG
00455       double d [2] ;
00456       d [0] = Lx [p] ;
00457       d [1] = Lx [q] ;
00458 #endif
00459       t = Lx [q+1] ;
00460 #ifdef LD
00461       y [0][0] = X [j  ][0] / d [0] ;
00462       y [0][1] = X [j  ][1] / d [0] ;
00463       y [0][2] = X [j  ][2] / d [0] ;
00464       y [1][0] = X [j-1][0] / d [1] ;
00465       y [1][1] = X [j-1][1] / d [1] ;
00466       y [1][2] = X [j-1][2] / d [1] ;
00467 #else
00468       y [0][0] = X [j  ][0] ;
00469       y [0][1] = X [j  ][1] ;
00470       y [0][2] = X [j  ][2] ;
00471       y [1][0] = X [j-1][0] ;
00472       y [1][1] = X [j-1][1] ;
00473       y [1][2] = X [j-1][2] ;
00474 #endif
00475       for (p++, q += 2 ; p < pend ; p++, q++)
00476       {
00477     Int i = Li [p] ;
00478     y [0][0] -= Lx [p] * X [i][0] ;
00479     y [0][1] -= Lx [p] * X [i][1] ;
00480     y [0][2] -= Lx [p] * X [i][2] ;
00481     y [1][0] -= Lx [q] * X [i][0] ;
00482     y [1][1] -= Lx [q] * X [i][1] ;
00483     y [1][2] -= Lx [q] * X [i][2] ;
00484       }
00485 #ifdef LL
00486       y [0][0] /= d [0] ;
00487       y [0][1] /= d [0] ;
00488       y [0][2] /= d [0] ;
00489       y [1][0] = (y [1][0] - t * y [0][0]) / d [1] ;
00490       y [1][1] = (y [1][1] - t * y [0][1]) / d [1] ;
00491       y [1][2] = (y [1][2] - t * y [0][2]) / d [1] ;
00492 #else
00493       y [1][0] -= t * y [0][0] ;
00494       y [1][1] -= t * y [0][1] ;
00495       y [1][2] -= t * y [0][2] ;
00496 #endif
00497       X [j  ][0] = y [0][0] ;
00498       X [j  ][1] = y [0][1] ;
00499       X [j  ][2] = y [0][2] ;
00500       X [j-1][0] = y [1][0] ;
00501       X [j-1][1] = y [1][1] ;
00502       X [j-1][2] = y [1][2] ;
00503       j -= 2 ;      /* advance to the next column of L */
00504
00505   }
00506   else
00507   {
00508
00509       /* -------------------------------------------------------------- */
00510       /* solve with a supernode of three columns of L */
00511       /* -------------------------------------------------------------- */
00512
00513       double y [3][3], t [3] ;
00514       Int q = Lp [j-1] ;
00515       Int r = Lp [j-2] ;
00516 #ifdef DIAG
00517       double d [3] ;
00518       d [0] = Lx [p] ;
00519       d [1] = Lx [q] ;
00520       d [2] = Lx [r] ;
00521 #endif
00522       t [0] = Lx [q+1] ;
00523       t [1] = Lx [r+1] ;
00524       t [2] = Lx [r+2] ;
00525 #ifdef LD
00526       y [0][0] = X [j  ][0] / d [0] ;
00527       y [0][1] = X [j  ][1] / d [0] ;
00528       y [0][2] = X [j  ][2] / d [0] ;
00529       y [1][0] = X [j-1][0] / d [1] ;
00530       y [1][1] = X [j-1][1] / d [1] ;
00531       y [1][2] = X [j-1][2] / d [1] ;
00532       y [2][0] = X [j-2][0] / d [2] ;
00533       y [2][1] = X [j-2][1] / d [2] ;
00534       y [2][2] = X [j-2][2] / d [2] ;
00535 #else
00536       y [0][0] = X [j  ][0] ;
00537       y [0][1] = X [j  ][1] ;
00538       y [0][2] = X [j  ][2] ;
00539       y [1][0] = X [j-1][0] ;
00540       y [1][1] = X [j-1][1] ;
00541       y [1][2] = X [j-1][2] ;
00542       y [2][0] = X [j-2][0] ;
00543       y [2][1] = X [j-2][1] ;
00544       y [2][2] = X [j-2][2] ;
00545 #endif
00546       for (p++, q += 2, r += 3 ; p < pend ; p++, q++, r++)
00547       {
00548     Int i = Li [p] ;
00549     y [0][0] -= Lx [p] * X [i][0] ;
00550     y [0][1] -= Lx [p] * X [i][1] ;
00551     y [0][2] -= Lx [p] * X [i][2] ;
00552     y [1][0] -= Lx [q] * X [i][0] ;
00553     y [1][1] -= Lx [q] * X [i][1] ;
00554     y [1][2] -= Lx [q] * X [i][2] ;
00555     y [2][0] -= Lx [r] * X [i][0] ;
00556     y [2][1] -= Lx [r] * X [i][1] ;
00557     y [2][2] -= Lx [r] * X [i][2] ;
00558       }
00559 #ifdef LL
00560       y [0][0] /= d [0] ;
00561       y [0][1] /= d [0] ;
00562       y [0][2] /= d [0] ;
00563       y [1][0] = (y [1][0] - t [0] * y [0][0]) / d [1] ;
00564       y [1][1] = (y [1][1] - t [0] * y [0][1]) / d [1] ;
00565       y [1][2] = (y [1][2] - t [0] * y [0][2]) / d [1] ;
00566       y [2][0] = (y [2][0] - t [2] * y [0][0] - t [1] * y [1][0]) / d [2];
00567       y [2][1] = (y [2][1] - t [2] * y [0][1] - t [1] * y [1][1]) / d [2];
00568       y [2][2] = (y [2][2] - t [2] * y [0][2] - t [1] * y [1][2]) / d [2];
00569 #else
00570       y [1][0] -= t [0] * y [0][0] ;
00571       y [1][1] -= t [0] * y [0][1] ;
00572       y [1][2] -= t [0] * y [0][2] ;
00573       y [2][0] -= t [2] * y [0][0] + t [1] * y [1][0] ;
00574       y [2][1] -= t [2] * y [0][1] + t [1] * y [1][1] ;
00575       y [2][2] -= t [2] * y [0][2] + t [1] * y [1][2] ;
00576 #endif
00577       X [j  ][0] = y [0][0] ;
00578       X [j  ][1] = y [0][1] ;
00579       X [j  ][2] = y [0][2] ;
00580       X [j-1][0] = y [1][0] ;
00581       X [j-1][1] = y [1][1] ;
00582       X [j-1][2] = y [1][2] ;
00583       X [j-2][0] = y [2][0] ;
00584       X [j-2][1] = y [2][1] ;
00585       X [j-2][2] = y [2][2] ;
00586       j -= 3 ;      /* advance to the next column of L */
00587   }
00588     }
00589 }
00590
00591
00592 /* ========================================================================== */
00593 /* === LSOLVE (4) =========================================================== */
00594 /* ========================================================================== */
00595
00596 /* Solve L'x=b, where b has 4 columns */
00597
00598 static void LSOLVE (PREFIX,4)
00599 (
00600     cholmod_factor *L,
00601     double X [ ][4]       /* n-by-4 in row form */
00602 )
00603 {
00604     double *Lx = L->x ;
00605     Int *Li = L->i ;
00606     Int *Lp = L->p ;
00607     Int *Lnz = L->nz ;
00608     Int j, n = L->n ;
00609
00610     for (j = n-1 ; j >= 0 ; )
00611     {
00612   /* get the start, end, and length of column j */
00613   Int p = Lp [j] ;
00614   Int lnz = Lnz [j] ;
00615   Int pend = p + lnz ;
00616
00617   /* find a chain of supernodes (up to j, j-1, and j-2) */
00618   if (j < 4 || lnz != Lnz [j-1] - 1 || Li [Lp [j-1]+1] != j)
00619   {
00620
00621       /* -------------------------------------------------------------- */
00622       /* solve with a single column of L */
00623       /* -------------------------------------------------------------- */
00624
00625       double y [4] ;
00626 #ifdef DIAG
00627       double d = Lx [p] ;
00628 #endif
00629 #ifdef LD
00630       y [0] = X [j][0] / d ;
00631       y [1] = X [j][1] / d ;
00632       y [2] = X [j][2] / d ;
00633       y [3] = X [j][3] / d ;
00634 #else
00635       y [0] = X [j][0] ;
00636       y [1] = X [j][1] ;
00637       y [2] = X [j][2] ;
00638       y [3] = X [j][3] ;
00639 #endif
00640       for (p++ ; p < pend ; p++)
00641       {
00642     Int i = Li [p] ;
00643     y [0] -= Lx [p] * X [i][0] ;
00644     y [1] -= Lx [p] * X [i][1] ;
00645     y [2] -= Lx [p] * X [i][2] ;
00646     y [3] -= Lx [p] * X [i][3] ;
00647       }
00648 #ifdef LL
00649       X [j][0] = y [0] / d ;
00650       X [j][1] = y [1] / d ;
00651       X [j][2] = y [2] / d ;
00652       X [j][3] = y [3] / d ;
00653 #else
00654       X [j][0] = y [0] ;
00655       X [j][1] = y [1] ;
00656       X [j][2] = y [2] ;
00657       X [j][3] = y [3] ;
00658 #endif
00659       j-- ; /* advance to the next column of L */
00660
00661   }
00662   else /* if (j == 1 || lnz != Lnz [j-2]-2 || Li [Lp [j-2]+2] != j) */
00663   {
00664
00665       /* -------------------------------------------------------------- */
00666       /* solve with a supernode of two columns of L */
00667       /* -------------------------------------------------------------- */
00668
00669       double y [2][4], t ;
00670       Int q = Lp [j-1] ;
00671 #ifdef DIAG
00672       double d [2] ;
00673       d [0] = Lx [p] ;
00674       d [1] = Lx [q] ;
00675 #endif
00676       t = Lx [q+1] ;
00677 #ifdef LD
00678       y [0][0] = X [j  ][0] / d [0] ;
00679       y [0][1] = X [j  ][1] / d [0] ;
00680       y [0][2] = X [j  ][2] / d [0] ;
00681       y [0][3] = X [j  ][3] / d [0] ;
00682       y [1][0] = X [j-1][0] / d [1] ;
00683       y [1][1] = X [j-1][1] / d [1] ;
00684       y [1][2] = X [j-1][2] / d [1] ;
00685       y [1][3] = X [j-1][3] / d [1] ;
00686 #else
00687       y [0][0] = X [j  ][0] ;
00688       y [0][1] = X [j  ][1] ;
00689       y [0][2] = X [j  ][2] ;
00690       y [0][3] = X [j  ][3] ;
00691       y [1][0] = X [j-1][0] ;
00692       y [1][1] = X [j-1][1] ;
00693       y [1][2] = X [j-1][2] ;
00694       y [1][3] = X [j-1][3] ;
00695 #endif
00696       for (p++, q += 2 ; p < pend ; p++, q++)
00697       {
00698     Int i = Li [p] ;
00699     y [0][0] -= Lx [p] * X [i][0] ;
00700     y [0][1] -= Lx [p] * X [i][1] ;
00701     y [0][2] -= Lx [p] * X [i][2] ;
00702     y [0][3] -= Lx [p] * X [i][3] ;
00703     y [1][0] -= Lx [q] * X [i][0] ;
00704     y [1][1] -= Lx [q] * X [i][1] ;
00705     y [1][2] -= Lx [q] * X [i][2] ;
00706     y [1][3] -= Lx [q] * X [i][3] ;
00707       }
00708 #ifdef LL
00709       y [0][0] /= d [0] ;
00710       y [0][1] /= d [0] ;
00711       y [0][2] /= d [0] ;
00712       y [0][3] /= d [0] ;
00713       y [1][0] = (y [1][0] - t * y [0][0]) / d [1] ;
00714       y [1][1] = (y [1][1] - t * y [0][1]) / d [1] ;
00715       y [1][2] = (y [1][2] - t * y [0][2]) / d [1] ;
00716       y [1][3] = (y [1][3] - t * y [0][3]) / d [1] ;
00717 #else
00718       y [1][0] -= t * y [0][0] ;
00719       y [1][1] -= t * y [0][1] ;
00720       y [1][2] -= t * y [0][2] ;
00721       y [1][3] -= t * y [0][3] ;
00722 #endif
00723       X [j  ][0] = y [0][0] ;
00724       X [j  ][1] = y [0][1] ;
00725       X [j  ][2] = y [0][2] ;
00726       X [j  ][3] = y [0][3] ;
00727       X [j-1][0] = y [1][0] ;
00728       X [j-1][1] = y [1][1] ;
00729       X [j-1][2] = y [1][2] ;
00730       X [j-1][3] = y [1][3] ;
00731       j -= 2 ;      /* advance to the next column of L */
00732   }
00733
00734   /* NOTE: with 4 right-hand-sides, it suffices to exploit dynamic
00735    * supernodes of just size 1 and 2.  3-column supernodes are not
00736    * needed. */
00737     }
00738 }
00739
00740 #endif
00741
00742 /* ========================================================================== */
00743 /* === LSOLVE (k) =========================================================== */
00744 /* ========================================================================== */
00745
00746 static void LSOLVE (PREFIX,k)
00747 (
00748     cholmod_factor *L,
00749     cholmod_dense *Y        /* nr-by-n where nr is 1 to 4 */
00750 )
00751 {
00752
00753 #ifndef REAL
00754 #ifdef DIAG
00755     double d [1] ;
00756 #endif
00757     double yx [2] ;
00758 #ifdef ZOMPLEX
00759     double yz [1] ;
00760     double *Lz = L->z ;
00761     double *Xz = Y->z ;
00762 #endif
00763     double *Lx = L->x ;
00764     double *Xx = Y->x ;
00765     Int *Li = L->i ;
00766     Int *Lp = L->p ;
00767     Int *Lnz = L->nz ;
00768     Int i, j, n = L->n ;
00769 #endif
00770
00771     ASSERT (L->xtype == Y->xtype) ; /* L and Y must have the same xtype */
00772     ASSERT (L->n == Y->ncol) ;      /* dimensions must match */
00773     ASSERT (Y->nrow == Y->d) ;      /* leading dimension of Y = # rows of Y */
00774     ASSERT (L->xtype != CHOLMOD_PATTERN) ;  /* L is not symbolic */
00775     ASSERT (!(L->is_super)) ;     /* L is simplicial LL' or LDL' */
00776
00777 #ifdef REAL
00778
00779     /* ---------------------------------------------------------------------- */
00780     /* solve a real linear system, with 1 to 4 RHS's and dynamic supernodes */
00781     /* ---------------------------------------------------------------------- */
00782
00783     ASSERT (Y->nrow <= 4) ;
00784     switch (Y->nrow)
00785     {
00786   case 1: LSOLVE (PREFIX,1) (L, Y->x) ; break ;
00787   case 2: LSOLVE (PREFIX,2) (L, Y->x) ; break ;
00788   case 3: LSOLVE (PREFIX,3) (L, Y->x) ; break ;
00789   case 4: LSOLVE (PREFIX,4) (L, Y->x) ; break ;
00790     }
00791
00792 #else
00793
00794     /* ---------------------------------------------------------------------- */
00795     /* solve a complex linear system, with just one right-hand-side */
00796     /* ---------------------------------------------------------------------- */
00797
00798     ASSERT (Y->nrow == 1) ;
00799
00800     for (j = n-1 ; j >= 0 ; j--)
00801     {
00802   /* get the start, end, and length of column j */
00803   Int p = Lp [j] ;
00804   Int lnz = Lnz [j] ;
00805   Int pend = p + lnz ;
00806
00807   /* y = X [j] ; */
00808   ASSIGN (yx,yz,0, Xx,Xz,j) ;
00809
00810 #ifdef DIAG
00811   /* d = Lx [p] ; */
00812   ASSIGN_REAL (d,0, Lx,p) ;
00813 #endif
00814 #ifdef LD
00815   /* y /= d ; */
00816   DIV_REAL (yx,yz,0, yx,yz,0, d,0) ;
00817 #endif
00818
00819   for (p++ ; p < pend ; p++)
00820   {
00821       /* y -= conj (Lx [p]) * X [Li [p]] ; */
00822       i = Li [p] ;
00823       MULTSUBCONJ (yx,yz,0, Lx,Lz,p, Xx,Xz,i) ;
00824   }
00825
00826 #ifdef LL
00827   /* X [j] = y / d ; */
00828   DIV_REAL (Xx,Xz,j, yx,yz,0, d,0) ;
00829 #else
00830   /* X [j] = y ; */
00831   ASSIGN (Xx,Xz,j, yx,yz,0) ;
00832 #endif
00833
00834     }
00835
00836 #endif
00837 }
00838
00839 /* prepare for the next inclusion of this file in cholmod_solve.c */
00840 #undef LL
00841 #undef LD
```