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