The R Project SVN R-packages

/* ========================================================================== */
/* === 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 blas_ok = TRUE ;

    /* 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)
    {
        blas_ok = CHOLMOD(super_lsolve) (L, Y, E, Common) ;    /* Y = L\Y */
        blas_ok = blas_ok &&
        CHOLMOD(super_ltsolve) (L, Y, E, Common) ;     /* Y = L'\Y*/
    }
    else if (sys == CHOLMOD_L || sys == CHOLMOD_LD)
    {
        blas_ok = CHOLMOD(super_lsolve) (L, Y, E, Common) ;    /* Y = L\Y */
    }
    else if (sys == CHOLMOD_Lt || sys == CHOLMOD_DLt)
    {
        blas_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 && !blas_ok)
    {
        /* Integer overflow in the BLAS.  This is probably impossible,
         * since the BLAS were used to create the supernodal factorization.
         * It might be possible for the calls to the BLAS to differ between
         * factorization and forward/backsolves, however.  This statement
         * is untested; it does not appear in the compiled code if
         * CHECK_BLAS_INT is true (when the same integer is used in CHOLMOD
         * and 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