Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_band.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === Core/cholmod_band ==================================================== */
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 = tril (triu (A,k1), k2)
00015  *
00016  * C is a matrix consisting of the diagonals of A from k1 to k2.
00017  *
00018  * k=0 is the main diagonal of A, k=1 is the superdiagonal, k=-1 is the
00019  * subdiagonal, and so on.  If A is m-by-n, then:
00020  *
00021  *  k1=-m       C = tril (A,k2)
00022  *  k2=n        C = triu (A,k1)
00023  *  k1=0 and k2=0     C = diag(A), except C is a matrix, not a vector
00024  *
00025  * Values of k1 and k2 less than -m are treated as -m, and values greater
00026  * than n are treated as n.
00027  *
00028  * A can be of any symmetry (upper, lower, or unsymmetric); C is returned in
00029  * the same form, and packed.  If A->stype > 0, entries in the lower
00030  * triangular part of A are ignored, and the opposite is true if
00031  * A->stype < 0.  If A has sorted columns, then so does C.
00032  * C has the same size as A.
00033  *
00034  * If inplace is TRUE, then the matrix A is modified in place.
00035  * Only packed matrices can be converted in place.
00036  *
00037  * C can be returned as a numerical valued matrix (if A has numerical values
00038  * and mode > 0), as a pattern-only (mode == 0), or as a pattern-only but with
00039  * the diagonal entries removed (mode < 0).
00040  *
00041  * workspace: none
00042  *
00043  * A can have an xtype of pattern or real.  Complex and zomplex cases supported
00044  * only if mode <= 0 (in which case the numerical values are ignored).
00045  */
00046 
00047 #include "amesos_cholmod_internal.h"
00048 #include "amesos_cholmod_core.h"
00049 
00050 static cholmod_sparse *band   /* returns C, or NULL if failure */
00051 (
00052     /* ---- input or in/out if inplace is TRUE --- */
00053     cholmod_sparse *A,
00054     /* ---- input ---- */
00055     UF_long k1,     /* ignore entries below the k1-st diagonal */
00056     UF_long k2,     /* ignore entries above the k2-nd diagonal */
00057     int mode,     /* >0: numerical, 0: pattern, <0: pattern (no diagonal) */
00058     int inplace,    /* if TRUE, then convert A in place */
00059     /* --------------- */
00060     cholmod_common *Common
00061 )
00062 {
00063     double *Ax, *Cx ;
00064     Int packed, nz, j, p, pend, i, ncol, nrow, jlo, jhi, ilo, ihi, sorted,
00065   values, diag ;
00066     Int *Ap, *Anz, *Ai, *Cp, *Ci ;
00067     cholmod_sparse *C ;
00068 
00069     /* ---------------------------------------------------------------------- */
00070     /* check inputs */
00071     /* ---------------------------------------------------------------------- */
00072 
00073     RETURN_IF_NULL_COMMON (NULL) ;
00074     RETURN_IF_NULL (A, NULL) ;
00075     values = (mode > 0) && (A->xtype != CHOLMOD_PATTERN) ;
00076     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN,
00077       values ? CHOLMOD_REAL : CHOLMOD_ZOMPLEX, NULL) ;
00078     packed = A->packed ;
00079     diag = (mode >= 0) ;
00080     if (inplace && !packed)
00081     {
00082   /* cannot operate on an unpacked matrix in place */
00083   ERROR (CHOLMOD_INVALID, "cannot operate on unpacked matrix in-place") ;
00084   return (NULL) ;
00085     }
00086     Common->status = CHOLMOD_OK ;
00087 
00088     /* ---------------------------------------------------------------------- */
00089     /* get inputs */
00090     /* ---------------------------------------------------------------------- */
00091 
00092 
00093     PRINT1 (("k1 %ld k2 %ld\n", k1, k2)) ;
00094     Ap  = A->p ;
00095     Anz = A->nz ;
00096     Ai  = A->i ;
00097     Ax  = A->x ;
00098     sorted = A->sorted ;
00099 
00100 
00101     if (A->stype > 0)
00102     {
00103   /* ignore any entries in strictly lower triangular part of A */
00104   k1 = MAX (k1, 0) ;
00105     }
00106     if (A->stype < 0)
00107     {
00108   /* ignore any entries in strictly upper triangular part of A */
00109   k2 = MIN (k2, 0) ;
00110     }
00111     ncol = A->ncol ;
00112     nrow = A->nrow ;
00113 
00114     /* ensure k1 and k2 are in the range -nrow to +ncol to
00115      * avoid possible integer overflow if k1 and k2 are huge */
00116     k1 = MAX (-nrow, k1) ;
00117     k1 = MIN (k1, ncol) ;
00118     k2 = MAX (-nrow, k2) ;
00119     k2 = MIN (k2, ncol) ;
00120 
00121     /* consider columns jlo to jhi.  columns outside this range are empty */
00122     jlo = MAX (k1, 0) ;
00123     jhi = MIN (k2+nrow, ncol) ;
00124 
00125     if (k1 > k2)
00126     {
00127   /* nothing to do */
00128   jlo = ncol ;
00129   jhi = ncol ;
00130     }
00131 
00132     /* ---------------------------------------------------------------------- */
00133     /* allocate C, or operate on A in place */
00134     /* ---------------------------------------------------------------------- */
00135 
00136     if (inplace)
00137     {
00138   /* convert A in place */
00139   C = A ;
00140     }
00141     else
00142     {
00143   /* count the number of entries in the result C */
00144   nz = 0 ;
00145   if (sorted)
00146   {
00147       for (j = jlo ; j < jhi ; j++)
00148       {
00149     ilo = j-k2 ;
00150     ihi = j-k1 ;
00151     p = Ap [j] ;
00152     pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
00153     for ( ; p < pend ; p++)
00154     {
00155         i = Ai [p] ;
00156         if (i > ihi)
00157         {
00158       break ;
00159         }
00160         if (i >= ilo && (diag || i != j))
00161         {
00162       nz++ ;
00163         }
00164     }
00165       }
00166   }
00167   else
00168   {
00169       for (j = jlo ; j < jhi ; j++)
00170       {
00171     ilo = j-k2 ;
00172     ihi = j-k1 ;
00173     p = Ap [j] ;
00174     pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
00175     for ( ; p < pend ; p++)
00176     {
00177         i = Ai [p] ;
00178         if (i >= ilo && i <= ihi && (diag || i != j))
00179         {
00180       nz++ ;
00181         }
00182     }
00183       }
00184   }
00185   /* allocate C; A will not be modified.  C is sorted if A is sorted */
00186   C = CHOLMOD(allocate_sparse) (A->nrow, ncol, nz, sorted, TRUE,
00187     A->stype, values ? A->xtype : CHOLMOD_PATTERN, Common) ;
00188   if (Common->status < CHOLMOD_OK)
00189   {
00190       return (NULL) ; /* out of memory */
00191   }
00192     }
00193 
00194     Cp = C->p ;
00195     Ci = C->i ;
00196     Cx = C->x ;
00197 
00198     /* ---------------------------------------------------------------------- */
00199     /* construct C */
00200     /* ---------------------------------------------------------------------- */
00201 
00202     /* columns 0 to jlo-1 are empty */
00203     for (j = 0 ; j < jlo ; j++)
00204     {
00205   Cp [j] = 0 ;
00206     }
00207 
00208     nz = 0 ;
00209     if (sorted)
00210     {
00211   if (values)
00212   {
00213       /* pattern and values */
00214       ASSERT (diag) ;
00215       for (j = jlo ; j < jhi ; j++)
00216       {
00217     ilo = j-k2 ;
00218     ihi = j-k1 ;
00219     p = Ap [j] ;
00220     pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
00221     Cp [j] = nz ;
00222     for ( ; p < pend ; p++)
00223     {
00224         i = Ai [p] ;
00225         if (i > ihi)
00226         {
00227       break ;
00228         }
00229         if (i >= ilo)
00230         {
00231       Ci [nz] = i ;
00232       Cx [nz] = Ax [p] ;
00233       nz++ ;
00234         }
00235     }
00236       }
00237   }
00238   else
00239   {
00240       /* pattern only, perhaps with no diagonal */
00241       for (j = jlo ; j < jhi ; j++)
00242       {
00243     ilo = j-k2 ;
00244     ihi = j-k1 ;
00245     p = Ap [j] ;
00246     pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
00247     Cp [j] = nz ;
00248     for ( ; p < pend ; p++)
00249     {
00250         i = Ai [p] ;
00251         if (i > ihi)
00252         {
00253       break ;
00254         }
00255         if (i >= ilo && (diag || i != j))
00256         {
00257       Ci [nz++] = i ;
00258         }
00259     }
00260       }
00261   }
00262     }
00263     else
00264     {
00265   if (values)
00266   {
00267       /* pattern and values */
00268       ASSERT (diag) ;
00269       for (j = jlo ; j < jhi ; j++)
00270       {
00271     ilo = j-k2 ;
00272     ihi = j-k1 ;
00273     p = Ap [j] ;
00274     pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
00275     Cp [j] = nz ;
00276     for ( ; p < pend ; p++)
00277     {
00278         i = Ai [p] ;
00279         if (i >= ilo && i <= ihi)
00280         {
00281       Ci [nz] = i ;
00282       Cx [nz] = Ax [p] ;
00283       nz++ ;
00284         }
00285     }
00286       }
00287   }
00288   else
00289   {
00290       /* pattern only, perhaps with no diagonal */
00291       for (j = jlo ; j < jhi ; j++)
00292       {
00293     ilo = j-k2 ;
00294     ihi = j-k1 ;
00295     p = Ap [j] ;
00296     pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
00297     Cp [j] = nz ;
00298     for ( ; p < pend ; p++)
00299     {
00300         i = Ai [p] ;
00301         if (i >= ilo && i <= ihi && (diag || i != j))
00302         {
00303       Ci [nz++] = i ;
00304         }
00305     }
00306       }
00307   }
00308     }
00309 
00310     /* columns jhi to ncol-1 are empty */
00311     for (j = jhi ; j <= ncol ; j++)
00312     {
00313   Cp [j] = nz ;
00314     }
00315 
00316     /* ---------------------------------------------------------------------- */
00317     /* reduce A in size if done in place */
00318     /* ---------------------------------------------------------------------- */
00319 
00320     if (inplace)
00321     {
00322   /* free the unused parts of A, and reduce A->i and A->x in size */
00323   ASSERT (MAX (1,nz) <= A->nzmax) ;
00324   CHOLMOD(reallocate_sparse) (nz, A, Common) ;
00325   ASSERT (Common->status >= CHOLMOD_OK) ;
00326     }
00327 
00328     /* ---------------------------------------------------------------------- */
00329     /* return the result C */
00330     /* ---------------------------------------------------------------------- */
00331 
00332     DEBUG (i = CHOLMOD(dump_sparse) (C, "band", Common)) ;
00333     ASSERT (IMPLIES (mode < 0, i == 0)) ;
00334     return (C) ;
00335 }
00336 
00337 
00338 /* ========================================================================== */
00339 /* === cholmod_band ========================================================= */
00340 /* ========================================================================== */
00341 
00342 cholmod_sparse *CHOLMOD(band)
00343 (
00344     /* ---- input ---- */
00345     cholmod_sparse *A,  /* matrix to extract band matrix from */
00346     UF_long k1,   /* ignore entries below the k1-st diagonal */
00347     UF_long k2,   /* ignore entries above the k2-nd diagonal */
00348     int mode,   /* >0: numerical, 0: pattern, <0: pattern (no diag) */
00349     /* --------------- */
00350     cholmod_common *Common
00351 )
00352 {
00353     return (band (A, k1, k2, mode, FALSE, Common)) ;
00354 }
00355 
00356 
00357 /* ========================================================================== */
00358 /* === cholmod_band_inplace ================================================= */
00359 /* ========================================================================== */
00360 
00361 int CHOLMOD(band_inplace)
00362 (
00363     /* ---- input ---- */
00364     UF_long k1,   /* ignore entries below the k1-st diagonal */
00365     UF_long k2,   /* ignore entries above the k2-nd diagonal */
00366     int mode,   /* >0: numerical, 0: pattern, <0: pattern (no diag) */
00367     /* ---- in/out --- */
00368     cholmod_sparse *A,  /* matrix from which entries not in band are removed */
00369     /* --------------- */
00370     cholmod_common *Common
00371 )
00372 {
00373     return (band (A, k1, k2, mode, TRUE, Common) != NULL) ;
00374 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines