Amesos Package Browser (Single Doxygen Collection) Development
amesos_klu_l.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  *
00036  * Copyright 2004-2007, Tim Davis.  All rights reserved.  See the README
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 /* This file should make the long int version of KLU */
00062 #define DLONG 1
00063 
00064 #include "amesos_klu_internal.h"
00065 
00066 size_t KLU_kernel_factor      /* 0 if failure, size of LU if OK */
00067 (
00068     /* inputs, not modified */
00069     Int n,      /* A is n-by-n. n must be > 0. */
00070     Int Ap [ ],     /* size n+1, column pointers for A */
00071     Int Ai [ ],     /* size nz = Ap [n], row indices for A */
00072     Entry Ax [ ],   /* size nz, values of A */
00073     Int Q [ ],      /* size n, optional column permutation */
00074     double Lsize,   /* estimate of number of nonzeros in L */
00075 
00076     /* outputs, not defined on input */
00077     Unit **p_LU,  /* row indices and values of L and U */
00078     Entry Udiag [ ],  /* size n, diagonal of U */
00079     Int Llen [ ], /* size n, column length of L */
00080     Int Ulen [ ], /* size n, column length of U */
00081     Int Lip [ ],  /* size n, column pointers for L */
00082     Int Uip [ ],  /* size n, column pointers for U */
00083     Int P [ ],    /* row permutation, size n */
00084     Int *lnz,   /* size of L */
00085     Int *unz,   /* size of U */
00086 
00087     /* workspace, undefined on input */
00088     Entry *X,     /* size n double's, zero on output */
00089     Int *Work,      /* size 5n Int's */
00090 
00091     /* inputs, not modified on output */
00092     Int k1,   /* the block of A is from k1 to k2-1 */
00093     Int PSinv [ ],    /* inverse of P from symbolic factorization */
00094     double Rs [ ],    /* scale factors for A */
00095 
00096     /* inputs, modified on output */
00097     Int Offp [ ],   /* off-diagonal matrix (modified by this routine) */
00098     Int Offi [ ],
00099     Entry Offx [ ],
00100     /* --------------- */
00101     KLU_common *Common
00102 )
00103 {
00104     double maxlnz, dunits ;
00105     Unit *LU ;
00106     Int *Pinv, *Lpend, *Stack, *Flag, *Ap_pos, *W ;
00107     Int lsize, usize, anz, ok ;
00108     size_t lusize ;
00109     ASSERT (Common != NULL) ;
00110 
00111     /* ---------------------------------------------------------------------- */
00112     /* get control parameters, or use defaults */
00113     /* ---------------------------------------------------------------------- */
00114 
00115     n = MAX (1, n) ;
00116     anz = Ap [n+k1] - Ap [k1] ;
00117 
00118     if (Lsize <= 0)
00119     {
00120   Lsize = -Lsize ;
00121   Lsize = MAX (Lsize, 1.0) ;
00122   lsize = Lsize * anz + n ;
00123     }
00124     else
00125     {
00126   lsize = Lsize ;
00127     }
00128 
00129     usize = lsize ;
00130 
00131     lsize  = MAX (n+1, lsize) ;
00132     usize  = MAX (n+1, usize) ;
00133 
00134     maxlnz = (((double) n) * ((double) n) + ((double) n)) / 2. ;
00135     maxlnz = MIN (maxlnz, ((double) INT_MAX)) ;
00136     lsize  = MIN (maxlnz, lsize) ;
00137     usize  = MIN (maxlnz, usize) ;
00138 
00139     PRINTF (("Welcome to klu: n %d anz %d k1 %d lsize %d usize %d maxlnz %g\n",
00140   n, anz, k1, lsize, usize, maxlnz)) ;
00141 
00142     /* ---------------------------------------------------------------------- */
00143     /* allocate workspace and outputs */
00144     /* ---------------------------------------------------------------------- */
00145 
00146     /* return arguments are not yet assigned */
00147     *p_LU = (Unit *) NULL ;
00148 
00149     /* these computations are safe from size_t overflow */
00150     W = Work ;
00151     Pinv = (Int *) W ;      W += n ;
00152     Stack = (Int *) W ;     W += n ;
00153     Flag = (Int *) W ;      W += n ;
00154     Lpend = (Int *) W ;     W += n ;
00155     Ap_pos = (Int *) W ;    W += n ;
00156 
00157     dunits = DUNITS (Int, lsize) + DUNITS (Entry, lsize) +
00158              DUNITS (Int, usize) + DUNITS (Entry, usize) ;
00159     lusize = (size_t) dunits ;
00160     ok = !INT_OVERFLOW (dunits) ; 
00161     LU = ok ? KLU_malloc (lusize, sizeof (Unit), Common) : NULL ;
00162     if (LU == NULL)
00163     {
00164   /* out of memory, or problem too large */
00165   Common->status = KLU_OUT_OF_MEMORY ;
00166   lusize = 0 ;
00167   return (lusize) ;
00168     }
00169 
00170     /* ---------------------------------------------------------------------- */
00171     /* factorize */
00172     /* ---------------------------------------------------------------------- */
00173 
00174     /* with pruning, and non-recursive depth-first-search */
00175     lusize = KLU_kernel (n, Ap, Ai, Ax, Q, lusize,
00176       Pinv, P, &LU, Udiag, Llen, Ulen, Lip, Uip, lnz, unz,
00177       X, Stack, Flag, Ap_pos, Lpend,
00178       k1, PSinv, Rs, Offp, Offi, Offx, Common) ;
00179 
00180     /* ---------------------------------------------------------------------- */
00181     /* return LU factors, or return nothing if an error occurred */
00182     /* ---------------------------------------------------------------------- */
00183 
00184     if (Common->status < KLU_OK)
00185     {
00186   LU = KLU_free (LU, lusize, sizeof (Unit), Common) ;
00187   lusize = 0 ;
00188     }
00189     *p_LU = LU ;
00190     PRINTF ((" in klu noffdiag %d\n", Common->noffdiag)) ;
00191     return (lusize) ;
00192 }
00193 
00194 
00195 /* ========================================================================== */
00196 /* === KLU_lsolve =========================================================== */
00197 /* ========================================================================== */
00198 
00199 /* Solve Lx=b.  Assumes L is unit lower triangular and where the unit diagonal
00200  * entry is NOT stored.  Overwrites B  with the solution X.  B is n-by-nrhs
00201  * and is stored in ROW form with row dimension nrhs.  nrhs must be in the
00202  * range 1 to 4. */
00203 
00204 void KLU_lsolve
00205 (
00206     /* inputs, not modified: */
00207     Int n,
00208     Int Lip [ ],
00209     Int Llen [ ],
00210     Unit LU [ ],
00211     Int nrhs,
00212     /* right-hand-side on input, solution to Lx=b on output */
00213     Entry X [ ]
00214 )
00215 {
00216     Entry x [4], lik ;
00217     Int *Li ;
00218     Entry *Lx ;
00219     Int k, p, len, i ;
00220 
00221     switch (nrhs)
00222     {
00223 
00224   case 1:
00225       for (k = 0 ; k < n ; k++)
00226       {
00227     x [0] = X [k] ;
00228     GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
00229     /* unit diagonal of L is not stored*/
00230     for (p = 0 ; p < len ; p++)
00231     {
00232         /* X [Li [p]] -= Lx [p] * x [0] ; */
00233         MULT_SUB (X [Li [p]], Lx [p], x [0]) ;
00234     }
00235       }
00236       break ;
00237 
00238   case 2:
00239 
00240       for (k = 0 ; k < n ; k++)
00241       {
00242     x [0] = X [2*k    ] ;
00243     x [1] = X [2*k + 1] ;
00244     GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
00245     for (p = 0 ; p < len ; p++)
00246     {
00247         i = Li [p] ;
00248         lik = Lx [p] ;
00249         MULT_SUB (X [2*i], lik, x [0]) ;
00250         MULT_SUB (X [2*i + 1], lik, x [1]) ;
00251     }
00252       }
00253       break ;
00254 
00255   case 3:
00256 
00257       for (k = 0 ; k < n ; k++)
00258       {
00259     x [0] = X [3*k    ] ;
00260     x [1] = X [3*k + 1] ;
00261     x [2] = X [3*k + 2] ;
00262     GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
00263     for (p = 0 ; p < len ; p++)
00264     {
00265         i = Li [p] ;
00266         lik = Lx [p] ;
00267         MULT_SUB (X [3*i], lik, x [0]) ;
00268         MULT_SUB (X [3*i + 1], lik, x [1]) ;
00269         MULT_SUB (X [3*i + 2], lik, x [2]) ;
00270     }
00271       }
00272       break ;
00273 
00274   case 4:
00275 
00276       for (k = 0 ; k < n ; k++)
00277       {
00278     x [0] = X [4*k    ] ;
00279     x [1] = X [4*k + 1] ;
00280     x [2] = X [4*k + 2] ;
00281     x [3] = X [4*k + 3] ;
00282     GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
00283     for (p = 0 ; p < len ; p++)
00284     {
00285         i = Li [p] ;
00286         lik = Lx [p] ;
00287         MULT_SUB (X [4*i], lik, x [0]) ;
00288         MULT_SUB (X [4*i + 1], lik, x [1]) ;
00289         MULT_SUB (X [4*i + 2], lik, x [2]) ;
00290         MULT_SUB (X [4*i + 3], lik, x [3]) ;
00291     }
00292       }
00293       break ;
00294 
00295     }
00296 }
00297 
00298 /* ========================================================================== */
00299 /* === KLU_usolve =========================================================== */
00300 /* ========================================================================== */
00301 
00302 /* Solve Ux=b.  Assumes U is non-unit upper triangular and where the diagonal
00303  * entry is NOT stored.  Overwrites B with the solution X.  B is n-by-nrhs
00304  * and is stored in ROW form with row dimension nrhs.  nrhs must be in the
00305  * range 1 to 4. */
00306 
00307 void KLU_usolve
00308 (
00309     /* inputs, not modified: */
00310     Int n,
00311     Int Uip [ ],
00312     Int Ulen [ ],
00313     Unit LU [ ],
00314     Entry Udiag [ ],
00315     Int nrhs,
00316     /* right-hand-side on input, solution to Ux=b on output */
00317     Entry X [ ]
00318 )
00319 {
00320     Entry x [4], uik, ukk ;
00321     Int *Ui ;
00322     Entry *Ux ;
00323     Int k, p, len, i ;
00324 
00325     switch (nrhs)
00326     {
00327 
00328   case 1:
00329 
00330       for (k = n-1 ; k >= 0 ; k--)
00331       {
00332     GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
00333     /* x [0] = X [k] / Udiag [k] ; */
00334     DIV (x [0], X [k], Udiag [k]) ;
00335     X [k] = x [0] ;
00336     for (p = 0 ; p < len ; p++)
00337     {
00338         /* X [Ui [p]] -= Ux [p] * x [0] ; */
00339         MULT_SUB (X [Ui [p]], Ux [p], x [0]) ;
00340 
00341     }
00342       }
00343 
00344       break ;
00345 
00346   case 2:
00347 
00348       for (k = n-1 ; k >= 0 ; k--)
00349       {
00350     GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
00351     ukk = Udiag [k] ;
00352     /* x [0] = X [2*k    ] / ukk ;
00353     x [1] = X [2*k + 1] / ukk ; */
00354     DIV (x [0], X [2*k], ukk) ;
00355     DIV (x [1], X [2*k + 1], ukk) ;
00356 
00357     X [2*k    ] = x [0] ;
00358     X [2*k + 1] = x [1] ;
00359     for (p = 0 ; p < len ; p++)
00360     {
00361         i = Ui [p] ;
00362         uik = Ux [p] ;
00363         /* X [2*i    ] -= uik * x [0] ;
00364         X [2*i + 1] -= uik * x [1] ; */
00365         MULT_SUB (X [2*i], uik, x [0]) ;
00366         MULT_SUB (X [2*i + 1], uik, x [1]) ;
00367     }
00368       }
00369 
00370       break ;
00371 
00372   case 3:
00373 
00374       for (k = n-1 ; k >= 0 ; k--)
00375       {
00376     GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
00377     ukk = Udiag [k] ;
00378 
00379     DIV (x [0], X [3*k], ukk) ;
00380     DIV (x [1], X [3*k + 1], ukk) ;
00381     DIV (x [2], X [3*k + 2], ukk) ;
00382 
00383     X [3*k    ] = x [0] ;
00384     X [3*k + 1] = x [1] ;
00385     X [3*k + 2] = x [2] ;
00386     for (p = 0 ; p < len ; p++)
00387     {
00388         i = Ui [p] ;
00389         uik = Ux [p] ;
00390         MULT_SUB (X [3*i], uik, x [0]) ;
00391         MULT_SUB (X [3*i + 1], uik, x [1]) ;
00392         MULT_SUB (X [3*i + 2], uik, x [2]) ;
00393     }
00394       }
00395 
00396       break ;
00397 
00398   case 4:
00399 
00400       for (k = n-1 ; k >= 0 ; k--)
00401       {
00402     GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
00403     ukk = Udiag [k] ;
00404 
00405     DIV (x [0], X [4*k], ukk) ;
00406     DIV (x [1], X [4*k + 1], ukk) ;
00407     DIV (x [2], X [4*k + 2], ukk) ;
00408     DIV (x [3], X [4*k + 3], ukk) ;
00409 
00410     X [4*k    ] = x [0] ;
00411     X [4*k + 1] = x [1] ;
00412     X [4*k + 2] = x [2] ;
00413     X [4*k + 3] = x [3] ;
00414     for (p = 0 ; p < len ; p++)
00415     {
00416         i = Ui [p] ;
00417         uik = Ux [p] ;
00418 
00419         MULT_SUB (X [4*i], uik, x [0]) ;
00420         MULT_SUB (X [4*i + 1], uik, x [1]) ;
00421         MULT_SUB (X [4*i + 2], uik, x [2]) ;
00422         MULT_SUB (X [4*i + 3], uik, x [3]) ;
00423     }
00424       }
00425 
00426       break ;
00427 
00428     }
00429 }
00430 
00431 
00432 /* ========================================================================== */
00433 /* === KLU_ltsolve ========================================================== */
00434 /* ========================================================================== */
00435 
00436 /* Solve L'x=b.  Assumes L is unit lower triangular and where the unit diagonal
00437  * entry is NOT stored.  Overwrites B with the solution X.  B is n-by-nrhs
00438  * and is stored in ROW form with row dimension nrhs.  nrhs must in the
00439  * range 1 to 4. */
00440 
00441 void KLU_ltsolve
00442 (
00443     /* inputs, not modified: */
00444     Int n,
00445     Int Lip [ ],
00446     Int Llen [ ],
00447     Unit LU [ ],
00448     Int nrhs,
00449 #ifdef COMPLEX
00450     Int conj_solve,
00451 #endif
00452     /* right-hand-side on input, solution to L'x=b on output */
00453     Entry X [ ]
00454 )
00455 {
00456     Entry x [4], lik ;
00457     Int *Li ;
00458     Entry *Lx ;
00459     Int k, p, len, i ;
00460 
00461     switch (nrhs)
00462     {
00463 
00464   case 1:
00465 
00466       for (k = n-1 ; k >= 0 ; k--)
00467       {
00468     GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
00469     x [0] = X [k] ;
00470     for (p = 0 ; p < len ; p++)
00471     {
00472 #ifdef COMPLEX
00473         if (conj_solve)
00474         {
00475       /* x [0] -= CONJ (Lx [p]) * X [Li [p]] ; */
00476       MULT_SUB_CONJ (x [0], X [Li [p]], Lx [p]) ;
00477         }
00478         else
00479 #endif
00480         {
00481       /*x [0] -= Lx [p] * X [Li [p]] ;*/
00482       MULT_SUB (x [0], Lx [p], X [Li [p]]) ;
00483         }
00484     }
00485     X [k] = x [0] ;
00486       }
00487       break ;
00488 
00489   case 2:
00490 
00491       for (k = n-1 ; k >= 0 ; k--)
00492       {
00493     x [0] = X [2*k    ] ;
00494     x [1] = X [2*k + 1] ;
00495     GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
00496     for (p = 0 ; p < len ; p++)
00497     {
00498         i = Li [p] ;
00499 #ifdef COMPLEX
00500         if (conj_solve)
00501         {
00502       CONJ (lik, Lx [p]) ;
00503         }
00504         else
00505 #endif
00506         {
00507       lik = Lx [p] ;
00508         }
00509         MULT_SUB (x [0], lik, X [2*i]) ;
00510         MULT_SUB (x [1], lik, X [2*i + 1]) ;
00511     }
00512     X [2*k    ] = x [0] ;
00513     X [2*k + 1] = x [1] ;
00514       }
00515       break ;
00516 
00517   case 3:
00518 
00519       for (k = n-1 ; k >= 0 ; k--)
00520       {
00521     x [0] = X [3*k    ] ;
00522     x [1] = X [3*k + 1] ;
00523     x [2] = X [3*k + 2] ;
00524     GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
00525     for (p = 0 ; p < len ; p++)
00526     {
00527         i = Li [p] ;
00528 #ifdef COMPLEX
00529         if (conj_solve)
00530         {
00531       CONJ (lik, Lx [p]) ;
00532         }
00533         else
00534 #endif
00535         {
00536       lik = Lx [p] ;
00537         }
00538         MULT_SUB (x [0], lik, X [3*i]) ;
00539         MULT_SUB (x [1], lik, X [3*i + 1]) ;
00540         MULT_SUB (x [2], lik, X [3*i + 2]) ;
00541     }
00542     X [3*k    ] = x [0] ;
00543     X [3*k + 1] = x [1] ;
00544     X [3*k + 2] = x [2] ;
00545       }
00546       break ;
00547 
00548   case 4:
00549 
00550       for (k = n-1 ; k >= 0 ; k--)
00551       {
00552     x [0] = X [4*k    ] ;
00553     x [1] = X [4*k + 1] ;
00554     x [2] = X [4*k + 2] ;
00555     x [3] = X [4*k + 3] ;
00556     GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
00557     for (p = 0 ; p < len ; p++)
00558     {
00559         i = Li [p] ;
00560 #ifdef COMPLEX
00561         if (conj_solve)
00562         {
00563       CONJ (lik, Lx [p]) ;
00564         }
00565         else
00566 #endif
00567         {
00568       lik = Lx [p] ;
00569         }
00570         MULT_SUB (x [0], lik, X [4*i]) ;
00571         MULT_SUB (x [1], lik, X [4*i + 1]) ;
00572         MULT_SUB (x [2], lik, X [4*i + 2]) ;
00573         MULT_SUB (x [3], lik, X [4*i + 3]) ;
00574     }
00575     X [4*k    ] = x [0] ;
00576     X [4*k + 1] = x [1] ;
00577     X [4*k + 2] = x [2] ;
00578     X [4*k + 3] = x [3] ;
00579       }
00580       break ;
00581     }
00582 }
00583 
00584 
00585 /* ========================================================================== */
00586 /* === KLU_utsolve ========================================================== */
00587 /* ========================================================================== */
00588 
00589 /* Solve U'x=b.  Assumes U is non-unit upper triangular and where the diagonal
00590  * entry is stored (and appears last in each column of U).  Overwrites B
00591  * with the solution X.  B is n-by-nrhs and is stored in ROW form with row
00592  * dimension nrhs.  nrhs must be in the range 1 to 4. */
00593 
00594 void KLU_utsolve
00595 (
00596     /* inputs, not modified: */
00597     Int n,
00598     Int Uip [ ],
00599     Int Ulen [ ],
00600     Unit LU [ ],
00601     Entry Udiag [ ],
00602     Int nrhs,
00603 #ifdef COMPLEX
00604     Int conj_solve,
00605 #endif
00606     /* right-hand-side on input, solution to Ux=b on output */
00607     Entry X [ ]
00608 )
00609 {
00610     Entry x [4], uik, ukk ;
00611     Int k, p, len, i ;
00612     Int *Ui ;
00613     Entry *Ux ;
00614 
00615     switch (nrhs)
00616     {
00617 
00618   case 1:
00619 
00620       for (k = 0 ; k < n ; k++)
00621       {
00622     GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
00623     x [0] = X [k] ;
00624     for (p = 0 ; p < len ; p++)
00625     {
00626 #ifdef COMPLEX
00627         if (conj_solve)
00628         {
00629       /* x [0] -= CONJ (Ux [p]) * X [Ui [p]] ; */
00630       MULT_SUB_CONJ (x [0], X [Ui [p]], Ux [p]) ;
00631         }
00632         else
00633 #endif
00634         {
00635       /* x [0] -= Ux [p] * X [Ui [p]] ; */
00636       MULT_SUB (x [0], Ux [p], X [Ui [p]]) ;
00637         }
00638     }
00639 #ifdef COMPLEX
00640     if (conj_solve)
00641     {
00642         CONJ (ukk, Udiag [k]) ;
00643     }
00644     else
00645 #endif
00646     {
00647         ukk = Udiag [k] ;
00648     }
00649     DIV (X [k], x [0], ukk) ;
00650       }
00651       break ;
00652 
00653   case 2:
00654 
00655       for (k = 0 ; k < n ; k++)
00656       {
00657     GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
00658     x [0] = X [2*k    ] ;
00659     x [1] = X [2*k + 1] ;
00660     for (p = 0 ; p < len ; p++)
00661     {
00662         i = Ui [p] ;
00663 #ifdef COMPLEX
00664         if (conj_solve)
00665         {
00666       CONJ (uik, Ux [p]) ;
00667         }
00668         else
00669 #endif
00670         {
00671       uik = Ux [p] ;
00672         }
00673         MULT_SUB (x [0], uik, X [2*i]) ;
00674         MULT_SUB (x [1], uik, X [2*i + 1]) ;
00675     }
00676 #ifdef COMPLEX
00677     if (conj_solve)
00678     {
00679         CONJ (ukk, Udiag [k]) ;
00680     }
00681     else
00682 #endif
00683     {
00684         ukk = Udiag [k] ;
00685     }
00686     DIV (X [2*k], x [0], ukk) ;
00687     DIV (X [2*k + 1], x [1], ukk) ;
00688       }
00689       break ;
00690 
00691   case 3:
00692 
00693       for (k = 0 ; k < n ; k++)
00694       {
00695     GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
00696     x [0] = X [3*k    ] ;
00697     x [1] = X [3*k + 1] ;
00698     x [2] = X [3*k + 2] ;
00699     for (p = 0 ; p < len ; p++)
00700     {
00701         i = Ui [p] ;
00702 #ifdef COMPLEX
00703         if (conj_solve)
00704         {
00705       CONJ (uik, Ux [p]) ;
00706         }
00707         else
00708 #endif
00709         {
00710       uik = Ux [p] ;
00711         }
00712         MULT_SUB (x [0], uik, X [3*i]) ;
00713         MULT_SUB (x [1], uik, X [3*i + 1]) ;
00714         MULT_SUB (x [2], uik, X [3*i + 2]) ;
00715     }
00716 #ifdef COMPLEX
00717     if (conj_solve)
00718     {
00719         CONJ (ukk, Udiag [k]) ;
00720     }
00721     else
00722 #endif
00723     {
00724         ukk = Udiag [k] ;
00725     }
00726     DIV (X [3*k], x [0], ukk) ;
00727     DIV (X [3*k + 1], x [1], ukk) ;
00728     DIV (X [3*k + 2], x [2], ukk) ;
00729       }
00730       break ;
00731 
00732   case 4:
00733 
00734       for (k = 0 ; k < n ; k++)
00735       {
00736     GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
00737     x [0] = X [4*k    ] ;
00738     x [1] = X [4*k + 1] ;
00739     x [2] = X [4*k + 2] ;
00740     x [3] = X [4*k + 3] ;
00741     for (p = 0 ; p < len ; p++)
00742     {
00743         i = Ui [p] ;
00744 #ifdef COMPLEX
00745         if (conj_solve)
00746         {
00747       CONJ (uik, Ux [p]) ;
00748         }
00749         else
00750 #endif
00751         {
00752       uik = Ux [p] ;
00753         }
00754         MULT_SUB (x [0], uik, X [4*i]) ;
00755         MULT_SUB (x [1], uik, X [4*i + 1]) ;
00756         MULT_SUB (x [2], uik, X [4*i + 2]) ;
00757         MULT_SUB (x [3], uik, X [4*i + 3]) ;
00758     }
00759 #ifdef COMPLEX
00760     if (conj_solve)
00761     {
00762         CONJ (ukk, Udiag [k]) ;
00763     }
00764     else
00765 #endif
00766     {
00767         ukk = Udiag [k] ;
00768     }
00769     DIV (X [4*k], x [0], ukk) ;
00770     DIV (X [4*k + 1], x [1], ukk) ;
00771     DIV (X [4*k + 2], x [2], ukk) ;
00772     DIV (X [4*k + 3], x [3], ukk) ;
00773       }
00774       break ;
00775     }
00776 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines