Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_add.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === Core/cholmod_add ===================================================== */
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 = alpha*A + beta*B, or spones(A+B).  Result is packed, with sorted or
00015  * unsorted columns.  This routine is much faster and takes less memory if C
00016  * is allowed to have unsorted columns.
00017  *
00018  * If A and B are both symmetric (in upper form) then C is the same.  Likewise,
00019  * if A and B are both symmetric (in lower form) then C is the same.
00020  * Otherwise, C is unsymmetric.  A and B must have the same dimension.
00021  *
00022  * workspace: Flag (nrow), W (nrow) if values, Iwork (max (nrow,ncol)).
00023  *  allocates temporary copies for A and B if they are symmetric.
00024  *  allocates temporary copy of C if it is to be returned sorted.
00025  *
00026  * A and B can have an xtype of pattern or real.  Complex or zomplex cases
00027  * are supported only if the "values" input parameter is FALSE.
00028  */
00029 
00030 #include "amesos_cholmod_internal.h"
00031 #include "amesos_cholmod_core.h"
00032 
00033 cholmod_sparse *CHOLMOD(add)
00034 (
00035     /* ---- input ---- */
00036     cholmod_sparse *A,      /* matrix to add */
00037     cholmod_sparse *B,      /* matrix to add */
00038     double alpha [2],     /* scale factor for A */
00039     double beta [2],      /* scale factor for B */
00040     int values,       /* if TRUE compute the numerical values of C */
00041     int sorted,       /* if TRUE, sort columns of C */
00042     /* --------------- */
00043     cholmod_common *Common
00044 )
00045 {
00046     double *Ax, *Bx, *Cx, *W ;
00047     Int apacked, up, lo, nrow, ncol, bpacked, nzmax, pa, paend, pb, pbend, i,
00048   j, p, mark, nz ;
00049     Int *Ap, *Ai, *Anz, *Bp, *Bi, *Bnz, *Flag, *Cp, *Ci ;
00050     cholmod_sparse *A2, *B2, *C ;
00051 
00052     /* ---------------------------------------------------------------------- */
00053     /* check inputs */
00054     /* ---------------------------------------------------------------------- */
00055 
00056     RETURN_IF_NULL_COMMON (NULL) ;
00057     RETURN_IF_NULL (A, NULL) ;
00058     RETURN_IF_NULL (B, NULL) ;
00059     values = values &&
00060   (A->xtype != CHOLMOD_PATTERN) && (B->xtype != CHOLMOD_PATTERN) ;
00061     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN,
00062       values ? CHOLMOD_REAL : CHOLMOD_ZOMPLEX, NULL) ;
00063     RETURN_IF_XTYPE_INVALID (B, CHOLMOD_PATTERN,
00064       values ? CHOLMOD_REAL : CHOLMOD_ZOMPLEX, NULL) ;
00065     if (A->nrow != B->nrow || A->ncol != B->ncol)
00066     {
00067   /* A and B must have the same dimensions */
00068   ERROR (CHOLMOD_INVALID, "A and B dimesions do not match") ;
00069   return (NULL) ;
00070     }
00071     /* A and B must have the same numerical type if values is TRUE (both must
00072      * be CHOLMOD_REAL, this is implicitly checked above) */
00073 
00074     Common->status = CHOLMOD_OK ;
00075 
00076     /* ---------------------------------------------------------------------- */
00077     /* allocate workspace */
00078     /* ---------------------------------------------------------------------- */
00079 
00080     nrow = A->nrow ;
00081     ncol = A->ncol ;
00082     CHOLMOD(allocate_work) (nrow, MAX (nrow,ncol), values ? nrow : 0, Common) ;
00083     if (Common->status < CHOLMOD_OK)
00084     {
00085   return (NULL) ;     /* out of memory */
00086     }
00087 
00088     /* ---------------------------------------------------------------------- */
00089     /* get inputs */
00090     /* ---------------------------------------------------------------------- */
00091 
00092     if (nrow <= 1)
00093     {
00094   /* C will be implicitly sorted, so no need to sort it here */
00095   sorted = FALSE ;
00096     }
00097 
00098     /* convert A or B to unsymmetric, if necessary */
00099     A2 = NULL ;
00100     B2 = NULL ;
00101 
00102     if (A->stype != B->stype)
00103     {
00104   if (A->stype)
00105   {
00106       /* workspace: Iwork (max (nrow,ncol)) */
00107       A2 = CHOLMOD(copy) (A, 0, values, Common) ;
00108       if (Common->status < CHOLMOD_OK)
00109       {
00110     return (NULL) ;     /* out of memory */
00111       }
00112       A = A2 ;
00113   }
00114   if (B->stype)
00115   {
00116       /* workspace: Iwork (max (nrow,ncol)) */
00117       B2 = CHOLMOD(copy) (B, 0, values, Common) ;
00118       if (Common->status < CHOLMOD_OK)
00119       {
00120     CHOLMOD(free_sparse) (&A2, Common) ;
00121     return (NULL) ;     /* out of memory */
00122       }
00123       B = B2 ;
00124   }
00125     }
00126 
00127     /* get the A matrix */
00128     ASSERT (A->stype == B->stype) ;
00129     up = (A->stype > 0) ;
00130     lo = (A->stype < 0) ;
00131 
00132     Ap  = A->p ;
00133     Anz = A->nz ;
00134     Ai  = A->i ;
00135     Ax  = A->x ;
00136     apacked = A->packed ;
00137 
00138     /* get the B matrix */
00139     Bp  = B->p ;
00140     Bnz = B->nz ;
00141     Bi  = B->i ;
00142     Bx  = B->x ;
00143     bpacked = B->packed ;
00144 
00145     /* get workspace */
00146     W = Common->Xwork ;     /* size nrow, used if values is TRUE */
00147     Flag = Common->Flag ;   /* size nrow, Flag [0..nrow-1] < mark on input */
00148 
00149     /* ---------------------------------------------------------------------- */
00150     /* allocate the result C */
00151     /* ---------------------------------------------------------------------- */
00152 
00153     /* If integer overflow occurs, nzmax < 0 and the allocate fails properly
00154      * (likewise in most other matrix manipulation routines). */
00155     nzmax = A->nzmax + B->nzmax ;
00156     C = CHOLMOD(allocate_sparse) (nrow, ncol, nzmax, FALSE, TRUE,
00157       SIGN (A->stype), values ? A->xtype : CHOLMOD_PATTERN, Common) ;
00158     if (Common->status < CHOLMOD_OK)
00159     {
00160   CHOLMOD(free_sparse) (&A2, Common) ;
00161   CHOLMOD(free_sparse) (&B2, Common) ;
00162   return (NULL) ;     /* out of memory */
00163     }
00164 
00165     Cp = C->p ;
00166     Ci = C->i ;
00167     Cx = C->x ;
00168 
00169     /* ---------------------------------------------------------------------- */
00170     /* compute C = alpha*A + beta*B */
00171     /* ---------------------------------------------------------------------- */
00172 
00173     nz = 0 ;
00174     for (j = 0 ; j < ncol ; j++)
00175     {
00176   Cp [j] = nz ;
00177 
00178   /* clear the Flag array */
00179   mark = CHOLMOD(clear_flag) (Common) ;
00180 
00181   /* scatter B into W */
00182   pb = Bp [j] ;
00183   pbend = (bpacked) ? (Bp [j+1]) : (pb + Bnz [j]) ;
00184   for (p = pb ; p < pbend ; p++)
00185   {
00186       i = Bi [p] ;
00187       if ((up && i > j) || (lo && i < j))
00188       {
00189     continue ;
00190       }
00191       Flag [i] = mark ;
00192       if (values)
00193       {
00194     W [i] = beta [0] * Bx [p] ;
00195       }
00196   }
00197 
00198   /* add A and gather from W into C(:,j) */
00199   pa = Ap [j] ;
00200   paend = (apacked) ? (Ap [j+1]) : (pa + Anz [j]) ;
00201   for (p = pa ; p < paend ; p++)
00202   {
00203       i = Ai [p] ;
00204       if ((up && i > j) || (lo && i < j))
00205       {
00206     continue ;
00207       }
00208       Flag [i] = EMPTY ;
00209       Ci [nz] = i ;
00210       if (values)
00211       {
00212     Cx [nz] = W [i] + alpha [0] * Ax [p] ;
00213     W [i] = 0 ;
00214       }
00215       nz++ ;
00216   }
00217 
00218   /* gather remaining entries into C(:,j), using pattern of B */
00219   for (p = pb ; p < pbend ; p++)
00220   {
00221       i = Bi [p] ;
00222       if ((up && i > j) || (lo && i < j))
00223       {
00224     continue ;
00225       }
00226       if (Flag [i] == mark)
00227       {
00228     Ci [nz] = i ;
00229     if (values)
00230     {
00231         Cx [nz] = W [i] ;
00232         W [i] = 0 ;
00233     }
00234     nz++ ;
00235       }
00236   }
00237     }
00238 
00239     Cp [ncol] = nz ;
00240 
00241     /* ---------------------------------------------------------------------- */
00242     /* reduce C in size and free temporary matrices */
00243     /* ---------------------------------------------------------------------- */
00244 
00245     ASSERT (MAX (1,nz) <= C->nzmax) ;
00246     CHOLMOD(reallocate_sparse) (nz, C, Common) ;
00247     ASSERT (Common->status >= CHOLMOD_OK) ;
00248 
00249     /* clear the Flag array */
00250     mark = CHOLMOD(clear_flag) (Common) ;
00251 
00252     CHOLMOD(free_sparse) (&A2, Common) ;
00253     CHOLMOD(free_sparse) (&B2, Common) ;
00254 
00255     /* ---------------------------------------------------------------------- */
00256     /* sort C, if requested */
00257     /* ---------------------------------------------------------------------- */
00258 
00259     if (sorted)
00260     {
00261   /* workspace: Iwork (max (nrow,ncol)) */
00262   if (!CHOLMOD(sort) (C, Common))
00263   {
00264       CHOLMOD(free_sparse) (&C, Common) ;
00265       if (Common->status < CHOLMOD_OK)
00266       {
00267     return (NULL) ;   /* out of memory */
00268       }
00269   }
00270     }
00271 
00272     /* ---------------------------------------------------------------------- */
00273     /* return result */
00274     /* ---------------------------------------------------------------------- */
00275 
00276     ASSERT (CHOLMOD(dump_sparse) (C, "add", Common) >= 0) ;
00277     return (C) ;
00278 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines