Amesos Package Browser (Single Doxygen Collection) Development
amesos_amd_l_order.c
Go to the documentation of this file.
00001 /* ========================================================================= */
00002 /* === AMD_order =========================================================== */
00003 /* ========================================================================= */
00004 
00005 /* ------------------------------------------------------------------------- */
00006 /* AMD, Copyright (c) Timothy A. Davis,              */
00007 /* Patrick R. Amestoy, and Iain S. Duff.  See ../README.txt for License.     */
00008 /* email: davis at cise.ufl.edu    CISE Department, Univ. of Florida.        */
00009 /* web: http://www.cise.ufl.edu/research/sparse/amd                          */
00010 /* ------------------------------------------------------------------------- */
00011 
00012 /* User-callable AMD minimum degree ordering routine.  See amd.h for
00013  * documentation.
00014  */
00015 
00016 /* This file should make the long int version of AMD */
00017 #define DLONG 1
00018 
00019 #include "amesos_amd_internal.h"
00020 
00021 /* ========================================================================= */
00022 /* === AMD_order =========================================================== */
00023 /* ========================================================================= */
00024 
00025 GLOBAL Int AMD_order
00026 (
00027     Int n,
00028     const Int Ap [ ],
00029     const Int Ai [ ],
00030     Int P [ ],
00031     double Control [ ],
00032     double Info [ ]
00033 )
00034 {
00035     Int *Len, *S, nz, i, *Pinv, info, status, *Rp, *Ri, *Cp, *Ci, ok ;
00036     size_t nzaat, slen ;
00037     double mem = 0 ;
00038 
00039 #ifndef NDEBUG
00040     AMD_debug_init ("amd") ;
00041 #endif
00042 
00043     /* clear the Info array, if it exists */
00044     info = Info != (double *) NULL ;
00045     if (info)
00046     {
00047   for (i = 0 ; i < AMD_INFO ; i++)
00048   {
00049       Info [i] = EMPTY ;
00050   }
00051   Info [AMD_N] = n ;
00052   Info [AMD_STATUS] = AMD_OK ;
00053     }
00054 
00055     /* make sure inputs exist and n is >= 0 */
00056     if (Ai == (Int *) NULL || Ap == (Int *) NULL || P == (Int *) NULL || n < 0)
00057     {
00058   if (info) Info [AMD_STATUS] = AMD_INVALID ;
00059   return (AMD_INVALID) ;      /* arguments are invalid */
00060     }
00061 
00062     if (n == 0)
00063     {
00064   return (AMD_OK) ;     /* n is 0 so there's nothing to do */
00065     }
00066 
00067     nz = Ap [n] ;
00068     if (info)
00069     {
00070   Info [AMD_NZ] = nz ;
00071     }
00072     if (nz < 0)
00073     {
00074   if (info) Info [AMD_STATUS] = AMD_INVALID ;
00075   return (AMD_INVALID) ;
00076     }
00077 
00078     /* check if n or nz will cause size_t overflow */
00079     if (((size_t) n) >= SIZE_T_MAX / sizeof (Int)
00080      || ((size_t) nz) >= SIZE_T_MAX / sizeof (Int))
00081     {
00082   if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ;
00083   return (AMD_OUT_OF_MEMORY) ;      /* problem too large */
00084     }
00085 
00086     /* check the input matrix:  AMD_OK, AMD_INVALID, or AMD_OK_BUT_JUMBLED */
00087     status = AMD_valid (n, n, Ap, Ai) ;
00088 
00089     if (status == AMD_INVALID)
00090     {
00091   if (info) Info [AMD_STATUS] = AMD_INVALID ;
00092   return (AMD_INVALID) ;      /* matrix is invalid */
00093     }
00094 
00095     /* allocate two size-n integer workspaces */
00096     Len = amesos_amd_malloc (n * sizeof (Int)) ;
00097     Pinv = amesos_amd_malloc (n * sizeof (Int)) ;
00098     mem += n ;
00099     mem += n ;
00100     if (!Len || !Pinv)
00101     {
00102   /* :: out of memory :: */
00103   amesos_amd_free (Len) ;
00104   amesos_amd_free (Pinv) ;
00105   if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ;
00106   return (AMD_OUT_OF_MEMORY) ;
00107     }
00108 
00109     if (status == AMD_OK_BUT_JUMBLED)
00110     {
00111   /* sort the input matrix and remove duplicate entries */
00112   AMD_DEBUG1 (("Matrix is jumbled\n")) ;
00113   Rp = amesos_amd_malloc ((n+1) * sizeof (Int)) ;
00114   Ri = amesos_amd_malloc (MAX (nz,1) * sizeof (Int)) ;
00115   mem += (n+1) ;
00116   mem += MAX (nz,1) ;
00117   if (!Rp || !Ri)
00118   {
00119       /* :: out of memory :: */
00120       amesos_amd_free (Rp) ;
00121       amesos_amd_free (Ri) ;
00122       amesos_amd_free (Len) ;
00123       amesos_amd_free (Pinv) ;
00124       if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ;
00125       return (AMD_OUT_OF_MEMORY) ;
00126   }
00127   /* use Len and Pinv as workspace to create R = A' */
00128   AMD_preprocess (n, Ap, Ai, Rp, Ri, Len, Pinv) ;
00129   Cp = Rp ;
00130   Ci = Ri ;
00131     }
00132     else
00133     {
00134   /* order the input matrix as-is.  No need to compute R = A' first */
00135   Rp = NULL ;
00136   Ri = NULL ;
00137   Cp = (Int *) Ap ;
00138   Ci = (Int *) Ai ;
00139     }
00140 
00141     /* --------------------------------------------------------------------- */
00142     /* determine the symmetry and count off-diagonal nonzeros in A+A' */
00143     /* --------------------------------------------------------------------- */
00144 
00145     nzaat = AMD_aat (n, Cp, Ci, Len, P, Info) ;
00146     AMD_DEBUG1 (("nzaat: %g\n", (double) nzaat)) ;
00147     ASSERT ((MAX (nz-n, 0) <= nzaat) && (nzaat <= 2 * (size_t) nz)) ;
00148 
00149     /* --------------------------------------------------------------------- */
00150     /* allocate workspace for matrix, elbow room, and 6 size-n vectors */
00151     /* --------------------------------------------------------------------- */
00152 
00153     S = NULL ;
00154     slen = nzaat ;      /* space for matrix */
00155     ok = ((slen + nzaat/5) >= slen) ; /* check for size_t overflow */
00156     slen += nzaat/5 ;     /* add elbow room */
00157     for (i = 0 ; ok && i < 7 ; i++)
00158     {
00159   ok = ((slen + n) > slen) ;  /* check for size_t overflow */
00160   slen += n ;     /* size-n elbow room, 6 size-n work */
00161     }
00162     mem += slen ;
00163     ok = ok && (slen < SIZE_T_MAX / sizeof (Int)) ; /* check for overflow */
00164     ok = ok && (slen < Int_MAX) ; /* S[i] for Int i must be OK */
00165     if (ok)
00166     {
00167   S = amesos_amd_malloc (slen * sizeof (Int)) ;
00168     }
00169     AMD_DEBUG1 (("slen %g\n", (double) slen)) ;
00170     if (!S)
00171     {
00172   /* :: out of memory :: (or problem too large) */
00173   amesos_amd_free (Rp) ;
00174   amesos_amd_free (Ri) ;
00175   amesos_amd_free (Len) ;
00176   amesos_amd_free (Pinv) ;
00177   if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ;
00178   return (AMD_OUT_OF_MEMORY) ;
00179     }
00180     if (info)
00181     {
00182   /* memory usage, in bytes. */
00183   Info [AMD_MEMORY] = mem * sizeof (Int) ;
00184     }
00185 
00186     /* --------------------------------------------------------------------- */
00187     /* order the matrix */
00188     /* --------------------------------------------------------------------- */
00189 
00190     AMD_1 (n, Cp, Ci, P, Pinv, Len, slen, S, Control, Info) ;
00191 
00192     /* --------------------------------------------------------------------- */
00193     /* free the workspace */
00194     /* --------------------------------------------------------------------- */
00195 
00196     amesos_amd_free (Rp) ;
00197     amesos_amd_free (Ri) ;
00198     amesos_amd_free (Len) ;
00199     amesos_amd_free (Pinv) ;
00200     amesos_amd_free (S) ;
00201     if (info) Info [AMD_STATUS] = status ;
00202     return (status) ;     /* successful ordering */
00203 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines