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