amesos_cholmod_solve.c

Go to the documentation of this file.

/* ========================================================================== */
/* === 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 "amesos_cholmod_internal.h"
#include "amesos_cholmod_cholesky.h"


/* ========================================================================== */
/* === TEMPLATE ============================================================= */
/* ========================================================================== */

#define REAL
#include "amesos_t_cholmod_solve.c"

#define COMPLEX
#include "amesos_t_cholmod_solve.c"

#define ZOMPLEX
#include "amesos_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 amesos_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 amesos_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 amesos_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 amesos_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 */
  /* ------------------------------------------------------------------ */

  amesos_perm (B, Perm, 0, nrhs, X) ;

    }
    else if (sys == CHOLMOD_Pt)
    {

  /* ------------------------------------------------------------------ */
  /* x = P'*b */
  /* ------------------------------------------------------------------ */

  amesos_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) ;
  }

  amesos_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) ;

  amesos_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 (:,...))' */
      /* -------------------------------------------------------------- */

      amesos_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' */
      /* -------------------------------------------------------------- */

      amesos_iptrans (Y, Perm, k1, ncols, X) ;
  }
    }

    CHOLMOD(free_dense) (&Y, Common) ;
    DEBUG (CHOLMOD(dump_dense) (X, "X result", Common)) ;
    return (X) ;
}
#endif