Amesos Package Browser (Single Doxygen Collection) Development
amesos_klu_analyze_given.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === klu_analyze_given ==================================================== */
00003 /* ========================================================================== */
00004 
00005 /* Given an input permutation P and Q, create the Symbolic object.  BTF can
00006  * be done to modify the user's P and Q (does not perform the max transversal;
00007  * just finds the strongly-connected components). */
00008 
00009 #include "amesos_klu_internal.h"
00010 
00011 /* ========================================================================== */
00012 /* === klu_alloc_symbolic =================================================== */
00013 /* ========================================================================== */
00014 
00015 /* Allocate Symbolic object, and check input matrix.  Not user callable. */
00016 
00017 KLU_symbolic *KLU_alloc_symbolic
00018 (
00019     Int n,
00020     Int *Ap,
00021     Int *Ai,
00022     KLU_common *Common
00023 )
00024 {
00025     KLU_symbolic *Symbolic ;
00026     Int *P, *Q, *R ;
00027     double *Lnz ;
00028     Int nz, i, j, p, pend ;
00029 
00030     if (Common == NULL)
00031     {
00032   return (NULL) ;
00033     }
00034     Common->status = KLU_OK ;
00035 
00036     /* A is n-by-n, with n > 0.  Ap [0] = 0 and nz = Ap [n] >= 0 required.
00037      * Ap [j] <= Ap [j+1] must hold for all j = 0 to n-1.  Row indices in Ai
00038      * must be in the range 0 to n-1, and no duplicate entries can be present.
00039      * The list of row indices in each column of A need not be sorted.
00040      */
00041 
00042     if (n <= 0 || Ap == NULL || Ai == NULL)
00043     {
00044   /* Ap and Ai must be present, and n must be > 0 */
00045   Common->status = KLU_INVALID ;
00046   return (NULL) ;
00047     }
00048 
00049     nz = Ap [n] ;
00050     if (Ap [0] != 0 || nz < 0)
00051     {
00052   /* nz must be >= 0 and Ap [0] must equal zero */
00053   Common->status = KLU_INVALID ;
00054   return (NULL) ;
00055     }
00056 
00057     for (j = 0 ; j < n ; j++)
00058     {
00059   if (Ap [j] > Ap [j+1])
00060   {
00061       /* column pointers must be non-decreasing */
00062       Common->status = KLU_INVALID ;
00063       return (NULL) ;
00064   }
00065     }
00066     P = KLU_malloc (n, sizeof (Int), Common) ;
00067     if (Common->status < KLU_OK)
00068     {
00069   /* out of memory */
00070   Common->status = KLU_OUT_OF_MEMORY ;
00071   return (NULL) ;
00072     }
00073     for (i = 0 ; i < n ; i++)
00074     {
00075   P [i] = EMPTY ;
00076     }
00077     for (j = 0 ; j < n ; j++)
00078     {
00079   pend = Ap [j+1] ;
00080   for (p = Ap [j] ; p < pend ; p++)
00081   {
00082       i = Ai [p] ;
00083       if (i < 0 || i >= n || P [i] == j)
00084       {
00085     /* row index out of range, or duplicate entry */
00086     KLU_free (P, n, sizeof (Int), Common) ;
00087     Common->status = KLU_INVALID ;
00088     return (NULL) ;
00089       }
00090       /* flag row i as appearing in column j */
00091       P [i] = j ;
00092   }
00093     }
00094 
00095     /* ---------------------------------------------------------------------- */
00096     /* allocate the Symbolic object */
00097     /* ---------------------------------------------------------------------- */
00098 
00099     Symbolic = KLU_malloc (sizeof (KLU_symbolic), 1, Common) ;
00100     if (Common->status < KLU_OK)
00101     {
00102   /* out of memory */
00103   KLU_free (P, n, sizeof (Int), Common) ;
00104   Common->status = KLU_OUT_OF_MEMORY ;
00105   return (NULL) ;
00106     }
00107 
00108     Q = KLU_malloc (n, sizeof (Int), Common) ;
00109     R = KLU_malloc (n+1, sizeof (Int), Common) ;
00110     Lnz = KLU_malloc (n, sizeof (double), Common) ;
00111 
00112     Symbolic->n = n ;
00113     Symbolic->nz = nz ;
00114     Symbolic->P = P ;
00115     Symbolic->Q = Q ;
00116     Symbolic->R = R ;
00117     Symbolic->Lnz = Lnz ;
00118 
00119     if (Common->status < KLU_OK)
00120     {
00121   /* out of memory */
00122   KLU_free_symbolic (&Symbolic, Common) ;
00123   Common->status = KLU_OUT_OF_MEMORY ;
00124   return (NULL) ;
00125     }
00126 
00127     return (Symbolic) ;
00128 }
00129 
00130 
00131 /* ========================================================================== */
00132 /* === KLU_analyze_given ==================================================== */
00133 /* ========================================================================== */
00134 
00135 KLU_symbolic *KLU_analyze_given     /* returns NULL if error, or a valid
00136                KLU_symbolic object if successful */
00137 (
00138     /* inputs, not modified */
00139     Int n,    /* A is n-by-n */
00140     Int Ap [ ],   /* size n+1, column pointers */
00141     Int Ai [ ],   /* size nz, row indices */
00142     Int Puser [ ],  /* size n, user's row permutation (may be NULL) */
00143     Int Quser [ ],  /* size n, user's column permutation (may be NULL) */
00144     /* -------------------- */
00145     KLU_common *Common
00146 )
00147 {
00148     KLU_symbolic *Symbolic ;
00149     double *Lnz ;
00150     Int nblocks, nz, block, maxblock, *P, *Q, *R, nzoff, p, pend, do_btf, k ;
00151 
00152     /* ---------------------------------------------------------------------- */
00153     /* determine if input matrix is valid, and get # of nonzeros */
00154     /* ---------------------------------------------------------------------- */
00155 
00156     Symbolic = KLU_alloc_symbolic (n, Ap, Ai, Common) ;
00157     if (Symbolic == NULL)
00158     {
00159   return (NULL) ;
00160     }
00161     P = Symbolic->P ;
00162     Q = Symbolic->Q ;
00163     R = Symbolic->R ;
00164     Lnz = Symbolic->Lnz ;
00165     nz = Symbolic->nz ;
00166 
00167     /* ---------------------------------------------------------------------- */
00168     /* Q = Quser, or identity if Quser is NULL */
00169     /* ---------------------------------------------------------------------- */
00170 
00171     if (Quser == (Int *) NULL)
00172     {
00173   for (k = 0 ; k < n ; k++)
00174   {
00175       Q [k] = k ;
00176   }
00177     }
00178     else
00179     {
00180   for (k = 0 ; k < n ; k++)
00181   {
00182       Q [k] = Quser [k] ;
00183   }
00184     }
00185 
00186     /* ---------------------------------------------------------------------- */
00187     /* get the control parameters for BTF and ordering method */
00188     /* ---------------------------------------------------------------------- */
00189 
00190     do_btf = Common->btf ;
00191     do_btf = (do_btf) ? TRUE : FALSE ;
00192     Symbolic->ordering = 2 ;
00193     Symbolic->do_btf = do_btf ;
00194 
00195     /* ---------------------------------------------------------------------- */
00196     /* find the block triangular form, if requested */
00197     /* ---------------------------------------------------------------------- */
00198 
00199     if (do_btf)
00200     {
00201 
00202   /* ------------------------------------------------------------------ */
00203   /* get workspace for BTF_strongcomp */
00204   /* ------------------------------------------------------------------ */
00205 
00206   Int *Pinv, *Work, *Bi, k1, k2, nk, oldcol ;
00207 
00208   Work = KLU_malloc (4*n, sizeof (Int), Common) ;
00209   Pinv = KLU_malloc (n, sizeof (Int), Common) ;
00210   if (Puser != (Int *) NULL)
00211   {
00212       Bi = KLU_malloc (nz+1, sizeof (Int), Common) ;
00213   }
00214   else
00215   {
00216       Bi = Ai ;
00217   }
00218 
00219   if (Common->status < KLU_OK)
00220   {
00221       /* out of memory */
00222       KLU_free (Work, 4*n, sizeof (Int), Common) ;
00223       KLU_free (Pinv, n, sizeof (Int), Common) ;
00224       if (Puser != (Int *) NULL)
00225       {
00226     KLU_free (Bi, nz+1, sizeof (Int), Common) ;
00227       }
00228       KLU_free_symbolic (&Symbolic, Common) ;
00229       Common->status = KLU_OUT_OF_MEMORY ;
00230       return (NULL) ;
00231   }
00232 
00233   /* ------------------------------------------------------------------ */
00234   /* B = Puser * A */
00235   /* ------------------------------------------------------------------ */
00236 
00237   if (Puser != (Int *) NULL)
00238   {
00239       for (k = 0 ; k < n ; k++)
00240       {
00241     Pinv [Puser [k]] = k ;
00242       }
00243       for (p = 0 ; p < nz ; p++)
00244       {
00245     Bi [p] = Pinv [Ai [p]] ;
00246       }
00247   }
00248 
00249   /* ------------------------------------------------------------------ */
00250   /* find the strongly-connected components */
00251   /* ------------------------------------------------------------------ */
00252 
00253   /* modifies Q, and determines P and R */
00254   nblocks = BTF_strongcomp (n, Ap, Bi, Q, P, R, Work) ;
00255 
00256   /* ------------------------------------------------------------------ */
00257   /* P = P * Puser */
00258   /* ------------------------------------------------------------------ */
00259 
00260   if (Puser != (Int *) NULL)
00261   {
00262       for (k = 0 ; k < n ; k++)
00263       {
00264     Work [k] = Puser [P [k]] ;
00265       }
00266       for (k = 0 ; k < n ; k++)
00267       {
00268     P [k] = Work [k] ;
00269       }
00270   }
00271 
00272   /* ------------------------------------------------------------------ */
00273   /* Pinv = inverse of P */
00274   /* ------------------------------------------------------------------ */
00275 
00276   for (k = 0 ; k < n ; k++)
00277   {
00278       Pinv [P [k]] = k ;
00279   }
00280 
00281   /* ------------------------------------------------------------------ */
00282   /* analyze each block */
00283   /* ------------------------------------------------------------------ */
00284 
00285   nzoff = 0 ;     /* nz in off-diagonal part */
00286   maxblock = 1 ;      /* size of the largest block */
00287 
00288   for (block = 0 ; block < nblocks ; block++)
00289   {
00290 
00291       /* -------------------------------------------------------------- */
00292       /* the block is from rows/columns k1 to k2-1 */
00293       /* -------------------------------------------------------------- */
00294 
00295       k1 = R [block] ;
00296       k2 = R [block+1] ;
00297       nk = k2 - k1 ;
00298       PRINTF (("BLOCK %d, k1 %d k2-1 %d nk %d\n", block, k1, k2-1, nk)) ;
00299       maxblock = MAX (maxblock, nk) ;
00300 
00301       /* -------------------------------------------------------------- */
00302       /* scan the kth block, C */
00303       /* -------------------------------------------------------------- */
00304 
00305       for (k = k1 ; k < k2 ; k++)
00306       {
00307     oldcol = Q [k] ;
00308     pend = Ap [oldcol+1] ;
00309     for (p = Ap [oldcol] ; p < pend ; p++)
00310     {
00311         if (Pinv [Ai [p]] < k1)
00312         {
00313       nzoff++ ;
00314         }
00315     }
00316       }
00317 
00318       /* fill-in not estimated */
00319       Lnz [block] = EMPTY ;
00320   }
00321 
00322   /* ------------------------------------------------------------------ */
00323   /* free all workspace */
00324   /* ------------------------------------------------------------------ */
00325 
00326   KLU_free (Work, 4*n, sizeof (Int), Common) ;
00327   KLU_free (Pinv, n, sizeof (Int), Common) ;
00328   if (Puser != (Int *) NULL)
00329   {
00330       KLU_free (Bi, nz+1, sizeof (Int), Common) ;
00331   }
00332 
00333     }
00334     else
00335     {
00336 
00337   /* ------------------------------------------------------------------ */
00338   /* BTF not requested */
00339   /* ------------------------------------------------------------------ */
00340 
00341   nzoff = 0 ;
00342   nblocks = 1 ;
00343   maxblock = n ;
00344   R [0] = 0 ;
00345   R [1] = n ;
00346   Lnz [0] = EMPTY ;
00347 
00348   /* ------------------------------------------------------------------ */
00349   /* P = Puser, or identity if Puser is NULL */
00350   /* ------------------------------------------------------------------ */
00351 
00352   for (k = 0 ; k < n ; k++)
00353   {
00354       P [k] = (Puser == NULL) ? k : Puser [k] ;
00355   }
00356     }
00357 
00358     /* ---------------------------------------------------------------------- */
00359     /* return the symbolic object */
00360     /* ---------------------------------------------------------------------- */
00361 
00362     Symbolic->nblocks = nblocks ;
00363     Symbolic->maxblock = maxblock ;
00364     Symbolic->lnz = EMPTY ;
00365     Symbolic->unz = EMPTY ;
00366     Symbolic->nzoff = nzoff ;
00367 
00368     return (Symbolic) ;
00369 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines