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