/* ========================================================================== */
/* === klu ================================================================== */
/* ========================================================================== */
/* KLU: factorizes P*A into L*U, using the Gilbert-Peierls algorithm, with
* optional symmetric pruning by Eisenstat and Liu. The code is by Tim Davis.
* This algorithm is what appears as the default sparse LU routine in MATLAB
* version 6.0, and still appears in MATLAB 6.5 as [L,U,P] = lu (A). Note
* that no column ordering is provided (see COLAMD or AMD for suitable
* orderings). SuperLU is based on this algorithm, except that it adds the
* use of dense matrix operations on "supernodes" (adjacent columns with
* identical). This code doesn't use supernodes, thus its name ("Kent" LU,
* as in "Clark Kent", in contrast with Super-LU...). This algorithm is slower
* than SuperLU and UMFPACK for large matrices with lots of nonzeros in their
* factors (such as for most finite-element problems). However, for matrices
* with very sparse LU factors, this algorithm is typically faster than both
* SuperLU and UMFPACK, since in this case there is little chance to exploit
* dense matrix kernels (the BLAS).
*
* NOTE: no error checking is done on the inputs. This version is not meant to
* be called directly by the user. Use klu_factor instead.
*
* No fill-reducing ordering is provided. The ordering quality of
* klu_kernel_factor is your responsibility. You must pre-permute A to reduce
* fill-in, or provide a fill-reducing input permutation Q.
*
* The input matrix A must be in compressed-column form, with either sorted
* or unsorted row indices. Row indices for column j of A is in
* Ai [Ap [j] ... Ap [j+1]-1] and the same range of indices in Ax holds the
* numerical values. No duplicate entries are allowed.
*
* Copyright August 2004, Tim Davis. All rights reserved. See the README
* file for details on permitted use. Note that no code from The MathWorks,
* Inc, or from SuperLU, or from any other source appears here. The code is
* written from scratch, from the algorithmic description in Gilbert & Peierls'
* and Eisenstat & Liu's journal papers.
*
* If an input permutation Q is provided, the factorization L*U = A (P,Q)
* is computed, where P is determined by partial pivoting, and Q is the input
* ordering. If the pivot tolerance is less than 1, the "diagonal" entry that
* KLU attempts to choose is the diagonal of A (Q,Q). In other words, the
* input permutation is applied symmetrically to the input matrix. The output
* permutation P includes both the partial pivoting ordering and the input
* permutation. If Q is (int *) NULL, then it is assumed to be the identity
* permutation. Q is not modified.
*
* TODO: Cite Gilbert, Eisenstat here. Write README file.
*/
/* ========================================================================== */
#include "klu_internal.h"
int KLU_kernel_factor
(
/* inputs, not modified */
int n, /* A is n-by-n. n must be > 0. */
int Ap [ ], /* size n+1, column pointers for A */
int Ai [ ], /* size nz = Ap [n], row indices for A */
Entry Ax [ ], /* size nz, values of A */
int Q [ ], /* size n, optional column permutation */
double Lsize,
/* outputs, not defined on input */
Unit **p_LU, /* row indices and values of L and U*/
Unit Udiag [ ], /* size n, diagonal of U */
int Llen [ ], /* size n, column length of L */
int Ulen [ ], /* size n, column length of U */
int Lip [ ], /* size n+1, column pointers for L*/
int Uip [ ], /* size n+1, column pointers for U*/
int P [ ], /* row permutation, size n */
int *lnz, /* size of L*/
int *unz, /* size of U*/
/* workspace, undefined on input */
Entry *X, /* size n double's, zero on output */
int *Work, /* size 5n int's */
/* ---- the following are only used in the BTF case --- */
/* inputs, not modified on output */
int k1, /* the block of A is from k1 to k2-1 */
int PSinv [ ], /* inverse of P from symbolic factorization */
double Rs [ ], /* scale factors for A */
/* inputs, modified on output */
int Offp [ ], /* off-diagonal matrix (modified by this routine) */
int Offi [ ],
Entry Offx [ ],
/* --------------- */
klu_common *Common
)
{
double maxlnz, dunits ;
Unit *LU ;
int *Pinv, *Lpend, *Stack, *Flag, *Ap_pos, *W ;
int lsize, usize, anz, no_btf ;
size_t lusize ;
/* ---------------------------------------------------------------------- */
/* get control parameters, or use defaults */
/* ---------------------------------------------------------------------- */
if (Common == NULL)
{
return (FALSE) ;
}
no_btf = (Offp == (int *) NULL) ;
n = MAX (1, n) ;
if (no_btf)
{
anz = Ap [n] ;
}
else
{
anz = Ap [n+k1] - Ap [k1] ;
}
if (Lsize <= 0)
{
Lsize = -Lsize ;
Lsize = MAX (Lsize, 1.0) ;
lsize = Lsize * anz + n ;
}
else
{
lsize = Lsize ;
}
usize = lsize ;
lsize = MAX (n+1, lsize) ;
usize = MAX (n+1, usize) ;
maxlnz = (((double) n) * ((double) n) + ((double) n)) / 2. ;
maxlnz = MIN (maxlnz, ((double) INT_MAX)) ;
lsize = MIN (maxlnz, lsize) ;
usize = MIN (maxlnz, usize) ;
PRINTF (("Welcome to klu: n %d anz %d k1 %d lsize %d usize %d maxlnz %g\n",
n, anz, k1, lsize, usize, maxlnz)) ;
/* ---------------------------------------------------------------------- */
/* allocate workspace and outputs */
/* ---------------------------------------------------------------------- */
/* return arguments are not yet assigned */
*p_LU = (Unit *) NULL ;
/* these computations are safe from size_t overflow */
W = Work ;
Pinv = (int *) W ; W += n ;
Stack = (int *) W ; W += n ;
Flag = (int *) W ; W += n ;
Lpend = (int *) W ; W += n ;
Ap_pos = (int *) W ; W += n ;
dunits = DUNITS (int, lsize) + DUNITS (Entry, lsize) +
DUNITS (int, usize) + DUNITS (Entry, usize) ;
if (INT_OVERFLOW (dunits))
{
/* problem is too large, lusize is larger than INT_MAX */
Common->status = KLU_TOO_LARGE ;
return (FALSE) ;
}
lusize = dunits ;
LU = klu_malloc (lusize, sizeof (Unit), Common) ;
if (Common->status < KLU_OK)
{
return (FALSE) ;
}
/* ---------------------------------------------------------------------- */
/* factorize */
/* ---------------------------------------------------------------------- */
/* with pruning, and non-recursive depth-first-search */
KLU_kernel (n, Ap, Ai, Ax, Q, lusize,
Pinv, P, &LU, Udiag, Llen, Ulen, Lip, Uip, lnz, unz,
X, Stack, Flag, Ap_pos, Lpend,
/* BTF and scaling case: */
no_btf, k1, PSinv, Rs, Offp, Offi, Offx, Common) ;
/* ---------------------------------------------------------------------- */
/* return LU factors, or return nothing if an error occurred */
/* ---------------------------------------------------------------------- */
if (Common->status < KLU_OK)
{
LU = klu_free (LU, Common) ;
}
*p_LU = LU ;
PRINTF ((" in klu noffdiag %d\n", Common->noffdiag)) ;
return (TRUE) ;
}
/* ========================================================================== */
/* === klu_lsolve =========================================================== */
/* ========================================================================== */
/* Solve Lx=b. Assumes L is unit lower triangular and where the unit diagonal
* entry is NOT stored. Overwrites B with the solution X. B is n-by-nrhs
* and is stored in ROW form with row dimension nrhs. nrhs must be in the
* range 1 to 4. */
void KLU_lsolve
(
/* inputs, not modified: */
int n,
int Lip [ ],
int Llen [ ],
Unit LU [ ],
int nrhs,
/* right-hand-side on input, solution to Lx=b on output */
Entry X [ ]
)
{
Entry x [4], lik ;
int *Li ;
Entry *Lx ;
int k, p, len, i ;
switch (nrhs)
{
case 1:
for (k = 0 ; k < n ; k++)
{
x [0] = X [k] ;
GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
/* unit diagonal of L is not stored*/
for (p = 0 ; p < len ; p++)
{
/* X [Li [p]] -= Lx [p] * x [0] ; */
MULT_SUB (X [Li [p]], Lx [p], x [0]) ;
}
}
break ;
case 2:
for (k = 0 ; k < n ; k++)
{
x [0] = X [2*k ] ;
x [1] = X [2*k + 1] ;
GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
for (p = 0 ; p < len ; p++)
{
i = Li [p] ;
lik = Lx [p] ;
MULT_SUB (X [2*i], lik, x [0]) ;
MULT_SUB (X [2*i + 1], lik, x [1]) ;
}
}
break ;
case 3:
for (k = 0 ; k < n ; k++)
{
x [0] = X [3*k ] ;
x [1] = X [3*k + 1] ;
x [2] = X [3*k + 2] ;
GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
for (p = 0 ; p < len ; p++)
{
i = Li [p] ;
lik = Lx [p] ;
MULT_SUB (X [3*i], lik, x [0]) ;
MULT_SUB (X [3*i + 1], lik, x [1]) ;
MULT_SUB (X [3*i + 2], lik, x [2]) ;
}
}
break ;
case 4:
for (k = 0 ; k < n ; k++)
{
x [0] = X [4*k ] ;
x [1] = X [4*k + 1] ;
x [2] = X [4*k + 2] ;
x [3] = X [4*k + 3] ;
GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
for (p = 0 ; p < len ; p++)
{
i = Li [p] ;
lik = Lx [p] ;
MULT_SUB (X [4*i], lik, x [0]) ;
MULT_SUB (X [4*i + 1], lik, x [1]) ;
MULT_SUB (X [4*i + 2], lik, x [2]) ;
MULT_SUB (X [4*i + 3], lik, x [3]) ;
}
}
break ;
}
}
/* ========================================================================== */
/* === klu_usolve =========================================================== */
/* ========================================================================== */
/* Solve Ux=b. Assumes U is non-unit upper triangular and where the diagonal
* entry is NOT stored. Overwrites B with the solution X. B is n-by-nrhs
* and is stored in ROW form with row dimension nrhs. nrhs must be in the
* range 1 to 4. */
void KLU_usolve
(
/* inputs, not modified: */
int n,
int Uip [ ],
int Ulen [ ],
Unit LU [ ],
Unit Udiag [ ],
int nrhs,
/* right-hand-side on input, solution to Ux=b on output */
Entry X [ ]
)
{
Entry x [4], uik, ukk ;
int *Ui ;
Entry *Ux, *Udiag_entry ;
int k, p, len, i ;
Udiag_entry = (Entry *) Udiag ;
switch (nrhs)
{
case 1:
for (k = n-1 ; k >= 0 ; k--)
{
GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
/* x [0] = X [k] / Udiag_entry [k] ; */
DIV (x [0], X [k], Udiag_entry [k]) ;
X [k] = x [0] ;
for (p = 0 ; p < len ; p++)
{
/* X [Ui [p]] -= Ux [p] * x [0] ; */
MULT_SUB (X [Ui [p]], Ux [p], x [0]) ;
}
}
break ;
case 2:
for (k = n-1 ; k >= 0 ; k--)
{
GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
ukk = Udiag_entry [k] ;
/* x [0] = X [2*k ] / ukk ;
x [1] = X [2*k + 1] / ukk ; */
DIV (x [0], X [2*k], ukk) ;
DIV (x [1], X [2*k + 1], ukk) ;
X [2*k ] = x [0] ;
X [2*k + 1] = x [1] ;
for (p = 0 ; p < len ; p++)
{
i = Ui [p] ;
uik = Ux [p] ;
/* X [2*i ] -= uik * x [0] ;
X [2*i + 1] -= uik * x [1] ; */
MULT_SUB (X [2*i], uik, x [0]) ;
MULT_SUB (X [2*i + 1], uik, x [1]) ;
}
}
break ;
case 3:
for (k = n-1 ; k >= 0 ; k--)
{
GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
ukk = Udiag_entry [k] ;
DIV (x [0], X [3*k], ukk) ;
DIV (x [1], X [3*k + 1], ukk) ;
DIV (x [2], X [3*k + 2], ukk) ;
X [3*k ] = x [0] ;
X [3*k + 1] = x [1] ;
X [3*k + 2] = x [2] ;
for (p = 0 ; p < len ; p++)
{
i = Ui [p] ;
uik = Ux [p] ;
MULT_SUB (X [3*i], uik, x [0]) ;
MULT_SUB (X [3*i + 1], uik, x [1]) ;
MULT_SUB (X [3*i + 2], uik, x [2]) ;
}
}
break ;
case 4:
for (k = n-1 ; k >= 0 ; k--)
{
GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
ukk = Udiag_entry [k] ;
DIV (x [0], X [4*k], ukk) ;
DIV (x [1], X [4*k + 1], ukk) ;
DIV (x [2], X [4*k + 2], ukk) ;
DIV (x [3], X [4*k + 3], ukk) ;
X [4*k ] = x [0] ;
X [4*k + 1] = x [1] ;
X [4*k + 2] = x [2] ;
X [4*k + 3] = x [3] ;
for (p = 0 ; p < len ; p++)
{
i = Ui [p] ;
uik = Ux [p] ;
MULT_SUB (X [4*i], uik, x [0]) ;
MULT_SUB (X [4*i + 1], uik, x [1]) ;
MULT_SUB (X [4*i + 2], uik, x [2]) ;
MULT_SUB (X [4*i + 3], uik, x [3]) ;
}
}
break ;
}
}
/* ========================================================================== */
/* === klu_ltsolve ========================================================== */
/* ========================================================================== */
/* Solve L'x=b. Assumes L is unit lower triangular and where the unit diagonal
* entry is NOT stored. Overwrites B with the solution X. B is n-by-nrhs
* and is stored in ROW form with row dimension nrhs. nrhs must in the
* range 1 to 4. */
void KLU_ltsolve
(
/* inputs, not modified: */
int n,
int Lip [ ],
int Llen [ ],
Unit LU [ ],
int nrhs,
#ifdef COMPLEX
int conj_solve,
#endif
/* right-hand-side on input, solution to L'x=b on output */
Entry X [ ]
)
{
Entry x [4], lik ;
int *Li ;
Entry *Lx ;
int k, p, len, i ;
switch (nrhs)
{
case 1:
for (k = n-1 ; k >= 0 ; k--)
{
GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
x [0] = X [k] ;
for (p = 0 ; p < len ; p++)
{
#ifdef COMPLEX
if (conj_solve)
{
/* x [0] -= CONJ (Lx [p]) * X [Li [p]] ; */
MULT_SUB_CONJ (x [0], X [Li [p]], Lx [p]) ;
}
else
#endif
{
/*x [0] -= Lx [p] * X [Li [p]] ;*/
MULT_SUB (x [0], Lx [p], X [Li [p]]) ;
}
}
X [k] = x [0] ;
}
break ;
case 2:
for (k = n-1 ; k >= 0 ; k--)
{
x [0] = X [2*k ] ;
x [1] = X [2*k + 1] ;
GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
for (p = 0 ; p < len ; p++)
{
i = Li [p] ;
#ifdef COMPLEX
if (conj_solve)
{
CONJ (lik, Lx [p]) ;
}
else
#endif
{
lik = Lx [p] ;
}
MULT_SUB (x [0], lik, X [2*i]) ;
MULT_SUB (x [1], lik, X [2*i + 1]) ;
}
X [2*k ] = x [0] ;
X [2*k + 1] = x [1] ;
}
break ;
case 3:
for (k = n-1 ; k >= 0 ; k--)
{
x [0] = X [3*k ] ;
x [1] = X [3*k + 1] ;
x [2] = X [3*k + 2] ;
GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
for (p = 0 ; p < len ; p++)
{
i = Li [p] ;
#ifdef COMPLEX
if (conj_solve)
{
CONJ (lik, Lx [p]) ;
}
else
#endif
{
lik = Lx [p] ;
}
MULT_SUB (x [0], lik, X [3*i]) ;
MULT_SUB (x [1], lik, X [3*i + 1]) ;
MULT_SUB (x [2], lik, X [3*i + 2]) ;
}
X [3*k ] = x [0] ;
X [3*k + 1] = x [1] ;
X [3*k + 2] = x [2] ;
}
break ;
case 4:
for (k = n-1 ; k >= 0 ; k--)
{
x [0] = X [4*k ] ;
x [1] = X [4*k + 1] ;
x [2] = X [4*k + 2] ;
x [3] = X [4*k + 3] ;
GET_POINTER (LU, Lip, Llen, Li, Lx, k, len) ;
for (p = 0 ; p < len ; p++)
{
i = Li [p] ;
#ifdef COMPLEX
if (conj_solve)
{
CONJ (lik, Lx [p]) ;
}
else
#endif
{
lik = Lx [p] ;
}
MULT_SUB (x [0], lik, X [4*i]) ;
MULT_SUB (x [1], lik, X [4*i + 1]) ;
MULT_SUB (x [2], lik, X [4*i + 2]) ;
MULT_SUB (x [3], lik, X [4*i + 3]) ;
}
X [4*k ] = x [0] ;
X [4*k + 1] = x [1] ;
X [4*k + 2] = x [2] ;
X [4*k + 3] = x [3] ;
}
break ;
}
}
/* ========================================================================== */
/* === klu_utsolve ========================================================== */
/* ========================================================================== */
/* Solve U'x=b. Assumes U is non-unit upper triangular and where the diagonal
* entry is stored (and appears last in each column of U). Overwrites B
* with the solution X. B is n-by-nrhs and is stored in ROW form with row
* dimension nrhs. nrhs must be in the range 1 to 4. */
void KLU_utsolve
(
/* inputs, not modified: */
int n,
int Uip [ ],
int Ulen [ ],
Unit LU [ ],
Unit Udiag [ ],
int nrhs,
#ifdef COMPLEX
int conj_solve,
#endif
/* right-hand-side on input, solution to Ux=b on output */
Entry X [ ]
)
{
Entry x [4], uik, ukk ;
int k, p, len, i ;
int *Ui ;
Entry *Ux, *Udiag_entry ;
Udiag_entry = (Entry *) Udiag ;
switch (nrhs)
{
case 1:
for (k = 0 ; k < n ; k++)
{
GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
x [0] = X [k] ;
for (p = 0 ; p < len ; p++)
{
#ifdef COMPLEX
if (conj_solve)
{
/* x [0] -= CONJ (Ux [p]) * X [Ui [p]] ; */
MULT_SUB_CONJ (x [0], X [Ui [p]], Ux [p]) ;
}
else
#endif
{
/* x [0] -= Ux [p] * X [Ui [p]] ; */
MULT_SUB (x [0], Ux [p], X [Ui [p]]) ;
}
}
#ifdef COMPLEX
if (conj_solve)
{
CONJ (ukk, Udiag_entry [k]) ;
}
else
#endif
{
ukk = Udiag_entry [k] ;
}
DIV (X [k], x [0], ukk) ;
}
break ;
case 2:
for (k = 0 ; k < n ; k++)
{
GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
x [0] = X [2*k ] ;
x [1] = X [2*k + 1] ;
for (p = 0 ; p < len ; p++)
{
i = Ui [p] ;
#ifdef COMPLEX
if (conj_solve)
{
CONJ (uik, Ux [p]) ;
}
else
#endif
{
uik = Ux [p] ;
}
MULT_SUB (x [0], uik, X [2*i]) ;
MULT_SUB (x [1], uik, X [2*i + 1]) ;
}
#ifdef COMPLEX
if (conj_solve)
{
CONJ (ukk, Udiag_entry [k]) ;
}
else
#endif
{
ukk = Udiag_entry [k] ;
}
DIV (X [2*k], x [0], ukk) ;
DIV (X [2*k + 1], x [1], ukk) ;
}
break ;
case 3:
for (k = 0 ; k < n ; k++)
{
GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
x [0] = X [3*k ] ;
x [1] = X [3*k + 1] ;
x [2] = X [3*k + 2] ;
for (p = 0 ; p < len ; p++)
{
i = Ui [p] ;
#ifdef COMPLEX
if (conj_solve)
{
CONJ (uik, Ux [p]) ;
}
else
#endif
{
uik = Ux [p] ;
}
MULT_SUB (x [0], uik, X [3*i]) ;
MULT_SUB (x [1], uik, X [3*i + 1]) ;
MULT_SUB (x [2], uik, X [3*i + 2]) ;
}
#ifdef COMPLEX
if (conj_solve)
{
CONJ (ukk, Udiag_entry [k]) ;
}
else
#endif
{
ukk = Udiag_entry [k] ;
}
DIV (X [3*k], x [0], ukk) ;
DIV (X [3*k + 1], x [1], ukk) ;
DIV (X [3*k + 2], x [2], ukk) ;
}
break ;
case 4:
for (k = 0 ; k < n ; k++)
{
GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, len) ;
x [0] = X [4*k ] ;
x [1] = X [4*k + 1] ;
x [2] = X [4*k + 2] ;
x [3] = X [4*k + 3] ;
for (p = 0 ; p < len ; p++)
{
i = Ui [p] ;
#ifdef COMPLEX
if (conj_solve)
{
CONJ (uik, Ux [p]) ;
}
else
#endif
{
uik = Ux [p] ;
}
MULT_SUB (x [0], uik, X [4*i]) ;
MULT_SUB (x [1], uik, X [4*i + 1]) ;
MULT_SUB (x [2], uik, X [4*i + 2]) ;
MULT_SUB (x [3], uik, X [4*i + 3]) ;
}
#ifdef COMPLEX
if (conj_solve)
{
CONJ (ukk, Udiag_entry [k]) ;
}
else
#endif
{
ukk = Udiag_entry [k] ;
}
DIV (X [4*k], x [0], ukk) ;
DIV (X [4*k + 1], x [1], ukk) ;
DIV (X [4*k + 2], x [2], ukk) ;
DIV (X [4*k + 3], x [3], ukk) ;
}
break ;
}
}
syntax highlighted by Code2HTML, v. 0.9.1