The R Project SVN
R-packages
(root)/branches/Matrix-nolmer/src/CHOLMOD/Cholesky/t_cholmod_lsolve.c – Rev 3548
Rev 1820 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
/* ========================================================================== *//* === Cholesky/t_cholmod_lsolve ============================================ *//* ========================================================================== *//* -----------------------------------------------------------------------------* 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 Lx=b with unit or non-unit diagonal, or* solve LDx=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/* -------------------------------------------------------------------------- *//* define the method *//* -------------------------------------------------------------------------- */#ifdef LL/* LL': solve Lx=b with non-unit diagonal */#define FORM_NAME(prefix,rank) prefix ## ll_lsolve_ ## rank#elif defined (LD)/* LDL': solve LDx=b */#define FORM_NAME(prefix,rank) prefix ## ldl_ldsolve_ ## rank#else/* LDL': solve Lx=b with unit diagonal */#define FORM_NAME(prefix,rank) prefix ## ldl_lsolve_ ## 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 Lx=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 ;double *X = Y->x ;Int *Li = L->i ;Int *Lp = L->p ;Int *Lnz = L->nz ;Int j, n = L->n ;for (j = 0 ; j < n ; ){/* 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 (lnz < 4 || lnz != Lnz [j+1] + 1 || Li [p+1] != j+1){/* -------------------------------------------------------------- *//* solve with a single column of L *//* -------------------------------------------------------------- */double y = X [j] ;#ifdef LLy /= Lx [p] ;X [j] = y ;#elif defined (LD)X [j] = y / Lx [p] ;#endiffor (p++ ; p < pend ; p++){X [Li [p]] -= Lx [p] * y ;}j++ ;}else if (lnz != Lnz [j+2] + 2 || Li [p+2] != j+2){/* -------------------------------------------------------------- *//* solve with a supernode of two columns of L *//* -------------------------------------------------------------- */double y [2] ;Int q = Lp [j+1] ;#ifdef LLy [0] = X [j] / Lx [p] ;y [1] = (X [j+1] - Lx [p+1] * y [0]) / Lx [q] ;X [j ] = y [0] ;X [j+1] = y [1] ;#elif defined (LD)y [0] = X [j] ;y [1] = X [j+1] - Lx [p+1] * y [0] ;X [j ] = y [0] / Lx [p] ;X [j+1] = y [1] / Lx [q] ;#elsey [0] = X [j] ;y [1] = X [j+1] - Lx [p+1] * y [0] ;X [j+1] = y [1] ;#endiffor (p += 2, q++ ; p < pend ; p++, q++){X [Li [p]] -= Lx [p] * y [0] + Lx [q] * y [1] ;}j += 2 ;}else{/* -------------------------------------------------------------- *//* solve with a supernode of three columns of L *//* -------------------------------------------------------------- */double y [3] ;Int q = Lp [j+1] ;Int r = Lp [j+2] ;#ifdef LLy [0] = X [j] / Lx [p] ;y [1] = (X [j+1] - Lx [p+1] * y [0]) / Lx [q] ;y [2] = (X [j+2] - Lx [p+2] * y [0] - Lx [q+1] * y [1]) / Lx [r] ;X [j ] = y [0] ;X [j+1] = y [1] ;X [j+2] = y [2] ;#elif defined (LD)y [0] = X [j] ;y [1] = X [j+1] - Lx [p+1] * y [0] ;y [2] = X [j+2] - Lx [p+2] * y [0] - Lx [q+1] * y [1] ;X [j ] = y [0] / Lx [p] ;X [j+1] = y [1] / Lx [q] ;X [j+2] = y [2] / Lx [r] ;#elsey [0] = X [j] ;y [1] = X [j+1] - Lx [p+1] * y [0] ;y [2] = X [j+2] - Lx [p+2] * y [0] - Lx [q+1] * y [1] ;X [j+1] = y [1] ;X [j+2] = y [2] ;#endiffor (p += 3, q += 2, r++ ; p < pend ; p++, q++, r++){X [Li [p]] -= Lx [p] * y [0] + Lx [q] * y [1] + Lx [r] * y [2] ;}j += 3 ;}}}/* ========================================================================== *//* === LSOLVE (2) =========================================================== *//* ========================================================================== *//* Solve Lx=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 ;double *X = Y->x ;Int *Li = L->i ;Int *Lp = L->p ;Int *Lnz = L->nz ;Int j, n = L->n ;for (j = 0 ; j < n ; ){/* 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 (lnz < 4 || lnz != Lnz [j+1] + 1 || Li [p+1] != j+1){/* -------------------------------------------------------------- *//* solve with a single column of L *//* -------------------------------------------------------------- */double y [2] ;y [0] = X [2*j ] ;y [1] = X [2*j+1] ;#ifdef LLy [0] /= Lx [p] ;y [1] /= Lx [p] ;X [2*j ] = y [0] ;X [2*j+1] = y [1] ;#elif defined (LD)X [2*j ] = y [0] / Lx [p] ;X [2*j+1] = y [1] / Lx [p] ;#endiffor (p++ ; p < pend ; p++){Int i = 2 * Li [p] ;X [i ] -= Lx [p] * y [0] ;X [i+1] -= Lx [p] * y [1] ;}j++ ;}else if (lnz != Lnz [j+2] + 2 || Li [p+2] != j+2){/* -------------------------------------------------------------- *//* solve with a supernode of two columns of L *//* -------------------------------------------------------------- */double y [2][2] ;Int q = Lp [j+1] ;y [0][0] = X [2*j ] ;y [0][1] = X [2*j+1] ;#ifdef LLy [0][0] /= Lx [p] ;y [0][1] /= Lx [p] ;y [1][0] = (X [2*j+2] - Lx [p+1] * y [0][0]) / Lx [q] ;y [1][1] = (X [2*j+3] - Lx [p+1] * y [0][1]) / Lx [q] ;X [2*j ] = y [0][0] ;X [2*j+1] = y [0][1] ;X [2*j+2] = y [1][0] ;X [2*j+3] = y [1][1] ;#elif defined (LD)y [1][0] = X [2*j+2] - Lx [p+1] * y [0][0] ;y [1][1] = X [2*j+3] - Lx [p+1] * y [0][1] ;X [2*j ] = y [0][0] / Lx [p] ;X [2*j+1] = y [0][1] / Lx [p] ;X [2*j+2] = y [1][0] / Lx [q] ;X [2*j+3] = y [1][1] / Lx [q] ;#elsey [1][0] = X [2*j+2] - Lx [p+1] * y [0][0] ;y [1][1] = X [2*j+3] - Lx [p+1] * y [0][1] ;X [2*j+2] = y [1][0] ;X [2*j+3] = y [1][1] ;#endiffor (p += 2, q++ ; p < pend ; p++, q++){Int i = 2 * Li [p] ;X [i ] -= Lx [p] * y [0][0] + Lx [q] * y [1][0] ;X [i+1] -= Lx [p] * y [0][1] + Lx [q] * y [1][1] ;}j += 2 ;}else{/* -------------------------------------------------------------- *//* solve with a supernode of three columns of L *//* -------------------------------------------------------------- */double y [3][2] ;Int q = Lp [j+1] ;Int r = Lp [j+2] ;y [0][0] = X [2*j ] ;y [0][1] = X [2*j+1] ;#ifdef LLy [0][0] /= Lx [p] ;y [0][1] /= Lx [p] ;y [1][0] = (X[2*j+2] - Lx[p+1] * y[0][0]) / Lx [q] ;y [1][1] = (X[2*j+3] - Lx[p+1] * y[0][1]) / Lx [q] ;y [2][0] = (X[2*j+4] - Lx[p+2] * y[0][0] - Lx[q+1] * y[1][0])/Lx[r];y [2][1] = (X[2*j+5] - Lx[p+2] * y[0][1] - Lx[q+1] * y[1][1])/Lx[r];X [2*j ] = y [0][0] ;X [2*j+1] = y [0][1] ;X [2*j+2] = y [1][0] ;X [2*j+3] = y [1][1] ;X [2*j+4] = y [2][0] ;X [2*j+5] = y [2][1] ;#elif defined (LD)y [1][0] = X [2*j+2] - Lx [p+1] * y [0][0] ;y [1][1] = X [2*j+3] - Lx [p+1] * y [0][1] ;y [2][0] = X [2*j+4] - Lx [p+2] * y [0][0] - Lx [q+1] * y [1][0] ;y [2][1] = X [2*j+5] - Lx [p+2] * y [0][1] - Lx [q+1] * y [1][1] ;X [2*j ] = y [0][0] / Lx [p] ;X [2*j+1] = y [0][1] / Lx [p] ;X [2*j+2] = y [1][0] / Lx [q] ;X [2*j+3] = y [1][1] / Lx [q] ;X [2*j+4] = y [2][0] / Lx [r] ;X [2*j+5] = y [2][1] / Lx [r] ;#elsey [1][0] = X [2*j+2] - Lx [p+1] * y [0][0] ;y [1][1] = X [2*j+3] - Lx [p+1] * y [0][1] ;y [2][0] = X [2*j+4] - Lx [p+2] * y [0][0] - Lx [q+1] * y [1][0] ;y [2][1] = X [2*j+5] - Lx [p+2] * y [0][1] - Lx [q+1] * y [1][1] ;X [2*j+2] = y [1][0] ;X [2*j+3] = y [1][1] ;X [2*j+4] = y [2][0] ;X [2*j+5] = y [2][1] ;#endiffor (p += 3, q += 2, r++ ; p < pend ; p++, q++, r++){Int i = 2 * Li [p] ;X[i ] -= Lx[p] * y[0][0] + Lx[q] * y[1][0] + Lx[r] * y[2][0] ;X[i+1] -= Lx[p] * y[0][1] + Lx[q] * y[1][1] + Lx[r] * y[2][1] ;}j += 3 ;}}}/* ========================================================================== *//* === LSOLVE (3) =========================================================== *//* ========================================================================== *//* Solve Lx=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 ;double *X = Y->x ;Int *Li = L->i ;Int *Lp = L->p ;Int *Lnz = L->nz ;Int j, n = L->n ;for (j = 0 ; j < n ; ){/* 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 (lnz < 4 || lnz != Lnz [j+1] + 1 || Li [p+1] != j+1){/* -------------------------------------------------------------- *//* solve with a single column of L *//* -------------------------------------------------------------- */double y [3] ;y [0] = X [3*j ] ;y [1] = X [3*j+1] ;y [2] = X [3*j+2] ;#ifdef LLy [0] /= Lx [p] ;y [1] /= Lx [p] ;y [2] /= Lx [p] ;X [3*j ] = y [0] ;X [3*j+1] = y [1] ;X [3*j+2] = y [2] ;#elif defined (LD)X [3*j ] = y [0] / Lx [p] ;X [3*j+1] = y [1] / Lx [p] ;X [3*j+2] = y [2] / Lx [p] ;#endiffor (p++ ; p < pend ; p++){Int i = 3 * Li [p] ;double lx = Lx [p] ;X [i ] -= lx * y [0] ;X [i+1] -= lx * y [1] ;X [i+2] -= lx * y [2] ;}j++ ;}else if (lnz != Lnz [j+2] + 2 || Li [p+2] != j+2){/* -------------------------------------------------------------- *//* solve with a supernode of two columns of L *//* -------------------------------------------------------------- */double y [2][3] ;Int q = Lp [j+1] ;y [0][0] = X [3*j ] ;y [0][1] = X [3*j+1] ;y [0][2] = X [3*j+2] ;#ifdef LLy [0][0] /= Lx [p] ;y [0][1] /= Lx [p] ;y [0][2] /= Lx [p] ;y [1][0] = (X [3*j+3] - Lx [p+1] * y [0][0]) / Lx [q] ;y [1][1] = (X [3*j+4] - Lx [p+1] * y [0][1]) / Lx [q] ;y [1][2] = (X [3*j+5] - Lx [p+1] * y [0][2]) / Lx [q] ;X [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+4] = y [1][1] ;X [3*j+5] = y [1][2] ;#elif defined (LD)y [1][0] = X [3*j+3] - Lx [p+1] * y [0][0] ;y [1][1] = X [3*j+4] - Lx [p+1] * y [0][1] ;y [1][2] = X [3*j+5] - Lx [p+1] * y [0][2] ;X [3*j ] = y [0][0] / Lx [p] ;X [3*j+1] = y [0][1] / Lx [p] ;X [3*j+2] = y [0][2] / Lx [p] ;X [3*j+3] = y [1][0] / Lx [q] ;X [3*j+4] = y [1][1] / Lx [q] ;X [3*j+5] = y [1][2] / Lx [q] ;#elsey [1][0] = X [3*j+3] - Lx [p+1] * y [0][0] ;y [1][1] = X [3*j+4] - Lx [p+1] * y [0][1] ;y [1][2] = X [3*j+5] - Lx [p+1] * y [0][2] ;X [3*j+3] = y [1][0] ;X [3*j+4] = y [1][1] ;X [3*j+5] = y [1][2] ;#endiffor (p += 2, q++ ; p < pend ; p++, q++){Int i = 3 * Li [p] ;double lx [2] ;lx [0] = Lx [p] ;lx [1] = Lx [q] ;X [i ] -= lx [0] * y [0][0] + lx [1] * y [1][0] ;X [i+1] -= lx [0] * y [0][1] + lx [1] * y [1][1] ;X [i+2] -= lx [0] * y [0][2] + lx [1] * y [1][2] ;}j += 2 ;}else{/* -------------------------------------------------------------- *//* solve with a supernode of three columns of L *//* -------------------------------------------------------------- */double y [3][3] ;Int q = Lp [j+1] ;Int r = Lp [j+2] ;y [0][0] = X [3*j ] ;y [0][1] = X [3*j+1] ;y [0][2] = X [3*j+2] ;#ifdef LLy [0][0] /= Lx [p] ;y [0][1] /= Lx [p] ;y [0][2] /= Lx [p] ;y [1][0] = (X[3*j+3] - Lx[p+1] * y[0][0]) / Lx [q] ;y [1][1] = (X[3*j+4] - Lx[p+1] * y[0][1]) / Lx [q] ;y [1][2] = (X[3*j+5] - Lx[p+1] * y[0][2]) / Lx [q] ;y [2][0] = (X[3*j+6] - Lx[p+2] * y[0][0] - Lx[q+1] * y[1][0])/Lx[r];y [2][1] = (X[3*j+7] - Lx[p+2] * y[0][1] - Lx[q+1] * y[1][1])/Lx[r];y [2][2] = (X[3*j+8] - Lx[p+2] * y[0][2] - Lx[q+1] * y[1][2])/Lx[r];X [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+4] = y [1][1] ;X [3*j+5] = y [1][2] ;X [3*j+6] = y [2][0] ;X [3*j+7] = y [2][1] ;X [3*j+8] = y [2][2] ;#elif defined (LD)y [1][0] = X [3*j+3] - Lx [p+1] * y [0][0] ;y [1][1] = X [3*j+4] - Lx [p+1] * y [0][1] ;y [1][2] = X [3*j+5] - Lx [p+1] * y [0][2] ;y [2][0] = X [3*j+6] - Lx [p+2] * y [0][0] - Lx [q+1] * y [1][0] ;y [2][1] = X [3*j+7] - Lx [p+2] * y [0][1] - Lx [q+1] * y [1][1] ;y [2][2] = X [3*j+8] - Lx [p+2] * y [0][2] - Lx [q+1] * y [1][2] ;X [3*j ] = y [0][0] / Lx [p] ;X [3*j+1] = y [0][1] / Lx [p] ;X [3*j+2] = y [0][2] / Lx [p] ;X [3*j+3] = y [1][0] / Lx [q] ;X [3*j+4] = y [1][1] / Lx [q] ;X [3*j+5] = y [1][2] / Lx [q] ;X [3*j+6] = y [2][0] / Lx [r] ;X [3*j+7] = y [2][1] / Lx [r] ;X [3*j+8] = y [2][2] / Lx [r] ;#elsey [1][0] = X [3*j+3] - Lx [p+1] * y [0][0] ;y [1][1] = X [3*j+4] - Lx [p+1] * y [0][1] ;y [1][2] = X [3*j+5] - Lx [p+1] * y [0][2] ;y [2][0] = X [3*j+6] - Lx [p+2] * y [0][0] - Lx [q+1] * y [1][0] ;y [2][1] = X [3*j+7] - Lx [p+2] * y [0][1] - Lx [q+1] * y [1][1] ;y [2][2] = X [3*j+8] - Lx [p+2] * y [0][2] - Lx [q+1] * y [1][2] ;X [3*j+3] = y [1][0] ;X [3*j+4] = y [1][1] ;X [3*j+5] = y [1][2] ;X [3*j+6] = y [2][0] ;X [3*j+7] = y [2][1] ;X [3*j+8] = y [2][2] ;#endiffor (p += 3, q += 2, r++ ; p < pend ; p++, q++, r++){Int i = 3 * Li [p] ;double lx [3] ;lx [0] = Lx [p] ;lx [1] = Lx [q] ;lx [2] = Lx [r] ;X [i ] -= lx[0] * y[0][0] + lx[1] * y[1][0] + lx[2] * y[2][0] ;X [i+1] -= lx[0] * y[0][1] + lx[1] * y[1][1] + lx[2] * y[2][1] ;X [i+2] -= lx[0] * y[0][2] + lx[1] * y[1][2] + lx[2] * y[2][2] ;}j += 3 ;}}}/* ========================================================================== *//* === LSOLVE (4) =========================================================== *//* ========================================================================== *//* Solve Lx=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 ;double *X = Y->x ;Int *Li = L->i ;Int *Lp = L->p ;Int *Lnz = L->nz ;Int j, n = L->n ;for (j = 0 ; j < n ; ){/* 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 (lnz < 4 || lnz != Lnz [j+1] + 1 || Li [p+1] != j+1){/* -------------------------------------------------------------- *//* solve with a single column of L *//* -------------------------------------------------------------- */double y [4] ;y [0] = X [4*j ] ;y [1] = X [4*j+1] ;y [2] = X [4*j+2] ;y [3] = X [4*j+3] ;#ifdef LLy [0] /= Lx [p] ;y [1] /= Lx [p] ;y [2] /= Lx [p] ;y [3] /= Lx [p] ;X [4*j ] = y [0] ;X [4*j+1] = y [1] ;X [4*j+2] = y [2] ;X [4*j+3] = y [3] ;#elif defined (LD)X [4*j ] = y [0] / Lx [p] ;X [4*j+1] = y [1] / Lx [p] ;X [4*j+2] = y [2] / Lx [p] ;X [4*j+3] = y [3] / Lx [p] ;#endiffor (p++ ; p < pend ; p++){Int i = 4 * Li [p] ;double lx = Lx [p] ;X [i ] -= lx * y [0] ;X [i+1] -= lx * y [1] ;X [i+2] -= lx * y [2] ;X [i+3] -= lx * y [3] ;}j++ ;}else if (lnz != Lnz [j+2] + 2 || Li [p+2] != j+2){/* -------------------------------------------------------------- *//* solve with a supernode of two columns of L *//* -------------------------------------------------------------- */double y [2][4] ;Int q = Lp [j+1] ;y [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] ;#ifdef LLy [0][0] /= Lx [p] ;y [0][1] /= Lx [p] ;y [0][2] /= Lx [p] ;y [0][3] /= Lx [p] ;y [1][0] = (X [4*j+4] - Lx [p+1] * y [0][0]) / Lx [q] ;y [1][1] = (X [4*j+5] - Lx [p+1] * y [0][1]) / Lx [q] ;y [1][2] = (X [4*j+6] - Lx [p+1] * y [0][2]) / Lx [q] ;y [1][3] = (X [4*j+7] - Lx [p+1] * y [0][3]) / Lx [q] ;X [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+5] = y [1][1] ;X [4*j+6] = y [1][2] ;X [4*j+7] = y [1][3] ;#elif defined (LD)y [1][0] = X [4*j+4] - Lx [p+1] * y [0][0] ;y [1][1] = X [4*j+5] - Lx [p+1] * y [0][1] ;y [1][2] = X [4*j+6] - Lx [p+1] * y [0][2] ;y [1][3] = X [4*j+7] - Lx [p+1] * y [0][3] ;X [4*j ] = y [0][0] / Lx [p] ;X [4*j+1] = y [0][1] / Lx [p] ;X [4*j+2] = y [0][2] / Lx [p] ;X [4*j+3] = y [0][3] / Lx [p] ;X [4*j+4] = y [1][0] / Lx [q] ;X [4*j+5] = y [1][1] / Lx [q] ;X [4*j+6] = y [1][2] / Lx [q] ;X [4*j+7] = y [1][3] / Lx [q] ;#elsey [1][0] = X [4*j+4] - Lx [p+1] * y [0][0] ;y [1][1] = X [4*j+5] - Lx [p+1] * y [0][1] ;y [1][2] = X [4*j+6] - Lx [p+1] * y [0][2] ;y [1][3] = X [4*j+7] - Lx [p+1] * y [0][3] ;X [4*j+4] = y [1][0] ;X [4*j+5] = y [1][1] ;X [4*j+6] = y [1][2] ;X [4*j+7] = y [1][3] ;#endiffor (p += 2, q++ ; p < pend ; p++, q++){Int i = 4 * Li [p] ;double lx [2] ;lx [0] = Lx [p] ;lx [1] = Lx [q] ;X [i ] -= lx [0] * y [0][0] + lx [1] * y [1][0] ;X [i+1] -= lx [0] * y [0][1] + lx [1] * y [1][1] ;X [i+2] -= lx [0] * y [0][2] + lx [1] * y [1][2] ;X [i+3] -= lx [0] * y [0][3] + lx [1] * y [1][3] ;}j += 2 ;}else{/* -------------------------------------------------------------- *//* solve with a supernode of three columns of L *//* -------------------------------------------------------------- */double y [3][4] ;Int q = Lp [j+1] ;Int r = Lp [j+2] ;y [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] ;#ifdef LLy [0][0] /= Lx [p] ;y [0][1] /= Lx [p] ;y [0][2] /= Lx [p] ;y [0][3] /= Lx [p] ;y [1][0] = (X[4*j+4] - Lx[p+1] * y[0][0]) / Lx [q] ;y [1][1] = (X[4*j+5] - Lx[p+1] * y[0][1]) / Lx [q] ;y [1][2] = (X[4*j+6] - Lx[p+1] * y[0][2]) / Lx [q] ;y [1][3] = (X[4*j+7] - Lx[p+1] * y[0][3]) / Lx [q] ;y [2][0] = (X[4*j+8] - Lx[p+2] * y[0][0] - Lx[q+1] * y[1][0])/Lx[r];y [2][1] = (X[4*j+9] - Lx[p+2] * y[0][1] - Lx[q+1] * y[1][1])/Lx[r];y [2][2] = (X[4*j+10]- Lx[p+2] * y[0][2] - Lx[q+1] * y[1][2])/Lx[r];y [2][3] = (X[4*j+11]- Lx[p+2] * y[0][3] - Lx[q+1] * y[1][3])/Lx[r];X [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+5] = y [1][1] ;X [4*j+6] = y [1][2] ;X [4*j+7] = y [1][3] ;X [4*j+8] = y [2][0] ;X [4*j+9] = y [2][1] ;X [4*j+10]= y [2][2] ;X [4*j+11]= y [2][3] ;#elif defined (LD)y [1][0] = X [4*j+4] - Lx [p+1] * y [0][0] ;y [1][1] = X [4*j+5] - Lx [p+1] * y [0][1] ;y [1][2] = X [4*j+6] - Lx [p+1] * y [0][2] ;y [1][3] = X [4*j+7] - Lx [p+1] * y [0][3] ;y [2][0] = X [4*j+8] - Lx [p+2] * y [0][0] - Lx [q+1] * y [1][0] ;y [2][1] = X [4*j+9] - Lx [p+2] * y [0][1] - Lx [q+1] * y [1][1] ;y [2][2] = X [4*j+10]- Lx [p+2] * y [0][2] - Lx [q+1] * y [1][2] ;y [2][3] = X [4*j+11]- Lx [p+2] * y [0][3] - Lx [q+1] * y [1][3] ;X [4*j ] = y [0][0] / Lx [p] ;X [4*j+1] = y [0][1] / Lx [p] ;X [4*j+2] = y [0][2] / Lx [p] ;X [4*j+3] = y [0][3] / Lx [p] ;X [4*j+4] = y [1][0] / Lx [q] ;X [4*j+5] = y [1][1] / Lx [q] ;X [4*j+6] = y [1][2] / Lx [q] ;X [4*j+7] = y [1][3] / Lx [q] ;X [4*j+8] = y [2][0] / Lx [r] ;X [4*j+9] = y [2][1] / Lx [r] ;X [4*j+10]= y [2][2] / Lx [r] ;X [4*j+11]= y [2][3] / Lx [r] ;#elsey [1][0] = X [4*j+4] - Lx [p+1] * y [0][0] ;y [1][1] = X [4*j+5] - Lx [p+1] * y [0][1] ;y [1][2] = X [4*j+6] - Lx [p+1] * y [0][2] ;y [1][3] = X [4*j+7] - Lx [p+1] * y [0][3] ;y [2][0] = X [4*j+8] - Lx [p+2] * y [0][0] - Lx [q+1] * y [1][0] ;y [2][1] = X [4*j+9] - Lx [p+2] * y [0][1] - Lx [q+1] * y [1][1] ;y [2][2] = X [4*j+10]- Lx [p+2] * y [0][2] - Lx [q+1] * y [1][2] ;y [2][3] = X [4*j+11]- Lx [p+2] * y [0][3] - Lx [q+1] * y [1][3] ;X [4*j+4] = y [1][0] ;X [4*j+5] = y [1][1] ;X [4*j+6] = y [1][2] ;X [4*j+7] = y [1][3] ;X [4*j+8] = y [2][0] ;X [4*j+9] = y [2][1] ;X [4*j+10]= y [2][2] ;X [4*j+11]= y [2][3] ;#endiffor (p += 3, q += 2, r++ ; p < pend ; p++, q++, r++){Int i = 4 * Li [p] ;double lx [3] ;lx [0] = Lx [p] ;lx [1] = Lx [q] ;lx [2] = Lx [r] ;X [i ] -= lx[0] * y[0][0] + lx[1] * y[1][0] + lx[2] * y[2][0] ;X [i+1] -= lx[0] * y[0][1] + lx[1] * y[1][1] + lx[2] * y[2][1] ;X [i+2] -= lx[0] * y[0][2] + lx[1] * y[1][2] + lx[2] * y[2][2] ;X [i+3] -= lx[0] * y[0][3] + lx[1] * y[1][3] + lx[2] * 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 REALdouble 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 = 0 ; j < n ; 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 LL/* y /= Lx [p] ; *//* X [j] = y ; */DIV_REAL (yx,yz,0, yx,yz,0, Lx,p) ;ASSIGN (Xx,Xz,j, yx,yz,0) ;#elif defined (LD)/* X [j] = y / Lx [p] ; */DIV_REAL (Xx,Xz,j, yx,yz,0, Lx,p) ;#endiffor (p++ ; p < pend ; p++){/* X [Li [p]] -= Lx [p] * y ; */i = Li [p] ;MULTSUB (Xx,Xz,i, Lx,Lz,p, yx,yz,0) ;}}#endif}/* prepare for the next inclusion of this file in cholmod_solve.c */#undef LL#undef LD