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