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