Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_copy.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === Core/cholmod_copy ==================================================== */
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, which allocates C and copies A into C, with possible change of
00015  * stype.  The diagonal can optionally be removed.  The numerical entries
00016  * can optionally be copied.  This routine differs from cholmod_copy_sparse,
00017  * which makes an exact copy of a sparse matrix.
00018  *
00019  * A can be of any type (packed/unpacked, upper/lower/unsymmetric).  C is
00020  * packed and can be of any stype (upper/lower/unsymmetric), except that if
00021  * A is rectangular C can only be unsymmetric.  If the stype of A and C
00022  * differ, then the appropriate conversion is made.
00023  *
00024  * Symmetry of A (A->stype):
00025  * <0: lower: assume A is symmetric with just tril(A); the rest of A is ignored
00026  *  0  unsym: assume A is unsymmetric; consider all entries in A
00027  * >0  upper: assume A is symmetric with just triu(A); the rest of A is ignored
00028  *
00029  * Symmetry of C (stype parameter):
00030  * <0  lower: return just tril(C)
00031  *  0  unsym: return all of C
00032  * >0  upper: return just triu(C)
00033  *
00034  * In MATLAB:             Using cholmod_copy:
00035  * ----------             ----------------------------
00036  * C = A ;              A unsymmetric, C unsymmetric
00037  * C = tril (A) ;           A unsymmetric, C lower
00038  * C = triu (A) ;           A unsymmetric, C upper
00039  * U = triu (A) ; L = tril (U',-1) ; C = L+U ;      A upper, C unsymmetric
00040  * C = triu (A)' ;            A upper, C lower
00041  * C = triu (A) ;           A upper, C upper
00042  * L = tril (A) ; U = triu (L',1) ; C = L+U ;     A lower, C unsymmetric
00043  * C = tril (A) ;           A lower, C lower
00044  * C = tril (A)' ;            A lower, C upper
00045  *
00046  * workspace: Iwork (max (nrow,ncol))
00047  *
00048  * A can have an xtype of pattern or real.  Complex and zomplex cases only
00049  * supported when mode <= 0 (in which case the numerical values are ignored).
00050  */
00051 
00052 #include "amesos_cholmod_internal.h"
00053 #include "amesos_cholmod_core.h"
00054 
00055 
00056 /* ========================================================================== */
00057 /* === copy_sym_to_unsym ==================================================== */
00058 /* ========================================================================== */
00059 
00060 /* Construct an unsymmetric copy of a symmetric sparse matrix.  This does the
00061  * work for as C = cholmod_copy (A, 0, mode, Common) when A is symmetric.
00062  * In this case, extra space can be added to C.
00063  */
00064 
00065 static cholmod_sparse *copy_sym_to_unsym
00066 (
00067     /* ---- input ---- */
00068     cholmod_sparse *A,  /* matrix to copy */
00069     int mode,   /* >0: numerical, 0: pattern, <0: pattern (no diag)
00070        * -2: pattern only, no diagonal, add 50% + n extra
00071        * space to C */
00072     /* --------------- */
00073     cholmod_common *Common
00074 )
00075 {
00076     double aij ;
00077     double *Ax, *Cx ;
00078     Int *Ap, *Ai, *Anz, *Cp, *Ci, *Wj, *Iwork ;
00079     cholmod_sparse *C ;
00080     Int nrow, ncol, nz, packed, j, p, pend, i, pc, up, lo, values, diag,
00081   astype, extra ;
00082 
00083     /* ---------------------------------------------------------------------- */
00084     /* get inputs */
00085     /* ---------------------------------------------------------------------- */
00086 
00087     nrow = A->nrow ;
00088     ncol = A->ncol ;
00089     Ap  = A->p ;
00090     Anz = A->nz ;
00091     Ai  = A->i ;
00092     Ax  = A->x ;
00093     packed = A->packed ;
00094     values = (mode > 0) && (A->xtype != CHOLMOD_PATTERN) ;
00095     diag = (mode >= 0) ;
00096 
00097     astype = SIGN (A->stype) ;
00098     up = (astype > 0) ;
00099     lo = (astype < 0) ;
00100     ASSERT (astype != 0) ;
00101 
00102     /* ---------------------------------------------------------------------- */
00103     /* create an unsymmetric copy of a symmetric matrix */
00104     /* ---------------------------------------------------------------------- */
00105 
00106     Iwork = Common->Iwork ;
00107     Wj = Iwork ;        /* size ncol (i/i/l) */
00108 
00109     /* In MATLAB notation, for converting a symmetric/upper matrix:
00110      *  U = triu (A) ;
00111      *  L = tril (U',-1) ;
00112      *  C = L + U ;
00113      *
00114      * For converting a symmetric/lower matrix to unsymmetric:
00115      *  L = tril (A) ;
00116      *  U = triu (L',1) ;
00117      *  C = L + U ;
00118      */
00119     ASSERT (up || lo) ;
00120     PRINT1 (("copy: convert symmetric to unsym\n")) ;
00121 
00122     /* count the number of entries in each column of C */
00123     for (j = 0 ; j < ncol ; j++)
00124     {
00125   Wj [j] = 0 ;
00126     }
00127     for (j = 0 ; j < ncol ; j++)
00128     {
00129   p = Ap [j] ;
00130   pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
00131   for ( ; p < pend ; p++)
00132   {
00133       i = Ai [p] ;
00134       if (i == j)
00135       {
00136     /* the diagonal entry A(i,i) will appear just once
00137      * (unless it is excluded with mode < 0) */
00138     if (diag)
00139     {
00140         Wj [j]++ ;
00141     }
00142       }
00143       else if ((up && i < j) || (lo && i > j))
00144       {
00145     /* upper case:  A(i,j) is in the strictly upper part;
00146      * A(j,i) will be added to the strictly lower part of C.
00147      * lower case is the opposite. */
00148     Wj [j]++ ;
00149     Wj [i]++ ;
00150       }
00151   }
00152     }
00153     nz = 0 ;
00154     for (j = 0 ; j < ncol ; j++)
00155     {
00156   nz += Wj [j] ;
00157     }
00158 
00159     extra = (mode == -2) ? (nz/2 + ncol) : 0 ;
00160 
00161     /* allocate C.  C is sorted if and only if A is sorted */
00162     C = CHOLMOD(allocate_sparse) (nrow, ncol, nz + extra, A->sorted, TRUE, 0,
00163       values ? A->xtype : CHOLMOD_PATTERN, Common) ;
00164     if (Common->status < CHOLMOD_OK)
00165     {
00166   return (NULL) ;
00167     }
00168 
00169     Cp = C->p ;
00170     Ci = C->i ;
00171     Cx = C->x ;
00172 
00173     /* construct the column pointers for C */
00174     p = 0 ;
00175     for (j = 0 ; j < ncol ; j++)
00176     {
00177   Cp [j] = p ;
00178   p += Wj [j] ;
00179     }
00180     Cp [ncol] = p ;
00181     for (j = 0 ; j < ncol ; j++)
00182     {
00183   Wj [j] = Cp [j] ;
00184     }
00185 
00186     /* construct C */
00187     if (values)
00188     {
00189 
00190   /* pattern and values */
00191   ASSERT (diag) ;
00192   for (j = 0 ; j < ncol ; j++)
00193   {
00194       p = Ap [j] ;
00195       pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
00196       for ( ; p < pend ; p++)
00197       {
00198     i = Ai [p] ;
00199     aij = Ax [p] ;
00200     if (i == j)
00201     {
00202         /* add diagonal entry A(i,i) to column i */
00203         pc = Wj [i]++ ;
00204         Ci [pc] = i ;
00205         Cx [pc] = aij ;
00206     }
00207     else if ((up && i < j) || (lo && i > j))
00208     {
00209         /* add A(i,j) to column j */
00210         pc = Wj [j]++ ;
00211         Ci [pc] = i ;
00212         Cx [pc] = aij ;
00213         /* add A(j,i) to column i */
00214         pc = Wj [i]++ ;
00215         Ci [pc] = j ;
00216         Cx [pc] = aij ;
00217     }
00218       }
00219   }
00220 
00221     }
00222     else
00223     {
00224 
00225   /* pattern only, possibly excluding the diagonal */
00226   for (j = 0 ; j < ncol ; j++)
00227   {
00228       p = Ap [j] ;
00229       pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
00230       for ( ; p < pend ; p++)
00231       {
00232     i = Ai [p] ;
00233     if (i == j)
00234     {
00235         /* add diagonal entry A(i,i) to column i
00236          * (unless it is excluded with mode < 0) */
00237         if (diag)
00238         {
00239       Ci [Wj [i]++] = i ;
00240         }
00241     }
00242     else if ((up && i < j) || (lo && i > j))
00243     {
00244         /* add A(i,j) to column j */
00245         Ci [Wj [j]++] = i ;
00246         /* add A(j,i) to column i */
00247         Ci [Wj [i]++] = j ;
00248     }
00249       }
00250   }
00251     }
00252 
00253     /* ---------------------------------------------------------------------- */
00254     /* return the result */
00255     /* ---------------------------------------------------------------------- */
00256 
00257     DEBUG (i = CHOLMOD(dump_sparse) (C, "copy_sym_to_unsym", Common)) ;
00258     PRINT1 (("mode %d nnzdiag "ID"\n", mode, i)) ;
00259     ASSERT (IMPLIES (mode < 0, i == 0)) ;
00260     return (C) ;
00261 }
00262 
00263 
00264 /* ========================================================================== */
00265 /* === cholmod_copy ========================================================= */
00266 /* ========================================================================== */
00267 
00268 cholmod_sparse *CHOLMOD(copy)
00269 (
00270     /* ---- input ---- */
00271     cholmod_sparse *A,  /* matrix to copy */
00272     int stype,    /* requested stype of C */
00273     int mode,   /* >0: numerical, 0: pattern, <0: pattern (no diag) */
00274     /* --------------- */
00275     cholmod_common *Common
00276 )
00277 {
00278     cholmod_sparse *C ;
00279     Int nrow, ncol, up, lo, values, diag, astype ;
00280 
00281     /* ---------------------------------------------------------------------- */
00282     /* check inputs */
00283     /* ---------------------------------------------------------------------- */
00284 
00285     RETURN_IF_NULL_COMMON (NULL) ;
00286     RETURN_IF_NULL (A, NULL) ;
00287     values = (mode > 0) && (A->xtype != CHOLMOD_PATTERN) ;
00288     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN,
00289       values ? CHOLMOD_REAL : CHOLMOD_ZOMPLEX, NULL) ;
00290     nrow = A->nrow ;
00291     ncol = A->ncol ;
00292     if ((stype || A->stype) && nrow != ncol)
00293     {
00294   /* inputs invalid */
00295   ERROR (CHOLMOD_INVALID, "matrix invalid") ;
00296   return (NULL) ;
00297     }
00298     Common->status = CHOLMOD_OK ;
00299 
00300     /* ---------------------------------------------------------------------- */
00301     /* allocate workspace */
00302     /* ---------------------------------------------------------------------- */
00303 
00304     CHOLMOD(allocate_work) (0, MAX (nrow,ncol), 0, Common) ;
00305     if (Common->status < CHOLMOD_OK)
00306     {
00307   /* out of memory */
00308   return (NULL) ;
00309     }
00310 
00311     /* ---------------------------------------------------------------------- */
00312     /* get inputs */
00313     /* ---------------------------------------------------------------------- */
00314 
00315     diag = (mode >= 0) ;
00316     astype = SIGN (A->stype) ;
00317     stype = SIGN (stype) ;
00318     up = (astype > 0) ;
00319     lo = (astype < 0) ;
00320 
00321     /* ---------------------------------------------------------------------- */
00322     /* copy the matrix */
00323     /* ---------------------------------------------------------------------- */
00324 
00325     if (astype == stype)
00326     {
00327 
00328   /* ------------------------------------------------------------------ */
00329   /* symmetry of A and C are the same */
00330   /* ------------------------------------------------------------------ */
00331 
00332   /* copy A into C, keeping the same symmetry.  If A is symmetric
00333    * entries in the ignored part of A are not copied into C */
00334   C = CHOLMOD(band) (A, -nrow, ncol, mode, Common) ;
00335 
00336     }
00337     else if (!astype)
00338     {
00339 
00340   /* ------------------------------------------------------------------ */
00341   /* convert unsymmetric matrix A into a symmetric matrix C */
00342   /* ------------------------------------------------------------------ */
00343 
00344   if (stype > 0)
00345   {
00346       /* C = triu (A) */
00347       C = CHOLMOD(band) (A, 0, ncol, mode, Common) ;
00348   }
00349   else
00350   {
00351       /* C = tril (A) */
00352       C = CHOLMOD(band) (A, -nrow, 0, mode, Common) ;
00353   }
00354   if (Common->status < CHOLMOD_OK)
00355   {
00356       /* out of memory */
00357       return (NULL) ;
00358   }
00359   C->stype = stype ;
00360 
00361     }
00362     else if (astype == -stype)
00363     {
00364 
00365   /* ------------------------------------------------------------------ */
00366   /* transpose a symmetric matrix */
00367   /* ------------------------------------------------------------------ */
00368 
00369   /* converting upper to lower or lower to upper */
00370   /* workspace: Iwork (nrow) */
00371   C = CHOLMOD(transpose) (A, values, Common) ;
00372   if (!diag)
00373   {
00374       /* remove diagonal, if requested */
00375       CHOLMOD(band_inplace) (-nrow, ncol, -1, C, Common) ;
00376   }
00377 
00378     }
00379     else
00380     {
00381 
00382   /* ------------------------------------------------------------------ */
00383   /* create an unsymmetric copy of a symmetric matrix */
00384   /* ------------------------------------------------------------------ */
00385 
00386   C = copy_sym_to_unsym (A, mode, Common) ;
00387     }
00388 
00389     /* ---------------------------------------------------------------------- */
00390     /* return if error */
00391     /* ---------------------------------------------------------------------- */
00392 
00393     if (Common->status < CHOLMOD_OK)
00394     {
00395   /* out of memory */
00396   return (NULL) ;
00397     }
00398 
00399     /* ---------------------------------------------------------------------- */
00400     /* return the result */
00401     /* ---------------------------------------------------------------------- */
00402 
00403     DEBUG (diag = CHOLMOD(dump_sparse) (C, "copy", Common)) ;
00404     PRINT1 (("mode %d nnzdiag "ID"\n", mode, diag)) ;
00405     ASSERT (IMPLIES (mode < 0, diag == 0)) ;
00406     return (C) ;
00407 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines