Amesos Package Browser (Single Doxygen Collection) Development
amesos_t_cholmod_transpose.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === Core/t_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 /* Template routine for cholmod_transpose.  All xtypes are supported.  For
00015  * complex matrices, either the array tranpose or complex conjugate transpose
00016  * can be computed. */
00017 
00018 #include "amesos_cholmod_template.h"
00019 
00020 /* ========================================================================== */
00021 /* === t_cholmod_transpose_unsym ============================================ */
00022 /* ========================================================================== */
00023 
00024 /* Compute F = A', A (:,f)', or A (p,f)', where A is unsymmetric and F is
00025  * already allocated.  The complex case performs either the array transpose
00026  * or complex conjugate transpose.
00027  *
00028  * workspace:
00029  * Iwork (MAX (nrow,ncol)) if fset is present
00030  * Iwork (nrow) if fset is NULL
00031  */
00032 
00033 static int TEMPLATE (cholmod_transpose_unsym)
00034 (
00035     /* ---- input ---- */
00036     cholmod_sparse *A,  /* matrix to transpose */
00037     Int *Perm,    /* size nrow, if present (can be NULL) */
00038     Int *fset,    /* subset of 0:(A->ncol)-1 */
00039     Int nf,   /* size of fset */
00040     /* ---- output --- */
00041     cholmod_sparse *F,  /* F = A', A(:,f)', or A(p,f)' */
00042     /* --------------- */
00043     cholmod_common *Common
00044 )
00045 {
00046     double *Ax, *Az, *Fx, *Fz ;
00047     Int *Ap, *Anz, *Ai, *Fp, *Fnz, *Fj, *Wi, *Iwork ;
00048     Int j, p, pend, nrow, ncol, Apacked, use_fset, fp, Fpacked, jj, permute ;
00049 
00050     /* ---------------------------------------------------------------------- */
00051     /* check inputs */
00052     /* ---------------------------------------------------------------------- */
00053 
00054     /* ensure the xtype of A and F match (ignored if this is pattern version) */
00055     if (!XTYPE_OK (A->xtype))
00056     {
00057   ERROR (CHOLMOD_INVALID, "real/complex mismatch") ;
00058   return (FALSE) ;
00059     }
00060 
00061     /* ---------------------------------------------------------------------- */
00062     /* get inputs */
00063     /* ---------------------------------------------------------------------- */
00064 
00065     use_fset = (fset != NULL) ;
00066     nrow = A->nrow ;
00067     ncol = A->ncol ;
00068 
00069     Ap = A->p ;   /* size A->ncol+1, column pointers of A */
00070     Ai = A->i ;   /* size nz = Ap [A->ncol], row indices of A */
00071     Ax = A->x ;   /* size nz, real values of A */
00072     Az = A->z ;   /* size nz, imag values of A */
00073     Anz = A->nz ;
00074     Apacked = A->packed ;
00075     ASSERT (IMPLIES (!Apacked, Anz != NULL)) ;
00076 
00077     permute = (Perm != NULL) ;
00078 
00079     Fp = F->p ;   /* size A->nrow+1, row pointers of F */
00080     Fj = F->i ;   /* size nz, column indices of F */
00081     Fx = F->x ;   /* size nz, real values of F */
00082     Fz = F->z ;   /* size nz, imag values of F */
00083     Fnz = F->nz ;
00084     Fpacked = F->packed ;
00085     ASSERT (IMPLIES (!Fpacked, Fnz != NULL)) ;
00086 
00087     nf = (use_fset) ? nf : ncol ;
00088 
00089     /* ---------------------------------------------------------------------- */
00090     /* get workspace */
00091     /* ---------------------------------------------------------------------- */
00092 
00093     Iwork = Common->Iwork ;
00094     Wi = Iwork ;    /* size nrow (i/l/l) */
00095 
00096     /* ---------------------------------------------------------------------- */
00097     /* construct the transpose */
00098     /* ---------------------------------------------------------------------- */
00099 
00100     for (jj = 0 ; jj < nf ; jj++)
00101     {
00102   j = (use_fset) ? (fset [jj]) : jj ;
00103   p = Ap [j] ;
00104   pend = (Apacked) ? (Ap [j+1]) : (p + Anz [j]) ;
00105   for ( ; p < pend ; p++)
00106   {
00107       fp = Wi [Ai [p]]++ ;
00108       Fj [fp] = j ;
00109 #ifdef NCONJUGATE
00110       ASSIGN (Fx, Fz, fp, Ax, Az, p) ;
00111 #else
00112       ASSIGN_CONJ (Fx, Fz, fp, Ax, Az, p) ;
00113 #endif
00114   }
00115     }
00116 
00117     return (TRUE) ;
00118 }
00119 
00120 
00121 /* ========================================================================== */
00122 /* === t_cholmod_transpose_sym ============================================== */
00123 /* ========================================================================== */
00124 
00125 /* Compute F = A' or A (p,p)', where A is symmetric and F is already allocated.
00126  * The complex case performs either the array transpose or complex conjugate
00127  * transpose.
00128  *
00129  * workspace:  Iwork (nrow) if Perm NULL, Iwork (2*nrow) if Perm non-NULL.
00130  */
00131 
00132 static int TEMPLATE (cholmod_transpose_sym)
00133 (
00134     /* ---- input ---- */
00135     cholmod_sparse *A,  /* matrix to transpose */
00136     Int *Perm,    /* size n, if present (can be NULL) */
00137     /* ---- output --- */
00138     cholmod_sparse *F,  /* F = A' or A(p,p)' */
00139     /* --------------- */
00140     cholmod_common *Common
00141 )
00142 {
00143     double *Ax, *Az, *Fx, *Fz ;
00144     Int *Ap, *Anz, *Ai, *Fp, *Fj, *Wi, *Pinv, *Iwork ;
00145     Int p, pend, packed, fp, upper, permute, jold, n, i, j, iold ;
00146 
00147     /* ---------------------------------------------------------------------- */
00148     /* check inputs */
00149     /* ---------------------------------------------------------------------- */
00150 
00151     /* ensure the xtype of A and F match (ignored if this is pattern version) */
00152     if (!XTYPE_OK (A->xtype))
00153     {
00154   ERROR (CHOLMOD_INVALID, "real/complex mismatch") ;
00155   return (FALSE) ;
00156     }
00157 
00158     /* ---------------------------------------------------------------------- */
00159     /* get inputs */
00160     /* ---------------------------------------------------------------------- */
00161 
00162     permute = (Perm != NULL) ;
00163     n = A->nrow ;
00164     Ap = A->p ;   /* size A->ncol+1, column pointers of A */
00165     Ai = A->i ;   /* size nz = Ap [A->ncol], row indices of A */
00166     Ax = A->x ;   /* size nz, real values of A */
00167     Az = A->z ;   /* size nz, imag values of A */
00168     Anz = A->nz ;
00169     packed = A->packed ;
00170     ASSERT (IMPLIES (!packed, Anz != NULL)) ;
00171     upper = (A->stype > 0) ;
00172 
00173     Fp = F->p ;   /* size A->nrow+1, row pointers of F */
00174     Fj = F->i ;   /* size nz, column indices of F */
00175     Fx = F->x ;   /* size nz, real values of F */
00176     Fz = F->z ;   /* size nz, imag values of F */
00177 
00178     /* ---------------------------------------------------------------------- */
00179     /* get workspace */
00180     /* ---------------------------------------------------------------------- */
00181 
00182     Iwork = Common->Iwork ;
00183     Wi = Iwork ;  /* size n (i/l/l) */
00184     Pinv = Iwork + n ;  /* size n (i/i/l) , unused if Perm NULL */
00185 
00186     /* ---------------------------------------------------------------------- */
00187     /* construct the transpose */
00188     /* ---------------------------------------------------------------------- */
00189 
00190     if (permute)
00191     {
00192   if (upper)
00193   {
00194       /* permuted, upper */
00195       for (j = 0 ; j < n ; j++)
00196       {
00197     jold = Perm [j] ;
00198     p = Ap [jold] ;
00199     pend = (packed) ? Ap [jold+1] : p + Anz [jold] ;
00200     for ( ; p < pend ; p++)
00201     {
00202         iold = Ai [p] ;
00203         if (iold <= jold)
00204         {
00205       i = Pinv [iold] ;
00206       if (i < j)
00207       {
00208           fp = Wi [i]++ ;
00209           Fj [fp] = j ;
00210 #ifdef NCONJUGATE
00211           ASSIGN (Fx, Fz, fp, Ax, Az, p) ;
00212 #else
00213           ASSIGN_CONJ (Fx, Fz, fp, Ax, Az, p) ;
00214 #endif
00215       }
00216       else
00217       {
00218           fp = Wi [j]++ ;
00219           Fj [fp] = i ;
00220           ASSIGN (Fx, Fz, fp, Ax, Az, p) ;
00221       }
00222         }
00223     }
00224       }
00225   }
00226   else
00227   {
00228       /* permuted, lower */
00229       for (j = 0 ; j < n ; j++)
00230       {
00231     jold = Perm [j] ;
00232     p = Ap [jold] ;
00233     pend = (packed) ? Ap [jold+1] : p + Anz [jold] ;
00234     for ( ; p < pend ; p++)
00235     {
00236         iold = Ai [p] ;
00237         if (iold >= jold)
00238         {
00239       i = Pinv [iold] ;
00240       if (i > j)
00241       {
00242           fp = Wi [i]++ ;
00243           Fj [fp] = j ;
00244 #ifdef NCONJUGATE
00245           ASSIGN (Fx, Fz, fp, Ax, Az, p) ;
00246 #else
00247           ASSIGN_CONJ (Fx, Fz, fp, Ax, Az, p) ;
00248 #endif
00249       }
00250       else
00251       {
00252           fp = Wi [j]++ ;
00253           Fj [fp] = i ;
00254           ASSIGN (Fx, Fz, fp, Ax, Az, p) ;
00255       }
00256         }
00257     }
00258       }
00259   }
00260     }
00261     else
00262     {
00263   if (upper)
00264   {
00265       /* unpermuted, upper */
00266       for (j = 0 ; j < n ; j++)
00267       {
00268     p = Ap [j] ;
00269     pend = (packed) ? Ap [j+1] : p + Anz [j] ;
00270     for ( ; p < pend ; p++)
00271     {
00272         i = Ai [p] ;
00273         if (i <= j)
00274         {
00275       fp = Wi [i]++ ;
00276       Fj [fp] = j ;
00277 #ifdef NCONJUGATE
00278       ASSIGN (Fx, Fz, fp, Ax, Az, p) ;
00279 #else
00280       ASSIGN_CONJ (Fx, Fz, fp, Ax, Az, p) ;
00281 #endif
00282         }
00283     }
00284       }
00285   }
00286   else
00287   {
00288       /* unpermuted, lower */
00289       for (j = 0 ; j < n ; j++)
00290       {
00291     p = Ap [j] ;
00292     pend = (packed) ? Ap [j+1] : p + Anz [j] ;
00293     for ( ; p < pend ; p++)
00294     {
00295         i = Ai [p] ;
00296         if (i >= j)
00297         {
00298       fp = Wi [i]++ ;
00299       Fj [fp] = j ;
00300 #ifdef NCONJUGATE
00301       ASSIGN (Fx, Fz, fp, Ax, Az, p) ;
00302 #else
00303       ASSIGN_CONJ (Fx, Fz, fp, Ax, Az, p) ;
00304 #endif
00305         }
00306     }
00307       }
00308   }
00309     }
00310 
00311     return (TRUE) ;
00312 }
00313 
00314 #undef PATTERN
00315 #undef REAL
00316 #undef COMPLEX
00317 #undef ZOMPLEX
00318 #undef NCONJUGATE
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines