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