/* ========================================================================== */
/* === Cholesky/cholmod_solve =============================================== */
/* ========================================================================== */
/* -----------------------------------------------------------------------------
* CHOLMOD/Cholesky Module. Copyright (C) 2005-2006, Timothy A. Davis
* The CHOLMOD/Cholesky Module is licensed under Version 2.1 of the GNU
* Lesser General Public License. See lesser.txt for a text of the license.
* CHOLMOD is also available under other licenses; contact authors for details.
* http://www.cise.ufl.edu/research/sparse
* -------------------------------------------------------------------------- */
/* Solve one of the following systems:
*
* Ax=b 0: CHOLMOD_A also applies the permutation L->Perm
* LDL'x=b 1: CHOLMOD_LDLt does not apply L->Perm
* LDx=b 2: CHOLMOD_LD
* DL'x=b 3: CHOLMOD_DLt
* Lx=b 4: CHOLMOD_L
* L'x=b 5: CHOLMOD_Lt
* Dx=b 6: CHOLMOD_D
* x=Pb 7: CHOLMOD_P apply a permutation (P is L->Perm)
* x=P'b 8: CHOLMOD_Pt apply an inverse permutation
*
* The factorization can be simplicial LDL', simplicial LL', or supernodal LL'.
* For an LL' factorization, D is the identity matrix. Thus CHOLMOD_LD and
* CHOLMOD_L solve the same system if an LL' factorization was performed,
* for example.
*
* The supernodal solver uses BLAS routines dtrsv, dgemv, dtrsm, and dgemm,
* or their complex counterparts ztrsv, zgemv, ztrsm, and zgemm.
*
* If both L and B are real, then X is returned real. If either is complex
* or zomplex, X is returned as either complex or zomplex, depending on the
* Common->prefer_zomplex parameter.
*
* Supports any numeric xtype (pattern-only matrices not supported).
*
* This routine does not check to see if the diagonal of L or D is zero,
* because sometimes a partial solve can be done with indefinite or singular
* matrix. If you wish to check in your own code, test L->minor. If
* L->minor == L->n, then the matrix has no zero diagonal entries.
* If k = L->minor < L->n, then L(k,k) is zero for an LL' factorization, or
* D(k,k) is zero for an LDL' factorization.
*
* This routine returns X as NULL only if it runs out of memory. If L is
* indefinite or singular, then X may contain Inf's or NaN's, but it will
* exist on output.
*/
#ifndef NCHOLESKY
#include "cholmod_internal.h"
#include "cholmod_cholesky.h"
#ifndef NSUPERNODAL
#include "cholmod_supernodal.h"
#endif
/* ========================================================================== */
/* === TEMPLATE ============================================================= */
/* ========================================================================== */
#define REAL
#include "t_cholmod_solve.c"
#define COMPLEX
#include "t_cholmod_solve.c"
#define ZOMPLEX
#include "t_cholmod_solve.c"
/* ========================================================================== */
/* === Permutation macro ==================================================== */
/* ========================================================================== */
/* If Perm is NULL, it is interpretted as the identity permutation */
#define P(k) ((Perm == NULL) ? (k) : Perm [k])
/* ========================================================================== */
/* === perm ================================================================= */
/* ========================================================================== */
/* Y = B (P (1:nrow), k1 : min (k1+ncols,ncol)-1) where B is nrow-by-ncol.
*
* Creates a permuted copy of a contiguous set of columns of B.
* Y is already allocated on input. Y must be of sufficient size. Let nk be
* the number of columns accessed in B. Y->xtype determines the complexity of
* the result.
*
* If B is real and Y is complex (or zomplex), only the real part of B is
* copied into Y. The imaginary part of Y is set to zero.
*
* If B is complex (or zomplex) and Y is real, both the real and imaginary and
* parts of B are returned in Y. Y is returned as nrow-by-2*nk. The even
* columns of Y contain the real part of B and the odd columns contain the
* imaginary part of B. Y->nzmax must be >= 2*nrow*nk. Otherise, Y is
* returned as nrow-by-nk with leading dimension nrow. Y->nzmax must be >=
* nrow*nk.
*
* The case where the input (B) is real and the output (Y) is zomplex is
* not used.
*/
static void perm
(
/* ---- input ---- */
cholmod_dense *B, /* input matrix B */
Int *Perm, /* optional input permutation (can be NULL) */
Int k1, /* first column of B to copy */
Int ncols, /* last column to copy is min(k1+ncols,B->ncol)-1 */
/* ---- in/out --- */
cholmod_dense *Y /* output matrix Y, already allocated */
)
{
double *Yx, *Yz, *Bx, *Bz ;
Int k2, nk, p, k, j, nrow, ncol, d, dual, dj, j2 ;
/* ---------------------------------------------------------------------- */
/* get inputs */
/* ---------------------------------------------------------------------- */
ncol = B->ncol ;
nrow = B->nrow ;
k2 = MIN (k1+ncols, ncol) ;
nk = MAX (k2 - k1, 0) ;
dual = (Y->xtype == CHOLMOD_REAL && B->xtype != CHOLMOD_REAL) ? 2 : 1 ;
d = B->d ;
Bx = B->x ;
Bz = B->z ;
Yx = Y->x ;
Yz = Y->z ;
Y->nrow = nrow ;
Y->ncol = dual*nk ;
Y->d = nrow ;
ASSERT (((Int) Y->nzmax) >= nrow*nk*dual) ;
/* ---------------------------------------------------------------------- */
/* Y = B (P (1:nrow), k1:k2-1) */
/* ---------------------------------------------------------------------- */
switch (Y->xtype)
{
case CHOLMOD_REAL:
switch (B->xtype)
{
case CHOLMOD_REAL:
/* Y real, B real */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = nrow * (j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Yx [k + j2] = Bx [p] ; /* real */
}
}
break ;
case CHOLMOD_COMPLEX:
/* Y real, B complex. Y is nrow-by-2*nk */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = nrow * 2 * (j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Yx [k + j2 ] = Bx [2*p ] ; /* real */
Yx [k + j2 + nrow] = Bx [2*p+1] ; /* imag */
}
}
break ;
case CHOLMOD_ZOMPLEX:
/* Y real, B zomplex. Y is nrow-by-2*nk */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = nrow * 2 * (j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Yx [k + j2 ] = Bx [p] ; /* real */
Yx [k + j2 + nrow] = Bz [p] ; /* imag */
}
}
break ;
}
break ;
case CHOLMOD_COMPLEX:
switch (B->xtype)
{
case CHOLMOD_REAL:
/* Y complex, B real */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = nrow * 2 * (j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Yx [2*k + j2] = Bx [p] ; /* real */
Yx [2*k+1 + j2] = 0 ; /* imag */
}
}
break ;
case CHOLMOD_COMPLEX:
/* Y complex, B complex */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = nrow * 2 * (j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Yx [2*k + j2] = Bx [2*p ] ; /* real */
Yx [2*k+1 + j2] = Bx [2*p+1] ; /* imag */
}
}
break ;
case CHOLMOD_ZOMPLEX:
/* Y complex, B zomplex */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = nrow * 2 * (j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Yx [2*k + j2] = Bx [p] ; /* real */
Yx [2*k+1 + j2] = Bz [p] ; /* imag */
}
}
break ;
}
break ;
case CHOLMOD_ZOMPLEX:
switch (B->xtype)
{
#if 0
case CHOLMOD_REAL:
/* this case is not used */
break ;
#endif
case CHOLMOD_COMPLEX:
/* Y zomplex, B complex */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = nrow * (j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Yx [k + j2] = Bx [2*p ] ; /* real */
Yz [k + j2] = Bx [2*p+1] ; /* imag */
}
}
break ;
case CHOLMOD_ZOMPLEX:
/* Y zomplex, B zomplex */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = nrow * (j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Yx [k + j2] = Bx [p] ; /* real */
Yz [k + j2] = Bz [p] ; /* imag */
}
}
break ;
}
break ;
}
}
/* ========================================================================== */
/* === iperm ================================================================ */
/* ========================================================================== */
/* X (P (1:nrow), k1 : min (k1+ncols,ncol)-1) = Y where X is nrow-by-ncol.
*
* Copies and permutes Y into a contiguous set of columns of X. X is already
* allocated on input. Y must be of sufficient size. Let nk be the number
* of columns accessed in X. X->xtype determines the complexity of the result.
*
* If X is real and Y is complex (or zomplex), only the real part of B is
* copied into X. The imaginary part of Y is ignored.
*
* If X is complex (or zomplex) and Y is real, both the real and imaginary and
* parts of Y are returned in X. Y is nrow-by-2*nk. The even
* columns of Y contain the real part of B and the odd columns contain the
* imaginary part of B. Y->nzmax must be >= 2*nrow*nk. Otherise, Y is
* nrow-by-nk with leading dimension nrow. Y->nzmax must be >= nrow*nk.
*
* The case where the input (Y) is complex and the output (X) is real,
* and the case where the input (Y) is zomplex and the output (X) is real,
* are not used.
*/
static void iperm
(
/* ---- input ---- */
cholmod_dense *Y, /* input matrix Y */
Int *Perm, /* optional input permutation (can be NULL) */
Int k1, /* first column of B to copy */
Int ncols, /* last column to copy is min(k1+ncols,B->ncol)-1 */
/* ---- in/out --- */
cholmod_dense *X /* output matrix X, already allocated */
)
{
double *Yx, *Yz, *Xx, *Xz ;
Int k2, nk, p, k, j, nrow, ncol, d, dj, j2 ;
/* ---------------------------------------------------------------------- */
/* get inputs */
/* ---------------------------------------------------------------------- */
ncol = X->ncol ;
nrow = X->nrow ;
k2 = MIN (k1+ncols, ncol) ;
nk = MAX (k2 - k1, 0) ;
d = X->d ;
Xx = X->x ;
Xz = X->z ;
Yx = Y->x ;
Yz = Y->z ;
ASSERT (((Int) Y->nzmax) >= nrow*nk*
((X->xtype != CHOLMOD_REAL && Y->xtype == CHOLMOD_REAL) ? 2:1)) ;
/* ---------------------------------------------------------------------- */
/* X (P (1:nrow), k1:k2-1) = Y */
/* ---------------------------------------------------------------------- */
switch (Y->xtype)
{
case CHOLMOD_REAL:
switch (X->xtype)
{
case CHOLMOD_REAL:
/* Y real, X real */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = nrow * (j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Xx [p] = Yx [k + j2] ; /* real */
}
}
break ;
case CHOLMOD_COMPLEX:
/* Y real, X complex. Y is nrow-by-2*nk */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = nrow * 2 * (j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Xx [2*p ] = Yx [k + j2 ] ; /* real */
Xx [2*p+1] = Yx [k + j2 + nrow] ; /* imag */
}
}
break ;
case CHOLMOD_ZOMPLEX:
/* Y real, X zomplex. Y is nrow-by-2*nk */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = nrow * 2 * (j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Xx [p] = Yx [k + j2 ] ; /* real */
Xz [p] = Yx [k + j2 + nrow] ; /* imag */
}
}
break ;
}
break ;
case CHOLMOD_COMPLEX:
switch (X->xtype)
{
#if 0
case CHOLMOD_REAL:
/* this case is not used */
break ;
#endif
case CHOLMOD_COMPLEX:
/* Y complex, X complex */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = nrow * 2 * (j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Xx [2*p ] = Yx [2*k + j2] ; /* real */
Xx [2*p+1] = Yx [2*k+1 + j2] ; /* imag */
}
}
break ;
case CHOLMOD_ZOMPLEX:
/* Y complex, X zomplex */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = nrow * 2 * (j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Xx [p] = Yx [2*k + j2] ; /* real */
Xz [p] = Yx [2*k+1 + j2] ; /* imag */
}
}
break ;
}
break ;
case CHOLMOD_ZOMPLEX:
switch (X->xtype)
{
#if 0
case CHOLMOD_REAL:
/* this case is not used */
break ;
#endif
case CHOLMOD_COMPLEX:
/* Y zomplex, X complex */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = nrow * (j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Xx [2*p ] = Yx [k + j2] ; /* real */
Xx [2*p+1] = Yz [k + j2] ; /* imag */
}
}
break ;
case CHOLMOD_ZOMPLEX:
/* Y zomplex, X zomplex */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = nrow * (j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Xx [p] = Yx [k + j2] ; /* real */
Xz [p] = Yz [k + j2] ; /* imag */
}
}
break ;
}
break ;
}
}
/* ========================================================================== */
/* === ptrans =============================================================== */
/* ========================================================================== */
/* Y = B (P (1:nrow), k1 : min (k1+ncols,ncol)-1)' where B is nrow-by-ncol.
*
* Creates a permuted and transposed copy of a contiguous set of columns of B.
* Y is already allocated on input. Y must be of sufficient size. Let nk be
* the number of columns accessed in B. Y->xtype determines the complexity of
* the result.
*
* If B is real and Y is complex (or zomplex), only the real part of B is
* copied into Y. The imaginary part of Y is set to zero.
*
* If B is complex (or zomplex) and Y is real, both the real and imaginary and
* parts of B are returned in Y. Y is returned as 2*nk-by-nrow. The even
* rows of Y contain the real part of B and the odd rows contain the
* imaginary part of B. Y->nzmax must be >= 2*nrow*nk. Otherise, Y is
* returned as nk-by-nrow with leading dimension nk. Y->nzmax must be >=
* nrow*nk.
*
* The array transpose is performed, not the complex conjugate transpose.
*/
static void ptrans
(
/* ---- input ---- */
cholmod_dense *B, /* input matrix B */
Int *Perm, /* optional input permutation (can be NULL) */
Int k1, /* first column of B to copy */
Int ncols, /* last column to copy is min(k1+ncols,B->ncol)-1 */
/* ---- in/out --- */
cholmod_dense *Y /* output matrix Y, already allocated */
)
{
double *Yx, *Yz, *Bx, *Bz ;
Int k2, nk, p, k, j, nrow, ncol, d, dual, dj, j2 ;
/* ---------------------------------------------------------------------- */
/* get inputs */
/* ---------------------------------------------------------------------- */
ncol = B->ncol ;
nrow = B->nrow ;
k2 = MIN (k1+ncols, ncol) ;
nk = MAX (k2 - k1, 0) ;
dual = (Y->xtype == CHOLMOD_REAL && B->xtype != CHOLMOD_REAL) ? 2 : 1 ;
d = B->d ;
Bx = B->x ;
Bz = B->z ;
Yx = Y->x ;
Yz = Y->z ;
Y->nrow = dual*nk ;
Y->ncol = nrow ;
Y->d = dual*nk ;
ASSERT (((Int) Y->nzmax) >= nrow*nk*dual) ;
/* ---------------------------------------------------------------------- */
/* Y = B (P (1:nrow), k1:k2-1)' */
/* ---------------------------------------------------------------------- */
switch (Y->xtype)
{
case CHOLMOD_REAL:
switch (B->xtype)
{
case CHOLMOD_REAL:
/* Y real, B real */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = j-k1 ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Yx [j2 + k*nk] = Bx [p] ; /* real */
}
}
break ;
case CHOLMOD_COMPLEX:
/* Y real, B complex. Y is 2*nk-by-nrow */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = 2*(j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Yx [j2 + k*2*nk] = Bx [2*p ] ; /* real */
Yx [j2+1 + k*2*nk] = Bx [2*p+1] ; /* imag */
}
}
break ;
case CHOLMOD_ZOMPLEX:
/* Y real, B zomplex. Y is 2*nk-by-nrow */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = 2*(j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Yx [j2 + k*2*nk] = Bx [p] ; /* real */
Yx [j2+1 + k*2*nk] = Bz [p] ; /* imag */
}
}
break ;
}
break ;
case CHOLMOD_COMPLEX:
switch (B->xtype)
{
case CHOLMOD_REAL:
/* Y complex, B real */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = 2*(j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Yx [j2 + k*2*nk] = Bx [p] ; /* real */
Yx [j2+1 + k*2*nk] = 0 ; /* imag */
}
}
break ;
case CHOLMOD_COMPLEX:
/* Y complex, B complex */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = 2*(j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Yx [j2 + k*2*nk] = Bx [2*p ] ; /* real */
Yx [j2+1 + k*2*nk] = Bx [2*p+1] ; /* imag */
}
}
break ;
case CHOLMOD_ZOMPLEX:
/* Y complex, B zomplex */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = 2*(j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Yx [j2 + k*2*nk] = Bx [p] ; /* real */
Yx [j2+1 + k*2*nk] = Bz [p] ; /* imag */
}
}
break ;
}
break ;
case CHOLMOD_ZOMPLEX:
switch (B->xtype)
{
case CHOLMOD_REAL:
/* Y zomplex, B real */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = j-k1 ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Yx [j2 + k*nk] = Bx [p] ; /* real */
Yz [j2 + k*nk] = 0 ; /* imag */
}
}
break ;
case CHOLMOD_COMPLEX:
/* Y zomplex, B complex */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = j-k1 ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Yx [j2 + k*nk] = Bx [2*p ] ; /* real */
Yz [j2 + k*nk] = Bx [2*p+1] ; /* imag */
}
}
break ;
case CHOLMOD_ZOMPLEX:
/* Y zomplex, B zomplex */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = j-k1 ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Yx [j2 + k*nk] = Bx [p] ; /* real */
Yz [j2 + k*nk] = Bz [p] ; /* imag */
}
}
break ;
}
break ;
}
}
/* ========================================================================== */
/* === iptrans ============================================================== */
/* ========================================================================== */
/* X (P (1:nrow), k1 : min (k1+ncols,ncol)-1) = Y' where X is nrow-by-ncol.
*
* Copies into a permuted and transposed contiguous set of columns of X.
* X is already allocated on input. Y must be of sufficient size. Let nk be
* the number of columns accessed in X. X->xtype determines the complexity of
* the result.
*
* If X is real and Y is complex (or zomplex), only the real part of Y is
* copied into X. The imaginary part of Y is ignored.
*
* If X is complex (or zomplex) and Y is real, both the real and imaginary and
* parts of X are returned in Y. Y is 2*nk-by-nrow. The even
* rows of Y contain the real part of X and the odd rows contain the
* imaginary part of X. Y->nzmax must be >= 2*nrow*nk. Otherise, Y is
* nk-by-nrow with leading dimension nk. Y->nzmax must be >= nrow*nk.
*
* The case where Y is complex or zomplex, and X is real, is not used.
*
* The array transpose is performed, not the complex conjugate transpose.
*/
static void iptrans
(
/* ---- input ---- */
cholmod_dense *Y, /* input matrix Y */
Int *Perm, /* optional input permutation (can be NULL) */
Int k1, /* first column of X to copy into */
Int ncols, /* last column to copy is min(k1+ncols,X->ncol)-1 */
/* ---- in/out --- */
cholmod_dense *X /* output matrix X, already allocated */
)
{
double *Yx, *Yz, *Xx, *Xz ;
Int k2, nk, p, k, j, nrow, ncol, d, dj, j2 ;
/* ---------------------------------------------------------------------- */
/* get inputs */
/* ---------------------------------------------------------------------- */
ncol = X->ncol ;
nrow = X->nrow ;
k2 = MIN (k1+ncols, ncol) ;
nk = MAX (k2 - k1, 0) ;
d = X->d ;
Xx = X->x ;
Xz = X->z ;
Yx = Y->x ;
Yz = Y->z ;
ASSERT (((Int) Y->nzmax) >= nrow*nk*
((X->xtype != CHOLMOD_REAL && Y->xtype == CHOLMOD_REAL) ? 2:1)) ;
/* ---------------------------------------------------------------------- */
/* X (P (1:nrow), k1:k2-1) = Y' */
/* ---------------------------------------------------------------------- */
switch (Y->xtype)
{
case CHOLMOD_REAL:
switch (X->xtype)
{
case CHOLMOD_REAL:
/* Y real, X real */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = j-k1 ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Xx [p] = Yx [j2 + k*nk] ; /* real */
}
}
break ;
case CHOLMOD_COMPLEX:
/* Y real, X complex. Y is 2*nk-by-nrow */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = 2*(j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Xx [2*p ] = Yx [j2 + k*2*nk] ; /* real */
Xx [2*p+1] = Yx [j2+1 + k*2*nk] ; /* imag */
}
}
break ;
case CHOLMOD_ZOMPLEX:
/* Y real, X zomplex. Y is 2*nk-by-nrow */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = 2*(j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Xx [p] = Yx [j2 + k*2*nk] ; /* real */
Xz [p] = Yx [j2+1 + k*2*nk] ; /* imag */
}
}
break ;
}
break ;
case CHOLMOD_COMPLEX:
switch (X->xtype)
{
#if 0
case CHOLMOD_REAL:
/* this case is not used */
break ;
#endif
case CHOLMOD_COMPLEX:
/* Y complex, X complex */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = 2*(j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Xx [2*p ] = Yx [j2 + k*2*nk] ; /* real */
Xx [2*p+1] = Yx [j2+1 + k*2*nk] ; /* imag */
}
}
break ;
case CHOLMOD_ZOMPLEX:
/* Y complex, X zomplex */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = 2*(j-k1) ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Xx [p] = Yx [j2 + k*2*nk] ; /* real */
Xz [p] = Yx [j2+1 + k*2*nk] ; /* imag */
}
}
break ;
}
break ;
case CHOLMOD_ZOMPLEX:
switch (X->xtype)
{
#if 0
case CHOLMOD_REAL:
/* this case is not used */
break ;
#endif
case CHOLMOD_COMPLEX:
/* Y zomplex, X complex */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = j-k1 ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Xx [2*p ] = Yx [j2 + k*nk] ; /* real */
Xx [2*p+1] = Yz [j2 + k*nk] ; /* imag */
}
}
break ;
case CHOLMOD_ZOMPLEX:
/* Y zomplex, X zomplex */
for (j = k1 ; j < k2 ; j++)
{
dj = d*j ;
j2 = j-k1 ;
for (k = 0 ; k < nrow ; k++)
{
p = P(k) + dj ;
Xx [p] = Yx [j2 + k*nk] ; /* real */
Xz [p] = Yz [j2 + k*nk] ; /* imag */
}
}
break ;
}
break ;
}
}
/* ========================================================================== */
/* === cholmod_solve ======================================================== */
/* ========================================================================== */
/* Solve a linear system.
*
* The factorization can be simplicial LDL', simplicial LL', or supernodal LL'.
* The Dx=b solve returns silently for the LL' factorizations (it is implicitly
* identity).
*/
cholmod_dense *CHOLMOD(solve)
(
/* ---- input ---- */
int sys, /* system to solve */
cholmod_factor *L, /* factorization to use */
cholmod_dense *B, /* right-hand-side */
/* --------------- */
cholmod_common *Common
)
{
cholmod_dense *Y = NULL, *X = NULL ;
Int *Perm ;
Int n, nrhs, ncols, ctype, xtype, k1, nr, ytype ;
/* ---------------------------------------------------------------------- */
/* check inputs */
/* ---------------------------------------------------------------------- */
RETURN_IF_NULL_COMMON (NULL) ;
RETURN_IF_NULL (L, NULL) ;
RETURN_IF_NULL (B, NULL) ;
RETURN_IF_XTYPE_INVALID (L, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, NULL) ;
RETURN_IF_XTYPE_INVALID (B, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, NULL) ;
if (sys < CHOLMOD_A || sys > CHOLMOD_Pt)
{
ERROR (CHOLMOD_INVALID, "invalid system") ;
return (NULL) ;
}
if (B->d < L->n || B->nrow != L->n)
{
ERROR (CHOLMOD_INVALID, "dimensions of L and B do not match") ;
return (NULL) ;
}
DEBUG (CHOLMOD(dump_factor) (L, "L", Common)) ;
DEBUG (CHOLMOD(dump_dense) (B, "B", Common)) ;
Common->status = CHOLMOD_OK ;
/* ---------------------------------------------------------------------- */
/* get inputs */
/* ---------------------------------------------------------------------- */
if ((sys == CHOLMOD_P || sys == CHOLMOD_Pt || sys == CHOLMOD_A)
&& L->ordering != CHOLMOD_NATURAL)
{
Perm = L->Perm ;
}
else
{
/* do not use L->Perm; use the identity permutation instead */
Perm = NULL ;
}
nrhs = B->ncol ;
n = L->n ;
/* ---------------------------------------------------------------------- */
/* allocate the result X */
/* ---------------------------------------------------------------------- */
ctype = (Common->prefer_zomplex) ? CHOLMOD_ZOMPLEX : CHOLMOD_COMPLEX ;
if (sys == CHOLMOD_P || sys == CHOLMOD_Pt)
{
/* x=Pb and x=P'b return X real if B is real; X is the preferred
* complex/zcomplex type if B is complex or zomplex */
xtype = (B->xtype == CHOLMOD_REAL) ? CHOLMOD_REAL : ctype ;
}
else if (L->xtype == CHOLMOD_REAL && B->xtype == CHOLMOD_REAL)
{
/* X is real if both L and B are real */
xtype = CHOLMOD_REAL ;
}
else
{
/* X is complex, use the preferred complex/zomplex type */
xtype = ctype ;
}
X = CHOLMOD(allocate_dense) (n, nrhs, n, xtype, Common) ;
if (Common->status < CHOLMOD_OK)
{
return (NULL) ;
}
/* ---------------------------------------------------------------------- */
/* solve using L, D, L', P, or some combination */
/* ---------------------------------------------------------------------- */
if (sys == CHOLMOD_P)
{
/* ------------------------------------------------------------------ */
/* x = P*b */
/* ------------------------------------------------------------------ */
perm (B, Perm, 0, nrhs, X) ;
}
else if (sys == CHOLMOD_Pt)
{
/* ------------------------------------------------------------------ */
/* x = P'*b */
/* ------------------------------------------------------------------ */
iperm (B, Perm, 0, nrhs, X) ;
}
else if (L->is_super)
{
/* ------------------------------------------------------------------ */
/* solve using a supernodal LL' factorization */
/* ------------------------------------------------------------------ */
#ifndef NSUPERNODAL
Int ok ;
/* allocate workspace */
cholmod_dense *E ;
Int dual ;
dual = (L->xtype == CHOLMOD_REAL && B->xtype != CHOLMOD_REAL) ? 2 : 1 ;
Y = CHOLMOD(allocate_dense) (n, dual*nrhs, n, L->xtype, Common) ;
E = CHOLMOD(allocate_dense) (dual*nrhs, L->maxesize, dual*nrhs,
L->xtype, Common) ;
if (Common->status < CHOLMOD_OK)
{
/* out of memory */
CHOLMOD(free_dense) (&X, Common) ;
CHOLMOD(free_dense) (&Y, Common) ;
CHOLMOD(free_dense) (&E, Common) ;
return (NULL) ;
}
perm (B, Perm, 0, nrhs, Y) ; /* Y = P*B */
if (sys == CHOLMOD_A || sys == CHOLMOD_LDLt)
{
ok = CHOLMOD(super_lsolve) (L, Y, E, Common) ; /* Y = L\Y */
ok = ok && CHOLMOD(super_ltsolve) (L, Y, E, Common) ; /* Y = L'\Y*/
}
else if (sys == CHOLMOD_L || sys == CHOLMOD_LD)
{
ok = CHOLMOD(super_lsolve) (L, Y, E, Common) ; /* Y = L\Y */
}
else if (sys == CHOLMOD_Lt || sys == CHOLMOD_DLt)
{
ok = CHOLMOD(super_ltsolve) (L, Y, E, Common) ; /* Y = L'\Y*/
}
CHOLMOD(free_dense) (&E, Common) ;
iperm (Y, Perm, 0, nrhs, X) ; /* X = P'*Y */
if (CHECK_BLAS_INT && !ok)
{
/* integer overflow in the BLAS */
CHOLMOD(free_dense) (&X, Common) ;
}
#else
/* CHOLMOD Supernodal module not installed */
ERROR (CHOLMOD_NOT_INSTALLED,"Supernodal module not installed") ;
#endif
}
else
{
/* ------------------------------------------------------------------ */
/* solve using a simplicial LL' or LDL' factorization */
/* ------------------------------------------------------------------ */
if (L->xtype == CHOLMOD_REAL && B->xtype == CHOLMOD_REAL)
{
/* L, B, and Y are all real */
/* solve with up to 4 columns of B at a time */
ncols = 4 ;
nr = MAX (4, nrhs) ;
ytype = CHOLMOD_REAL ;
}
else if (L->xtype == CHOLMOD_REAL)
{
/* solve with one column of B (real/imag), at a time */
ncols = 1 ;
nr = 2 ;
ytype = CHOLMOD_REAL ;
}
else
{
/* L is complex or zomplex, B is real/complex/zomplex, Y has the
* same complexity as L. Solve with one column of B at a time. */
ncols = 1 ;
nr = 1 ;
ytype = L->xtype ;
}
Y = CHOLMOD(allocate_dense) (nr, n, nr, ytype, Common) ;
if (Common->status < CHOLMOD_OK)
{
/* out of memory */
CHOLMOD(free_dense) (&X, Common) ;
CHOLMOD(free_dense) (&Y, Common) ;
return (NULL) ;
}
for (k1 = 0 ; k1 < nrhs ; k1 += ncols)
{
/* -------------------------------------------------------------- */
/* Y = B (P, k1:k1+ncols-1)' = (P * B (:,...))' */
/* -------------------------------------------------------------- */
ptrans (B, Perm, k1, ncols, Y) ;
/* -------------------------------------------------------------- */
/* solve Y = (L' \ (L \ Y'))', or other system, with template */
/* -------------------------------------------------------------- */
switch (L->xtype)
{
case CHOLMOD_REAL:
r_simplicial_solver (sys, L, Y) ;
break ;
case CHOLMOD_COMPLEX:
c_simplicial_solver (sys, L, Y) ;
break ;
case CHOLMOD_ZOMPLEX:
z_simplicial_solver (sys, L, Y) ;
break ;
}
/* -------------------------------------------------------------- */
/* X (P, k1:k2+ncols-1) = Y' */
/* -------------------------------------------------------------- */
iptrans (Y, Perm, k1, ncols, X) ;
}
}
CHOLMOD(free_dense) (&Y, Common) ;
DEBUG (CHOLMOD(dump_dense) (X, "X result", Common)) ;
return (X) ;
}
#endif
syntax highlighted by Code2HTML, v. 0.9.1