The R Project SVN
R-packages
(root)/branches/Matrix-nolmer/src/CHOLMOD/Cholesky/t_cholmod_ltsolve.c – Rev 3548
Rev 1820 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
/* ========================================================================== *//* === Cholesky/t_cholmod_ltsolve =========================================== *//* ========================================================================== *//* -----------------------------------------------------------------------------* CHOLMOD/Cholesky Module. Version 0.6. Copyright (C) 2005, 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* -------------------------------------------------------------------------- *//* Template routine to solve L'x=b with unit or non-unit diagonal, or* solve DL'x=b.** The numeric xtype of L and Y must match. Y contains b on input and x on* output, stored in row-form. Y is nrow-by-n, where nrow must equal 1 for the* complex or zomplex cases, and nrow <= 4 for the real case.** This file is not compiled separately. It is included in t_cholmod_solve.c* instead. It contains no user-callable routines.** workspace: none** Supports real, complex, and zomplex factors.*//* undefine all prior definitions */#undef FORM_NAME#undef LSOLVE#undef DIAG/* -------------------------------------------------------------------------- *//* define the method *//* -------------------------------------------------------------------------- */#ifdef LL/* LL': solve Lx=b with non-unit diagonal */#define FORM_NAME(prefix,rank) prefix ## ll_ltsolve_ ## rank#define DIAG#elif defined (LD)/* LDL': solve LDx=b */#define FORM_NAME(prefix,rank) prefix ## ldl_dltsolve_ ## rank#define DIAG#else/* LDL': solve Lx=b with unit diagonal */#define FORM_NAME(prefix,rank) prefix ## ldl_ltsolve_ ## rank#endif/* LSOLVE(k) defines the name of a routine for an n-by-k right-hand-side. */#define LSOLVE(prefix,rank) FORM_NAME(prefix,rank)#ifdef REAL/* ========================================================================== *//* === LSOLVE (1) =========================================================== *//* ========================================================================== *//* Solve L'x=b, where b has 1 column */static void LSOLVE (PREFIX,1)(cholmod_factor *L,cholmod_dense *Y /* n-by-1 in row form */){double *Lx = L->x, *X = Y->x ;Int *Li = L->i ;Int *Lp = L->p ;Int *Lnz = L->nz ;Int j, n = L->n ;for (j = n-1 ; j >= 0 ; ){/* get the start, end, and length of column j */Int p = Lp [j] ;Int lnz = Lnz [j] ;Int pend = p + lnz ;/* find a chain of supernodes (up to j, j-1, and j-2) */if (j == 0 || lnz != Lnz [j-1] - 1 || Li [Lp [j-1]+1] != j){/* -------------------------------------------------------------- *//* solve with a single column of L *//* -------------------------------------------------------------- */double y = X [j] ;#ifdef DIAGdouble d = Lx [p] ;#endif#ifdef LDy /= d ;#endiffor (p++ ; p < pend ; p++){y -= Lx [p] * X [Li [p]] ;}#ifdef LLX [j] = y / d ;#elseX [j] = y ;#endifj-- ;}else if (j == 1 || lnz != Lnz [j-2]-2 || Li [Lp [j-2]+2] != j){/* -------------------------------------------------------------- *//* solve with a supernode of two columns of L *//* -------------------------------------------------------------- */double y [2], t ;Int q = Lp [j-1] ;#ifdef DIAGdouble d [2] ;d [0] = Lx [p] ;d [1] = Lx [q] ;#endift = Lx [q+1] ;#ifdef LDy [0] = X [j ] / d [0] ;y [1] = X [j-1] / d [1] ;#elsey [0] = X [j ] ;y [1] = X [j-1] ;#endiffor (p++, q += 2 ; p < pend ; p++, q++){Int i = Li [p] ;y [0] -= Lx [p] * X [i] ;y [1] -= Lx [q] * X [i] ;}#ifdef LLy [0] /= d [0] ;y [1] = (y [1] - t * y [0]) / d [1] ;#elsey [1] -= t * y [0] ;#endifX [j ] = y [0] ;X [j-1] = y [1] ;j -= 2 ;}else{/* -------------------------------------------------------------- *//* solve with a supernode of three columns of L *//* -------------------------------------------------------------- */double y [3], t [3] ;Int q = Lp [j-1] ;Int r = Lp [j-2] ;#ifdef DIAGdouble d [3] ;d [0] = Lx [p] ;d [1] = Lx [q] ;d [2] = Lx [r] ;#endift [0] = Lx [q+1] ;t [1] = Lx [r+1] ;t [2] = Lx [r+2] ;#ifdef LDy [0] = X [j] / d [0] ;y [1] = X [j-1] / d [1] ;y [2] = X [j-2] / d [2] ;#elsey [0] = X [j] ;y [1] = X [j-1] ;y [2] = X [j-2] ;#endiffor (p++, q += 2, r += 3 ; p < pend ; p++, q++, r++){Int i = Li [p] ;y [0] -= Lx [p] * X [i] ;y [1] -= Lx [q] * X [i] ;y [2] -= Lx [r] * X [i] ;}q = Lp [j-1] ;r = Lp [j-2] ;#ifdef LLy [0] /= d [0] ;y [1] = (y [1] - t [0] * y [0]) / d [1] ;y [2] = (y [2] - t [2] * y [0] - t [1] * y [1]) / d [2] ;#elsey [1] -= t [0] * y [0] ;y [2] -= t [2] * y [0] + t [1] * y [1] ;#endifX [j-2] = y [2] ;X [j-1] = y [1] ;X [j ] = y [0] ;j -= 3 ;}}}/* ========================================================================== *//* === LSOLVE (2) =========================================================== *//* ========================================================================== *//* Solve L'x=b, where b has 2 columns */static void LSOLVE (PREFIX,2)(cholmod_factor *L,cholmod_dense *Y /* n-by-2 in row form */){double *Lx = L->x, *X = Y->x ;Int *Li = L->i ;Int *Lp = L->p ;Int *Lnz = L->nz ;Int j, n = L->n ;for (j = n-1 ; j >= 0 ; ){/* get the start, end, and length of column j */Int p = Lp [j] ;Int lnz = Lnz [j] ;Int pend = p + lnz ;/* find a chain of supernodes (up to j, j-1, and j-2) */if (j == 0 || lnz != Lnz [j-1] - 1 || Li [Lp [j-1]+1] != j){/* -------------------------------------------------------------- *//* solve with a single column of L *//* -------------------------------------------------------------- */double y [2] ;#ifdef DIAGdouble d = Lx [p] ;#endif#ifdef LDy [0] = X [2*j ] / d ;y [1] = X [2*j+1] / d ;#elsey [0] = X [2*j ] ;y [1] = X [2*j+1] ;#endiffor (p++ ; p < pend ; p++){Int i = 2 * Li [p] ;y [0] -= Lx [p] * X [i ] ;y [1] -= Lx [p] * X [i+1] ;}#ifdef LLX [2*j ] = y [0] / d ;X [2*j+1] = y [1] / d ;#elseX [2*j ] = y [0] ;X [2*j+1] = y [1] ;#endifj-- ;}else if (j == 1 || lnz != Lnz [j-2]-2 || Li [Lp [j-2]+2] != j){/* -------------------------------------------------------------- *//* solve with a supernode of two columns of L *//* -------------------------------------------------------------- */double y [2][2], t ;Int q = Lp [j-1] ;#ifdef DIAGdouble d [2] ;d [0] = Lx [p] ;d [1] = Lx [q] ;#endift = Lx [q+1] ;#ifdef LDy [0][0] = X [2*j ] / d [0] ;y [0][1] = X [2*j+1] / d [0] ;y [1][0] = X [2*j-2] / d [1] ;y [1][1] = X [2*j-1] / d [1] ;#elsey [0][0] = X [2*j ] ;y [0][1] = X [2*j+1] ;y [1][0] = X [2*j-2] ;y [1][1] = X [2*j-1] ;#endiffor (p++, q += 2 ; p < pend ; p++, q++){Int i = 2 * Li [p] ;y [0][0] -= Lx [p] * X [i] ;y [0][1] -= Lx [p] * X [i+1] ;y [1][0] -= Lx [q] * X [i] ;y [1][1] -= Lx [q] * X [i+1] ;}#ifdef LLy [0][0] /= d [0] ;y [0][1] /= d [0] ;y [1][0] = (y [1][0] - t * y [0][0]) / d [1] ;y [1][1] = (y [1][1] - t * y [0][1]) / d [1] ;#elsey [1][0] -= t * y [0][0] ;y [1][1] -= t * y [0][1] ;#endifX [2*j ] = y [0][0] ;X [2*j+1] = y [0][1] ;X [2*j-2] = y [1][0] ;X [2*j-1] = y [1][1] ;j -= 2 ;}else{/* -------------------------------------------------------------- *//* solve with a supernode of three columns of L *//* -------------------------------------------------------------- */double y [3][2], t [3] ;Int q = Lp [j-1] ;Int r = Lp [j-2] ;#ifdef DIAGdouble d [3] ;d [0] = Lx [p] ;d [1] = Lx [q] ;d [2] = Lx [r] ;#endift [0] = Lx [q+1] ;t [1] = Lx [r+1] ;t [2] = Lx [r+2] ;#ifdef LDy [0][0] = X [2*j ] / d [0] ;y [0][1] = X [2*j+1] / d [0] ;y [1][0] = X [2*j-2] / d [1] ;y [1][1] = X [2*j-1] / d [1] ;y [2][0] = X [2*j-4] / d [2] ;y [2][1] = X [2*j-3] / d [2] ;#elsey [0][0] = X [2*j ] ;y [0][1] = X [2*j+1] ;y [1][0] = X [2*j-2] ;y [1][1] = X [2*j-1] ;y [2][0] = X [2*j-4] ;y [2][1] = X [2*j-3] ;#endiffor (p++, q += 2, r += 3 ; p < pend ; p++, q++, r++){Int i = 2 * Li [p] ;y [0][0] -= Lx [p] * X [i] ; y [0][1] -= Lx [p] * X [i+1] ;y [1][0] -= Lx [q] * X [i] ; y [1][1] -= Lx [q] * X [i+1] ;y [2][0] -= Lx [r] * X [i] ; y [2][1] -= Lx [r] * X [i+1] ;}#ifdef LLy [0][0] /= d [0] ;y [0][1] /= d [0] ;y [1][0] = (y [1][0] - t [0] * y [0][0]) / d [1] ;y [1][1] = (y [1][1] - t [0] * y [0][1]) / d [1] ;y [2][0] = (y [2][0] - t [2] * y [0][0] - t [1] * y [1][0]) / d [2];y [2][1] = (y [2][1] - t [2] * y [0][1] - t [1] * y [1][1]) / d [2];#elsey [1][0] -= t [0] * y [0][0] ;y [1][1] -= t [0] * y [0][1] ;y [2][0] -= t [2] * y [0][0] + t [1] * y [1][0] ;y [2][1] -= t [2] * y [0][1] + t [1] * y [1][1] ;#endifX [2*j ] = y [0][0] ;X [2*j+1] = y [0][1] ;X [2*j-2] = y [1][0] ;X [2*j-1] = y [1][1] ;X [2*j-4] = y [2][0] ;X [2*j-3] = y [2][1] ;j -= 3 ;}}}/* ========================================================================== *//* === LSOLVE (3) =========================================================== *//* ========================================================================== *//* Solve L'x=b, where b has 3 columns */static void LSOLVE (PREFIX,3)(cholmod_factor *L,cholmod_dense *Y /* n-by-3 in row form */){double *Lx = L->x, *X = Y->x ;Int *Li = L->i ;Int *Lp = L->p ;Int *Lnz = L->nz ;Int j, n = L->n ;for (j = n-1 ; j >= 0 ; ){/* get the start, end, and length of column j */Int p = Lp [j] ;Int lnz = Lnz [j] ;Int pend = p + lnz ;/* find a chain of supernodes (up to j, j-1, and j-2) */if (j == 0 || lnz != Lnz [j-1] - 1 || Li [Lp [j-1]+1] != j){/* -------------------------------------------------------------- *//* solve with a single column of L *//* -------------------------------------------------------------- */double y [3] ;#ifdef DIAGdouble d = Lx [p] ;#endif#ifdef LDy [0] = X [3*j ] / d ;y [1] = X [3*j+1] / d ;y [2] = X [3*j+2] / d ;#elsey [0] = X [3*j ] ;y [1] = X [3*j+1] ;y [2] = X [3*j+2] ;#endiffor (p++ ; p < pend ; p++){Int i = 3 * Li [p] ;y [0] -= Lx [p] * X [i ] ;y [1] -= Lx [p] * X [i+1] ;y [2] -= Lx [p] * X [i+2] ;}#ifdef LLX [3*j ] = y [0] / d ;X [3*j+1] = y [1] / d ;X [3*j+2] = y [2] / d ;#elseX [3*j ] = y [0] ;X [3*j+1] = y [1] ;X [3*j+2] = y [2] ;#endifj-- ;}else if (j == 1 || lnz != Lnz [j-2]-2 || Li [Lp [j-2]+2] != j){/* -------------------------------------------------------------- *//* solve with a supernode of two columns of L *//* -------------------------------------------------------------- */double y [2][3], t ;Int q = Lp [j-1] ;#ifdef DIAGdouble d [2] ;d [0] = Lx [p] ;d [1] = Lx [q] ;#endift = Lx [q+1] ;#ifdef LDy [0][0] = X [3*j ] / d [0] ;y [0][1] = X [3*j+1] / d [0] ;y [0][2] = X [3*j+2] / d [0] ;y [1][0] = X [3*j-3] / d [1] ;y [1][1] = X [3*j-2] / d [1] ;y [1][2] = X [3*j-1] / d [1] ;#elsey [0][0] = X [3*j ] ;y [0][1] = X [3*j+1] ;y [0][2] = X [3*j+2] ;y [1][0] = X [3*j-3] ;y [1][1] = X [3*j-2] ;y [1][2] = X [3*j-1] ;#endiffor (p++, q += 2 ; p < pend ; p++, q++){Int i = 3 * Li [p] ;y [0][0] -= Lx [p] * X [i] ;y [0][1] -= Lx [p] * X [i+1] ;y [0][2] -= Lx [p] * X [i+2] ;y [1][0] -= Lx [q] * X [i] ;y [1][1] -= Lx [q] * X [i+1] ;y [1][2] -= Lx [q] * X [i+2] ;}q = Lp [j-1] ;#ifdef LLy [0][0] /= d [0] ;y [0][1] /= d [0] ;y [0][2] /= d [0] ;y [1][0] = (y [1][0] - t * y [0][0]) / d [1] ;y [1][1] = (y [1][1] - t * y [0][1]) / d [1] ;y [1][2] = (y [1][2] - t * y [0][2]) / d [1] ;#elsey [1][0] -= t * y [0][0] ;y [1][1] -= t * y [0][1] ;y [1][2] -= t * y [0][2] ;#endifX [3*j ] = y [0][0] ;X [3*j+1] = y [0][1] ;X [3*j+2] = y [0][2] ;X [3*j-3] = y [1][0] ;X [3*j-2] = y [1][1] ;X [3*j-1] = y [1][2] ;j -= 2 ;}else{/* -------------------------------------------------------------- *//* solve with a supernode of three columns of L *//* -------------------------------------------------------------- */double y [3][3], t [3] ;Int q = Lp [j-1] ;Int r = Lp [j-2] ;#ifdef DIAGdouble d [3] ;d [0] = Lx [p] ;d [1] = Lx [q] ;d [2] = Lx [r] ;#endift [0] = Lx [q+1] ;t [1] = Lx [r+1] ;t [2] = Lx [r+2] ;#ifdef LDy [0][0] = X [3*j ] / d [0] ;y [0][1] = X [3*j+1] / d [0] ;y [0][2] = X [3*j+2] / d [0] ;y [1][0] = X [3*j-3] / d [1] ;y [1][1] = X [3*j-2] / d [1] ;y [1][2] = X [3*j-1] / d [1] ;y [2][0] = X [3*j-6] / d [2] ;y [2][1] = X [3*j-5] / d [2] ;y [2][2] = X [3*j-4] / d [2] ;#elsey [0][0] = X [3*j ] ;y [0][1] = X [3*j+1] ;y [0][2] = X [3*j+2] ;y [1][0] = X [3*j-3] ;y [1][1] = X [3*j-2] ;y [1][2] = X [3*j-1] ;y [2][0] = X [3*j-6] ;y [2][1] = X [3*j-5] ;y [2][2] = X [3*j-4] ;#endiffor (p++, q += 2, r += 3 ; p < pend ; p++, q++, r++){Int i = 3 * Li [p] ;y [0][0] -= Lx [p] * X [i] ;y [0][1] -= Lx [p] * X [i+1] ;y [0][2] -= Lx [p] * X [i+2] ;y [1][0] -= Lx [q] * X [i] ;y [1][1] -= Lx [q] * X [i+1] ;y [1][2] -= Lx [q] * X [i+2] ;y [2][0] -= Lx [r] * X [i] ;y [2][1] -= Lx [r] * X [i+1] ;y [2][2] -= Lx [r] * X [i+2] ;}#ifdef LLy [0][0] /= d [0] ;y [0][1] /= d [0] ;y [0][2] /= d [0] ;y [1][0] = (y [1][0] - t [0] * y [0][0]) / d [1] ;y [1][1] = (y [1][1] - t [0] * y [0][1]) / d [1] ;y [1][2] = (y [1][2] - t [0] * y [0][2]) / d [1] ;y [2][0] = (y [2][0] - t [2] * y [0][0] - t [1] * y [1][0]) / d [2];y [2][1] = (y [2][1] - t [2] * y [0][1] - t [1] * y [1][1]) / d [2];y [2][2] = (y [2][2] - t [2] * y [0][2] - t [1] * y [1][2]) / d [2];#elsey [1][0] -= t [0] * y [0][0] ;y [1][1] -= t [0] * y [0][1] ;y [1][2] -= t [0] * y [0][2] ;y [2][0] -= t [2] * y [0][0] + t [1] * y [1][0] ;y [2][1] -= t [2] * y [0][1] + t [1] * y [1][1] ;y [2][2] -= t [2] * y [0][2] + t [1] * y [1][2] ;#endifX [3*j ] = y [0][0] ;X [3*j+1] = y [0][1] ;X [3*j+2] = y [0][2] ;X [3*j-3] = y [1][0] ;X [3*j-2] = y [1][1] ;X [3*j-1] = y [1][2] ;X [3*j-6] = y [2][0] ;X [3*j-5] = y [2][1] ;X [3*j-4] = y [2][2] ;j -= 3 ;}}}/* ========================================================================== *//* === LSOLVE (4) =========================================================== *//* ========================================================================== *//* Solve L'x=b, where b has 4 columns */static void LSOLVE (PREFIX,4)(cholmod_factor *L,cholmod_dense *Y /* n-by-4 in row form */){double *Lx = L->x, *X = Y->x ;Int *Li = L->i ;Int *Lp = L->p ;Int *Lnz = L->nz ;Int j, n = L->n ;for (j = n-1 ; j >= 0 ; ){/* get the start, end, and length of column j */Int p = Lp [j] ;Int lnz = Lnz [j] ;Int pend = p + lnz ;/* find a chain of supernodes (up to j, j-1, and j-2) */if (j == 0 || lnz != Lnz [j-1] - 1 || Li [Lp [j-1]+1] != j){/* -------------------------------------------------------------- *//* solve with a single column of L *//* -------------------------------------------------------------- */double y [4] ;#ifdef DIAGdouble d = Lx [p] ;#endif#ifdef LDy [0] = X [4*j ] / d ;y [1] = X [4*j+1] / d ;y [2] = X [4*j+2] / d ;y [3] = X [4*j+3] / d ;#elsey [0] = X [4*j ] ;y [1] = X [4*j+1] ;y [2] = X [4*j+2] ;y [3] = X [4*j+3] ;#endiffor (p++ ; p < pend ; p++){Int i = 4 * Li [p] ;y [0] -= Lx [p] * X [i ] ;y [1] -= Lx [p] * X [i+1] ;y [2] -= Lx [p] * X [i+2] ;y [3] -= Lx [p] * X [i+3] ;}#ifdef LLX [4*j ] = y [0] / d ;X [4*j+1] = y [1] / d ;X [4*j+2] = y [2] / d ;X [4*j+3] = y [3] / d ;#elseX [4*j ] = y [0] ;X [4*j+1] = y [1] ;X [4*j+2] = y [2] ;X [4*j+3] = y [3] ;#endifj-- ;}else if (j == 1 || lnz != Lnz [j-2]-2 || Li [Lp [j-2]+2] != j){/* -------------------------------------------------------------- *//* solve with a supernode of two columns of L *//* -------------------------------------------------------------- */double y [2][4], t ;Int q = Lp [j-1] ;#ifdef DIAGdouble d [2] ;d [0] = Lx [p] ;d [1] = Lx [q] ;#endift = Lx [q+1] ;#ifdef LDy [0][0] = X [4*j ] / d [0] ;y [0][1] = X [4*j+1] / d [0] ;y [0][2] = X [4*j+2] / d [0] ;y [0][3] = X [4*j+3] / d [0] ;y [1][0] = X [4*j-4] / d [1] ;y [1][1] = X [4*j-3] / d [1] ;y [1][2] = X [4*j-2] / d [1] ;y [1][3] = X [4*j-1] / d [1] ;#elsey [0][0] = X [4*j ] ;y [0][1] = X [4*j+1] ;y [0][2] = X [4*j+2] ;y [0][3] = X [4*j+3] ;y [1][0] = X [4*j-4] ;y [1][1] = X [4*j-3] ;y [1][2] = X [4*j-2] ;y [1][3] = X [4*j-1] ;#endiffor (p++, q += 2 ; p < pend ; p++, q++){Int i = 4 * Li [p] ;y [0][0] -= Lx [p] * X [i] ;y [0][1] -= Lx [p] * X [i+1] ;y [0][2] -= Lx [p] * X [i+2] ;y [0][3] -= Lx [p] * X [i+3] ;y [1][0] -= Lx [q] * X [i] ;y [1][1] -= Lx [q] * X [i+1] ;y [1][2] -= Lx [q] * X [i+2] ;y [1][3] -= Lx [q] * X [i+3] ;}#ifdef LLy [0][0] /= d [0] ;y [0][1] /= d [0] ;y [0][2] /= d [0] ;y [0][3] /= d [0] ;y [1][0] = (y [1][0] - t * y [0][0]) / d [1] ;y [1][1] = (y [1][1] - t * y [0][1]) / d [1] ;y [1][2] = (y [1][2] - t * y [0][2]) / d [1] ;y [1][3] = (y [1][3] - t * y [0][3]) / d [1] ;#elsey [1][0] -= t * y [0][0] ;y [1][1] -= t * y [0][1] ;y [1][2] -= t * y [0][2] ;y [1][3] -= t * y [0][3] ;#endifX [4*j ] = y [0][0] ;X [4*j+1] = y [0][1] ;X [4*j+2] = y [0][2] ;X [4*j+3] = y [0][3] ;X [4*j-4] = y [1][0] ;X [4*j-3] = y [1][1] ;X [4*j-2] = y [1][2] ;X [4*j-1] = y [1][3] ;j -= 2 ;}else{/* -------------------------------------------------------------- *//* solve with a supernode of three columns of L *//* -------------------------------------------------------------- */double y [3][4], t [3] ;Int q = Lp [j-1] ;Int r = Lp [j-2] ;#ifdef DIAGdouble d [3] ;d [0] = Lx [p] ;d [1] = Lx [q] ;d [2] = Lx [r] ;#endift [0] = Lx [q+1] ;t [1] = Lx [r+1] ;t [2] = Lx [r+2] ;#ifdef LDy [0][0] = X [4*j ] / d [0] ;y [0][1] = X [4*j+1] / d [0] ;y [0][2] = X [4*j+2] / d [0] ;y [0][3] = X [4*j+3] / d [0] ;y [1][0] = X [4*j-4] / d [1] ;y [1][1] = X [4*j-3] / d [1] ;y [1][2] = X [4*j-2] / d [1] ;y [1][3] = X [4*j-1] / d [1] ;y [2][0] = X [4*j-8] / d [2] ;y [2][1] = X [4*j-7] / d [2] ;y [2][2] = X [4*j-6] / d [2] ;y [2][3] = X [4*j-5] / d [2] ;#elsey [0][0] = X [4*j ] ;y [0][1] = X [4*j+1] ;y [0][2] = X [4*j+2] ;y [0][3] = X [4*j+3] ;y [1][0] = X [4*j-4] ;y [1][1] = X [4*j-3] ;y [1][2] = X [4*j-2] ;y [1][3] = X [4*j-1] ;y [2][0] = X [4*j-8] ;y [2][1] = X [4*j-7] ;y [2][2] = X [4*j-6] ;y [2][3] = X [4*j-5] ;#endiffor (p++, q += 2, r += 3 ; p < pend ; p++, q++, r++){Int i = 4 * Li [p] ;y [0][0] -= Lx [p] * X [i] ;y [0][1] -= Lx [p] * X [i+1] ;y [0][2] -= Lx [p] * X [i+2] ;y [0][3] -= Lx [p] * X [i+3] ;y [1][0] -= Lx [q] * X [i] ;y [1][1] -= Lx [q] * X [i+1] ;y [1][2] -= Lx [q] * X [i+2] ;y [1][3] -= Lx [q] * X [i+3] ;y [2][0] -= Lx [r] * X [i] ;y [2][1] -= Lx [r] * X [i+1] ;y [2][2] -= Lx [r] * X [i+2] ;y [2][3] -= Lx [r] * X [i+3] ;}#ifdef LLy [0][0] /= d [0] ;y [0][1] /= d [0] ;y [0][2] /= d [0] ;y [0][3] /= d [0] ;y [1][0] = (y [1][0] - t [0] * y [0][0]) / d [1] ;y [1][1] = (y [1][1] - t [0] * y [0][1]) / d [1] ;y [1][2] = (y [1][2] - t [0] * y [0][2]) / d [1] ;y [1][3] = (y [1][3] - t [0] * y [0][3]) / d [1] ;y [2][0] = (y [2][0] - t [2] * y [0][0] - t [1] * y [1][0]) / d [2];y [2][1] = (y [2][1] - t [2] * y [0][1] - t [1] * y [1][1]) / d [2];y [2][2] = (y [2][2] - t [2] * y [0][2] - t [1] * y [1][2]) / d [2];y [2][3] = (y [2][3] - t [2] * y [0][3] - t [1] * y [1][3]) / d [2];#elsey [1][0] -= t [0] * y [0][0] ;y [1][1] -= t [0] * y [0][1] ;y [1][2] -= t [0] * y [0][2] ;y [1][3] -= t [0] * y [0][3] ;y [2][0] -= t [2] * y [0][0] + t [1] * y [1][0] ;y [2][1] -= t [2] * y [0][1] + t [1] * y [1][1] ;y [2][2] -= t [2] * y [0][2] + t [1] * y [1][2] ;y [2][3] -= t [2] * y [0][3] + t [1] * y [1][3] ;#endifX [4*j ] = y [0][0] ;X [4*j+1] = y [0][1] ;X [4*j+2] = y [0][2] ;X [4*j+3] = y [0][3] ;X [4*j-4] = y [1][0] ;X [4*j-3] = y [1][1] ;X [4*j-2] = y [1][2] ;X [4*j-1] = y [1][3] ;X [4*j-8] = y [2][0] ;X [4*j-7] = y [2][1] ;X [4*j-6] = y [2][2] ;X [4*j-5] = y [2][3] ;j -= 3 ;}}}#endif/* ========================================================================== *//* === LSOLVE (k) =========================================================== *//* ========================================================================== */static void LSOLVE (PREFIX,k)(cholmod_factor *L,cholmod_dense *Y /* nr-by-n where nr is 1 to 4 */){#ifndef REAL#ifdef DIAGdouble d [1] ;#endifdouble yx [2] ;#ifdef ZOMPLEXdouble yz [1] ;double *Lz = L->z ;double *Xz = Y->z ;#endifdouble *Lx = L->x ;double *Xx = Y->x ;Int *Li = L->i ;Int *Lp = L->p ;Int *Lnz = L->nz ;Int i, j, n = L->n ;#endifASSERT (L->xtype == Y->xtype) ; /* L and Y must have the same xtype */ASSERT (L->n == Y->ncol) ; /* dimensions must match */ASSERT (Y->nrow == Y->d) ; /* leading dimension of Y = # rows of Y */ASSERT (L->xtype != CHOLMOD_PATTERN) ; /* L is not symbolic */ASSERT (!(L->is_super)) ; /* L is simplicial LL' or LDL' */#ifdef REAL/* ---------------------------------------------------------------------- *//* solve a real linear system, with 1 to 4 RHS's and dynamic supernodes *//* ---------------------------------------------------------------------- */ASSERT (Y->nrow <= 4) ;switch (Y->nrow){case 1: LSOLVE (PREFIX,1) (L, Y) ; break ;case 2: LSOLVE (PREFIX,2) (L, Y) ; break ;case 3: LSOLVE (PREFIX,3) (L, Y) ; break ;case 4: LSOLVE (PREFIX,4) (L, Y) ; break ;}#else/* ---------------------------------------------------------------------- *//* solve a complex linear system, with just one right-hand-side *//* ---------------------------------------------------------------------- */ASSERT (Y->nrow == 1) ;for (j = n-1 ; j >= 0 ; j--){/* get the start, end, and length of column j */Int p = Lp [j] ;Int lnz = Lnz [j] ;Int pend = p + lnz ;/* y = X [j] ; */ASSIGN (yx,yz,0, Xx,Xz,j) ;#ifdef DIAG/* d = Lx [p] ; */ASSIGN_REAL (d,0, Lx,p) ;#endif#ifdef LD/* y /= d ; */DIV_REAL (yx,yz,0, yx,yz,0, d,0) ;#endiffor (p++ ; p < pend ; p++){/* y -= conj (Lx [p]) * X [Li [p]] ; */i = Li [p] ;MULTSUBCONJ (yx,yz,0, Lx,Lz,p, Xx,Xz,i) ;}#ifdef LL/* X [j] = y / d ; */DIV_REAL (Xx,Xz,j, yx,yz,0, d,0) ;#else/* X [j] = y ; */ASSIGN (Xx,Xz,j, yx,yz,0) ;#endif}#endif}/* prepare for the next inclusion of this file in cholmod_solve.c */#undef LL#undef LD