Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_aat.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === Core/cholmod_aat ===================================================== */
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 /* C = A*A' or C = A(:,f)*A(:,f)'
00015  *
00016  * A can be packed or unpacked, sorted or unsorted, but must be stored with
00017  * both upper and lower parts (A->stype of zero).  C is returned as packed,
00018  * C->stype of zero (both upper and lower parts present), and unsorted.  See
00019  * cholmod_ssmult in the MatrixOps Module for a more general matrix-matrix
00020  * multiply.
00021  *
00022  * You can trivially convert C into a symmetric upper/lower matrix by
00023  * changing C->stype = 1 or -1 after calling this routine.
00024  *
00025  * workspace:
00026  *  Flag (A->nrow),
00027  *  Iwork (max (A->nrow, A->ncol)) if fset present,
00028  *  Iwork (A->nrow) if no fset,
00029  *  W (A->nrow) if mode > 0,
00030  *  allocates temporary copy for A'.
00031  *
00032  * A can be pattern or real.  Complex or zomplex cases are supported only
00033  * if the mode is <= 0 (in which case the numerical values are ignored).
00034  */
00035 
00036 #include "amesos_cholmod_internal.h"
00037 #include "amesos_cholmod_core.h"
00038 
00039 cholmod_sparse *CHOLMOD(aat)
00040 (
00041     /* ---- input ---- */
00042     cholmod_sparse *A,  /* input matrix; C=A*A' is constructed */
00043     Int *fset,    /* subset of 0:(A->ncol)-1 */
00044     size_t fsize, /* size of fset */
00045     int mode,   /* >0: numerical, 0: pattern, <0: pattern (no diag)
00046        * -2: pattern only, no diagonal, add 50% + n extra
00047        * space to C */
00048     /* --------------- */
00049     cholmod_common *Common
00050 )
00051 {
00052     double fjt ;
00053     double *Ax, *Fx, *Cx, *W ;
00054     Int *Ap, *Anz, *Ai, *Fp, *Fi, *Cp, *Ci, *Flag ;
00055     cholmod_sparse *C, *F ;
00056     Int packed, j, i, pa, paend, pf, pfend, n, mark, cnz, t, p, values, diag,
00057   extra ;
00058 
00059     /* ---------------------------------------------------------------------- */
00060     /* check inputs */
00061     /* ---------------------------------------------------------------------- */
00062 
00063     RETURN_IF_NULL_COMMON (NULL) ;
00064     RETURN_IF_NULL (A, NULL) ;
00065     values = (mode > 0) && (A->xtype != CHOLMOD_PATTERN) ;
00066     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN,
00067       values ? CHOLMOD_REAL : CHOLMOD_ZOMPLEX, NULL) ;
00068     if (A->stype)
00069     {
00070   ERROR (CHOLMOD_INVALID, "matrix cannot be symmetric") ;
00071   return (NULL) ;
00072     }
00073     Common->status = CHOLMOD_OK ;
00074 
00075     /* ---------------------------------------------------------------------- */
00076     /* allocate workspace */
00077     /* ---------------------------------------------------------------------- */
00078 
00079     diag = (mode >= 0) ;
00080     n = A->nrow ;
00081     CHOLMOD(allocate_work) (n, MAX (A->ncol, A->nrow), values ? n : 0, Common) ;
00082     if (Common->status < CHOLMOD_OK)
00083     {
00084   return (NULL) ;     /* out of memory */
00085     }
00086     ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, values ? n : 0, Common)) ;
00087 
00088     /* ---------------------------------------------------------------------- */
00089     /* get inputs */
00090     /* ---------------------------------------------------------------------- */
00091 
00092     ASSERT (CHOLMOD(dump_sparse) (A, "A", Common) >= 0) ;
00093 
00094     /* get the A matrix */
00095     Ap  = A->p ;
00096     Anz = A->nz ;
00097     Ai  = A->i ;
00098     Ax  = A->x ;
00099     packed = A->packed ;
00100 
00101     /* get workspace */
00102     W = Common->Xwork ;   /* size n, unused if values is FALSE */
00103     Flag = Common->Flag ; /* size n, Flag [0..n-1] < mark on input*/
00104 
00105     /* ---------------------------------------------------------------------- */
00106     /* F = A' or A(:,f)' */
00107     /* ---------------------------------------------------------------------- */
00108 
00109     /* workspace: Iwork (nrow if no fset; MAX (nrow,ncol) if fset)*/
00110     F = CHOLMOD(ptranspose) (A, values, NULL, fset, fsize, Common) ;
00111     if (Common->status < CHOLMOD_OK)
00112     {
00113   return (NULL) ;     /* out of memory */
00114     }
00115 
00116     Fp = F->p ;
00117     Fi = F->i ;
00118     Fx = F->x ;
00119 
00120     /* ---------------------------------------------------------------------- */
00121     /* count the number of entries in the result C */
00122     /* ---------------------------------------------------------------------- */
00123 
00124     cnz = 0 ;
00125     for (j = 0 ; j < n ; j++)
00126     {
00127   /* clear the Flag array */
00128   mark = CHOLMOD(clear_flag) (Common) ;
00129 
00130   /* exclude the diagonal, if requested */
00131   if (!diag)
00132   {
00133       Flag [j] = mark ;
00134   }
00135 
00136   /* for each nonzero F(t,j) in column j, do: */
00137   pfend = Fp [j+1] ;
00138   for (pf = Fp [j] ; pf < pfend ; pf++)
00139   {
00140       /* F(t,j) is nonzero */
00141       t = Fi [pf] ;
00142 
00143       /* add the nonzero pattern of A(:,t) to the pattern of C(:,j) */
00144       pa = Ap [t] ;
00145       paend = (packed) ? (Ap [t+1]) : (pa + Anz [t]) ;
00146       for ( ; pa < paend ; pa++)
00147       {
00148     i = Ai [pa] ;
00149     if (Flag [i] != mark)
00150     {
00151         Flag [i] = mark ;
00152         cnz++ ;
00153     }
00154       }
00155   }
00156   if (cnz < 0)
00157   {
00158       break ;     /* integer overflow case */
00159   }
00160     }
00161 
00162     extra = (mode == -2) ? (cnz/2 + n) : 0 ;
00163 
00164     mark = CHOLMOD(clear_flag) (Common) ;
00165 
00166     /* ---------------------------------------------------------------------- */
00167     /* check for integer overflow */
00168     /* ---------------------------------------------------------------------- */
00169 
00170     if (cnz < 0 || (cnz + extra) < 0)
00171     {
00172   ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
00173   CHOLMOD(clear_flag) (Common) ;
00174   CHOLMOD(free_sparse) (&F, Common) ;
00175   return (NULL) ;     /* problem too large */
00176     }
00177 
00178     /* ---------------------------------------------------------------------- */
00179     /* allocate C */
00180     /* ---------------------------------------------------------------------- */
00181 
00182     C = CHOLMOD(allocate_sparse) (n, n, cnz + extra, FALSE, TRUE, 0,
00183       values ? A->xtype : CHOLMOD_PATTERN, Common) ;
00184     if (Common->status < CHOLMOD_OK)
00185     {
00186   CHOLMOD(free_sparse) (&F, Common) ;
00187   return (NULL) ;     /* out of memory */
00188     }
00189 
00190     Cp = C->p ;
00191     Ci = C->i ;
00192     Cx = C->x ;
00193 
00194     /* ---------------------------------------------------------------------- */
00195     /* C = A*A' */
00196     /* ---------------------------------------------------------------------- */
00197 
00198     cnz = 0 ;
00199 
00200     if (values)
00201     {
00202 
00203   /* pattern and values */
00204   for (j = 0 ; j < n ; j++)
00205   {
00206       /* clear the Flag array */
00207       mark = CHOLMOD(clear_flag) (Common) ;
00208 
00209       /* start column j of C */
00210       Cp [j] = cnz ;
00211 
00212       /* for each nonzero F(t,j) in column j, do: */
00213       pfend = Fp [j+1] ;
00214       for (pf = Fp [j] ; pf < pfend ; pf++)
00215       {
00216     /* F(t,j) is nonzero */
00217     t = Fi [pf] ;
00218     fjt = Fx [pf] ;
00219 
00220     /* add the nonzero pattern of A(:,t) to the pattern of C(:,j)
00221      * and scatter the values into W */
00222     pa = Ap [t] ;
00223     paend = (packed) ? (Ap [t+1]) : (pa + Anz [t]) ;
00224     for ( ; pa < paend ; pa++)
00225     {
00226         i = Ai [pa] ;
00227         if (Flag [i] != mark)
00228         {
00229       Flag [i] = mark ;
00230       Ci [cnz++] = i ;
00231         }
00232         W [i] += Ax [pa] * fjt ;
00233     }
00234       }
00235 
00236       /* gather the values into C(:,j) */
00237       for (p = Cp [j] ; p < cnz ; p++)
00238       {
00239     i = Ci [p] ;
00240     Cx [p] = W [i] ;
00241     W [i] = 0 ;
00242       }
00243   }
00244 
00245     }
00246     else
00247     {
00248 
00249   /* pattern only */
00250   for (j = 0 ; j < n ; j++)
00251   {
00252       /* clear the Flag array */
00253       mark = CHOLMOD(clear_flag) (Common) ;
00254 
00255       /* exclude the diagonal, if requested */
00256       if (!diag)
00257       {
00258     Flag [j] = mark ;
00259       }
00260 
00261       /* start column j of C */
00262       Cp [j] = cnz ;
00263 
00264       /* for each nonzero F(t,j) in column j, do: */
00265       pfend = Fp [j+1] ;
00266       for (pf = Fp [j] ; pf < pfend ; pf++)
00267       {
00268     /* F(t,j) is nonzero */
00269     t = Fi [pf] ;
00270 
00271     /* add the nonzero pattern of A(:,t) to the pattern of C(:,j) */
00272     pa = Ap [t] ;
00273     paend = (packed) ? (Ap [t+1]) : (pa + Anz [t]) ;
00274     for ( ; pa < paend ; pa++)
00275     {
00276         i = Ai [pa] ;
00277         if (Flag [i] != mark)
00278         {
00279       Flag [i] = mark ;
00280       Ci [cnz++] = i ;
00281         }
00282     }
00283       }
00284   }
00285     }
00286 
00287     Cp [n] = cnz ;
00288     ASSERT (IMPLIES (mode != -2, MAX (1,cnz) == C->nzmax)) ;
00289 
00290     /* ---------------------------------------------------------------------- */
00291     /* clear workspace and free temporary matrices and return result */
00292     /* ---------------------------------------------------------------------- */
00293 
00294     CHOLMOD(free_sparse) (&F, Common) ;
00295     CHOLMOD(clear_flag) (Common) ;
00296     ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, values ? n : 0, Common)) ;
00297     DEBUG (i = CHOLMOD(dump_sparse) (C, "aat", Common)) ;
00298     ASSERT (IMPLIES (mode < 0, i == 0)) ;
00299     return (C) ;
00300 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines