Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_spsolve.c
Go to the documentation of this file.
```00001 /* ========================================================================== */
00002 /* === Cholesky/cholmod_spsolve ============================================= */
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 /* Given an LL' or LDL' factorization of A, 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  * where b and x are sparse.  If L and b are real, then x is real.  Otherwise,
00026  * x is complex or zomplex, depending on the Common->prefer_zomplex parameter.
00027  * All xtypes of x and b are supported (real, complex, and zomplex).
00028  */
00029
00030 #ifndef NCHOLESKY
00031
00032 #include "amesos_cholmod_internal.h"
00033 #include "amesos_cholmod_cholesky.h"
00034
00035 /* ========================================================================== */
00036 /* === EXPAND_AS_NEEDED ===================================================== */
00037 /* ========================================================================== */
00038
00039 /* Double the size of the sparse matrix X, if we have run out of space. */
00040
00041 #define EXPAND_AS_NEEDED \
00042 if (xnz >= nzmax) \
00043 { \
00044     nzmax *= 2 ; \
00045     CHOLMOD(reallocate_sparse) (nzmax, X, Common) ; \
00046     if (Common->status < CHOLMOD_OK) \
00047     { \
00048   CHOLMOD(free_sparse) (&X, Common) ; \
00049   CHOLMOD(free_dense) (&X4, Common) ; \
00050   CHOLMOD(free_dense) (&B4, Common) ; \
00051   return (NULL) ; \
00052     } \
00053     Xi = X->i ; \
00054     Xx = X->x ; \
00055     Xz = X->z ; \
00056 }
00057
00058
00059 /* ========================================================================== */
00060 /* === cholmod_spolve ======================================================= */
00061 /* ========================================================================== */
00062
00063 cholmod_sparse *CHOLMOD(spsolve)      /* returns the sparse solution X */
00064 (
00065     /* ---- input ---- */
00066     int sys,    /* system to solve */
00067     cholmod_factor *L,  /* factorization to use */
00068     cholmod_sparse *B,  /* right-hand-side */
00069     /* --------------- */
00070     cholmod_common *Common
00071 )
00072 {
00073     double x, z ;
00074     cholmod_dense *X4, *B4 ;
00075     cholmod_sparse *X ;
00076     double *Bx, *Bz, *Xx, *Xz, *B4x, *B4z, *X4x, *X4z ;
00077     Int *Bi, *Bp, *Xp, *Xi, *Bnz ;
00078     Int n, nrhs, q, p, i, j, j1, j2, packed, block, pend, jn, xtype ;
00079     size_t xnz, nzmax ;
00080
00081     /* ---------------------------------------------------------------------- */
00082     /* check inputs */
00083     /* ---------------------------------------------------------------------- */
00084
00085     RETURN_IF_NULL_COMMON (NULL) ;
00086     RETURN_IF_NULL (L, NULL) ;
00087     RETURN_IF_NULL (B, NULL) ;
00088     RETURN_IF_XTYPE_INVALID (L, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, NULL) ;
00089     RETURN_IF_XTYPE_INVALID (B, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, NULL) ;
00090     if (L->n != B->nrow)
00091     {
00092   ERROR (CHOLMOD_INVALID, "dimensions of L and B do not match") ;
00093   return (NULL) ;
00094     }
00095     if (B->stype)
00096     {
00097   ERROR (CHOLMOD_INVALID, "B cannot be stored in symmetric mode") ;
00098   return (NULL) ;
00099     }
00100     Common->status = CHOLMOD_OK ;
00101
00102     /* ---------------------------------------------------------------------- */
00103     /* allocate workspace B4 and initial result X */
00104     /* ---------------------------------------------------------------------- */
00105
00106     n = L->n ;
00107     nrhs = B->ncol ;
00108
00109     /* X is real if both L and B are real, complex/zomplex otherwise */
00110     xtype = (L->xtype == CHOLMOD_REAL && B->xtype == CHOLMOD_REAL) ?
00111   CHOLMOD_REAL :
00112   (Common->prefer_zomplex ? CHOLMOD_ZOMPLEX : CHOLMOD_COMPLEX) ;
00113
00114     /* solve up to 4 columns at a time */
00115     block = MIN (nrhs, 4) ;
00116
00117     /* initial size of X is at most 4*n */
00118     nzmax = n*block ;
00119
00120     X = CHOLMOD(spzeros) (n, nrhs, nzmax, xtype, Common) ;
00121     B4 = CHOLMOD(zeros) (n, block, B->xtype, Common) ;
00122     if (Common->status < CHOLMOD_OK)
00123     {
00124   CHOLMOD(free_sparse) (&X, Common) ;
00125   CHOLMOD(free_dense) (&B4, Common) ;
00126   return (NULL) ;
00127     }
00128
00129     Bp = B->p ;
00130     Bi = B->i ;
00131     Bx = B->x ;
00132     Bz = B->z ;
00133     Bnz = B->nz ;
00134     packed = B->packed ;
00135
00136     Xp = X->p ;
00137     Xi = X->i ;
00138     Xx = X->x ;
00139     Xz = X->z ;
00140
00141     xnz = 0 ;
00142
00143     B4x = B4->x ;
00144     B4z = B4->z ;
00145
00146     /* ---------------------------------------------------------------------- */
00147     /* solve in chunks of 4 columns at a time */
00148     /* ---------------------------------------------------------------------- */
00149
00150     for (j1 = 0 ; j1 < nrhs ; j1 += block)
00151     {
00152
00153   /* ------------------------------------------------------------------ */
00154   /* adjust the number of columns of B4 */
00155   /* ------------------------------------------------------------------ */
00156
00157   j2 = MIN (nrhs, j1 + block) ;
00158   B4->ncol = j2 - j1 ;
00159
00160   /* ------------------------------------------------------------------ */
00161   /* scatter B(j1:j2-1) into B4 */
00162   /* ------------------------------------------------------------------ */
00163
00164   for (j = j1 ; j < j2 ; j++)
00165   {
00166       p = Bp [j] ;
00167       pend = (packed) ? (Bp [j+1]) : (p + Bnz [j]) ;
00168       jn = (j-j1)*n ;
00169
00170       switch (B->xtype)
00171       {
00172
00173     case CHOLMOD_REAL:
00174         for ( ; p < pend ; p++)
00175         {
00176       B4x [Bi [p] + jn] = Bx [p] ;
00177         }
00178         break ;
00179
00180     case CHOLMOD_COMPLEX:
00181         for ( ; p < pend ; p++)
00182         {
00183       q = Bi [p] + jn ;
00184       B4x [2*q  ] = Bx [2*p  ] ;
00185       B4x [2*q+1] = Bx [2*p+1] ;
00186         }
00187         break ;
00188
00189     case CHOLMOD_ZOMPLEX:
00190         for ( ; p < pend ; p++)
00191         {
00192       q = Bi [p] + jn ;
00193       B4x [q] = Bx [p] ;
00194       B4z [q] = Bz [p] ;
00195         }
00196         break ;
00197       }
00198   }
00199
00200   /* ------------------------------------------------------------------ */
00201   /* solve the system (X4 = A\B4 or other system) */
00202   /* ------------------------------------------------------------------ */
00203
00204   X4 = CHOLMOD(solve) (sys, L, B4, Common) ;
00205   if (Common->status < CHOLMOD_OK)
00206   {
00207       CHOLMOD(free_sparse) (&X, Common) ;
00208       CHOLMOD(free_dense) (&B4, Common) ;
00209       CHOLMOD(free_dense) (&X4, Common) ;
00210       return (NULL) ;
00211   }
00212   ASSERT (X4->xtype == xtype) ;
00213   X4x = X4->x ;
00214   X4z = X4->z ;
00215
00216   /* ------------------------------------------------------------------ */
00217   /* append the solution onto X */
00218   /* ------------------------------------------------------------------ */
00219
00220   for (j = j1 ; j < j2 ; j++)
00221   {
00222       Xp [j] = xnz ;
00223       jn = (j-j1)*n ;
00224       if ( xnz + n <= nzmax)
00225       {
00226
00227     /* ---------------------------------------------------------- */
00228     /* X is guaranteed to be large enough */
00229     /* ---------------------------------------------------------- */
00230
00231     switch (xtype)
00232     {
00233
00234         case CHOLMOD_REAL:
00235       for (i = 0 ; i < n ; i++)
00236       {
00237           x = X4x [i + jn] ;
00238           if (IS_NONZERO (x))
00239           {
00240         Xi [xnz] = i ;
00241         Xx [xnz] = x ;
00242         xnz++ ;
00243           }
00244       }
00245       break ;
00246
00247         case CHOLMOD_COMPLEX:
00248       for (i = 0 ; i < n ; i++)
00249       {
00250           x = X4x [2*(i + jn)  ] ;
00251           z = X4x [2*(i + jn)+1] ;
00252           if (IS_NONZERO (x) || IS_NONZERO (z))
00253           {
00254         Xi [xnz] = i ;
00255         Xx [2*xnz  ] = x ;
00256         Xx [2*xnz+1] = z ;
00257         xnz++ ;
00258           }
00259       }
00260       break ;
00261
00262         case CHOLMOD_ZOMPLEX:
00263       for (i = 0 ; i < n ; i++)
00264       {
00265           x = X4x [i + jn] ;
00266           z = X4z [i + jn] ;
00267           if (IS_NONZERO (x) || IS_NONZERO (z))
00268           {
00269         Xi [xnz] = i ;
00270         Xx [xnz] = x ;
00271         Xz [xnz] = z ;
00272         xnz++ ;
00273           }
00274       }
00275       break ;
00276     }
00277
00278       }
00279       else
00280       {
00281
00282     /* ---------------------------------------------------------- */
00283     /* X may need to increase in size */
00284     /* ---------------------------------------------------------- */
00285
00286     switch (xtype)
00287     {
00288
00289         case CHOLMOD_REAL:
00290       for (i = 0 ; i < n ; i++)
00291       {
00292           x = X4x [i + jn] ;
00293           if (IS_NONZERO (x))
00294           {
00295         EXPAND_AS_NEEDED ;
00296         Xi [xnz] = i ;
00297         Xx [xnz] = x ;
00298         xnz++ ;
00299           }
00300       }
00301       break ;
00302
00303         case CHOLMOD_COMPLEX:
00304       for (i = 0 ; i < n ; i++)
00305       {
00306           x = X4x [2*(i + jn)  ] ;
00307           z = X4x [2*(i + jn)+1] ;
00308           if (IS_NONZERO (x) || IS_NONZERO (z))
00309           {
00310         EXPAND_AS_NEEDED ;
00311         Xi [xnz] = i ;
00312         Xx [2*xnz  ] = x ;
00313         Xx [2*xnz+1] = z ;
00314         xnz++ ;
00315           }
00316       }
00317       break ;
00318
00319         case CHOLMOD_ZOMPLEX:
00320       for (i = 0 ; i < n ; i++)
00321       {
00322           x = X4x [i + jn] ;
00323           z = X4z [i + jn] ;
00324           if (IS_NONZERO (x) || IS_NONZERO (z))
00325           {
00326         EXPAND_AS_NEEDED ;
00327         Xi [xnz] = i ;
00328         Xx [xnz] = x ;
00329         Xz [xnz] = z ;
00330         xnz++ ;
00331           }
00332       }
00333       break ;
00334     }
00335
00336       }
00337   }
00338   CHOLMOD(free_dense) (&X4, Common) ;
00339
00340   /* ------------------------------------------------------------------ */
00341   /* clear B4 for next iteration */
00342   /* ------------------------------------------------------------------ */
00343
00344   if (j2 < nrhs)
00345   {
00346
00347       for (j = j1 ; j < j2 ; j++)
00348       {
00349     p = Bp [j] ;
00350     pend = (packed) ? (Bp [j+1]) : (p + Bnz [j]) ;
00351     jn = (j-j1)*n ;
00352
00353     switch (B->xtype)
00354     {
00355
00356         case CHOLMOD_REAL:
00357       for ( ; p < pend ; p++)
00358       {
00359           B4x [Bi [p] + jn] = 0 ;
00360       }
00361       break ;
00362
00363         case CHOLMOD_COMPLEX:
00364       for ( ; p < pend ; p++)
00365       {
00366           q = Bi [p] + jn ;
00367           B4x [2*q  ] = 0 ;
00368           B4x [2*q+1] = 0 ;
00369       }
00370       break ;
00371
00372         case CHOLMOD_ZOMPLEX:
00373       for ( ; p < pend ; p++)
00374       {
00375           q = Bi [p] + jn ;
00376           B4x [q] = 0 ;
00377           B4z [q] = 0 ;
00378       }
00379       break ;
00380     }
00381       }
00382   }
00383     }
00384
00385     Xp [nrhs] = xnz ;
00386
00387     /* ---------------------------------------------------------------------- */
00388     /* reduce X in size, free workspace, and return result */
00389     /* ---------------------------------------------------------------------- */
00390
00391     ASSERT (xnz <= X->nzmax) ;
00392     CHOLMOD(reallocate_sparse) (xnz, X, Common) ;
00393     ASSERT (Common->status == CHOLMOD_OK) ;
00394     CHOLMOD(free_dense) (&B4, Common) ;
00395     return (X) ;
00396 }
00397 #endif
```