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