Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_l_transpose.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === Core/cholmod_transpose =============================================== */
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 to
00015  * compute the transpose or permuted transpose of a matrix:
00016  *
00017  * Primary routines:
00018  * -----------------
00019  * cholmod_transpose    transpose sparse matrix
00020  * cholmod_ptranspose   transpose and permute sparse matrix
00021  * cholmod_sort     sort row indices in each column of sparse matrix
00022  *
00023  * Secondary routines:
00024  * -------------------
00025  * cholmod_transpose_unsym  transpose unsymmetric sparse matrix
00026  * cholmod_transpose_sym  transpose symmetric sparse matrix
00027  *
00028  * All xtypes (pattern, real, complex, and zomplex) are supported.
00029  *
00030  * ---------------------------------------
00031  * Unsymmetric case: A->stype is zero.
00032  * ---------------------------------------
00033  *
00034  * Computes F = A', F = A (:,f)' or F = A (p,f)', except that the indexing by
00035  * f does not work the same as the MATLAB notation (see below).  A->stype
00036  * is zero, which denotes that both the upper and lower triangular parts of
00037  * A are present (and used).  A may in fact be symmetric in pattern and/or
00038  * value; A->stype just denotes which part of A are stored.  A may be
00039  * rectangular.
00040  *
00041  * p is a permutation of 0:m-1, and f is a subset of 0:n-1, where A is m-by-n.
00042  * There can be no duplicate entries in p or f.
00043  *
00044  * The set f is held in fset and fsize.
00045  *  fset = NULL means ":" in MATLAB. fsize is ignored.
00046  *  fset != NULL means f = fset [0..fsize-1].
00047  *  fset != NULL and fsize = 0 means f is the empty set.
00048  *
00049  * Columns not in the set f are considered to be zero.  That is,
00050  * if A is 5-by-10 then F = A (:,[3 4])' is not 2-by-5, but 10-by-5, and rows
00051  * 3 and 4 of F are equal to columns 3 and 4 of A (the other rows of F are
00052  * zero).  More precisely, in MATLAB notation:
00053  *
00054  *  [m n] = size (A) ;
00055  *  F = A ;
00056  *  notf = ones (1,n) ;
00057  *  notf (f) = 0 ;
00058  *  F (:, find (notf)) = 0
00059  *  F = F'
00060  *
00061  * If you want the MATLAB equivalent F=A(p,f) operation, use cholmod_submatrix
00062  * instead (which does not compute the transpose).
00063  *
00064  * F->nzmax must be large enough to hold the matrix F.  It is not modified.
00065  * If F->nz is present then F->nz [j] = # of entries in column j of F.
00066  *
00067  * A can be sorted or unsorted, with packed or unpacked columns.
00068  *
00069  * If f is present and not sorted in ascending order, then F is unsorted
00070  * (that is, it may contain columns whose row indices do not appear in
00071  * ascending order).  Otherwise, F is sorted (the row indices in each
00072  * column of F appear in strictly ascending order).
00073  *
00074  * F is returned in packed or unpacked form, depending on F->packed on input.
00075  * If F->packed is false, then F is returned in unpacked form (F->nz must be
00076  * present).  Each row i of F is large enough to hold all the entries in row i
00077  * of A, even if f is provided.  That is, F->i and
00078  * F->x [F->p [i] .. F->p [i] + F->nz [i] - 1] contain all entries in A (i,f),
00079  * but F->p [i+1] - F->p [i] is equal to the number of nonzeros in A (i,:),
00080  * not just A (i,f).
00081  *
00082  * The cholmod_transpose_unsym routine is the only operation in CHOLMOD that
00083  * can produce an unpacked matrix.
00084  *
00085  * ---------------------------------------
00086  * Symmetric case: A->stype is nonzero.
00087  * ---------------------------------------
00088  *
00089  * Computes F = A' or F = A(p,p)', the transpose or permuted transpose, where
00090  * A->stype is nonzero.
00091  *
00092  * If A->stype > 0, then A is a symmetric matrix where just the upper part
00093  * of the matrix is stored.  Entries in the lower triangular part may be
00094  * present, but are ignored.  A must be square.  If F=A', then F is returned
00095  * sorted; otherwise F is unsorted for the F=A(p,p)' case.
00096  *
00097  * There can be no duplicate entries in p.
00098  * The fset and fsize parameters are not used.
00099  *
00100  * Three kinds of transposes are available, depending on the "values" parameter:
00101  * 0: do not transpose the numerical values; create a CHOLMOD_PATTERN matrix
00102  * 1: array transpose
00103  * 2: complex conjugate transpose (same as 2 if input is real or pattern)
00104  *
00105  * -----------------------------------------------------------------------------
00106  *
00107  * For cholmod_transpose_unsym and cholmod_transpose_sym, the output matrix
00108  * F must already be pre-allocated by the caller, with the correct dimensions.
00109  * If F is not valid or has the wrong dimensions, it is not modified.
00110  * Otherwise, if F is too small, the transpose is not computed; the contents
00111  * of F->p contain the column pointers of the resulting matrix, where
00112  * F->p [F->ncol] > F->nzmax.  In this case, the remaining contents of F are
00113  * not modified.  F can still be properly free'd with cholmod_free_sparse.
00114  */
00115 
00116 /* This file should make the long int version of CHOLMOD */
00117 #define DLONG 1
00118 
00119 #include "amesos_cholmod_internal.h"
00120 #include "amesos_cholmod_core.h"
00121 
00122 
00123 /* ========================================================================== */
00124 /* === TEMPLATE ============================================================= */
00125 /* ========================================================================== */
00126 
00127 #define PATTERN
00128 #include "amesos_t_cholmod_transpose.c"
00129 #define REAL
00130 #include "amesos_t_cholmod_transpose.c"
00131 #define COMPLEX
00132 #include "amesos_t_cholmod_transpose.c"
00133 #define COMPLEX
00134 #define NCONJUGATE
00135 #include "amesos_t_cholmod_transpose.c"
00136 #define ZOMPLEX
00137 #include "amesos_t_cholmod_transpose.c"
00138 #define ZOMPLEX
00139 #define NCONJUGATE
00140 #include "amesos_t_cholmod_transpose.c"
00141 
00142 
00143 /* ========================================================================== */
00144 /* === cholmod_transpose_unsym ============================================== */
00145 /* ========================================================================== */
00146 
00147 /* Compute F = A', A (:,f)', or A (p,f)', where A is unsymmetric and F is
00148  * already allocated.  See cholmod_transpose for a simpler routine.
00149  *
00150  * workspace:
00151  * Iwork (MAX (nrow,ncol)) if fset is present
00152  * Iwork (nrow) if fset is NULL
00153  *
00154  * The xtype of A and F must match, unless values is zero or F->xtype is
00155  * CHOLMOD_PATTERN (in which case only the pattern of A is transpose into F).
00156  */
00157 
00158 int CHOLMOD(transpose_unsym)
00159 (
00160     /* ---- input ---- */
00161     cholmod_sparse *A,  /* matrix to transpose */
00162     int values,   /* 2: complex conj. transpose, 1: array transpose,
00163          0: do not transpose the numerical values */
00164     Int *Perm,    /* size nrow, if present (can be NULL) */
00165     Int *fset,    /* subset of 0:(A->ncol)-1 */
00166     size_t fsize, /* size of fset */
00167     /* ---- output --- */
00168     cholmod_sparse *F,  /* F = A', A(:,f)', or A(p,f)' */
00169     /* --------------- */
00170     cholmod_common *Common
00171 )
00172 {
00173     Int *Fp, *Fnz, *Ap, *Ai, *Anz, *Wi ;
00174     Int nrow, ncol, permute, use_fset, Apacked, Fpacked, p, pend,
00175   i, j, k, Fsorted, nf, jj, jlast ;
00176     size_t s ;
00177     int ok = TRUE ;
00178 
00179     /* ---------------------------------------------------------------------- */
00180     /* check inputs */
00181     /* ---------------------------------------------------------------------- */
00182 
00183     RETURN_IF_NULL_COMMON (FALSE) ;
00184     RETURN_IF_NULL (A, FALSE) ;
00185     RETURN_IF_NULL (F, FALSE) ;
00186     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00187     RETURN_IF_XTYPE_INVALID (F, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00188     if (A->nrow != F->ncol || A->ncol != F->nrow)
00189     {
00190   ERROR (CHOLMOD_INVALID, "F has the wrong dimensions") ;
00191   return (FALSE) ;
00192     }
00193     Common->status = CHOLMOD_OK ;
00194 
00195     /* ---------------------------------------------------------------------- */
00196     /* get inputs */
00197     /* ---------------------------------------------------------------------- */
00198 
00199     nf = fsize ;
00200     use_fset = (fset != NULL) ;
00201     nrow = A->nrow ;
00202     ncol = A->ncol ;
00203 
00204     Ap = A->p ;   /* size A->ncol+1, column pointers of A */
00205     Ai = A->i ;   /* size nz = Ap [A->ncol], row indices of A */
00206     Anz = A->nz ;
00207     Apacked = A->packed ;
00208     ASSERT (IMPLIES (!Apacked, Anz != NULL)) ;
00209 
00210     permute = (Perm != NULL) ;
00211 
00212     Fp = F->p ;   /* size A->nrow+1, row pointers of F */
00213     Fnz = F->nz ;
00214     Fpacked = F->packed ;
00215     ASSERT (IMPLIES (!Fpacked, Fnz != NULL)) ;
00216 
00217     nf = (use_fset) ? nf : ncol ;
00218 
00219     /* ---------------------------------------------------------------------- */
00220     /* allocate workspace */
00221     /* ---------------------------------------------------------------------- */
00222 
00223     /* s = nrow + ((fset != NULL) ? ncol : 0) */
00224     s = CHOLMOD(add_size_t) (nrow, ((fset != NULL) ? ncol : 0), &ok) ;
00225     if (!ok)
00226     {
00227   ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
00228   return (FALSE) ;
00229     }
00230 
00231     CHOLMOD(allocate_work) (0, s, 0, Common) ;
00232     if (Common->status < CHOLMOD_OK)
00233     {
00234   return (FALSE) ;  /* out of memory */
00235     }
00236 
00237     Wi = Common->Iwork ;  /* size nrow (i/l/l) */
00238 
00239     /* ---------------------------------------------------------------------- */
00240     /* check Perm and fset */
00241     /* ---------------------------------------------------------------------- */
00242 
00243     if (permute)
00244     {
00245   for (i = 0 ; i < nrow ; i++)
00246   {
00247       Wi [i] = 1 ;
00248   }
00249   for (k = 0 ; k < nrow ; k++)
00250   {
00251       i = Perm [k] ;
00252       if (i < 0 || i > nrow || Wi [i] == 0)
00253       {
00254     ERROR (CHOLMOD_INVALID, "invalid permutation") ;
00255     return (FALSE) ;
00256       }
00257       Wi [i] = 0 ;
00258   }
00259     }
00260 
00261     if (use_fset)
00262     {
00263   for (j = 0 ; j < ncol ; j++)
00264   {
00265       Wi [j] = 1 ;
00266   }
00267   for (k = 0 ; k < nf ; k++)
00268   {
00269       j = fset [k] ;
00270       if (j < 0 || j > ncol || Wi [j] == 0)
00271       {
00272     ERROR (CHOLMOD_INVALID, "invalid fset") ;
00273     return (FALSE) ;
00274       }
00275       Wi [j] = 0 ;
00276   }
00277     }
00278 
00279     /* Perm and fset are now valid */
00280     ASSERT (CHOLMOD(dump_perm) (Perm, nrow, nrow, "Perm", Common)) ;
00281     ASSERT (CHOLMOD(dump_perm) (fset, nf, ncol, "fset", Common)) ;
00282 
00283     /* ---------------------------------------------------------------------- */
00284     /* count the entries in each row of A or A(:,f) */
00285     /* ---------------------------------------------------------------------- */
00286 
00287     for (i = 0 ; i < nrow ; i++)
00288     {
00289   Wi [i] = 0 ;
00290     }
00291 
00292     jlast = EMPTY ;
00293     Fsorted = TRUE ;
00294 
00295     if (use_fset)
00296     {
00297   /* count entries in each row of A(:,f) */
00298   for (jj = 0 ; jj < nf ; jj++)
00299   {
00300       j = fset [jj] ;
00301       if (j <= jlast)
00302       {
00303     Fsorted = FALSE ;
00304       }
00305       p = Ap [j] ;
00306       pend = (Apacked) ? (Ap [j+1]) : (p + Anz [j]) ;
00307       for ( ; p < pend ; p++)
00308       {
00309     Wi [Ai [p]]++ ;
00310       }
00311       jlast = j ;
00312   }
00313 
00314   /* save the nz counts if F is unpacked, and recount all of A */
00315   if (!Fpacked)
00316   {
00317       if (permute)
00318       {
00319     for (i = 0 ; i < nrow ; i++)
00320     {
00321         Fnz [i] = Wi [Perm [i]] ;
00322     }
00323       }
00324       else
00325       {
00326     for (i = 0 ; i < nrow ; i++)
00327     {
00328         Fnz [i] = Wi [i] ;
00329     }
00330       }
00331       for (i = 0 ; i < nrow ; i++)
00332       {
00333     Wi [i] = 0 ;
00334       }
00335 
00336       /* count entries in each row of A */
00337       for (j = 0 ; j < ncol ; j++)
00338       {
00339     p = Ap [j] ;
00340     pend = (Apacked) ? (Ap [j+1]) : (p + Anz [j]) ;
00341     for ( ; p < pend ; p++)
00342     {
00343         Wi [Ai [p]]++ ;
00344     }
00345       }
00346   }
00347 
00348     }
00349     else
00350     {
00351 
00352   /* count entries in each row of A */
00353   for (j = 0 ; j < ncol ; j++)
00354   {
00355       p = Ap [j] ;
00356       pend = (Apacked) ? (Ap [j+1]) : (p + Anz [j]) ;
00357       for ( ; p < pend ; p++)
00358       {
00359     Wi [Ai [p]]++ ;
00360       }
00361   }
00362 
00363   /* save the nz counts if F is unpacked */
00364   if (!Fpacked)
00365   {
00366       if (permute)
00367       {
00368     for (i = 0 ; i < nrow ; i++)
00369     {
00370         Fnz [i] = Wi [Perm [i]] ;
00371     }
00372       }
00373       else
00374       {
00375     for (i = 0 ; i < nrow ; i++)
00376     {
00377         Fnz [i] = Wi [i] ;
00378     }
00379       }
00380   }
00381     }
00382 
00383     /* ---------------------------------------------------------------------- */
00384     /* compute the row pointers */
00385     /* ---------------------------------------------------------------------- */
00386 
00387     p = 0 ;
00388     if (permute)
00389     {
00390   for (i = 0 ; i < nrow ; i++)
00391   {
00392       Fp [i] = p ;
00393       p += Wi [Perm [i]] ;
00394   }
00395   for (i = 0 ; i < nrow ; i++)
00396   {
00397       Wi [Perm [i]] = Fp [i] ;
00398   }
00399     }
00400     else
00401     {
00402   for (i = 0 ; i < nrow ; i++)
00403   {
00404       Fp [i] = p ;
00405       p += Wi [i] ;
00406   }
00407   for (i = 0 ; i < nrow ; i++)
00408   {
00409       Wi [i] = Fp [i] ;
00410   }
00411     }
00412     Fp [nrow] = p ;
00413 
00414     if (p > (Int) (F->nzmax))
00415     {
00416   ERROR (CHOLMOD_INVALID, "F is too small") ;
00417   return (FALSE) ;
00418     }
00419 
00420     /* ---------------------------------------------------------------------- */
00421     /* transpose matrix, using template routine */
00422     /* ---------------------------------------------------------------------- */
00423 
00424     ok = FALSE ;
00425     if (values == 0 || F->xtype == CHOLMOD_PATTERN)
00426     {
00427   ok = amesos_p_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
00428     }
00429     else if (F->xtype == CHOLMOD_REAL)
00430     {
00431   ok = amesos_r_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
00432     }
00433     else if (F->xtype == CHOLMOD_COMPLEX)
00434     {
00435   if (values == 1)
00436   {
00437       /* array transpose */
00438       ok = amesos_ct_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
00439   }
00440   else
00441   {
00442       /* complex conjugate transpose */
00443       ok = amesos_c_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
00444   }
00445     }
00446     else if (F->xtype == CHOLMOD_ZOMPLEX)
00447     {
00448   if (values == 1)
00449   {
00450       /* array transpose */
00451       ok = amesos_zt_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
00452   }
00453   else
00454   {
00455       /* complex conjugate transpose */
00456       ok = amesos_z_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
00457   }
00458     }
00459 
00460     /* ---------------------------------------------------------------------- */
00461     /* finalize result F */
00462     /* ---------------------------------------------------------------------- */
00463 
00464     if (ok)
00465     {
00466   F->sorted = Fsorted ;
00467     }
00468     ASSERT (CHOLMOD(dump_sparse) (F, "output F unsym", Common) >= 0) ;
00469     return (ok) ;
00470 }
00471 
00472 
00473 /* ========================================================================== */
00474 /* === cholmod_transpose_sym ================================================ */
00475 /* ========================================================================== */
00476 
00477 /* Compute F = A' or A (p,p)', where A is symmetric and F is already allocated.
00478  * See cholmod_transpose for a simpler routine.
00479  *
00480  * workspace:  Iwork (nrow) if Perm NULL, Iwork (2*nrow) if Perm non-NULL.
00481  */
00482 
00483 int CHOLMOD(transpose_sym)
00484 (
00485     /* ---- input ---- */
00486     cholmod_sparse *A,  /* matrix to transpose */
00487     int values,   /* 2: complex conj. transpose, 1: array transpose,
00488          0: do not transpose the numerical values */
00489     Int *Perm,    /* size nrow, if present (can be NULL) */
00490     /* ---- output --- */
00491     cholmod_sparse *F,  /* F = A' or A(p,p)' */
00492     /* --------------- */
00493     cholmod_common *Common
00494 )
00495 {
00496     Int *Ap, *Anz, *Ai, *Fp, *Wi, *Pinv, *Iwork ;
00497     Int p, pend, packed, upper, permute, jold, n, i, j, k, iold ;
00498     size_t s ;
00499     int ok = TRUE ;
00500 
00501     /* ---------------------------------------------------------------------- */
00502     /* check inputs */
00503     /* ---------------------------------------------------------------------- */
00504 
00505     RETURN_IF_NULL_COMMON (FALSE) ;
00506     RETURN_IF_NULL (A, FALSE) ;
00507     RETURN_IF_NULL (F, FALSE) ;
00508     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00509     RETURN_IF_XTYPE_INVALID (F, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00510     if (A->nrow != A->ncol || A->stype == 0)
00511     {
00512   /* this routine handles square symmetric matrices only */
00513   ERROR (CHOLMOD_INVALID, "matrix must be symmetric") ;
00514   return (FALSE) ;
00515     }
00516     if (A->nrow != F->ncol || A->ncol != F->nrow)
00517     {
00518   ERROR (CHOLMOD_INVALID, "F has the wrong dimensions") ;
00519   return (FALSE) ;
00520     }
00521     Common->status = CHOLMOD_OK ;
00522 
00523     /* ---------------------------------------------------------------------- */
00524     /* get inputs */
00525     /* ---------------------------------------------------------------------- */
00526 
00527     permute = (Perm != NULL) ;
00528     n = A->nrow ;
00529     Ap = A->p ;   /* size A->ncol+1, column pointers of A */
00530     Ai = A->i ;   /* size nz = Ap [A->ncol], row indices of A */
00531     Anz = A->nz ;
00532     packed = A->packed ;
00533     ASSERT (IMPLIES (!packed, Anz != NULL)) ;
00534     upper = (A->stype > 0) ;
00535 
00536     Fp = F->p ;   /* size A->nrow+1, row pointers of F */
00537 
00538     /* ---------------------------------------------------------------------- */
00539     /* allocate workspace */
00540     /* ---------------------------------------------------------------------- */
00541 
00542     /* s = (Perm != NULL) ? 2*n : n */
00543     s = CHOLMOD(add_size_t) (n, ((Perm != NULL) ? n : 0), &ok) ;
00544     if (!ok)
00545     {
00546   ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
00547   return (FALSE) ;
00548     }
00549 
00550     CHOLMOD(allocate_work) (0, s, 0, Common) ;
00551     if (Common->status < CHOLMOD_OK)
00552     {
00553   return (FALSE) ;  /* out of memory */
00554     }
00555 
00556     /* ---------------------------------------------------------------------- */
00557     /* get workspace */
00558     /* ---------------------------------------------------------------------- */
00559 
00560     Iwork = Common->Iwork ;
00561     Wi   = Iwork ;      /* size n (i/l/l) */
00562     Pinv = Iwork + n ;      /* size n (i/i/l) , unused if Perm NULL */
00563 
00564     /* ---------------------------------------------------------------------- */
00565     /* check Perm and construct inverse permutation */
00566     /* ---------------------------------------------------------------------- */
00567 
00568     if (permute)
00569     {
00570   for (i = 0 ; i < n ; i++)
00571   {
00572       Pinv [i] = EMPTY ;
00573   }
00574   for (k = 0 ; k < n ; k++)
00575   {
00576       i = Perm [k] ;
00577       if (i < 0 || i > n || Pinv [i] != EMPTY)
00578       {
00579     ERROR (CHOLMOD_INVALID, "invalid permutation") ;
00580     return (FALSE) ;
00581       }
00582       Pinv [i] = k ;
00583   }
00584     }
00585 
00586     /* Perm is now valid */
00587     ASSERT (CHOLMOD(dump_perm) (Perm, n, n, "Perm", Common)) ;
00588 
00589     /* ---------------------------------------------------------------------- */
00590     /* count the entries in each row of F */
00591     /* ---------------------------------------------------------------------- */
00592 
00593     for (i = 0 ; i < n ; i++)
00594     {
00595   Wi [i] = 0 ;
00596     }
00597 
00598     if (packed)
00599     {
00600   if (permute)
00601   {
00602       if (upper)
00603       {
00604     /* packed, permuted, upper */
00605     for (j = 0 ; j < n ; j++)
00606     {
00607         jold = Perm [j] ;
00608         pend = Ap [jold+1] ;
00609         for (p = Ap [jold] ; p < pend ; p++)
00610         {
00611       iold = Ai [p] ;
00612       if (iold <= jold)
00613       {
00614           i = Pinv [iold] ;
00615           Wi [MIN (i, j)]++ ;
00616       }
00617         }
00618     }
00619       }
00620       else
00621       {
00622     /* packed, permuted, lower */
00623     for (j = 0 ; j < n ; j++)
00624     {
00625         jold = Perm [j] ;
00626         pend = Ap [jold+1] ;
00627         for (p = Ap [jold] ; p < pend ; p++)
00628         {
00629       iold = Ai [p] ;
00630       if (iold >= jold)
00631       {
00632           i = Pinv [iold] ;
00633           Wi [MAX (i, j)]++ ;
00634       }
00635         }
00636     }
00637       }
00638   }
00639   else
00640   {
00641       if (upper)
00642       {
00643     /* packed, unpermuted, upper */
00644     for (j = 0 ; j < n ; j++)
00645     {
00646         pend = Ap [j+1] ;
00647         for (p = Ap [j] ; p < pend ; p++)
00648         {
00649       i = Ai [p] ;
00650       if (i <= j)
00651       {
00652           Wi [i]++ ;
00653       }
00654         }
00655     }
00656       }
00657       else
00658       {
00659     /* packed, unpermuted, lower */
00660     for (j = 0 ; j < n ; j++)
00661     {
00662         pend = Ap [j+1] ;
00663         for (p = Ap [j] ; p < pend ; p++)
00664         {
00665       i = Ai [p] ;
00666       if (i >= j)
00667       {
00668           Wi [i]++ ;
00669       }
00670         }
00671     }
00672       }
00673   }
00674     }
00675     else
00676     {
00677   if (permute)
00678   {
00679       if (upper)
00680       {
00681     /* unpacked, permuted, upper */
00682     for (j = 0 ; j < n ; j++)
00683     {
00684         jold = Perm [j] ;
00685         p = Ap [jold] ;
00686         pend = p + Anz [jold] ;
00687         for ( ; p < pend ; p++)
00688         {
00689       iold = Ai [p] ;
00690       if (iold <= jold)
00691       {
00692           i = Pinv [iold] ;
00693           Wi [MIN (i, j)]++ ;
00694       }
00695         }
00696     }
00697       }
00698       else
00699       {
00700     /* unpacked, permuted, lower */
00701     for (j = 0 ; j < n ; j++)
00702     {
00703         jold = Perm [j] ;
00704         p = Ap [jold] ;
00705         pend = p + Anz [jold] ;
00706         for ( ; p < pend ; p++)
00707         {
00708       iold = Ai [p] ;
00709       if (iold >= jold)
00710       {
00711           i = Pinv [iold] ;
00712           Wi [MAX (i, j)]++ ;
00713       }
00714         }
00715     }
00716       }
00717   }
00718   else
00719   {
00720       if (upper)
00721       {
00722     /* unpacked, unpermuted, upper */
00723     for (j = 0 ; j < n ; j++)
00724     {
00725         p = Ap [j] ;
00726         pend = p + Anz [j] ;
00727         for ( ; p < pend ; p++)
00728         {
00729       i = Ai [p] ;
00730       if (i <= j)
00731       {
00732           Wi [i]++ ;
00733       }
00734         }
00735     }
00736       }
00737       else
00738       {
00739     /* unpacked, unpermuted, lower */
00740     for (j = 0 ; j < n ; j++)
00741     {
00742         p = Ap [j] ;
00743         pend = p + Anz [j] ;
00744         for ( ; p < pend ; p++)
00745         {
00746       i = Ai [p] ;
00747       if (i >= j)
00748       {
00749           Wi [i]++ ;
00750       }
00751         }
00752     }
00753       }
00754   }
00755     }
00756 
00757     /* ---------------------------------------------------------------------- */
00758     /* compute the row pointers */
00759     /* ---------------------------------------------------------------------- */
00760 
00761     p = 0 ;
00762     for (i = 0 ; i < n ; i++)
00763     {
00764   Fp [i] = p ;
00765   p += Wi [i] ;
00766     }
00767     Fp [n] = p ;
00768     for (i = 0 ; i < n ; i++)
00769     {
00770   Wi [i] = Fp [i] ;
00771     }
00772 
00773     if (p > (Int) (F->nzmax))
00774     {
00775   ERROR (CHOLMOD_INVALID, "F is too small") ;
00776   return (FALSE) ;
00777     }
00778 
00779     /* ---------------------------------------------------------------------- */
00780     /* transpose matrix, using template routine */
00781     /* ---------------------------------------------------------------------- */
00782 
00783     ok = FALSE ;
00784     if (values == 0 || F->xtype == CHOLMOD_PATTERN)
00785     {
00786   PRINT2 (("\n:::: p_transpose_sym Perm %p\n", Perm)) ;
00787   ok = amesos_p_cholmod_transpose_sym (A, Perm, F, Common) ;
00788     }
00789     else if (F->xtype == CHOLMOD_REAL)
00790     {
00791   PRINT2 (("\n:::: r_transpose_sym Perm %p\n", Perm)) ;
00792   ok = amesos_r_cholmod_transpose_sym (A, Perm, F, Common) ;
00793     }
00794     else if (F->xtype == CHOLMOD_COMPLEX)
00795     {
00796   if (values == 1)
00797   {
00798       /* array transpose */
00799       PRINT2 (("\n:::: ct_transpose_sym Perm %p\n", Perm)) ;
00800       ok = amesos_ct_cholmod_transpose_sym (A, Perm, F, Common) ;
00801   }
00802   else
00803   {
00804       /* complex conjugate transpose */
00805       PRINT2 (("\n:::: c_transpose_sym Perm %p\n", Perm)) ;
00806       ok = amesos_c_cholmod_transpose_sym (A, Perm, F, Common) ;
00807   }
00808     }
00809     else if (F->xtype == CHOLMOD_ZOMPLEX)
00810     {
00811   if (values == 1)
00812   {
00813       /* array transpose */
00814       PRINT2 (("\n:::: zt_transpose_sym Perm %p\n", Perm)) ;
00815       ok = amesos_zt_cholmod_transpose_sym (A, Perm, F, Common) ;
00816   }
00817   else
00818   {
00819       /* complex conjugate transpose */
00820       PRINT2 (("\n:::: z_transpose_sym Perm %p\n", Perm)) ;
00821       ok = amesos_z_cholmod_transpose_sym (A, Perm, F, Common) ;
00822   }
00823     }
00824 
00825     /* ---------------------------------------------------------------------- */
00826     /* finalize result F */
00827     /* ---------------------------------------------------------------------- */
00828 
00829     /* F is sorted if there is no permutation vector */
00830     if (ok)
00831     {
00832   F->sorted = !permute ;
00833   F->packed = TRUE ;
00834   F->stype = - SIGN (A->stype) ;  /* flip the stype */
00835   ASSERT (CHOLMOD(dump_sparse) (F, "output F sym", Common) >= 0) ;
00836     }
00837     return (ok) ;
00838 }
00839 
00840 
00841 /* ========================================================================== */
00842 /* === cholmod_transpose ==================================================== */
00843 /* ========================================================================== */
00844 
00845 /* Returns A'.  See also cholmod_ptranspose below. */
00846 
00847 cholmod_sparse *CHOLMOD(transpose)
00848 (
00849     /* ---- input ---- */
00850     cholmod_sparse *A,  /* matrix to transpose */
00851     int values,   /* 2: complex conj. transpose, 1: array transpose,
00852          0: do not transpose the numerical values
00853          (returns its result as CHOLMOD_PATTERN) */
00854     /* --------------- */
00855     cholmod_common *Common
00856 )
00857 {
00858     return (CHOLMOD(ptranspose) (A, values, NULL, NULL, 0, Common)) ;
00859 }
00860 
00861 
00862 /* ========================================================================== */
00863 /* === cholmod_ptranspose =================================================== */
00864 /* ========================================================================== */
00865 
00866 /* Return A' or A(p,p)' if A is symmetric.  Return A', A(:,f)', or A(p,f)' if
00867  * A is unsymmetric.
00868  *
00869  * workspace:
00870  * Iwork (MAX (nrow,ncol)) if unsymmetric and fset is non-NULL
00871  * Iwork (nrow) if unsymmetric and fset is NULL
00872  * Iwork (2*nrow) if symmetric and Perm is non-NULL.
00873  * Iwork (nrow) if symmetric and Perm is NULL.
00874  *
00875  * A simple worst-case upper bound on the workspace is nrow+ncol.
00876  */
00877 
00878 cholmod_sparse *CHOLMOD(ptranspose)
00879 (
00880     /* ---- input ---- */
00881     cholmod_sparse *A,  /* matrix to transpose */
00882     int values,   /* 2: complex conj. transpose, 1: array transpose,
00883          0: do not transpose the numerical values */
00884     Int *Perm,    /* if non-NULL, F = A(p,f) or A(p,p) */
00885     Int *fset,    /* subset of 0:(A->ncol)-1 */
00886     size_t fsize, /* size of fset */
00887     /* --------------- */
00888     cholmod_common *Common
00889 )
00890 {
00891     Int *Ap, *Anz ;
00892     cholmod_sparse *F ;
00893     Int nrow, ncol, use_fset, j, jj, fnz, packed, stype, nf, xtype ;
00894     size_t ineed ;
00895     int ok = TRUE ;
00896 
00897     nf = fsize ;
00898 
00899     /* ---------------------------------------------------------------------- */
00900     /* check inputs */
00901     /* ---------------------------------------------------------------------- */
00902 
00903     RETURN_IF_NULL_COMMON (NULL) ;
00904     RETURN_IF_NULL (A, FALSE) ;
00905     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, NULL) ;
00906     stype = A->stype ;
00907     Common->status = CHOLMOD_OK ;
00908 
00909     /* ---------------------------------------------------------------------- */
00910     /* allocate workspace */
00911     /* ---------------------------------------------------------------------- */
00912 
00913     nrow = A->nrow ;
00914     ncol = A->ncol ;
00915 
00916     if (stype != 0)
00917     {
00918   use_fset = FALSE ;
00919   if (Perm != NULL)
00920   {
00921       ineed = CHOLMOD(mult_size_t) (A->nrow, 2, &ok) ;
00922   }
00923   else
00924   {
00925       ineed = A->nrow ;
00926   }
00927     }
00928     else
00929     {
00930   use_fset = (fset != NULL) ;
00931   if (use_fset)
00932   {
00933       ineed = MAX (A->nrow, A->ncol) ;
00934   }
00935   else
00936   {
00937       ineed = A->nrow ;
00938   }
00939     }
00940 
00941     if (!ok)
00942     {
00943   ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
00944   return (NULL) ;
00945     }
00946 
00947     CHOLMOD(allocate_work) (0, ineed, 0, Common) ;
00948     if (Common->status < CHOLMOD_OK)
00949     {
00950   return (NULL) ;     /* out of memory */
00951     }
00952 
00953     /* ---------------------------------------------------------------------- */
00954     /* get inputs */
00955     /* ---------------------------------------------------------------------- */
00956 
00957     Ap = A->p ;
00958     Anz = A->nz ;
00959     packed = A->packed ;
00960     ASSERT (IMPLIES (!packed, Anz != NULL)) ;
00961     xtype = values ? A->xtype : CHOLMOD_PATTERN ;
00962 
00963     /* ---------------------------------------------------------------------- */
00964     /* allocate F */
00965     /* ---------------------------------------------------------------------- */
00966 
00967     /* determine # of nonzeros in F */
00968     if (stype != 0)
00969     {
00970   /* F=A' or F=A(p,p)', fset is ignored */
00971   fnz = CHOLMOD(nnz) (A, Common) ;
00972     }
00973     else
00974     {
00975   nf = (use_fset) ? nf : ncol ;
00976   if (use_fset)
00977   {
00978       fnz = 0 ;
00979       /* F=A(:,f)' or F=A(p,f)' */
00980       for (jj = 0 ; jj < nf ; jj++)
00981       {
00982     /* The fset is not yet checked; it will be thoroughly checked
00983      * in cholmod_transpose_unsym.  For now, just make sure we don't
00984      * access Ap and Anz out of bounds. */
00985     j = fset [jj] ;
00986     if (j >= 0 && j < ncol)
00987     {
00988         fnz += packed ? (Ap [j+1] - Ap [j]) : MAX (0, Anz [j]) ;
00989     }
00990       }
00991   }
00992   else
00993   {
00994       /* F=A' or F=A(p,:)' */
00995       fnz = CHOLMOD(nnz) (A, Common) ;
00996   }
00997     }
00998 
00999     /* F is ncol-by-nrow, fnz nonzeros, sorted unless f is present and unsorted,
01000      * packed, of opposite stype as A, and with/without numerical values */
01001     F = CHOLMOD(allocate_sparse) (ncol, nrow, fnz, TRUE, TRUE, -SIGN(stype),
01002       xtype, Common) ;
01003     if (Common->status < CHOLMOD_OK)
01004     {
01005   return (NULL) ;     /* out of memory */
01006     }
01007 
01008     /* ---------------------------------------------------------------------- */
01009     /* transpose and optionally permute the matrix A */
01010     /* ---------------------------------------------------------------------- */
01011 
01012     if (stype != 0)
01013     {
01014   /* F = A (p,p)', using upper or lower triangular part of A only */
01015   ok = CHOLMOD(transpose_sym) (A, values, Perm, F, Common) ;
01016     }
01017     else
01018     {
01019   /* F = A (p,f)' */
01020   ok = CHOLMOD(transpose_unsym) (A, values, Perm, fset, nf, F, Common) ;
01021     }
01022 
01023     /* ---------------------------------------------------------------------- */
01024     /* return the matrix F, or NULL if an error occured */
01025     /* ---------------------------------------------------------------------- */
01026 
01027     if (!ok)
01028     {
01029   CHOLMOD(free_sparse) (&F, Common) ;
01030     }
01031     return (F) ;
01032 }
01033 
01034 
01035 /* ========================================================================== */
01036 /* === cholmod_sort ========================================================= */
01037 /* ========================================================================== */
01038 
01039 /* Sort the columns of A, in place.  Returns A in packed form, even if it
01040  * starts as unpacked.  Removes entries in the ignored part of a symmetric
01041  * matrix.
01042  *
01043  * workspace: Iwork (max (nrow,ncol)).  Allocates additional workspace for a
01044  * temporary copy of A'.
01045  */
01046 
01047 int CHOLMOD(sort)
01048 (
01049     /* ---- in/out --- */
01050     cholmod_sparse *A,  /* matrix to sort */
01051     /* --------------- */
01052     cholmod_common *Common
01053 )
01054 {
01055     Int *Ap ;
01056     cholmod_sparse *F ;
01057     Int anz, ncol, nrow, stype ;
01058 
01059     /* ---------------------------------------------------------------------- */
01060     /* check inputs */
01061     /* ---------------------------------------------------------------------- */
01062 
01063     RETURN_IF_NULL_COMMON (FALSE) ;
01064     RETURN_IF_NULL (A, FALSE) ;
01065     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
01066     Common->status = CHOLMOD_OK ;
01067     nrow = A->nrow ;
01068     if (nrow <= 1)
01069     {
01070   /* a 1-by-n sparse matrix must be sorted */
01071   A->sorted = TRUE ;
01072   return (TRUE) ;
01073     }
01074 
01075     /* ---------------------------------------------------------------------- */
01076     /* allocate workspace */
01077     /* ---------------------------------------------------------------------- */
01078 
01079     ncol = A->ncol ;
01080     CHOLMOD(allocate_work) (0, MAX (nrow, ncol), 0, Common) ;
01081     if (Common->status < CHOLMOD_OK)
01082     {
01083   return (FALSE) ;  /* out of memory */
01084     }
01085 
01086     /* ---------------------------------------------------------------------- */
01087     /* get inputs */
01088     /* ---------------------------------------------------------------------- */
01089 
01090     anz = CHOLMOD(nnz) (A, Common) ;
01091     stype = A->stype ;
01092 
01093     /* ---------------------------------------------------------------------- */
01094     /* sort the columns of the matrix */
01095     /* ---------------------------------------------------------------------- */
01096 
01097     /* allocate workspace for transpose: ncol-by-nrow, same # of nonzeros as A,
01098      * sorted, packed, same stype as A, and of the same numeric type as A. */
01099     F = CHOLMOD(allocate_sparse) (ncol, nrow, anz, TRUE, TRUE, stype,
01100       A->xtype, Common) ;
01101     if (Common->status < CHOLMOD_OK)
01102     {
01103   return (FALSE) ;  /* out of memory */
01104     }
01105 
01106     if (stype != 0)
01107     {
01108   /* F = A', upper or lower triangular part only */
01109   CHOLMOD(transpose_sym) (A, 1, NULL, F, Common) ;
01110   A->packed = TRUE ;
01111   /* A = F' */
01112   CHOLMOD(transpose_sym) (F, 1, NULL, A, Common) ;
01113     }
01114     else
01115     {
01116   /* F = A' */
01117   CHOLMOD(transpose_unsym) (A, 1, NULL, NULL, 0, F, Common) ;
01118   A->packed = TRUE ;
01119   /* A = F' */
01120   CHOLMOD(transpose_unsym) (F, 1, NULL, NULL, 0, A, Common) ;
01121     }
01122 
01123     ASSERT (A->sorted && A->packed) ;
01124     ASSERT (CHOLMOD(dump_sparse) (A, "Asorted", Common) >= 0) ;
01125 
01126     /* ---------------------------------------------------------------------- */
01127     /* reduce A in size, if needed.  This must succeed. */
01128     /* ---------------------------------------------------------------------- */
01129 
01130     Ap = A->p ;
01131     anz = Ap [ncol] ;
01132     ASSERT ((size_t) anz <= A->nzmax) ;
01133     CHOLMOD(reallocate_sparse) (anz, A, Common) ;
01134     ASSERT (Common->status >= CHOLMOD_OK) ;
01135 
01136     /* ---------------------------------------------------------------------- */
01137     /* free workspace */
01138     /* ---------------------------------------------------------------------- */
01139 
01140     CHOLMOD(free_sparse) (&F, Common) ;
01141     return (TRUE) ;
01142 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines