Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_triplet.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === Core/cholmod_triplet ================================================= */
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_triplet object:
00015  *
00016  * A sparse matrix held in triplet form is the simplest one for a user to
00017  * create.  It consists of a list of nz entries in arbitrary order, held in
00018  * three arrays: i, j, and x, each of length nk.  The kth entry is in row i[k],
00019  * column j[k], with value x[k].  There may be duplicate values; if A(i,j)
00020  * appears more than once, its value is the sum of the entries with those row
00021  * and column indices.
00022  *
00023  * Primary routines:
00024  * -----------------
00025  * cholmod_allocate_triplet allocate a triplet matrix
00026  * cholmod_free_triplet   free a triplet matrix
00027  *
00028  * Secondary routines:
00029  * -------------------
00030  * cholmod_reallocate_triplet reallocate a triplet matrix
00031  * cholmod_sparse_to_triplet  create a triplet matrix copy of a sparse matrix
00032  * cholmod_triplet_to_sparse  create a sparse matrix copy of a triplet matrix
00033  * cholmod_copy_triplet   create a copy of a triplet matrix
00034  *
00035  * The relationship between an m-by-n cholmod_sparse matrix A and a
00036  * cholmod_triplet matrix (i, j, and x) is identical to how they are used in
00037  * the MATLAB "sparse" and "find" functions:
00038  *
00039  *  [i j x] = find (A)
00040  *  [m n] = size (A)
00041  *  A = sparse (i,j,x,m,n)
00042  *
00043  * with the exception that the cholmod_sparse matrix may be "unpacked", may
00044  * have either sorted or unsorted columns (depending on the option selected),
00045  * and may be symmetric with just the upper or lower triangular part stored.
00046  * Likewise, the cholmod_triplet matrix may contain just the entries in the
00047  * upper or lower triangular part of a symmetric matrix.
00048  *
00049  * MATLAB sparse matrices are always "packed", always have sorted columns,
00050  * and always store both parts of a symmetric matrix.  In some cases, MATLAB
00051  * behaves like CHOLMOD by ignoring entries in the upper or lower triangular
00052  * part of a matrix that is otherwise assumed to be symmetric (such as the
00053  * input to chol).  In CHOLMOD, that option is a characteristic of the object.
00054  * In MATLAB, that option is based on how a matrix is used as the input to
00055  * a function.
00056  *
00057  * The triplet matrix is provided to give the user a simple way of constructing
00058  * a sparse matrix.  There are very few operations supported for triplet
00059  * matrices.  The assumption is that they will be converted to cholmod_sparse
00060  * matrix form first.
00061  *
00062  * Adding two triplet matrices simply involves concatenating the contents of
00063  * the three arrays (i, j, and x).   To permute a triplet matrix, just replace
00064  * the row and column indices with their permuted values.  For example, if
00065  * P is a permutation vector, then P [k] = j means row/column j is the kth
00066  * row/column in C=P*A*P'.  In MATLAB notation, C=A(p,p).  If Pinv is an array
00067  * of size n and T is the triplet form of A, then:
00068  *
00069  *  Ti = T->i ;
00070  *  Tj = T->j ;
00071  *  for (k = 0 ; k < n  ; k++) Pinv [P [k]] = k ;
00072  *  for (k = 0 ; k < nz ; k++) Ti [k] = Pinv [Ti [k]] ;
00073  *  for (k = 0 ; k < nz ; k++) Tj [k] = Pinv [Tj [k]] ;
00074  *
00075  * overwrites T with the triplet form of C=P*A*P'.  The conversion
00076  *
00077  *  C = cholmod_triplet_to_sparse (T, 0, &Common) ;
00078  *
00079  * will then return the matrix C = P*A*P'.
00080  *
00081  * Note that T->stype > 0 means that entries in the lower triangular part of
00082  * T are transposed into the upper triangular part when T is converted to
00083  * sparse matrix (cholmod_sparse) form with cholmod_triplet_to_sparse.  The
00084  * opposite is true for T->stype < 0.
00085  *
00086  * Since the triplet matrix T is so simple to generate, it's quite easy
00087  * to remove entries that you do not want, prior to converting T to the
00088  * cholmod_sparse form.  So if you include these entries in T, CHOLMOD
00089  * assumes that there must be a reason (such as the one above).  Thus,
00090  * no entry in a triplet matrix is ever ignored.
00091  *
00092  * Other operations, such as extacting a submatrix, horizontal and vertical
00093  * concatenation, multiply a triplet matrix times a dense matrix, are also
00094  * simple.  Multiplying two triplet matrices is not trivial; the simplest
00095  * method is to convert them to cholmod_sparse matrices first.
00096  *
00097  * Supports all xtypes (pattern, real, complex, and zomplex).
00098  */
00099  
00100 #include "amesos_cholmod_internal.h"
00101 #include "amesos_cholmod_core.h"
00102 
00103 
00104 /* ========================================================================== */
00105 /* === TEMPLATE ============================================================= */
00106 /* ========================================================================== */
00107 
00108 #define PATTERN
00109 #include "amesos_t_cholmod_triplet.c"
00110 #define REAL
00111 #include "amesos_t_cholmod_triplet.c"
00112 #define COMPLEX
00113 #include "amesos_t_cholmod_triplet.c"
00114 #define ZOMPLEX
00115 #include "amesos_t_cholmod_triplet.c"
00116 
00117 
00118 /* ========================================================================== */
00119 /* === cholmod_allocate_triplet ============================================= */
00120 /* ========================================================================== */
00121 
00122 /* allocate space for a triplet matrix
00123  *
00124  * workspace: none
00125  */
00126 
00127 cholmod_triplet *CHOLMOD(allocate_triplet)
00128 (
00129     /* ---- input ---- */
00130     size_t nrow,  /* # of rows of T */
00131     size_t ncol,  /* # of columns of T */
00132     size_t nzmax, /* max # of nonzeros of T */
00133     int stype,    /* stype of T */
00134     int xtype,    /* CHOLMOD_PATTERN, _REAL, _COMPLEX, or _ZOMPLEX */
00135     /* --------------- */
00136     cholmod_common *Common
00137 )
00138 {
00139     cholmod_triplet *T ;
00140     size_t nzmax0 ;
00141     int ok = TRUE ;
00142 
00143     /* ---------------------------------------------------------------------- */
00144     /* get inputs */
00145     /* ---------------------------------------------------------------------- */
00146 
00147     RETURN_IF_NULL_COMMON (NULL) ;
00148     if (xtype < CHOLMOD_PATTERN || xtype > CHOLMOD_ZOMPLEX)
00149     {
00150   ERROR (CHOLMOD_INVALID, "xtype invalid") ;
00151   return (NULL) ;
00152     }
00153     /* ensure the dimensions do not cause integer overflow */
00154     (void) CHOLMOD(add_size_t) (ncol, 2, &ok) ;
00155     if (!ok || nrow > Int_max || ncol > Int_max || nzmax > Int_max)
00156     {
00157   ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
00158   return (NULL) ;
00159     }
00160 
00161     Common->status = CHOLMOD_OK ;
00162 
00163     /* ---------------------------------------------------------------------- */
00164     /* allocate header */
00165     /* ---------------------------------------------------------------------- */
00166 
00167     T = CHOLMOD(malloc) (sizeof (cholmod_triplet), 1, Common) ;
00168     if (Common->status < CHOLMOD_OK)
00169     {
00170   return (NULL) ;     /* out of memory */
00171     }
00172 
00173     PRINT1 (("cholmod_allocate_triplet %d-by-%d nzmax %d xtype %d\n",
00174     nrow, ncol, nzmax, xtype)) ;
00175 
00176     nzmax = MAX (1, nzmax) ;
00177 
00178     T->nrow = nrow ;
00179     T->ncol = ncol ;
00180     T->nzmax = nzmax ;
00181     T->nnz = 0 ;
00182     T->stype = stype ;
00183     T->itype = ITYPE ;
00184     T->xtype = xtype ;
00185     T->dtype = DTYPE ;
00186 
00187     T->j = NULL ;
00188     T->i = NULL ;
00189     T->x = NULL ;
00190     T->z = NULL ;
00191 
00192     /* ---------------------------------------------------------------------- */
00193     /* allocate the matrix itself */
00194     /* ---------------------------------------------------------------------- */
00195 
00196     nzmax0 = 0 ;
00197     CHOLMOD(realloc_multiple) (nzmax, 2, xtype, &(T->i), &(T->j),
00198     &(T->x), &(T->z), &nzmax0, Common) ;
00199 
00200     if (Common->status < CHOLMOD_OK)
00201     {
00202   CHOLMOD(free_triplet) (&T, Common) ;
00203   return (NULL) ;     /* out of memory */
00204     }
00205 
00206     return (T) ;
00207 }
00208 
00209 
00210 /* ========================================================================== */
00211 /* === cholmod_free_triplet ================================================= */
00212 /* ========================================================================== */
00213 
00214 /* free a triplet matrix
00215  *
00216  * workspace: none
00217  */
00218 
00219 int CHOLMOD(free_triplet)
00220 (
00221     /* ---- in/out --- */
00222     cholmod_triplet **THandle,    /* matrix to deallocate, NULL on output */
00223     /* --------------- */
00224     cholmod_common *Common
00225 )
00226 {
00227     Int nz ;
00228     cholmod_triplet *T ;
00229 
00230     RETURN_IF_NULL_COMMON (FALSE) ;
00231 
00232     if (THandle == NULL)
00233     {
00234   /* nothing to do */
00235   return (TRUE) ;
00236     }
00237     T = *THandle ;
00238     if (T == NULL)
00239     {
00240   /* nothing to do */
00241   return (TRUE) ;
00242     }
00243     nz = T->nzmax ;
00244     T->j = CHOLMOD(free) (nz, sizeof (Int), T->j, Common) ;
00245     T->i = CHOLMOD(free) (nz, sizeof (Int), T->i, Common) ;
00246     if (T->xtype == CHOLMOD_REAL)
00247     {
00248   T->x = CHOLMOD(free) (nz, sizeof (double), T->x, Common) ;
00249     }
00250     else if (T->xtype == CHOLMOD_COMPLEX)
00251     {
00252   T->x = CHOLMOD(free) (nz, 2*sizeof (double), T->x, Common) ;
00253     }
00254     else if (T->xtype == CHOLMOD_ZOMPLEX)
00255     {
00256   T->x = CHOLMOD(free) (nz, sizeof (double), T->x, Common) ;
00257   T->z = CHOLMOD(free) (nz, sizeof (double), T->z, Common) ;
00258     }
00259     *THandle = CHOLMOD(free) (1, sizeof (cholmod_triplet), (*THandle), Common) ;
00260     return (TRUE) ;
00261 }
00262 
00263 
00264 /* ========================================================================== */
00265 /* === cholmod_reallocate_triplet =========================================== */
00266 /* ========================================================================== */
00267 
00268 /* Change the size of T->i, T->j, and T->x, or allocate them if their current
00269  * size is zero.  T->x is not modified if T->xtype is CHOLMOD_PATTERN.
00270  *
00271  * workspace: none
00272  */
00273 
00274 int CHOLMOD(reallocate_triplet)
00275 (
00276     /* ---- input ---- */
00277     size_t nznew, /* new # of entries in T */
00278     /* ---- in/out --- */
00279     cholmod_triplet *T, /* triplet matrix to modify */
00280     /* --------------- */
00281     cholmod_common *Common
00282 )
00283 {
00284 
00285     /* ---------------------------------------------------------------------- */
00286     /* get inputs */
00287     /* ---------------------------------------------------------------------- */
00288 
00289     RETURN_IF_NULL_COMMON (FALSE) ;
00290     RETURN_IF_NULL (T, FALSE) ;
00291     RETURN_IF_XTYPE_INVALID (T, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00292     Common->status = CHOLMOD_OK ;
00293     PRINT1 (("realloc triplet %d to %d, xtype: %d\n",
00294     T->nzmax, nznew, T->xtype)) ;
00295 
00296     /* ---------------------------------------------------------------------- */
00297     /* resize the matrix */
00298     /* ---------------------------------------------------------------------- */
00299 
00300     CHOLMOD(realloc_multiple) (MAX (1,nznew), 2, T->xtype, &(T->i), &(T->j),
00301       &(T->x), &(T->z), &(T->nzmax), Common) ;
00302 
00303     return (Common->status == CHOLMOD_OK) ;
00304 }
00305 
00306 
00307 /* ========================================================================== */
00308 /* === cholmod_triplet_to_sparse ============================================ */
00309 /* ========================================================================== */
00310 
00311 /* Convert a set of triplets into a cholmod_sparse matrix.  In MATLAB notation,
00312  * for unsymmetric matrices:
00313  *
00314  *  A = sparse (Ti, Tj, Tx, nrow, ncol, nzmax) ;
00315  *
00316  * For the symmetric upper case:
00317  *
00318  *  A = sparse (min(Ti,Tj), max(Ti,Tj), Tx, nrow, ncol, nzmax) ;
00319  *
00320  * For the symmetric lower case:
00321  *
00322  *  A = sparse (max(Ti,Tj), min(Ti,Tj), Tx, nrow, ncol, nzmax) ;
00323  *
00324  * If Tx is NULL, then A->x is not allocated, and only the pattern of A is
00325  * computed.  A is returned in packed form, and can be of any stype
00326  * (upper/lower/unsymmetric).  It has enough space to hold the values in T,
00327  * or nzmax, whichever is larger.
00328  *
00329  * workspace: Iwork (max (nrow,ncol))
00330  *  allocates a temporary copy of its output matrix.
00331  *
00332  * The resulting sparse matrix has the same xtype as the input triplet matrix.
00333  */
00334 
00335 cholmod_sparse *CHOLMOD(triplet_to_sparse)
00336 (
00337     /* ---- input ---- */
00338     cholmod_triplet *T, /* matrix to copy */
00339     size_t nzmax, /* allocate at least this much space in output matrix */
00340     /* --------------- */
00341     cholmod_common *Common
00342 )
00343 {
00344     cholmod_sparse *R, *A = NULL ;
00345     Int *Wj, *Rp, *Ri, *Rnz, *Ti, *Tj ;
00346     Int i, j, p, k, stype, nrow, ncol, nz, ok ;
00347     size_t anz = 0 ;
00348 
00349     /* ---------------------------------------------------------------------- */
00350     /* check inputs */
00351     /* ---------------------------------------------------------------------- */
00352 
00353     RETURN_IF_NULL_COMMON (NULL) ;
00354     RETURN_IF_NULL (T, NULL) ;
00355     Ti = T->i ;
00356     Tj = T->j ;
00357     RETURN_IF_NULL (Ti, NULL) ;
00358     RETURN_IF_NULL (Tj, NULL) ;
00359     RETURN_IF_XTYPE_INVALID (T, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, NULL) ;
00360     stype = SIGN (T->stype) ;
00361     if (stype && T->nrow != T->ncol)
00362     {
00363   /* inputs invalid */
00364   ERROR (CHOLMOD_INVALID, "matrix invalid") ;
00365   return (NULL) ;
00366     }
00367     Common->status = CHOLMOD_OK ;
00368     DEBUG (CHOLMOD(dump_triplet) (T, "T", Common)) ;
00369 
00370     /* ---------------------------------------------------------------------- */
00371     /* get inputs */
00372     /* ---------------------------------------------------------------------- */
00373 
00374     nrow = T->nrow ;
00375     ncol = T->ncol ;
00376     nz = T->nnz ;
00377 
00378     /* ---------------------------------------------------------------------- */
00379     /* allocate workspace */
00380     /* ---------------------------------------------------------------------- */
00381 
00382     CHOLMOD(allocate_work) (0, MAX (nrow, ncol), 0, Common) ;
00383     if (Common->status < CHOLMOD_OK)
00384     {
00385   return (NULL) ;     /* out of memory */
00386     }
00387 
00388     /* ---------------------------------------------------------------------- */
00389     /* allocate temporary matrix R */
00390     /* ---------------------------------------------------------------------- */
00391 
00392     R = CHOLMOD(allocate_sparse) (ncol, nrow, nz, FALSE, FALSE, -stype,
00393       T->xtype, Common) ;
00394 
00395     if (Common->status < CHOLMOD_OK)
00396     {
00397   return (NULL) ;     /* out of memory */
00398     }
00399 
00400     Rp = R->p ;
00401     Ri = R->i ;
00402     Rnz = R->nz ;
00403 
00404     /* ---------------------------------------------------------------------- */
00405     /* count the entries in each row of A (also counting duplicates) */
00406     /* ---------------------------------------------------------------------- */
00407 
00408     for (i = 0 ; i < nrow ; i++)
00409     {
00410   Rnz [i] = 0 ; 
00411     }
00412 
00413     if (stype > 0)
00414     {
00415   for (k = 0 ; k < nz ; k++)
00416   {
00417       i = Ti [k] ;
00418       j = Tj [k] ;
00419       if (i < 0 || i >= nrow || j < 0 || j >= ncol)
00420       {
00421     ERROR (CHOLMOD_INVALID, "index out of range") ;
00422     break ;
00423       }
00424       /* A will be symmetric with just the upper triangular part stored.
00425        * Create a matrix R that is lower triangular.  Entries in the
00426        * upper part of R are transposed to the lower part. */
00427       Rnz [MIN (i,j)]++ ;
00428   }
00429     }
00430     else if (stype < 0)
00431     {
00432   for (k = 0 ; k < nz ; k++)
00433   {
00434       i = Ti [k] ;
00435       j = Tj [k] ;
00436       if (i < 0 || i >= nrow || j < 0 || j >= ncol)
00437       {
00438     ERROR (CHOLMOD_INVALID, "index out of range") ;
00439     break ;
00440       }
00441       /* A will be symmetric with just the lower triangular part stored.
00442        * Create a matrix R that is upper triangular.  Entries in the
00443        * lower part of R are transposed to the upper part. */
00444       Rnz [MAX (i,j)]++ ;
00445   }
00446     }
00447     else
00448     {
00449   for (k = 0 ; k < nz ; k++)
00450   {
00451       i = Ti [k] ;
00452       j = Tj [k] ;
00453       if (i < 0 || i >= nrow || j < 0 || j >= ncol)
00454       {
00455     ERROR (CHOLMOD_INVALID, "index out of range") ;
00456     break ;
00457       }
00458       /* constructing an unsymmetric matrix */
00459       Rnz [i]++ ;
00460   }
00461     }
00462 
00463     if (Common->status < CHOLMOD_OK)
00464     {
00465   /* triplet matrix is invalid */
00466   CHOLMOD(free_sparse) (&R, Common) ;
00467   return (NULL) ;
00468     }
00469 
00470     /* ---------------------------------------------------------------------- */
00471     /* construct the row pointers */
00472     /* ---------------------------------------------------------------------- */
00473 
00474     p = 0 ;
00475     for (i = 0 ; i < nrow ; i++)
00476     {
00477   Rp [i] = p ;
00478   p += Rnz [i] ;
00479     }
00480     Rp [nrow] = p ;
00481 
00482     /* use Wj (i/l/l) as temporary row pointers */
00483     Wj = Common->Iwork ;  /* size MAX (nrow,ncol) FUTURE WORK: (i/l/l) */
00484     for (i = 0 ; i < nrow ; i++)
00485     {
00486   Wj [i] = Rp [i] ;
00487     }
00488 
00489     /* ---------------------------------------------------------------------- */
00490     /* construct triplet matrix, using template routine */
00491     /* ---------------------------------------------------------------------- */
00492 
00493     switch (T->xtype)
00494     {
00495   case CHOLMOD_PATTERN:
00496       anz = amesos_p_cholmod_triplet_to_sparse (T, R, Common) ;
00497       break ;
00498 
00499   case CHOLMOD_REAL:
00500       anz = amesos_r_cholmod_triplet_to_sparse (T, R, Common) ;
00501       break ;
00502 
00503   case CHOLMOD_COMPLEX:
00504       anz = amesos_c_cholmod_triplet_to_sparse (T, R, Common) ;
00505       break ;
00506 
00507   case CHOLMOD_ZOMPLEX:
00508       anz = amesos_z_cholmod_triplet_to_sparse (T, R, Common) ;
00509       break ;
00510     }
00511 
00512     /* ---------------------------------------------------------------------- */
00513     /* A = R' (array transpose, not complex conjugate transpose) */
00514     /* ---------------------------------------------------------------------- */
00515 
00516     /* workspace: Iwork (R->nrow), which is A->ncol */
00517 
00518     ASSERT (CHOLMOD(dump_sparse) (R, "R", Common) >= 0) ;
00519 
00520     A = CHOLMOD(allocate_sparse) (nrow, ncol, MAX (anz, nzmax), TRUE, TRUE,
00521   stype, T->xtype, Common) ;
00522 
00523     if (stype)
00524     {
00525   ok = CHOLMOD(transpose_sym) (R, 1, NULL, A, Common) ;
00526     }
00527     else
00528     {
00529   ok = CHOLMOD(transpose_unsym) (R, 1, NULL, NULL, 0, A, Common) ; 
00530     }
00531 
00532     CHOLMOD(free_sparse) (&R, Common) ;
00533     if (Common->status < CHOLMOD_OK)
00534     {
00535   CHOLMOD(free_sparse) (&A, Common) ;
00536     }
00537 
00538     /* ---------------------------------------------------------------------- */
00539     /* return result */
00540     /* ---------------------------------------------------------------------- */
00541 
00542     ASSERT (CHOLMOD(dump_sparse) (A, "A = triplet(T) result", Common) >= 0) ;
00543     return (A) ;
00544 }
00545 
00546 
00547 /* ========================================================================== */
00548 /* === cholmod_sparse_to_triplet ============================================ */
00549 /* ========================================================================== */
00550 
00551 /* Converts a sparse column-oriented matrix to triplet form.
00552  * The resulting triplet matrix has the same xtype as the sparse matrix.
00553  *
00554  * workspace: none
00555  */
00556 
00557 cholmod_triplet *CHOLMOD(sparse_to_triplet)
00558 (
00559     /* ---- input ---- */
00560     cholmod_sparse *A,  /* matrix to copy */
00561     /* --------------- */
00562     cholmod_common *Common
00563 )
00564 {
00565     double *Ax, *Az, *Tx, *Tz ;
00566     Int *Ap, *Ai, *Ti, *Tj, *Anz ;
00567     cholmod_triplet *T ;
00568     Int i, xtype, p, pend, k, j, nrow, ncol, nz, stype, packed, up, lo,
00569   both ;
00570 
00571     /* ---------------------------------------------------------------------- */
00572     /* check inputs */
00573     /* ---------------------------------------------------------------------- */
00574 
00575     RETURN_IF_NULL_COMMON (NULL) ;
00576     RETURN_IF_NULL (A, NULL) ;
00577     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, NULL) ;
00578     stype = SIGN (A->stype) ;
00579     nrow = A->nrow ;
00580     ncol = A->ncol ;
00581     if (stype && nrow != ncol)
00582     {
00583   /* inputs invalid */
00584   ERROR (CHOLMOD_INVALID, "matrix invalid") ;
00585   return (NULL) ;
00586     }
00587     Ax = A->x ;
00588     Az = A->z ;
00589     xtype = A->xtype ;
00590     Common->status = CHOLMOD_OK ;
00591 
00592     ASSERT (CHOLMOD(dump_sparse) (A, "A", Common) >= 0) ;
00593 
00594     /* ---------------------------------------------------------------------- */
00595     /* allocate triplet matrix */
00596     /* ---------------------------------------------------------------------- */
00597 
00598     nz = CHOLMOD(nnz) (A, Common) ;
00599     T = CHOLMOD(allocate_triplet) (nrow, ncol, nz, A->stype, A->xtype, Common) ;
00600 
00601     if (Common->status < CHOLMOD_OK)
00602     {
00603   return (NULL) ;     /* out of memory */
00604     }
00605 
00606     /* ---------------------------------------------------------------------- */
00607     /* convert to a sparse matrix */
00608     /* ---------------------------------------------------------------------- */
00609 
00610     Ap = A->p ;
00611     Ai = A->i ;
00612     Anz = A->nz ;
00613     packed = A->packed ;
00614 
00615     Ti = T->i ;
00616     Tj = T->j ;
00617     Tx = T->x ;
00618     Tz = T->z ;
00619     T->stype = A->stype ;
00620 
00621     both = (A->stype == 0) ;
00622     up = (A->stype > 0) ;
00623     lo = (A->stype < 0) ;
00624 
00625     k = 0 ;
00626 
00627     for (j = 0 ; j < ncol ; j++)
00628     {
00629   p = Ap [j] ;
00630   pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
00631   for ( ; p < pend ; p++)
00632   {
00633       i = Ai [p] ;
00634       if (both || (up && i <= j) || (lo && i >= j))
00635       {
00636     Ti [k] = Ai [p] ;
00637     Tj [k] = j ;
00638 
00639     if (xtype == CHOLMOD_REAL)
00640     {
00641         Tx [k] = Ax [p] ;
00642     }
00643     else if (xtype == CHOLMOD_COMPLEX)
00644     {
00645         Tx [2*k  ] = Ax [2*p  ] ;
00646         Tx [2*k+1] = Ax [2*p+1] ;
00647     }
00648     else if (xtype == CHOLMOD_ZOMPLEX)
00649     {
00650         Tx [k] = Ax [p] ;
00651         Tz [k] = Az [p] ;
00652     }
00653 
00654     k++ ;
00655     ASSERT (k <= nz) ;
00656       }
00657   }
00658     }
00659 
00660     T->nnz = k ;
00661 
00662     /* ---------------------------------------------------------------------- */
00663     /* return result */
00664     /* ---------------------------------------------------------------------- */
00665 
00666     ASSERT (CHOLMOD(dump_triplet) (T, "T", Common)) ;
00667     return (T) ;
00668 }
00669 
00670 
00671 /* ========================================================================== */
00672 /* === cholmod_copy_triplet ================================================= */
00673 /* ========================================================================== */
00674 
00675 /* Create an exact copy of a triplet matrix, except that entries in unused
00676  * space are not copied (they might not be initialized, and copying them would
00677  * cause program checkers such as purify and valgrind to complain).
00678  * The output triplet matrix has the same xtype as the input triplet matrix.
00679  */
00680 
00681 cholmod_triplet *CHOLMOD(copy_triplet)
00682 (
00683     /* ---- input ---- */
00684     cholmod_triplet *T, /* matrix to copy */
00685     /* --------------- */
00686     cholmod_common *Common
00687 )
00688 {
00689     double *Tx, *Tz, *Cx, *Cz ;
00690     Int *Ci, *Cj, *Ti, *Tj ;
00691     cholmod_triplet *C ;
00692     Int xtype, k, nz ;
00693 
00694     /* ---------------------------------------------------------------------- */
00695     /* check inputs */
00696     /* ---------------------------------------------------------------------- */
00697 
00698     RETURN_IF_NULL_COMMON (NULL) ;
00699     RETURN_IF_NULL (T, NULL) ;
00700     RETURN_IF_XTYPE_INVALID (T, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, NULL) ;
00701     nz = T->nnz ;
00702     Ti = T->i ;
00703     Tj = T->j ;
00704     Tx = T->x ;
00705     Tz = T->z ;
00706     xtype = T->xtype ;
00707     RETURN_IF_NULL (Ti, NULL) ;
00708     RETURN_IF_NULL (Tj, NULL) ;
00709     Common->status = CHOLMOD_OK ;
00710     DEBUG (CHOLMOD(dump_triplet) (T, "T input", Common)) ;
00711 
00712     /* ---------------------------------------------------------------------- */
00713     /* allocate copy */
00714     /* ---------------------------------------------------------------------- */
00715 
00716     C = CHOLMOD(allocate_triplet) (T->nrow, T->ncol, T->nzmax, T->stype,
00717       xtype, Common) ;
00718 
00719     if (Common->status < CHOLMOD_OK)
00720     {
00721   return (NULL) ;     /* out of memory */
00722     }
00723 
00724     /* ---------------------------------------------------------------------- */
00725     /* copy the triplet matrix */
00726     /* ---------------------------------------------------------------------- */
00727 
00728     Ci = C->i ;
00729     Cj = C->j ;
00730     Cx = C->x ;
00731     Cz = C->z ;
00732     C->nnz = nz ;
00733 
00734     for (k = 0 ; k < nz ; k++)
00735     {
00736   Ci [k] = Ti [k] ;
00737     }
00738     for (k = 0 ; k < nz ; k++)
00739     {
00740   Cj [k] = Tj [k] ;
00741     }
00742 
00743     if (xtype == CHOLMOD_REAL)
00744     {
00745   for (k = 0 ; k < nz ; k++)
00746   {
00747       Cx [k] = Tx [k] ;
00748   }
00749     }
00750     else if (xtype == CHOLMOD_COMPLEX)
00751     {
00752   for (k = 0 ; k < nz ; k++)
00753   {
00754       Cx [2*k  ] = Tx [2*k  ] ;
00755       Cx [2*k+1] = Tx [2*k+1] ;
00756   }
00757     }
00758     else if (xtype == CHOLMOD_ZOMPLEX)
00759     {
00760   for (k = 0 ; k < nz ; k++)
00761   {
00762       Cx [k] = Tx [k] ;
00763       Cz [k] = Tz [k] ;
00764   }
00765     }
00766 
00767     /* ---------------------------------------------------------------------- */
00768     /* return the result */
00769     /* ---------------------------------------------------------------------- */
00770 
00771     ASSERT (CHOLMOD(dump_triplet) (C, "C triplet copy", Common)) ;
00772     return (C) ;
00773 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines