Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_solve.c
Go to the documentation of this file.
```00001 /* ========================================================================== */
00002 /* === Cholesky/cholmod_solve =============================================== */
00003 /* ========================================================================== */
00004
00005 /* -----------------------------------------------------------------------------
00006  * CHOLMOD/Cholesky Module.  Copyright (C) 2005-2006, Timothy A. Davis
00008  * Lesser General Public License.  See lesser.txt for a text of the license.
00009  * CHOLMOD is also available under other licenses; contact authors for details.
00010  * http://www.cise.ufl.edu/research/sparse
00011  * -------------------------------------------------------------------------- */
00012
00013 /* Solve one of the following systems:
00014  *
00015  *  Ax=b      0: CHOLMOD_A  also applies the permutation L->Perm
00016  *  LDL'x=b     1: CHOLMOD_LDLt does not apply L->Perm
00017  *  LDx=b     2: CHOLMOD_LD
00018  *  DL'x=b      3: CHOLMOD_DLt
00019  *  Lx=b      4: CHOLMOD_L
00020  *  L'x=b     5: CHOLMOD_Lt
00021  *  Dx=b      6: CHOLMOD_D
00022  *  x=Pb      7: CHOLMOD_P  apply a permutation (P is L->Perm)
00023  *  x=P'b     8: CHOLMOD_Pt apply an inverse permutation
00024  *
00025  * The factorization can be simplicial LDL', simplicial LL', or supernodal LL'.
00026  * For an LL' factorization, D is the identity matrix.  Thus CHOLMOD_LD and
00027  * CHOLMOD_L solve the same system if an LL' factorization was performed,
00028  * for example.
00029  *
00030  * The supernodal solver uses BLAS routines dtrsv, dgemv, dtrsm, and dgemm,
00031  * or their complex counterparts ztrsv, zgemv, ztrsm, and zgemm.
00032  *
00033  * If both L and B are real, then X is returned real.  If either is complex
00034  * or zomplex, X is returned as either complex or zomplex, depending on the
00035  * Common->prefer_zomplex parameter.
00036  *
00037  * Supports any numeric xtype (pattern-only matrices not supported).
00038  *
00039  * This routine does not check to see if the diagonal of L or D is zero,
00040  * because sometimes a partial solve can be done with indefinite or singular
00041  * matrix.  If you wish to check in your own code, test L->minor.  If
00042  * L->minor == L->n, then the matrix has no zero diagonal entries.
00043  * If k = L->minor < L->n, then L(k,k) is zero for an LL' factorization, or
00044  * D(k,k) is zero for an LDL' factorization.
00045  *
00046  * This routine returns X as NULL only if it runs out of memory.  If L is
00047  * indefinite or singular, then X may contain Inf's or NaN's, but it will
00048  * exist on output.
00049  */
00050
00051 #ifndef NCHOLESKY
00052
00053 #include "amesos_cholmod_internal.h"
00054 #include "amesos_cholmod_cholesky.h"
00055
00056
00057 /* ========================================================================== */
00058 /* === TEMPLATE ============================================================= */
00059 /* ========================================================================== */
00060
00061 #define REAL
00062 #include "amesos_t_cholmod_solve.c"
00063
00064 #define COMPLEX
00065 #include "amesos_t_cholmod_solve.c"
00066
00067 #define ZOMPLEX
00068 #include "amesos_t_cholmod_solve.c"
00069
00070 /* ========================================================================== */
00071 /* === Permutation macro ==================================================== */
00072 /* ========================================================================== */
00073
00074 /* If Perm is NULL, it is interpretted as the identity permutation */
00075
00076 #define P(k) ((Perm == NULL) ? (k) : Perm [k])
00077
00078
00079 /* ========================================================================== */
00080 /* === perm ================================================================= */
00081 /* ========================================================================== */
00082
00083 /* Y = B (P (1:nrow), k1 : min (k1+ncols,ncol)-1) where B is nrow-by-ncol.
00084  *
00085  * Creates a permuted copy of a contiguous set of columns of B.
00086  * Y is already allocated on input.  Y must be of sufficient size.  Let nk be
00087  * the number of columns accessed in B.  Y->xtype determines the complexity of
00088  * the result.
00089  *
00090  * If B is real and Y is complex (or zomplex), only the real part of B is
00091  * copied into Y.  The imaginary part of Y is set to zero.
00092  *
00093  * If B is complex (or zomplex) and Y is real, both the real and imaginary and
00094  * parts of B are returned in Y.  Y is returned as nrow-by-2*nk. The even
00095  * columns of Y contain the real part of B and the odd columns contain the
00096  * imaginary part of B.  Y->nzmax must be >= 2*nrow*nk.  Otherise, Y is
00097  * returned as nrow-by-nk with leading dimension nrow.  Y->nzmax must be >=
00098  * nrow*nk.
00099  *
00100  * The case where the input (B) is real and the output (Y) is zomplex is
00101  * not used.
00102  */
00103
00104 static void amesos_perm
00105 (
00106     /* ---- input ---- */
00107     cholmod_dense *B, /* input matrix B */
00108     Int *Perm,    /* optional input permutation (can be NULL) */
00109     Int k1,   /* first column of B to copy */
00110     Int ncols,    /* last column to copy is min(k1+ncols,B->ncol)-1 */
00111     /* ---- in/out --- */
00112     cholmod_dense *Y  /* output matrix Y, already allocated */
00113 )
00114 {
00115     double *Yx, *Yz, *Bx, *Bz ;
00116     Int k2, nk, p, k, j, nrow, ncol, d, dual, dj, j2 ;
00117
00118     /* ---------------------------------------------------------------------- */
00119     /* get inputs */
00120     /* ---------------------------------------------------------------------- */
00121
00122     ncol = B->ncol ;
00123     nrow = B->nrow ;
00124     k2 = MIN (k1+ncols, ncol) ;
00125     nk = MAX (k2 - k1, 0) ;
00126     dual = (Y->xtype == CHOLMOD_REAL && B->xtype != CHOLMOD_REAL) ? 2 : 1 ;
00127     d = B->d ;
00128     Bx = B->x ;
00129     Bz = B->z ;
00130     Yx = Y->x ;
00131     Yz = Y->z ;
00132     Y->nrow = nrow ;
00133     Y->ncol = dual*nk ;
00134     Y->d = nrow ;
00135     ASSERT (((Int) Y->nzmax) >= nrow*nk*dual) ;
00136
00137     /* ---------------------------------------------------------------------- */
00138     /* Y = B (P (1:nrow), k1:k2-1) */
00139     /* ---------------------------------------------------------------------- */
00140
00141     switch (Y->xtype)
00142     {
00143
00144   case CHOLMOD_REAL:
00145
00146       switch (B->xtype)
00147       {
00148
00149     case CHOLMOD_REAL:
00150         /* Y real, B real */
00151         for (j = k1 ; j < k2 ; j++)
00152         {
00153       dj = d*j ;
00154       j2 = nrow * (j-k1) ;
00155       for (k = 0 ; k < nrow ; k++)
00156       {
00157           p = P(k) + dj ;
00158           Yx [k + j2] = Bx [p] ;    /* real */
00159       }
00160         }
00161         break ;
00162
00163     case CHOLMOD_COMPLEX:
00164         /* Y real, B complex. Y is nrow-by-2*nk */
00165         for (j = k1 ; j < k2 ; j++)
00166         {
00167       dj = d*j ;
00168       j2 = nrow * 2 * (j-k1) ;
00169       for (k = 0 ; k < nrow ; k++)
00170       {
00171           p = P(k) + dj ;
00172           Yx [k + j2       ] = Bx [2*p  ] ; /* real */
00173           Yx [k + j2 + nrow] = Bx [2*p+1] ; /* imag */
00174       }
00175         }
00176         break ;
00177
00178     case CHOLMOD_ZOMPLEX:
00179         /* Y real, B zomplex. Y is nrow-by-2*nk */
00180         for (j = k1 ; j < k2 ; j++)
00181         {
00182       dj = d*j ;
00183       j2 = nrow * 2 * (j-k1) ;
00184       for (k = 0 ; k < nrow ; k++)
00185       {
00186           p = P(k) + dj ;
00187           Yx [k + j2       ] = Bx [p] ; /* real */
00188           Yx [k + j2 + nrow] = Bz [p] ; /* imag */
00189       }
00190         }
00191         break ;
00192
00193       }
00194       break ;
00195
00196   case CHOLMOD_COMPLEX:
00197
00198       switch (B->xtype)
00199       {
00200
00201     case CHOLMOD_REAL:
00202         /* Y complex, B real */
00203         for (j = k1 ; j < k2 ; j++)
00204         {
00205       dj = d*j ;
00206       j2 = nrow * 2 * (j-k1) ;
00207       for (k = 0 ; k < nrow ; k++)
00208       {
00209           p = P(k) + dj ;
00210           Yx [2*k   + j2] = Bx [p] ;    /* real */
00211           Yx [2*k+1 + j2] = 0 ;   /* imag */
00212       }
00213         }
00214         break ;
00215
00216     case CHOLMOD_COMPLEX:
00217         /* Y complex, B complex */
00218         for (j = k1 ; j < k2 ; j++)
00219         {
00220       dj = d*j ;
00221       j2 = nrow * 2 * (j-k1) ;
00222       for (k = 0 ; k < nrow ; k++)
00223       {
00224           p = P(k) + dj ;
00225           Yx [2*k   + j2] = Bx [2*p  ] ;  /* real */
00226           Yx [2*k+1 + j2] = Bx [2*p+1] ;  /* imag */
00227       }
00228         }
00229         break ;
00230
00231     case CHOLMOD_ZOMPLEX:
00232         /* Y complex, B zomplex */
00233         for (j = k1 ; j < k2 ; j++)
00234         {
00235       dj = d*j ;
00236       j2 = nrow * 2 * (j-k1) ;
00237       for (k = 0 ; k < nrow ; k++)
00238       {
00239           p = P(k) + dj ;
00240           Yx [2*k   + j2] = Bx [p] ;    /* real */
00241           Yx [2*k+1 + j2] = Bz [p] ;    /* imag */
00242       }
00243         }
00244         break ;
00245
00246       }
00247       break ;
00248
00249   case CHOLMOD_ZOMPLEX:
00250
00251       switch (B->xtype)
00252       {
00253
00254 #if 0
00255     case CHOLMOD_REAL:
00256         /* this case is not used */
00257         break ;
00258 #endif
00259
00260     case CHOLMOD_COMPLEX:
00261         /* Y zomplex, B complex */
00262         for (j = k1 ; j < k2 ; j++)
00263         {
00264       dj = d*j ;
00265       j2 = nrow * (j-k1) ;
00266       for (k = 0 ; k < nrow ; k++)
00267       {
00268           p = P(k) + dj ;
00269           Yx [k + j2] = Bx [2*p  ] ;    /* real */
00270           Yz [k + j2] = Bx [2*p+1] ;    /* imag */
00271       }
00272         }
00273         break ;
00274
00275     case CHOLMOD_ZOMPLEX:
00276         /* Y zomplex, B zomplex */
00277         for (j = k1 ; j < k2 ; j++)
00278         {
00279       dj = d*j ;
00280       j2 = nrow * (j-k1) ;
00281       for (k = 0 ; k < nrow ; k++)
00282       {
00283           p = P(k) + dj ;
00284           Yx [k + j2] = Bx [p] ;    /* real */
00285           Yz [k + j2] = Bz [p] ;    /* imag */
00286       }
00287         }
00288         break ;
00289
00290       }
00291       break ;
00292
00293     }
00294 }
00295
00296
00297 /* ========================================================================== */
00298 /* === iperm ================================================================ */
00299 /* ========================================================================== */
00300
00301 /* X (P (1:nrow), k1 : min (k1+ncols,ncol)-1) = Y where X is nrow-by-ncol.
00302  *
00303  * Copies and permutes Y into a contiguous set of columns of X.  X is already
00304  * allocated on input.  Y must be of sufficient size.  Let nk be the number
00305  * of columns accessed in X.  X->xtype determines the complexity of the result.
00306  *
00307  * If X is real and Y is complex (or zomplex), only the real part of B is
00308  * copied into X.  The imaginary part of Y is ignored.
00309  *
00310  * If X is complex (or zomplex) and Y is real, both the real and imaginary and
00311  * parts of Y are returned in X.  Y is nrow-by-2*nk. The even
00312  * columns of Y contain the real part of B and the odd columns contain the
00313  * imaginary part of B.  Y->nzmax must be >= 2*nrow*nk.  Otherise, Y is
00314  * nrow-by-nk with leading dimension nrow.  Y->nzmax must be >= nrow*nk.
00315  *
00316  * The case where the input (Y) is complex and the output (X) is real,
00317  * and the case where the input (Y) is zomplex and the output (X) is real,
00318  * are not used.
00319  */
00320
00321 static void amesos_iperm
00322 (
00323     /* ---- input ---- */
00324     cholmod_dense *Y, /* input matrix Y */
00325     Int *Perm,    /* optional input permutation (can be NULL) */
00326     Int k1,   /* first column of B to copy */
00327     Int ncols,    /* last column to copy is min(k1+ncols,B->ncol)-1 */
00328     /* ---- in/out --- */
00329     cholmod_dense *X  /* output matrix X, already allocated */
00330 )
00331 {
00332     double *Yx, *Yz, *Xx, *Xz ;
00333     Int k2, nk, p, k, j, nrow, ncol, d, dj, j2 ;
00334
00335     /* ---------------------------------------------------------------------- */
00336     /* get inputs */
00337     /* ---------------------------------------------------------------------- */
00338
00339     ncol = X->ncol ;
00340     nrow = X->nrow ;
00341     k2 = MIN (k1+ncols, ncol) ;
00342     nk = MAX (k2 - k1, 0) ;
00343     d = X->d ;
00344     Xx = X->x ;
00345     Xz = X->z ;
00346     Yx = Y->x ;
00347     Yz = Y->z ;
00348     ASSERT (((Int) Y->nzmax) >= nrow*nk*
00349       ((X->xtype != CHOLMOD_REAL && Y->xtype == CHOLMOD_REAL) ? 2:1)) ;
00350
00351     /* ---------------------------------------------------------------------- */
00352     /* X (P (1:nrow), k1:k2-1) = Y */
00353     /* ---------------------------------------------------------------------- */
00354
00355     switch (Y->xtype)
00356     {
00357
00358   case CHOLMOD_REAL:
00359
00360       switch (X->xtype)
00361       {
00362
00363     case CHOLMOD_REAL:
00364         /* Y real, X real */
00365         for (j = k1 ; j < k2 ; j++)
00366         {
00367       dj = d*j ;
00368       j2 = nrow * (j-k1) ;
00369       for (k = 0 ; k < nrow ; k++)
00370       {
00371           p = P(k) + dj ;
00372           Xx [p] = Yx [k + j2] ;    /* real */
00373       }
00374         }
00375         break ;
00376
00377     case CHOLMOD_COMPLEX:
00378         /* Y real, X complex. Y is nrow-by-2*nk */
00379         for (j = k1 ; j < k2 ; j++)
00380         {
00381       dj = d*j ;
00382       j2 = nrow * 2 * (j-k1) ;
00383       for (k = 0 ; k < nrow ; k++)
00384       {
00385           p = P(k) + dj ;
00386           Xx [2*p  ] = Yx [k + j2       ] ; /* real */
00387           Xx [2*p+1] = Yx [k + j2 + nrow] ; /* imag */
00388       }
00389         }
00390         break ;
00391
00392     case CHOLMOD_ZOMPLEX:
00393         /* Y real, X zomplex. Y is nrow-by-2*nk */
00394         for (j = k1 ; j < k2 ; j++)
00395         {
00396       dj = d*j ;
00397       j2 = nrow * 2 * (j-k1) ;
00398       for (k = 0 ; k < nrow ; k++)
00399       {
00400           p = P(k) + dj ;
00401           Xx [p] = Yx [k + j2       ] ; /* real */
00402           Xz [p] = Yx [k + j2 + nrow] ; /* imag */
00403       }
00404         }
00405         break ;
00406
00407       }
00408       break ;
00409
00410   case CHOLMOD_COMPLEX:
00411
00412       switch (X->xtype)
00413       {
00414
00415 #if 0
00416     case CHOLMOD_REAL:
00417         /* this case is not used */
00418         break ;
00419 #endif
00420
00421     case CHOLMOD_COMPLEX:
00422         /* Y complex, X complex */
00423         for (j = k1 ; j < k2 ; j++)
00424         {
00425       dj = d*j ;
00426       j2 = nrow * 2 * (j-k1) ;
00427       for (k = 0 ; k < nrow ; k++)
00428       {
00429           p = P(k) + dj ;
00430           Xx [2*p  ] = Yx [2*k   + j2] ;  /* real */
00431           Xx [2*p+1] = Yx [2*k+1 + j2] ;  /* imag */
00432       }
00433         }
00434         break ;
00435
00436     case CHOLMOD_ZOMPLEX:
00437         /* Y complex, X zomplex */
00438         for (j = k1 ; j < k2 ; j++)
00439         {
00440       dj = d*j ;
00441       j2 = nrow * 2 * (j-k1) ;
00442       for (k = 0 ; k < nrow ; k++)
00443       {
00444           p = P(k) + dj ;
00445           Xx [p] = Yx [2*k   + j2] ;    /* real */
00446           Xz [p] = Yx [2*k+1 + j2] ;    /* imag */
00447       }
00448         }
00449         break ;
00450
00451       }
00452       break ;
00453
00454   case CHOLMOD_ZOMPLEX:
00455
00456       switch (X->xtype)
00457       {
00458
00459 #if 0
00460     case CHOLMOD_REAL:
00461         /* this case is not used */
00462         break ;
00463 #endif
00464
00465     case CHOLMOD_COMPLEX:
00466         /* Y zomplex, X complex */
00467         for (j = k1 ; j < k2 ; j++)
00468         {
00469       dj = d*j ;
00470       j2 = nrow * (j-k1) ;
00471       for (k = 0 ; k < nrow ; k++)
00472       {
00473           p = P(k) + dj ;
00474           Xx [2*p  ] = Yx [k + j2] ;    /* real */
00475           Xx [2*p+1] = Yz [k + j2] ;    /* imag */
00476       }
00477         }
00478         break ;
00479
00480     case CHOLMOD_ZOMPLEX:
00481         /* Y zomplex, X zomplex */
00482         for (j = k1 ; j < k2 ; j++)
00483         {
00484       dj = d*j ;
00485       j2 = nrow * (j-k1) ;
00486       for (k = 0 ; k < nrow ; k++)
00487       {
00488           p = P(k) + dj ;
00489           Xx [p] = Yx [k + j2] ;    /* real */
00490           Xz [p] = Yz [k + j2] ;    /* imag */
00491       }
00492         }
00493         break ;
00494
00495       }
00496       break ;
00497
00498     }
00499 }
00500
00501
00502 /* ========================================================================== */
00503 /* === ptrans =============================================================== */
00504 /* ========================================================================== */
00505
00506 /* Y = B (P (1:nrow), k1 : min (k1+ncols,ncol)-1)' where B is nrow-by-ncol.
00507  *
00508  * Creates a permuted and transposed copy of a contiguous set of columns of B.
00509  * Y is already allocated on input.  Y must be of sufficient size.  Let nk be
00510  * the number of columns accessed in B.  Y->xtype determines the complexity of
00511  * the result.
00512  *
00513  * If B is real and Y is complex (or zomplex), only the real part of B is
00514  * copied into Y.  The imaginary part of Y is set to zero.
00515  *
00516  * If B is complex (or zomplex) and Y is real, both the real and imaginary and
00517  * parts of B are returned in Y.  Y is returned as 2*nk-by-nrow. The even
00518  * rows of Y contain the real part of B and the odd rows contain the
00519  * imaginary part of B.  Y->nzmax must be >= 2*nrow*nk.  Otherise, Y is
00520  * returned as nk-by-nrow with leading dimension nk.  Y->nzmax must be >=
00521  * nrow*nk.
00522  *
00523  * The array transpose is performed, not the complex conjugate transpose.
00524  */
00525
00526 static void amesos_ptrans
00527 (
00528     /* ---- input ---- */
00529     cholmod_dense *B, /* input matrix B */
00530     Int *Perm,    /* optional input permutation (can be NULL) */
00531     Int k1,   /* first column of B to copy */
00532     Int ncols,    /* last column to copy is min(k1+ncols,B->ncol)-1 */
00533     /* ---- in/out --- */
00534     cholmod_dense *Y  /* output matrix Y, already allocated */
00535 )
00536 {
00537     double *Yx, *Yz, *Bx, *Bz ;
00538     Int k2, nk, p, k, j, nrow, ncol, d, dual, dj, j2 ;
00539
00540     /* ---------------------------------------------------------------------- */
00541     /* get inputs */
00542     /* ---------------------------------------------------------------------- */
00543
00544     ncol = B->ncol ;
00545     nrow = B->nrow ;
00546     k2 = MIN (k1+ncols, ncol) ;
00547     nk = MAX (k2 - k1, 0) ;
00548     dual = (Y->xtype == CHOLMOD_REAL && B->xtype != CHOLMOD_REAL) ? 2 : 1 ;
00549     d = B->d ;
00550     Bx = B->x ;
00551     Bz = B->z ;
00552     Yx = Y->x ;
00553     Yz = Y->z ;
00554     Y->nrow = dual*nk ;
00555     Y->ncol = nrow ;
00556     Y->d = dual*nk ;
00557     ASSERT (((Int) Y->nzmax) >= nrow*nk*dual) ;
00558
00559     /* ---------------------------------------------------------------------- */
00560     /* Y = B (P (1:nrow), k1:k2-1)' */
00561     /* ---------------------------------------------------------------------- */
00562
00563     switch (Y->xtype)
00564     {
00565
00566   case CHOLMOD_REAL:
00567
00568       switch (B->xtype)
00569       {
00570
00571     case CHOLMOD_REAL:
00572         /* Y real, B real  */
00573         for (j = k1 ; j < k2 ; j++)
00574         {
00575       dj = d*j ;
00576       j2 = j-k1 ;
00577       for (k = 0 ; k < nrow ; k++)
00578       {
00579           p = P(k) + dj ;
00580           Yx [j2 + k*nk] = Bx [p] ;   /* real */
00581       }
00582         }
00583         break ;
00584
00585     case CHOLMOD_COMPLEX:
00586         /* Y real, B complex. Y is 2*nk-by-nrow */
00587         for (j = k1 ; j < k2 ; j++)
00588         {
00589       dj = d*j ;
00590       j2 = 2*(j-k1) ;
00591       for (k = 0 ; k < nrow ; k++)
00592       {
00593           p = P(k) + dj ;
00594           Yx [j2   + k*2*nk] = Bx [2*p  ] ; /* real */
00595           Yx [j2+1 + k*2*nk] = Bx [2*p+1] ; /* imag */
00596       }
00597         }
00598         break ;
00599
00600     case CHOLMOD_ZOMPLEX:
00601         /* Y real, B zomplex. Y is 2*nk-by-nrow */
00602         for (j = k1 ; j < k2 ; j++)
00603         {
00604       dj = d*j ;
00605       j2 = 2*(j-k1) ;
00606       for (k = 0 ; k < nrow ; k++)
00607       {
00608           p = P(k) + dj ;
00609           Yx [j2   + k*2*nk] = Bx [p] ; /* real */
00610           Yx [j2+1 + k*2*nk] = Bz [p] ; /* imag */
00611       }
00612         }
00613         break ;
00614
00615       }
00616       break ;
00617
00618   case CHOLMOD_COMPLEX:
00619
00620       switch (B->xtype)
00621       {
00622
00623     case CHOLMOD_REAL:
00624         /* Y complex, B real  */
00625         for (j = k1 ; j < k2 ; j++)
00626         {
00627       dj = d*j ;
00628       j2 = 2*(j-k1) ;
00629       for (k = 0 ; k < nrow ; k++)
00630       {
00631           p = P(k) + dj ;
00632           Yx [j2   + k*2*nk] = Bx [p] ; /* real */
00633           Yx [j2+1 + k*2*nk] = 0 ;    /* imag */
00634       }
00635         }
00636         break ;
00637
00638     case CHOLMOD_COMPLEX:
00639         /* Y complex, B complex  */
00640         for (j = k1 ; j < k2 ; j++)
00641         {
00642       dj = d*j ;
00643       j2 = 2*(j-k1) ;
00644       for (k = 0 ; k < nrow ; k++)
00645       {
00646           p = P(k) + dj ;
00647           Yx [j2   + k*2*nk] = Bx [2*p  ] ; /* real */
00648           Yx [j2+1 + k*2*nk] = Bx [2*p+1] ; /* imag */
00649       }
00650         }
00651         break ;
00652
00653     case CHOLMOD_ZOMPLEX:
00654         /* Y complex, B zomplex  */
00655         for (j = k1 ; j < k2 ; j++)
00656         {
00657       dj = d*j ;
00658       j2 = 2*(j-k1) ;
00659       for (k = 0 ; k < nrow ; k++)
00660       {
00661           p = P(k) + dj ;
00662           Yx [j2   + k*2*nk] = Bx [p] ; /* real */
00663           Yx [j2+1 + k*2*nk] = Bz [p] ; /* imag */
00664       }
00665         }
00666         break ;
00667
00668       }
00669       break ;
00670
00671   case CHOLMOD_ZOMPLEX:
00672
00673       switch (B->xtype)
00674       {
00675
00676     case CHOLMOD_REAL:
00677         /* Y zomplex, B real  */
00678         for (j = k1 ; j < k2 ; j++)
00679         {
00680       dj = d*j ;
00681       j2 = j-k1 ;
00682       for (k = 0 ; k < nrow ; k++)
00683       {
00684           p = P(k) + dj ;
00685           Yx [j2 + k*nk] = Bx [p] ;   /* real */
00686           Yz [j2 + k*nk] = 0 ;    /* imag */
00687       }
00688         }
00689         break ;
00690
00691     case CHOLMOD_COMPLEX:
00692         /* Y zomplex, B complex  */
00693         for (j = k1 ; j < k2 ; j++)
00694         {
00695       dj = d*j ;
00696       j2 = j-k1 ;
00697       for (k = 0 ; k < nrow ; k++)
00698       {
00699           p = P(k) + dj ;
00700           Yx [j2 + k*nk] = Bx [2*p  ] ; /* real */
00701           Yz [j2 + k*nk] = Bx [2*p+1] ; /* imag */
00702       }
00703         }
00704         break ;
00705
00706     case CHOLMOD_ZOMPLEX:
00707         /* Y zomplex, B zomplex */
00708         for (j = k1 ; j < k2 ; j++)
00709         {
00710       dj = d*j ;
00711       j2 = j-k1 ;
00712       for (k = 0 ; k < nrow ; k++)
00713       {
00714           p = P(k) + dj ;
00715           Yx [j2 + k*nk] = Bx [p] ;   /* real */
00716           Yz [j2 + k*nk] = Bz [p] ;   /* imag */
00717       }
00718         }
00719         break ;
00720
00721       }
00722       break ;
00723
00724     }
00725 }
00726
00727
00728 /* ========================================================================== */
00729 /* === iptrans ============================================================== */
00730 /* ========================================================================== */
00731
00732 /* X (P (1:nrow), k1 : min (k1+ncols,ncol)-1) = Y' where X is nrow-by-ncol.
00733  *
00734  * Copies into a permuted and transposed contiguous set of columns of X.
00735  * X is already allocated on input.  Y must be of sufficient size.  Let nk be
00736  * the number of columns accessed in X.  X->xtype determines the complexity of
00737  * the result.
00738  *
00739  * If X is real and Y is complex (or zomplex), only the real part of Y is
00740  * copied into X.  The imaginary part of Y is ignored.
00741  *
00742  * If X is complex (or zomplex) and Y is real, both the real and imaginary and
00743  * parts of X are returned in Y.  Y is 2*nk-by-nrow. The even
00744  * rows of Y contain the real part of X and the odd rows contain the
00745  * imaginary part of X.  Y->nzmax must be >= 2*nrow*nk.  Otherise, Y is
00746  * nk-by-nrow with leading dimension nk.  Y->nzmax must be >= nrow*nk.
00747  *
00748  * The case where Y is complex or zomplex, and X is real, is not used.
00749  *
00750  * The array transpose is performed, not the complex conjugate transpose.
00751  */
00752
00753 static void amesos_iptrans
00754 (
00755     /* ---- input ---- */
00756     cholmod_dense *Y, /* input matrix Y */
00757     Int *Perm,    /* optional input permutation (can be NULL) */
00758     Int k1,   /* first column of X to copy into */
00759     Int ncols,    /* last column to copy is min(k1+ncols,X->ncol)-1 */
00760     /* ---- in/out --- */
00761     cholmod_dense *X  /* output matrix X, already allocated */
00762 )
00763 {
00764     double *Yx, *Yz, *Xx, *Xz ;
00765     Int k2, nk, p, k, j, nrow, ncol, d, dj, j2 ;
00766
00767     /* ---------------------------------------------------------------------- */
00768     /* get inputs */
00769     /* ---------------------------------------------------------------------- */
00770
00771     ncol = X->ncol ;
00772     nrow = X->nrow ;
00773     k2 = MIN (k1+ncols, ncol) ;
00774     nk = MAX (k2 - k1, 0) ;
00775     d = X->d ;
00776     Xx = X->x ;
00777     Xz = X->z ;
00778     Yx = Y->x ;
00779     Yz = Y->z ;
00780     ASSERT (((Int) Y->nzmax) >= nrow*nk*
00781       ((X->xtype != CHOLMOD_REAL && Y->xtype == CHOLMOD_REAL) ? 2:1)) ;
00782
00783     /* ---------------------------------------------------------------------- */
00784     /* X (P (1:nrow), k1:k2-1) = Y' */
00785     /* ---------------------------------------------------------------------- */
00786
00787     switch (Y->xtype)
00788     {
00789
00790   case CHOLMOD_REAL:
00791
00792       switch (X->xtype)
00793       {
00794
00795     case CHOLMOD_REAL:
00796         /* Y real, X real  */
00797         for (j = k1 ; j < k2 ; j++)
00798         {
00799       dj = d*j ;
00800       j2 = j-k1 ;
00801       for (k = 0 ; k < nrow ; k++)
00802       {
00803           p = P(k) + dj ;
00804           Xx [p] = Yx [j2 + k*nk] ;   /* real */
00805       }
00806         }
00807         break ;
00808
00809     case CHOLMOD_COMPLEX:
00810         /* Y real, X complex. Y is 2*nk-by-nrow */
00811         for (j = k1 ; j < k2 ; j++)
00812         {
00813       dj = d*j ;
00814       j2 = 2*(j-k1) ;
00815       for (k = 0 ; k < nrow ; k++)
00816       {
00817           p = P(k) + dj ;
00818           Xx [2*p  ] = Yx [j2   + k*2*nk] ; /* real */
00819           Xx [2*p+1] = Yx [j2+1 + k*2*nk] ; /* imag */
00820       }
00821         }
00822         break ;
00823
00824     case CHOLMOD_ZOMPLEX:
00825         /* Y real, X zomplex. Y is 2*nk-by-nrow */
00826         for (j = k1 ; j < k2 ; j++)
00827         {
00828       dj = d*j ;
00829       j2 = 2*(j-k1) ;
00830       for (k = 0 ; k < nrow ; k++)
00831       {
00832           p = P(k) + dj ;
00833           Xx [p] = Yx [j2   + k*2*nk] ; /* real */
00834           Xz [p] = Yx [j2+1 + k*2*nk] ; /* imag */
00835       }
00836         }
00837         break ;
00838
00839       }
00840       break ;
00841
00842   case CHOLMOD_COMPLEX:
00843
00844       switch (X->xtype)
00845       {
00846
00847 #if 0
00848     case CHOLMOD_REAL:
00849         /* this case is not used */
00850         break ;
00851 #endif
00852
00853     case CHOLMOD_COMPLEX:
00854         /* Y complex, X complex  */
00855         for (j = k1 ; j < k2 ; j++)
00856         {
00857       dj = d*j ;
00858       j2 = 2*(j-k1) ;
00859       for (k = 0 ; k < nrow ; k++)
00860       {
00861           p = P(k) + dj ;
00862           Xx [2*p  ] = Yx [j2   + k*2*nk] ; /* real */
00863           Xx [2*p+1] = Yx [j2+1 + k*2*nk] ; /* imag */
00864       }
00865         }
00866         break ;
00867
00868     case CHOLMOD_ZOMPLEX:
00869         /* Y complex, X zomplex  */
00870         for (j = k1 ; j < k2 ; j++)
00871         {
00872       dj = d*j ;
00873       j2 = 2*(j-k1) ;
00874       for (k = 0 ; k < nrow ; k++)
00875       {
00876           p = P(k) + dj ;
00877           Xx [p] = Yx [j2   + k*2*nk] ; /* real */
00878           Xz [p] = Yx [j2+1 + k*2*nk] ; /* imag */
00879       }
00880         }
00881         break ;
00882
00883       }
00884       break ;
00885
00886   case CHOLMOD_ZOMPLEX:
00887
00888       switch (X->xtype)
00889       {
00890
00891 #if 0
00892     case CHOLMOD_REAL:
00893         /* this case is not used */
00894         break ;
00895 #endif
00896
00897     case CHOLMOD_COMPLEX:
00898         /* Y zomplex, X complex  */
00899         for (j = k1 ; j < k2 ; j++)
00900         {
00901       dj = d*j ;
00902       j2 = j-k1 ;
00903       for (k = 0 ; k < nrow ; k++)
00904       {
00905           p = P(k) + dj ;
00906           Xx [2*p  ] = Yx [j2 + k*nk] ; /* real */
00907           Xx [2*p+1] = Yz [j2 + k*nk] ; /* imag */
00908       }
00909         }
00910         break ;
00911
00912     case CHOLMOD_ZOMPLEX:
00913         /* Y zomplex, X zomplex */
00914         for (j = k1 ; j < k2 ; j++)
00915         {
00916       dj = d*j ;
00917       j2 = j-k1 ;
00918       for (k = 0 ; k < nrow ; k++)
00919       {
00920           p = P(k) + dj ;
00921           Xx [p] = Yx [j2 + k*nk] ;   /* real */
00922           Xz [p] = Yz [j2 + k*nk] ;   /* imag */
00923       }
00924         }
00925         break ;
00926
00927       }
00928       break ;
00929
00930     }
00931 }
00932
00933
00934 /* ========================================================================== */
00935 /* === cholmod_solve ======================================================== */
00936 /* ========================================================================== */
00937
00938 /* Solve a linear system.
00939  *
00940  * The factorization can be simplicial LDL', simplicial LL', or supernodal LL'.
00941  * The Dx=b solve returns silently for the LL' factorizations (it is implicitly
00942  * identity).
00943  */
00944
00945 cholmod_dense *CHOLMOD(solve)
00946 (
00947     /* ---- input ---- */
00948     int sys,    /* system to solve */
00949     cholmod_factor *L,  /* factorization to use */
00950     cholmod_dense *B, /* right-hand-side */
00951     /* --------------- */
00952     cholmod_common *Common
00953 )
00954 {
00955     cholmod_dense *Y = NULL, *X = NULL ;
00956     Int *Perm ;
00957     Int n, nrhs, ncols, ctype, xtype, k1, nr, ytype ;
00958
00959     /* ---------------------------------------------------------------------- */
00960     /* check inputs */
00961     /* ---------------------------------------------------------------------- */
00962
00963     RETURN_IF_NULL_COMMON (NULL) ;
00964     RETURN_IF_NULL (L, NULL) ;
00965     RETURN_IF_NULL (B, NULL) ;
00966     RETURN_IF_XTYPE_INVALID (L, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, NULL) ;
00967     RETURN_IF_XTYPE_INVALID (B, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, NULL) ;
00968     if (sys < CHOLMOD_A || sys > CHOLMOD_Pt)
00969     {
00970   ERROR (CHOLMOD_INVALID, "invalid system") ;
00971   return (NULL) ;
00972     }
00973     if (B->d < L->n || B->nrow != L->n)
00974     {
00975   ERROR (CHOLMOD_INVALID, "dimensions of L and B do not match") ;
00976   return (NULL) ;
00977     }
00978     DEBUG (CHOLMOD(dump_factor) (L, "L", Common)) ;
00979     DEBUG (CHOLMOD(dump_dense) (B, "B", Common)) ;
00980     Common->status = CHOLMOD_OK ;
00981
00982     /* ---------------------------------------------------------------------- */
00983     /* get inputs */
00984     /* ---------------------------------------------------------------------- */
00985
00986     if ((sys == CHOLMOD_P || sys == CHOLMOD_Pt || sys == CHOLMOD_A)
00987       && L->ordering != CHOLMOD_NATURAL)
00988     {
00989   Perm = L->Perm ;
00990     }
00991     else
00992     {
00993   /* do not use L->Perm; use the identity permutation instead */
00994   Perm = NULL ;
00995     }
00996
00997     nrhs = B->ncol ;
00998     n = L->n ;
00999
01000     /* ---------------------------------------------------------------------- */
01001     /* allocate the result X */
01002     /* ---------------------------------------------------------------------- */
01003
01004     ctype = (Common->prefer_zomplex) ? CHOLMOD_ZOMPLEX : CHOLMOD_COMPLEX ;
01005
01006     if (sys == CHOLMOD_P || sys == CHOLMOD_Pt)
01007     {
01008   /* x=Pb and x=P'b return X real if B is real; X is the preferred
01009    * complex/zcomplex type if B is complex or zomplex */
01010   xtype = (B->xtype == CHOLMOD_REAL) ? CHOLMOD_REAL : ctype ;
01011     }
01012     else if (L->xtype == CHOLMOD_REAL && B->xtype == CHOLMOD_REAL)
01013     {
01014   /* X is real if both L and B are real */
01015   xtype = CHOLMOD_REAL ;
01016     }
01017     else
01018     {
01019   /* X is complex, use the preferred complex/zomplex type */
01020   xtype = ctype ;
01021     }
01022
01023     X = CHOLMOD(allocate_dense) (n, nrhs, n, xtype, Common) ;
01024     if (Common->status < CHOLMOD_OK)
01025     {
01026   return (NULL) ;
01027     }
01028
01029     /* ---------------------------------------------------------------------- */
01030     /* solve using L, D, L', P, or some combination */
01031     /* ---------------------------------------------------------------------- */
01032
01033     if (sys == CHOLMOD_P)
01034     {
01035
01036   /* ------------------------------------------------------------------ */
01037   /* x = P*b */
01038   /* ------------------------------------------------------------------ */
01039
01040   amesos_perm (B, Perm, 0, nrhs, X) ;
01041
01042     }
01043     else if (sys == CHOLMOD_Pt)
01044     {
01045
01046   /* ------------------------------------------------------------------ */
01047   /* x = P'*b */
01048   /* ------------------------------------------------------------------ */
01049
01050   amesos_iperm (B, Perm, 0, nrhs, X) ;
01051
01052     }
01053     else if (L->is_super)
01054     {
01055
01056   /* ------------------------------------------------------------------ */
01057   /* solve using a supernodal LL' factorization */
01058   /* ------------------------------------------------------------------ */
01059
01060 #ifndef NSUPERNODAL
01061   Int blas_ok = TRUE ;
01062
01063   /* allocate workspace */
01064   cholmod_dense *E ;
01065   Int dual ;
01066   dual = (L->xtype == CHOLMOD_REAL && B->xtype != CHOLMOD_REAL) ? 2 : 1 ;
01067   Y = CHOLMOD(allocate_dense) (n, dual*nrhs, n, L->xtype, Common) ;
01068   E = CHOLMOD(allocate_dense) (dual*nrhs, L->maxesize, dual*nrhs,
01069     L->xtype, Common) ;
01070
01071   if (Common->status < CHOLMOD_OK)
01072   {
01073       /* out of memory */
01074       CHOLMOD(free_dense) (&X, Common) ;
01075       CHOLMOD(free_dense) (&Y, Common) ;
01076       CHOLMOD(free_dense) (&E, Common) ;
01077       return (NULL) ;
01078   }
01079
01080   amesos_perm (B, Perm, 0, nrhs, Y) ;         /* Y = P*B */
01081
01082   if (sys == CHOLMOD_A || sys == CHOLMOD_LDLt)
01083   {
01084       blas_ok = CHOLMOD(super_lsolve) (L, Y, E, Common) ;    /* Y = L\Y */
01085       blas_ok = blas_ok &&
01086     CHOLMOD(super_ltsolve) (L, Y, E, Common) ;     /* Y = L'\Y*/
01087   }
01088   else if (sys == CHOLMOD_L || sys == CHOLMOD_LD)
01089   {
01090       blas_ok = CHOLMOD(super_lsolve) (L, Y, E, Common) ;    /* Y = L\Y */
01091   }
01092   else if (sys == CHOLMOD_Lt || sys == CHOLMOD_DLt)
01093   {
01094       blas_ok = CHOLMOD(super_ltsolve) (L, Y, E, Common) ;   /* Y = L'\Y*/
01095   }
01096   CHOLMOD(free_dense) (&E, Common) ;
01097
01098   amesos_iperm (Y, Perm, 0, nrhs, X) ;          /* X = P'*Y */
01099
01100   if (CHECK_BLAS_INT && !blas_ok)
01101   {
01102       /* Integer overflow in the BLAS.  This is probably impossible,
01103        * since the BLAS were used to create the supernodal factorization.
01104        * It might be possible for the calls to the BLAS to differ between
01105        * factorization and forward/backsolves, however.  This statement
01106        * is untested; it does not appear in the compiled code if
01107        * CHECK_BLAS_INT is true (when the same integer is used in CHOLMOD
01108        * and the BLAS. */
01109       CHOLMOD(free_dense) (&X, Common) ;
01110   }
01111
01112 #else
01113   /* CHOLMOD Supernodal module not installed */
01114   ERROR (CHOLMOD_NOT_INSTALLED,"Supernodal module not installed") ;
01115 #endif
01116
01117     }
01118     else
01119     {
01120
01121   /* ------------------------------------------------------------------ */
01122   /* solve using a simplicial LL' or LDL' factorization */
01123   /* ------------------------------------------------------------------ */
01124
01125   if (L->xtype == CHOLMOD_REAL && B->xtype == CHOLMOD_REAL)
01126   {
01127       /* L, B, and Y are all real */
01128       /* solve with up to 4 columns of B at a time */
01129       ncols = 4 ;
01130       nr = MAX (4, nrhs) ;
01131       ytype = CHOLMOD_REAL ;
01132   }
01133   else if (L->xtype == CHOLMOD_REAL)
01134   {
01135       /* solve with one column of B (real/imag), at a time */
01136       ncols = 1 ;
01137       nr = 2 ;
01138       ytype = CHOLMOD_REAL ;
01139   }
01140   else
01141   {
01142       /* L is complex or zomplex, B is real/complex/zomplex, Y has the
01143        * same complexity as L.  Solve with one column of B at a time. */
01144       ncols = 1 ;
01145       nr = 1 ;
01146       ytype = L->xtype ;
01147   }
01148
01149   Y = CHOLMOD(allocate_dense) (nr, n, nr, ytype, Common) ;
01150
01151   if (Common->status < CHOLMOD_OK)
01152   {
01153       /* out of memory */
01154       CHOLMOD(free_dense) (&X, Common) ;
01155       CHOLMOD(free_dense) (&Y, Common) ;
01156       return (NULL) ;
01157   }
01158
01159   for (k1 = 0 ; k1 < nrhs ; k1 += ncols)
01160   {
01161       /* -------------------------------------------------------------- */
01162       /* Y = B (P, k1:k1+ncols-1)' = (P * B (:,...))' */
01163       /* -------------------------------------------------------------- */
01164
01165       amesos_ptrans (B, Perm, k1, ncols, Y) ;
01166
01167       /* -------------------------------------------------------------- */
01168       /* solve Y = (L' \ (L \ Y'))', or other system, with template */
01169       /* -------------------------------------------------------------- */
01170
01171       switch (L->xtype)
01172       {
01173     case CHOLMOD_REAL:
01174         r_simplicial_solver (sys, L, Y) ;
01175         break ;
01176
01177     case CHOLMOD_COMPLEX:
01178         c_simplicial_solver (sys, L, Y) ;
01179         break ;
01180
01181     case CHOLMOD_ZOMPLEX:
01182         z_simplicial_solver (sys, L, Y) ;
01183         break ;
01184       }
01185
01186       /* -------------------------------------------------------------- */
01187       /* X (P, k1:k2+ncols-1) = Y' */
01188       /* -------------------------------------------------------------- */
01189
01190       amesos_iptrans (Y, Perm, k1, ncols, X) ;
01191   }
01192     }
01193
01194     CHOLMOD(free_dense) (&Y, Common) ;
01195     DEBUG (CHOLMOD(dump_dense) (X, "X result", Common)) ;
01196     return (X) ;
01197 }
01198 #endif
```