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