Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_sparse.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === Core/cholmod_sparse ================================================== */
00003 /* ========================================================================== */
00004 
00005 /* -----------------------------------------------------------------------------
00006  * CHOLMOD/Core Module.  Copyright (C) 2005-2006,
00007  * Univ. of Florida.  Author: Timothy A. Davis
00008  * The CHOLMOD/Core Module is licensed under Version 2.1 of the GNU
00009  * Lesser General Public License.  See lesser.txt for a text of the license.
00010  * CHOLMOD is also available under other licenses; contact authors for details.
00011  * http://www.cise.ufl.edu/research/sparse
00012  * -------------------------------------------------------------------------- */
00013 
00014 /* Core utility routines for the cholmod_sparse object:
00015  *
00016  * A sparse matrix is held in compressed column form.  In the basic type
00017  * ("packed", which corresponds to a MATLAB sparse matrix), an n-by-n matrix
00018  * with nz entries is held in three arrays: p of size n+1, i of size nz, and x
00019  * of size nz.  Row indices of column j are held in i [p [j] ... p [j+1]-1] and
00020  * in the same locations in x.  There may be no duplicate entries in a column.
00021  * Row indices in each column may be sorted or unsorted (CHOLMOD keeps track).
00022  *
00023  * Primary routines:
00024  * -----------------
00025  * cholmod_allocate_sparse  allocate a sparse matrix
00026  * cholmod_free_sparse    free a sparse matrix
00027  *
00028  * Secondary routines:
00029  * -------------------
00030  * cholmod_reallocate_sparse  change the size (# entries) of sparse matrix
00031  * cholmod_nnz      number of nonzeros in a sparse matrix
00032  * cholmod_speye    sparse identity matrix
00033  * cholmod_spzeros    sparse zero matrix
00034  * cholmod_copy_sparse    create a copy of a sparse matrix
00035  *
00036  * All xtypes are supported (pattern, real, complex, and zomplex)
00037  */
00038 
00039 #include "amesos_cholmod_internal.h"
00040 #include "amesos_cholmod_core.h"
00041 
00042 
00043 /* ========================================================================== */
00044 /* === cholmod_allocate_sparse ============================================== */
00045 /* ========================================================================== */
00046 
00047 /* Allocate space for a matrix.  A->i and A->x are not initialized.  A->p
00048  * (and A->nz if A is not packed) are set to zero, so a matrix containing no
00049  * entries (all zero) is returned.  See also cholmod_spzeros.
00050  *
00051  * workspace: none
00052  */
00053 
00054 cholmod_sparse *CHOLMOD(allocate_sparse)
00055 (
00056     /* ---- input ---- */
00057     size_t nrow,  /* # of rows of A */
00058     size_t ncol,  /* # of columns of A */
00059     size_t nzmax, /* max # of nonzeros of A */
00060     int sorted,   /* TRUE if columns of A sorted, FALSE otherwise */
00061     int packed,   /* TRUE if A will be packed, FALSE otherwise */
00062     int stype,    /* stype of A */
00063     int xtype,    /* CHOLMOD_PATTERN, _REAL, _COMPLEX, or _ZOMPLEX */
00064     /* --------------- */
00065     cholmod_common *Common
00066 )
00067 {
00068     cholmod_sparse *A ;
00069     Int *Ap, *Anz ;
00070     size_t nzmax0 ;
00071     Int j ;
00072     int ok = TRUE ;
00073 
00074     /* ---------------------------------------------------------------------- */
00075     /* get inputs */
00076     /* ---------------------------------------------------------------------- */
00077 
00078     RETURN_IF_NULL_COMMON (NULL) ;
00079     if (stype != 0 && nrow != ncol)
00080     {
00081   ERROR (CHOLMOD_INVALID, "rectangular matrix with stype != 0 invalid") ;
00082   return (NULL) ;
00083     }
00084     if (xtype < CHOLMOD_PATTERN || xtype > CHOLMOD_ZOMPLEX)
00085     {
00086   ERROR (CHOLMOD_INVALID, "xtype invalid") ;
00087   return (NULL) ;
00088     }
00089     /* ensure the dimensions do not cause integer overflow */
00090     (void) CHOLMOD(add_size_t) (ncol, 2, &ok) ;
00091     if (!ok || nrow > Int_max || ncol > Int_max || nzmax > Int_max)
00092     {
00093   ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
00094   return (NULL) ;
00095     }
00096     Common->status = CHOLMOD_OK ;
00097 
00098     /* ---------------------------------------------------------------------- */
00099     /* allocate header */
00100     /* ---------------------------------------------------------------------- */
00101 
00102     A = CHOLMOD(malloc) (sizeof (cholmod_sparse), 1, Common) ;
00103     if (Common->status < CHOLMOD_OK)
00104     {
00105   return (NULL) ;     /* out of memory */
00106     }
00107     PRINT1 (("cholmod_allocate_sparse %d-by-%d nzmax %d sorted %d packed %d"
00108     " xtype %d\n", nrow, ncol, nzmax, sorted, packed, xtype)) ;
00109 
00110     nzmax = MAX (1, nzmax) ;
00111 
00112     A->nrow = nrow ;
00113     A->ncol = ncol ;
00114     A->nzmax = nzmax ;
00115     A->packed = packed ;    /* default is packed (A->nz not present) */
00116     A->stype = stype ;
00117     A->itype = ITYPE ;
00118     A->xtype = xtype ;
00119     A->dtype = DTYPE ;
00120 
00121     A->nz = NULL ;
00122     A->p = NULL ;
00123     A->i = NULL ;
00124     A->x = NULL ;
00125     A->z = NULL ;
00126 
00127     /* A 1-by-m matrix always has sorted columns */
00128     A->sorted = (nrow <= 1) ? TRUE : sorted ;
00129 
00130     /* ---------------------------------------------------------------------- */
00131     /* allocate the matrix itself */
00132     /* ---------------------------------------------------------------------- */
00133 
00134     /* allocate O(ncol) space */
00135     A->p = CHOLMOD(malloc) (((size_t) ncol)+1, sizeof (Int), Common) ;
00136     if (!packed)
00137     {
00138   A->nz = CHOLMOD(malloc) (ncol, sizeof (Int), Common) ;
00139     }
00140 
00141     /* allocate O(nz) space */
00142     nzmax0 = 0 ;
00143     CHOLMOD(realloc_multiple) (nzmax, 1, xtype, &(A->i), NULL, &(A->x), &(A->z),
00144       &nzmax0, Common) ;
00145 
00146     if (Common->status < CHOLMOD_OK)
00147     {
00148   CHOLMOD(free_sparse) (&A, Common) ;
00149   return (NULL) ;     /* out of memory */
00150     }
00151 
00152     /* ---------------------------------------------------------------------- */
00153     /* initialize A->p and A->nz so that A is an empty matrix */
00154     /* ---------------------------------------------------------------------- */
00155 
00156     Ap = A->p ;
00157     for (j = 0 ; j <= (Int) ncol ; j++)
00158     {
00159   Ap [j] = 0 ;
00160     }
00161     if (!packed)
00162     {
00163   Anz = A->nz ;
00164   for (j = 0 ; j < (Int) ncol ; j++)
00165   {
00166       Anz [j] = 0 ;
00167   }
00168     }
00169     return (A) ;
00170 }
00171 
00172 
00173 /* ========================================================================== */
00174 /* === cholmod_free_sparse ================================================== */
00175 /* ========================================================================== */
00176 
00177 /* free a sparse matrix
00178  *
00179  * workspace: none
00180  */
00181 
00182 int CHOLMOD(free_sparse)
00183 (
00184     /* ---- in/out --- */
00185     cholmod_sparse **AHandle, /* matrix to deallocate, NULL on output */
00186     /* --------------- */
00187     cholmod_common *Common
00188 )
00189 {
00190     Int n, nz ;
00191     cholmod_sparse *A ;
00192 
00193     RETURN_IF_NULL_COMMON (FALSE) ;
00194 
00195     if (AHandle == NULL)
00196     {
00197   /* nothing to do */
00198   return (TRUE) ;
00199     }
00200     A = *AHandle ;
00201     if (A == NULL)
00202     {
00203   /* nothing to do */
00204   return (TRUE) ;
00205     }
00206     n = A->ncol ;
00207     nz = A->nzmax ;
00208     A->p  = CHOLMOD(free) (n+1, sizeof (Int), A->p,  Common) ;
00209     A->i  = CHOLMOD(free) (nz,  sizeof (Int), A->i,  Common) ;
00210     A->nz = CHOLMOD(free) (n,   sizeof (Int), A->nz, Common) ;
00211 
00212     switch (A->xtype)
00213     {
00214   case CHOLMOD_REAL:
00215       A->x = CHOLMOD(free) (nz, sizeof (double), A->x,  Common) ;
00216       break ;
00217 
00218   case CHOLMOD_COMPLEX:
00219       A->x = CHOLMOD(free) (nz, 2*sizeof (double), A->x,  Common) ;
00220       break ;
00221 
00222   case CHOLMOD_ZOMPLEX:
00223       A->x = CHOLMOD(free) (nz, sizeof (double), A->x,  Common) ;
00224       A->z = CHOLMOD(free) (nz, sizeof (double), A->z,  Common) ;
00225       break ;
00226     }
00227 
00228     *AHandle = CHOLMOD(free) (1, sizeof (cholmod_sparse), (*AHandle), Common) ;
00229     return (TRUE) ;
00230 }
00231 
00232 
00233 /* ========================================================================== */
00234 /* === cholmod_reallocate_sparse ============================================ */
00235 /* ========================================================================== */
00236 
00237 /* Change the size of A->i, A->x, and A->z, or allocate them if their current
00238  * size is zero.  A->x and A->z are not modified if A->xtype is CHOLMOD_PATTERN.
00239  * A->z is not modified unless A->xtype is CHOLMOD_ZOMPLEX.
00240  * 
00241  * workspace: none
00242  */
00243 
00244 int CHOLMOD(reallocate_sparse)
00245 (
00246     /* ---- input ---- */
00247     size_t nznew, /* new # of entries in A */
00248     /* ---- in/out --- */
00249     cholmod_sparse *A,  /* matrix to reallocate */
00250     /* --------------- */
00251     cholmod_common *Common
00252 )
00253 {
00254 
00255     /* ---------------------------------------------------------------------- */
00256     /* get inputs */
00257     /* ---------------------------------------------------------------------- */
00258 
00259     RETURN_IF_NULL_COMMON (FALSE) ;
00260     RETURN_IF_NULL (A, FALSE) ;
00261     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00262     Common->status = CHOLMOD_OK ;
00263     PRINT1 (("realloc matrix %d to %d, xtype: %d\n",
00264     A->nzmax, nznew, A->xtype)) ;
00265 
00266     /* ---------------------------------------------------------------------- */
00267     /* resize the matrix */
00268     /* ---------------------------------------------------------------------- */
00269 
00270     CHOLMOD(realloc_multiple) (MAX (1,nznew), 1, A->xtype, &(A->i), NULL,
00271       &(A->x), &(A->z), &(A->nzmax), Common) ;
00272 
00273     return (Common->status == CHOLMOD_OK) ;
00274 }
00275 
00276 
00277 /* ========================================================================== */
00278 /* === cholmod_speye ======================================================== */
00279 /* ========================================================================== */
00280 
00281 /* Return a sparse identity matrix. */
00282 
00283 cholmod_sparse *CHOLMOD(speye)
00284 (
00285     /* ---- input ---- */
00286     size_t nrow,  /* # of rows of A */
00287     size_t ncol,  /* # of columns of A */
00288     int xtype,    /* CHOLMOD_PATTERN, _REAL, _COMPLEX, or _ZOMPLEX */
00289     /* --------------- */
00290     cholmod_common *Common
00291 )
00292 {
00293     double *Ax, *Az ;
00294     cholmod_sparse *A ;
00295     Int *Ap, *Ai ;
00296     Int j, n ;
00297 
00298     /* ---------------------------------------------------------------------- */
00299     /* get inputs */
00300     /* ---------------------------------------------------------------------- */
00301 
00302     RETURN_IF_NULL_COMMON (NULL) ;
00303     Common->status = CHOLMOD_OK ;
00304 
00305     /* ---------------------------------------------------------------------- */
00306     /* allocate the matrix */
00307     /* ---------------------------------------------------------------------- */
00308 
00309     n = MIN (nrow, ncol) ;
00310     A = CHOLMOD(allocate_sparse) (nrow, ncol, n, TRUE, TRUE, 0, xtype,
00311       Common) ;
00312 
00313     if (Common->status < CHOLMOD_OK)
00314     {
00315   return (NULL) ;     /* out of memory or inputs invalid */
00316     }
00317 
00318     /* ---------------------------------------------------------------------- */
00319     /* create the identity matrix */
00320     /* ---------------------------------------------------------------------- */
00321 
00322     Ap = A->p ;
00323     Ai = A->i ;
00324     Ax = A->x ;
00325     Az = A->z ;
00326 
00327     for (j = 0 ; j < n ; j++)
00328     {
00329   Ap [j] = j ;
00330     }
00331     for (j = n ; j <= ((Int) ncol) ; j++)
00332     {
00333   Ap [j] = n ;
00334     }
00335     for (j = 0 ; j < n ; j++)
00336     {
00337   Ai [j] = j ;
00338     }
00339 
00340     switch (xtype)
00341     {
00342   case CHOLMOD_REAL:
00343       for (j = 0 ; j < n ; j++)
00344       {
00345     Ax [j] = 1 ;
00346       }
00347       break ;
00348 
00349   case CHOLMOD_COMPLEX:
00350       for (j = 0 ; j < n ; j++)
00351       {
00352     Ax [2*j  ] = 1 ;
00353     Ax [2*j+1] = 0 ;
00354       }
00355       break ;
00356 
00357   case CHOLMOD_ZOMPLEX:
00358       for (j = 0 ; j < n ; j++)
00359       {
00360     Ax [j] = 1 ;
00361       }
00362       for (j = 0 ; j < n ; j++)
00363       {
00364     Az [j] = 0 ;
00365       }
00366       break ;
00367     }
00368 
00369     return (A) ;
00370 }
00371 
00372 
00373 /* ========================================================================== */
00374 /* === cholmod_spzeros ====================================================== */
00375 /* ========================================================================== */
00376 
00377 /* Return a sparse zero matrix. */
00378 
00379 cholmod_sparse *CHOLMOD(spzeros)
00380 (
00381     /* ---- input ---- */
00382     size_t nrow,  /* # of rows of A */
00383     size_t ncol,  /* # of columns of A */
00384     size_t nzmax, /* max # of nonzeros of A */
00385     int xtype,    /* CHOLMOD_PATTERN, _REAL, _COMPLEX, or _ZOMPLEX */
00386     /* --------------- */
00387     cholmod_common *Common
00388 )
00389 {
00390 
00391     /* ---------------------------------------------------------------------- */
00392     /* get inputs */
00393     /* ---------------------------------------------------------------------- */
00394 
00395     RETURN_IF_NULL_COMMON (NULL) ;
00396     Common->status = CHOLMOD_OK ;
00397 
00398     /* ---------------------------------------------------------------------- */
00399     /* allocate the matrix */
00400     /* ---------------------------------------------------------------------- */
00401 
00402     return (CHOLMOD(allocate_sparse) (nrow, ncol, nzmax, TRUE, TRUE, 0, xtype,
00403       Common)) ;
00404 }
00405 
00406 
00407 /* ========================================================================== */
00408 /* === cholmod_nnz ========================================================== */
00409 /* ========================================================================== */
00410 
00411 /* Return the number of entries in a sparse matrix.
00412  *
00413  * workspace: none
00414  * integer overflow cannot occur, since the matrix is already allocated.
00415  */
00416 
00417 UF_long CHOLMOD(nnz)
00418 (
00419     /* ---- input ---- */
00420     cholmod_sparse *A,
00421     /* --------------- */
00422     cholmod_common *Common
00423 )
00424 {
00425     Int *Ap, *Anz ;
00426     size_t nz ;
00427     Int j, ncol ;
00428 
00429     /* ---------------------------------------------------------------------- */
00430     /* get inputs */
00431     /* ---------------------------------------------------------------------- */
00432 
00433     RETURN_IF_NULL_COMMON (EMPTY) ;
00434     RETURN_IF_NULL (A, EMPTY) ;
00435     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, EMPTY) ;
00436     Common->status = CHOLMOD_OK ;
00437 
00438     /* ---------------------------------------------------------------------- */
00439     /* return nnz (A) */
00440     /* ---------------------------------------------------------------------- */
00441 
00442     ncol = A->ncol ;
00443     if (A->packed)
00444     {
00445   Ap = A->p ;
00446   RETURN_IF_NULL (Ap, EMPTY) ;
00447   nz = Ap [ncol] ;
00448     }
00449     else
00450     {
00451   Anz = A->nz ;
00452   RETURN_IF_NULL (Anz, EMPTY) ;
00453   nz = 0 ;
00454   for (j = 0 ; j < ncol ; j++)
00455   {
00456       nz += MAX (0, Anz [j]) ;
00457   }
00458     }
00459     return (nz) ;
00460 }
00461 
00462 
00463 /* ========================================================================== */
00464 /* === cholmod_copy_sparse ================================================== */
00465 /* ========================================================================== */
00466 
00467 /* C = A.  Create an exact copy of a sparse matrix, with one exception.
00468  * Entries in unused space are not copied (they might not be initialized,
00469  * and copying them would cause program checkers such as purify and
00470  * valgrind to complain).  The xtype of the resulting matrix C is the same as
00471  * the xtype of the input matrix A.
00472  *
00473  * See also Core/cholmod_copy, which copies a matrix with possible changes
00474  * in stype, presence of diagonal entries, pattern vs. numerical values,
00475  * real and/or imaginary parts, and so on.
00476  */
00477 
00478 cholmod_sparse *CHOLMOD(copy_sparse)
00479 (
00480     /* ---- input ---- */
00481     cholmod_sparse *A,  /* matrix to copy */
00482     /* --------------- */
00483     cholmod_common *Common
00484 )
00485 {
00486     double *Ax, *Cx, *Az, *Cz ;
00487     Int *Ap, *Ai, *Anz, *Cp, *Ci, *Cnz ;
00488     cholmod_sparse *C ;
00489     Int p, pend, j, ncol, packed, nzmax, nz, xtype ;
00490 
00491     /* ---------------------------------------------------------------------- */
00492     /* check inputs */
00493     /* ---------------------------------------------------------------------- */
00494 
00495     RETURN_IF_NULL_COMMON (NULL) ;
00496     RETURN_IF_NULL (A, NULL) ;
00497     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, NULL) ;
00498     if (A->stype != 0 && A->nrow != A->ncol)
00499     {
00500   ERROR (CHOLMOD_INVALID, "rectangular matrix with stype != 0 invalid") ;
00501   return (NULL) ;
00502     }
00503     Common->status = CHOLMOD_OK ;
00504     ASSERT (CHOLMOD(dump_sparse) (A, "A original", Common) >= 0) ;
00505 
00506     /* ---------------------------------------------------------------------- */
00507     /* get inputs */
00508     /* ---------------------------------------------------------------------- */
00509 
00510     ncol = A->ncol ;
00511     nzmax = A->nzmax ;
00512     packed = A->packed ;
00513     Ap = A->p ;
00514     Ai = A->i ;
00515     Ax = A->x ;
00516     Az = A->z ;
00517     Anz = A->nz ;
00518     xtype = A->xtype ;
00519 
00520     /* ---------------------------------------------------------------------- */
00521     /* allocate the copy */
00522     /* ---------------------------------------------------------------------- */
00523 
00524     C = CHOLMOD(allocate_sparse) (A->nrow, A->ncol, A->nzmax, A->sorted,
00525       A->packed, A->stype, A->xtype, Common) ;
00526 
00527     if (Common->status < CHOLMOD_OK)
00528     {
00529   return (NULL) ;     /* out of memory */
00530     }
00531 
00532     Cp = C->p ;
00533     Ci = C->i ;
00534     Cx = C->x ;
00535     Cz = C->z ;
00536     Cnz = C->nz ;
00537 
00538     /* ---------------------------------------------------------------------- */
00539     /* copy the matrix */
00540     /* ---------------------------------------------------------------------- */
00541 
00542     for (j = 0 ; j <= ncol ; j++)
00543     {
00544   Cp [j] = Ap [j] ;
00545     }
00546 
00547     if (packed)
00548     {
00549   nz = Ap [ncol] ;
00550   for (p = 0 ; p < nz ; p++)
00551   {
00552       Ci [p] = Ai [p] ;
00553   }
00554 
00555   switch (xtype)
00556   {
00557       case CHOLMOD_REAL:
00558     for (p = 0 ; p < nz ; p++)
00559     {
00560         Cx [p] = Ax [p] ;
00561     }
00562     break ;
00563 
00564       case CHOLMOD_COMPLEX:
00565     for (p = 0 ; p < 2*nz ; p++)
00566     {
00567         Cx [p] = Ax [p] ;
00568     }
00569     break ;
00570 
00571       case CHOLMOD_ZOMPLEX:
00572     for (p = 0 ; p < nz ; p++)
00573     {
00574         Cx [p] = Ax [p] ;
00575         Cz [p] = Az [p] ;
00576     }
00577     break ;
00578   }
00579 
00580     }
00581     else
00582     {
00583 
00584   for (j = 0 ; j < ncol ; j++)
00585   {
00586       Cnz [j] = Anz [j] ;
00587   }
00588 
00589   switch (xtype)
00590   {
00591       case CHOLMOD_PATTERN:
00592     for (j = 0 ; j < ncol ; j++)
00593     {
00594         p = Ap [j] ;
00595         pend = p + Anz [j] ;
00596         for ( ; p < pend ; p++)
00597         {
00598       Ci [p] = Ai [p] ;
00599         }
00600     }
00601     break ;
00602 
00603       case CHOLMOD_REAL:
00604     for (j = 0 ; j < ncol ; j++)
00605     {
00606         p = Ap [j] ;
00607         pend = p + Anz [j] ;
00608         for ( ; p < pend ; p++)
00609         {
00610       Ci [p] = Ai [p] ;
00611       Cx [p] = Ax [p] ;
00612         }
00613     }
00614     break ;
00615 
00616       case CHOLMOD_COMPLEX:
00617     for (j = 0 ; j < ncol ; j++)
00618     {
00619         p = Ap [j] ;
00620         pend = p + Anz [j] ;
00621         for ( ; p < pend ; p++)
00622         {
00623       Ci [p] = Ai [p] ;
00624       Cx [2*p  ] = Ax [2*p  ] ;
00625       Cx [2*p+1] = Ax [2*p+1] ;
00626         }
00627     }
00628     break ;
00629 
00630       case CHOLMOD_ZOMPLEX:
00631     for (j = 0 ; j < ncol ; j++)
00632     {
00633         p = Ap [j] ;
00634         pend = p + Anz [j] ;
00635         for ( ; p < pend ; p++)
00636         {
00637       Ci [p] = Ai [p] ;
00638       Cx [p] = Ax [p] ;
00639       Cz [p] = Az [p] ;
00640         }
00641     }
00642     break ;
00643   }
00644     }
00645 
00646     /* ---------------------------------------------------------------------- */
00647     /* return the result */
00648     /* ---------------------------------------------------------------------- */
00649 
00650     ASSERT (CHOLMOD(dump_sparse) (C, "C copy", Common) >= 0) ;
00651     return (C) ;
00652 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines