Amesos Package Browser (Single Doxygen Collection) Development
amesos_klu_impl.c
Go to the documentation of this file.
```00001 /* ========================================================================== */
00002 /* === klu ================================================================== */
00003 /* ========================================================================== */
00004
00005 /* KLU: factorizes P*A into L*U, using the Gilbert-Peierls algorithm [1], with
00006  * optional symmetric pruning by Eisenstat and Liu [2].  The code is by Tim
00007  * Davis.  This algorithm is what appears as the default sparse LU routine in
00008  * MATLAB version 6.0, and still appears in MATLAB 6.5 as [L,U,P] = lu (A).
00009  * Note that no column ordering is provided (see COLAMD or AMD for suitable
00010  * orderings).  SuperLU is based on this algorithm, except that it adds the
00011  * use of dense matrix operations on "supernodes" (adjacent columns with
00012  * identical).  This code doesn't use supernodes, thus its name ("Kent" LU,
00013  * as in "Clark Kent", in contrast with Super-LU...).  This algorithm is slower
00014  * than SuperLU and UMFPACK for large matrices with lots of nonzeros in their
00015  * factors (such as for most finite-element problems).  However, for matrices
00016  * with very sparse LU factors, this algorithm is typically faster than both
00017  * SuperLU and UMFPACK, since in this case there is little chance to exploit
00018  * dense matrix kernels (the BLAS).
00019  *
00020  * Only one block of A is factorized, in the BTF form.  The input n is the
00021  * size of the block; k1 is the first row and column in the block.
00022  *
00023  * NOTE: no error checking is done on the inputs.  This version is not meant to
00024  * be called directly by the user.  Use klu_factor instead.
00025  *
00026  * No fill-reducing ordering is provided.  The ordering quality of
00027  * klu_kernel_factor is the responsibility of the caller.  The input A must
00028  * pre-permuted to reduce fill-in, or fill-reducing input permutation Q must
00029  * be provided.
00030  *
00031  * The input matrix A must be in compressed-column form, with either sorted
00032  * or unsorted row indices.  Row indices for column j of A is in
00033  * Ai [Ap [j] ... Ap [j+1]-1] and the same range of indices in Ax holds the
00034  * numerical values.  No duplicate entries are allowed.
00035  *
00037  * file for details on permitted use.  Note that no code from The MathWorks,
00038  * Inc, or from SuperLU, or from any other source appears here.  The code is
00039  * written from scratch, from the algorithmic description in Gilbert & Peierls'
00040  * and Eisenstat & Liu's journal papers [1,2].
00041  *
00042  * If an input permutation Q is provided, the factorization L*U = A (P,Q)
00043  * is computed, where P is determined by partial pivoting, and Q is the input
00044  * ordering.  If the pivot tolerance is less than 1, the "diagonal" entry that
00045  * KLU attempts to choose is the diagonal of A (Q,Q).  In other words, the
00046  * input permutation is applied symmetrically to the input matrix.  The output
00047  * permutation P includes both the partial pivoting ordering and the input
00048  * permutation.  If Q is NULL, then it is assumed to be the identity
00049  * permutation.  Q is not modified.
00050  *
00051  * [1] Gilbert, J. R. and Peierls, T., "Sparse Partial Pivoting in Time
00052  *  Proportional to Arithmetic Operations," SIAM J. Sci. Stat. Comp.,
00053  *  vol 9, pp.  862-874, 1988.
00054  * [2] Eisenstat, S. C. and Liu, J. W. H., "Exploiting Structural Symmetry in
00055  *  Unsymmetric Sparse Symbolic Factorization," SIAM J. Matrix Analysis &
00056  *  Applic., vol 13, pp.  202-211, 1992.
00057  */
00058
00059 /* ========================================================================== */
00060
00061 #include "amesos_klu_internal.h"
00062
00063 size_t KLU_kernel_factor      /* 0 if failure, size of LU if OK */
00064 (
00065     /* inputs, not modified */
00066     Int n,      /* A is n-by-n. n must be > 0. */
00067     Int Ap [ ],     /* size n+1, column pointers for A */
00068     Int Ai [ ],     /* size nz = Ap [n], row indices for A */
00069     Entry Ax [ ],   /* size nz, values of A */
00070     Int Q [ ],      /* size n, optional column permutation */
00071     double Lsize,   /* estimate of number of nonzeros in L */
00072
00073     /* outputs, not defined on input */
00074     Unit **p_LU,  /* row indices and values of L and U */
00075     Entry Udiag [ ],  /* size n, diagonal of U */
00076     Int Llen [ ], /* size n, column length of L */
00077     Int Ulen [ ], /* size n, column length of U */
00078     Int Lip [ ],  /* size n, column pointers for L */
00079     Int Uip [ ],  /* size n, column pointers for U */
00080     Int P [ ],    /* row permutation, size n */
00081     Int *lnz,   /* size of L */
00082     Int *unz,   /* size of U */
00083
00084     /* workspace, undefined on input */
00085     Entry *X,     /* size n double's, zero on output */
00086     Int *Work,      /* size 5n Int's */
00087
00088     /* inputs, not modified on output */
00089     Int k1,   /* the block of A is from k1 to k2-1 */
00090     Int PSinv [ ],    /* inverse of P from symbolic factorization */
00091     double Rs [ ],    /* scale factors for A */
00092
00093     /* inputs, modified on output */
00094     Int Offp [ ],   /* off-diagonal matrix (modified by this routine) */
00095     Int Offi [ ],
00096     Entry Offx [ ],
00097     /* --------------- */
00098     KLU_common *Common
00099 )
00100 {
00101     double maxlnz, dunits ;
00102     Unit *LU ;
00103     Int *Pinv, *Lpend, *Stack, *Flag, *Ap_pos, *W ;
00104     Int lsize, usize, anz, ok ;
00105     size_t lusize ;
00106     ASSERT (Common != NULL) ;
00107
00108     /* ---------------------------------------------------------------------- */
00109     /* get control parameters, or use defaults */
00110     /* ---------------------------------------------------------------------- */
00111
00112     n = MAX (1, n) ;
00113     anz = Ap [n+k1] - Ap [k1] ;
00114
00115     if (Lsize <= 0)
00116     {
00117   Lsize = -Lsize ;
00118   Lsize = MAX (Lsize, 1.0) ;
00119   lsize = (int) ( Lsize * anz + n ) ;
00120     }
00121     else
00122     {
00123   lsize = (int) Lsize ;
00124     }
00125
00126     usize = lsize ;
00127
00128     lsize  = MAX (n+1, lsize) ;
00129     usize  = MAX (n+1, usize) ;
00130
00131     maxlnz = (((double) n) * ((double) n) + ((double) n)) / 2. ;
00132     maxlnz = MIN (maxlnz, ((double) INT_MAX)) ;
00133     lsize  = (int) ( MIN (maxlnz, lsize) ) ;
00134     usize  = (int) ( MIN (maxlnz, usize) ) ;
00135
00136     PRINTF (("Welcome to klu: n %d anz %d k1 %d lsize %d usize %d maxlnz %g\n",
00137   n, anz, k1, lsize, usize, maxlnz)) ;
00138
00139     /* ---------------------------------------------------------------------- */
00140     /* allocate workspace and outputs */
00141     /* ---------------------------------------------------------------------- */
00142
00143     /* return arguments are not yet assigned */
00144     *p_LU = (Unit *) NULL ;
00145
00146     /* these computations are safe from size_t overflow */
00147     W = Work ;
00148     Pinv = (Int *) W ;      W += n ;
00149     Stack = (Int *) W ;     W += n ;
00150     Flag = (Int *) W ;      W += n ;
00151     Lpend = (Int *) W ;     W += n ;
00152     Ap_pos = (Int *) W ;    W += n ;
00153
00154     dunits = DUNITS (Int, lsize) + DUNITS (Entry, lsize) +
00155              DUNITS (Int, usize) + DUNITS (Entry, usize) ;
00156     lusize = (size_t) dunits ;
00157     ok = !INT_OVERFLOW (dunits) ;
00158     LU = ok ? KLU_malloc (lusize, sizeof (Unit), Common) : NULL ;
00159     if (LU == NULL)
00160     {
00161   /* out of memory, or problem too large */
00162   Common->status = KLU_OUT_OF_MEMORY ;
00163   lusize = 0 ;
00164   return (lusize) ;
00165     }
00166
00167     /* ---------------------------------------------------------------------- */
00168     /* factorize */
00169     /* ---------------------------------------------------------------------- */
00170
00171     /* with pruning, and non-recursive depth-first-search */
00172     lusize = KLU_kernel (n, Ap, Ai, Ax, Q, lusize,
00173       Pinv, P, &LU, Udiag, Llen, Ulen, Lip, Uip, lnz, unz,
00174       X, Stack, Flag, Ap_pos, Lpend,
00175       k1, PSinv, Rs, Offp, Offi, Offx, Common) ;
00176
00177     /* ---------------------------------------------------------------------- */
00178     /* return LU factors, or return nothing if an error occurred */
00179     /* ---------------------------------------------------------------------- */
00180
00181     if (Common->status < KLU_OK)
00182     {
00183   LU = KLU_free (LU, lusize, sizeof (Unit), Common) ;
00184   lusize = 0 ;
00185     }
00186     *p_LU = LU ;
00187     PRINTF ((" in klu noffdiag %d\n", Common->noffdiag)) ;
00188     return (lusize) ;
00189 }
00190
00191
00192 /* ========================================================================== */
00193 /* === KLU_lsolve =========================================================== */
00194 /* ========================================================================== */
00195
00196 /* Solve Lx=b.  Assumes L is unit lower triangular and where the unit diagonal
00197  * entry is NOT stored.  Overwrites B  with the solution X.  B is n-by-nrhs
00198  * and is stored in ROW form with row dimension nrhs.  nrhs must be in the
00199  * range 1 to 4. */
00200
00201 void KLU_lsolve
00202 (
00203     /* inputs, not modified: */
00204     Int n,
00205     Int Lip [ ],
00206     Int Llen [ ],
00207     Unit LU [ ],
00208     Int nrhs,
00209     /* right-hand-side on input, solution to Lx=b on output */
00210     Entry X [ ]
00211 )
00212 {
00213     Entry x [4], lik ;
00214     Int *Li ;
00215     Entry *Lx ;
00216     Int k, p, len, i ;
00217
00218     switch (nrhs)
00219     {
00220
00221   case 1:
00222       for (k = 0 ; k < n ; k++)
00223       {
00224     x [0] = X [k] ;
00225     GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
00226     /* unit diagonal of L is not stored*/
00227     for (p = 0 ; p < len ; p++)
00228     {
00229         /* X [Li [p]] -= Lx [p] * x [0] ; */
00230         MULT_SUB (X [Li [p]], Lx [p], x [0]) ;
00231     }
00232       }
00233       break ;
00234
00235   case 2:
00236
00237       for (k = 0 ; k < n ; k++)
00238       {
00239     x [0] = X [2*k    ] ;
00240     x [1] = X [2*k + 1] ;
00241     GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
00242     for (p = 0 ; p < len ; p++)
00243     {
00244         i = Li [p] ;
00245         lik = Lx [p] ;
00246         MULT_SUB (X [2*i], lik, x [0]) ;
00247         MULT_SUB (X [2*i + 1], lik, x [1]) ;
00248     }
00249       }
00250       break ;
00251
00252   case 3:
00253
00254       for (k = 0 ; k < n ; k++)
00255       {
00256     x [0] = X [3*k    ] ;
00257     x [1] = X [3*k + 1] ;
00258     x [2] = X [3*k + 2] ;
00259     GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
00260     for (p = 0 ; p < len ; p++)
00261     {
00262         i = Li [p] ;
00263         lik = Lx [p] ;
00264         MULT_SUB (X [3*i], lik, x [0]) ;
00265         MULT_SUB (X [3*i + 1], lik, x [1]) ;
00266         MULT_SUB (X [3*i + 2], lik, x [2]) ;
00267     }
00268       }
00269       break ;
00270
00271   case 4:
00272
00273       for (k = 0 ; k < n ; k++)
00274       {
00275     x [0] = X [4*k    ] ;
00276     x [1] = X [4*k + 1] ;
00277     x [2] = X [4*k + 2] ;
00278     x [3] = X [4*k + 3] ;
00279     GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
00280     for (p = 0 ; p < len ; p++)
00281     {
00282         i = Li [p] ;
00283         lik = Lx [p] ;
00284         MULT_SUB (X [4*i], lik, x [0]) ;
00285         MULT_SUB (X [4*i + 1], lik, x [1]) ;
00286         MULT_SUB (X [4*i + 2], lik, x [2]) ;
00287         MULT_SUB (X [4*i + 3], lik, x [3]) ;
00288     }
00289       }
00290       break ;
00291
00292     }
00293 }
00294
00295 /* ========================================================================== */
00296 /* === KLU_usolve =========================================================== */
00297 /* ========================================================================== */
00298
00299 /* Solve Ux=b.  Assumes U is non-unit upper triangular and where the diagonal
00300  * entry is NOT stored.  Overwrites B with the solution X.  B is n-by-nrhs
00301  * and is stored in ROW form with row dimension nrhs.  nrhs must be in the
00302  * range 1 to 4. */
00303
00304 void KLU_usolve
00305 (
00306     /* inputs, not modified: */
00307     Int n,
00308     Int Uip [ ],
00309     Int Ulen [ ],
00310     Unit LU [ ],
00311     Entry Udiag [ ],
00312     Int nrhs,
00313     /* right-hand-side on input, solution to Ux=b on output */
00314     Entry X [ ]
00315 )
00316 {
00317     Entry x [4], uik, ukk ;
00318     Int *Ui ;
00319     Entry *Ux ;
00320     Int k, p, len, i ;
00321
00322     switch (nrhs)
00323     {
00324
00325   case 1:
00326
00327       for (k = n-1 ; k >= 0 ; k--)
00328       {
00329     GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
00330     /* x [0] = X [k] / Udiag [k] ; */
00331     DIV (x [0], X [k], Udiag [k]) ;
00332     X [k] = x [0] ;
00333     for (p = 0 ; p < len ; p++)
00334     {
00335         /* X [Ui [p]] -= Ux [p] * x [0] ; */
00336         MULT_SUB (X [Ui [p]], Ux [p], x [0]) ;
00337
00338     }
00339       }
00340
00341       break ;
00342
00343   case 2:
00344
00345       for (k = n-1 ; k >= 0 ; k--)
00346       {
00347     GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
00348     ukk = Udiag [k] ;
00349     /* x [0] = X [2*k    ] / ukk ;
00350     x [1] = X [2*k + 1] / ukk ; */
00351     DIV (x [0], X [2*k], ukk) ;
00352     DIV (x [1], X [2*k + 1], ukk) ;
00353
00354     X [2*k    ] = x [0] ;
00355     X [2*k + 1] = x [1] ;
00356     for (p = 0 ; p < len ; p++)
00357     {
00358         i = Ui [p] ;
00359         uik = Ux [p] ;
00360         /* X [2*i    ] -= uik * x [0] ;
00361         X [2*i + 1] -= uik * x [1] ; */
00362         MULT_SUB (X [2*i], uik, x [0]) ;
00363         MULT_SUB (X [2*i + 1], uik, x [1]) ;
00364     }
00365       }
00366
00367       break ;
00368
00369   case 3:
00370
00371       for (k = n-1 ; k >= 0 ; k--)
00372       {
00373     GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
00374     ukk = Udiag [k] ;
00375
00376     DIV (x [0], X [3*k], ukk) ;
00377     DIV (x [1], X [3*k + 1], ukk) ;
00378     DIV (x [2], X [3*k + 2], ukk) ;
00379
00380     X [3*k    ] = x [0] ;
00381     X [3*k + 1] = x [1] ;
00382     X [3*k + 2] = x [2] ;
00383     for (p = 0 ; p < len ; p++)
00384     {
00385         i = Ui [p] ;
00386         uik = Ux [p] ;
00387         MULT_SUB (X [3*i], uik, x [0]) ;
00388         MULT_SUB (X [3*i + 1], uik, x [1]) ;
00389         MULT_SUB (X [3*i + 2], uik, x [2]) ;
00390     }
00391       }
00392
00393       break ;
00394
00395   case 4:
00396
00397       for (k = n-1 ; k >= 0 ; k--)
00398       {
00399     GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
00400     ukk = Udiag [k] ;
00401
00402     DIV (x [0], X [4*k], ukk) ;
00403     DIV (x [1], X [4*k + 1], ukk) ;
00404     DIV (x [2], X [4*k + 2], ukk) ;
00405     DIV (x [3], X [4*k + 3], ukk) ;
00406
00407     X [4*k    ] = x [0] ;
00408     X [4*k + 1] = x [1] ;
00409     X [4*k + 2] = x [2] ;
00410     X [4*k + 3] = x [3] ;
00411     for (p = 0 ; p < len ; p++)
00412     {
00413         i = Ui [p] ;
00414         uik = Ux [p] ;
00415
00416         MULT_SUB (X [4*i], uik, x [0]) ;
00417         MULT_SUB (X [4*i + 1], uik, x [1]) ;
00418         MULT_SUB (X [4*i + 2], uik, x [2]) ;
00419         MULT_SUB (X [4*i + 3], uik, x [3]) ;
00420     }
00421       }
00422
00423       break ;
00424
00425     }
00426 }
00427
00428
00429 /* ========================================================================== */
00430 /* === KLU_ltsolve ========================================================== */
00431 /* ========================================================================== */
00432
00433 /* Solve L'x=b.  Assumes L is unit lower triangular and where the unit diagonal
00434  * entry is NOT stored.  Overwrites B with the solution X.  B is n-by-nrhs
00435  * and is stored in ROW form with row dimension nrhs.  nrhs must in the
00436  * range 1 to 4. */
00437
00438 void KLU_ltsolve
00439 (
00440     /* inputs, not modified: */
00441     Int n,
00442     Int Lip [ ],
00443     Int Llen [ ],
00444     Unit LU [ ],
00445     Int nrhs,
00446 #ifdef COMPLEX
00447     Int conj_solve,
00448 #endif
00449     /* right-hand-side on input, solution to L'x=b on output */
00450     Entry X [ ]
00451 )
00452 {
00453     Entry x [4], lik ;
00454     Int *Li ;
00455     Entry *Lx ;
00456     Int k, p, len, i ;
00457
00458     switch (nrhs)
00459     {
00460
00461   case 1:
00462
00463       for (k = n-1 ; k >= 0 ; k--)
00464       {
00465     GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
00466     x [0] = X [k] ;
00467     for (p = 0 ; p < len ; p++)
00468     {
00469 #ifdef COMPLEX
00470         if (conj_solve)
00471         {
00472       /* x [0] -= CONJ (Lx [p]) * X [Li [p]] ; */
00473       MULT_SUB_CONJ (x [0], X [Li [p]], Lx [p]) ;
00474         }
00475         else
00476 #endif
00477         {
00478       /*x [0] -= Lx [p] * X [Li [p]] ;*/
00479       MULT_SUB (x [0], Lx [p], X [Li [p]]) ;
00480         }
00481     }
00482     X [k] = x [0] ;
00483       }
00484       break ;
00485
00486   case 2:
00487
00488       for (k = n-1 ; k >= 0 ; k--)
00489       {
00490     x [0] = X [2*k    ] ;
00491     x [1] = X [2*k + 1] ;
00492     GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
00493     for (p = 0 ; p < len ; p++)
00494     {
00495         i = Li [p] ;
00496 #ifdef COMPLEX
00497         if (conj_solve)
00498         {
00499       CONJ (lik, Lx [p]) ;
00500         }
00501         else
00502 #endif
00503         {
00504       lik = Lx [p] ;
00505         }
00506         MULT_SUB (x [0], lik, X [2*i]) ;
00507         MULT_SUB (x [1], lik, X [2*i + 1]) ;
00508     }
00509     X [2*k    ] = x [0] ;
00510     X [2*k + 1] = x [1] ;
00511       }
00512       break ;
00513
00514   case 3:
00515
00516       for (k = n-1 ; k >= 0 ; k--)
00517       {
00518     x [0] = X [3*k    ] ;
00519     x [1] = X [3*k + 1] ;
00520     x [2] = X [3*k + 2] ;
00521     GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
00522     for (p = 0 ; p < len ; p++)
00523     {
00524         i = Li [p] ;
00525 #ifdef COMPLEX
00526         if (conj_solve)
00527         {
00528       CONJ (lik, Lx [p]) ;
00529         }
00530         else
00531 #endif
00532         {
00533       lik = Lx [p] ;
00534         }
00535         MULT_SUB (x [0], lik, X [3*i]) ;
00536         MULT_SUB (x [1], lik, X [3*i + 1]) ;
00537         MULT_SUB (x [2], lik, X [3*i + 2]) ;
00538     }
00539     X [3*k    ] = x [0] ;
00540     X [3*k + 1] = x [1] ;
00541     X [3*k + 2] = x [2] ;
00542       }
00543       break ;
00544
00545   case 4:
00546
00547       for (k = n-1 ; k >= 0 ; k--)
00548       {
00549     x [0] = X [4*k    ] ;
00550     x [1] = X [4*k + 1] ;
00551     x [2] = X [4*k + 2] ;
00552     x [3] = X [4*k + 3] ;
00553     GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
00554     for (p = 0 ; p < len ; p++)
00555     {
00556         i = Li [p] ;
00557 #ifdef COMPLEX
00558         if (conj_solve)
00559         {
00560       CONJ (lik, Lx [p]) ;
00561         }
00562         else
00563 #endif
00564         {
00565       lik = Lx [p] ;
00566         }
00567         MULT_SUB (x [0], lik, X [4*i]) ;
00568         MULT_SUB (x [1], lik, X [4*i + 1]) ;
00569         MULT_SUB (x [2], lik, X [4*i + 2]) ;
00570         MULT_SUB (x [3], lik, X [4*i + 3]) ;
00571     }
00572     X [4*k    ] = x [0] ;
00573     X [4*k + 1] = x [1] ;
00574     X [4*k + 2] = x [2] ;
00575     X [4*k + 3] = x [3] ;
00576       }
00577       break ;
00578     }
00579 }
00580
00581
00582 /* ========================================================================== */
00583 /* === KLU_utsolve ========================================================== */
00584 /* ========================================================================== */
00585
00586 /* Solve U'x=b.  Assumes U is non-unit upper triangular and where the diagonal
00587  * entry is stored (and appears last in each column of U).  Overwrites B
00588  * with the solution X.  B is n-by-nrhs and is stored in ROW form with row
00589  * dimension nrhs.  nrhs must be in the range 1 to 4. */
00590
00591 void KLU_utsolve
00592 (
00593     /* inputs, not modified: */
00594     Int n,
00595     Int Uip [ ],
00596     Int Ulen [ ],
00597     Unit LU [ ],
00598     Entry Udiag [ ],
00599     Int nrhs,
00600 #ifdef COMPLEX
00601     Int conj_solve,
00602 #endif
00603     /* right-hand-side on input, solution to Ux=b on output */
00604     Entry X [ ]
00605 )
00606 {
00607     Entry x [4], uik, ukk ;
00608     Int k, p, len, i ;
00609     Int *Ui ;
00610     Entry *Ux ;
00611
00612     switch (nrhs)
00613     {
00614
00615   case 1:
00616
00617       for (k = 0 ; k < n ; k++)
00618       {
00619     GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
00620     x [0] = X [k] ;
00621     for (p = 0 ; p < len ; p++)
00622     {
00623 #ifdef COMPLEX
00624         if (conj_solve)
00625         {
00626       /* x [0] -= CONJ (Ux [p]) * X [Ui [p]] ; */
00627       MULT_SUB_CONJ (x [0], X [Ui [p]], Ux [p]) ;
00628         }
00629         else
00630 #endif
00631         {
00632       /* x [0] -= Ux [p] * X [Ui [p]] ; */
00633       MULT_SUB (x [0], Ux [p], X [Ui [p]]) ;
00634         }
00635     }
00636 #ifdef COMPLEX
00637     if (conj_solve)
00638     {
00639         CONJ (ukk, Udiag [k]) ;
00640     }
00641     else
00642 #endif
00643     {
00644         ukk = Udiag [k] ;
00645     }
00646     DIV (X [k], x [0], ukk) ;
00647       }
00648       break ;
00649
00650   case 2:
00651
00652       for (k = 0 ; k < n ; k++)
00653       {
00654     GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
00655     x [0] = X [2*k    ] ;
00656     x [1] = X [2*k + 1] ;
00657     for (p = 0 ; p < len ; p++)
00658     {
00659         i = Ui [p] ;
00660 #ifdef COMPLEX
00661         if (conj_solve)
00662         {
00663       CONJ (uik, Ux [p]) ;
00664         }
00665         else
00666 #endif
00667         {
00668       uik = Ux [p] ;
00669         }
00670         MULT_SUB (x [0], uik, X [2*i]) ;
00671         MULT_SUB (x [1], uik, X [2*i + 1]) ;
00672     }
00673 #ifdef COMPLEX
00674     if (conj_solve)
00675     {
00676         CONJ (ukk, Udiag [k]) ;
00677     }
00678     else
00679 #endif
00680     {
00681         ukk = Udiag [k] ;
00682     }
00683     DIV (X [2*k], x [0], ukk) ;
00684     DIV (X [2*k + 1], x [1], ukk) ;
00685       }
00686       break ;
00687
00688   case 3:
00689
00690       for (k = 0 ; k < n ; k++)
00691       {
00692     GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
00693     x [0] = X [3*k    ] ;
00694     x [1] = X [3*k + 1] ;
00695     x [2] = X [3*k + 2] ;
00696     for (p = 0 ; p < len ; p++)
00697     {
00698         i = Ui [p] ;
00699 #ifdef COMPLEX
00700         if (conj_solve)
00701         {
00702       CONJ (uik, Ux [p]) ;
00703         }
00704         else
00705 #endif
00706         {
00707       uik = Ux [p] ;
00708         }
00709         MULT_SUB (x [0], uik, X [3*i]) ;
00710         MULT_SUB (x [1], uik, X [3*i + 1]) ;
00711         MULT_SUB (x [2], uik, X [3*i + 2]) ;
00712     }
00713 #ifdef COMPLEX
00714     if (conj_solve)
00715     {
00716         CONJ (ukk, Udiag [k]) ;
00717     }
00718     else
00719 #endif
00720     {
00721         ukk = Udiag [k] ;
00722     }
00723     DIV (X [3*k], x [0], ukk) ;
00724     DIV (X [3*k + 1], x [1], ukk) ;
00725     DIV (X [3*k + 2], x [2], ukk) ;
00726       }
00727       break ;
00728
00729   case 4:
00730
00731       for (k = 0 ; k < n ; k++)
00732       {
00733     GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
00734     x [0] = X [4*k    ] ;
00735     x [1] = X [4*k + 1] ;
00736     x [2] = X [4*k + 2] ;
00737     x [3] = X [4*k + 3] ;
00738     for (p = 0 ; p < len ; p++)
00739     {
00740         i = Ui [p] ;
00741 #ifdef COMPLEX
00742         if (conj_solve)
00743         {
00744       CONJ (uik, Ux [p]) ;
00745         }
00746         else
00747 #endif
00748         {
00749       uik = Ux [p] ;
00750         }
00751         MULT_SUB (x [0], uik, X [4*i]) ;
00752         MULT_SUB (x [1], uik, X [4*i + 1]) ;
00753         MULT_SUB (x [2], uik, X [4*i + 2]) ;
00754         MULT_SUB (x [3], uik, X [4*i + 3]) ;
00755     }
00756 #ifdef COMPLEX
00757     if (conj_solve)
00758     {
00759         CONJ (ukk, Udiag [k]) ;
00760     }
00761     else
00762 #endif
00763     {
00764         ukk = Udiag [k] ;
00765     }
00766     DIV (X [4*k], x [0], ukk) ;
00767     DIV (X [4*k + 1], x [1], ukk) ;
00768     DIV (X [4*k + 2], x [2], ukk) ;
00769     DIV (X [4*k + 3], x [3], ukk) ;
00770       }
00771       break ;
00772     }
00773 }
```