Amesos Package Browser (Single Doxygen Collection) Development
amesos_klu_analyze.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === klu_analyze ========================================================== */
00003 /* ========================================================================== */
00004 
00005 /* Order the matrix using BTF (or not), and then AMD, COLAMD, the natural
00006  * ordering, or the user-provided-function on the blocks.  Does not support
00007  * using a given ordering (use klu_analyze_given for that case). */
00008 
00009 #include "amesos_klu_internal.h"
00010 
00011 /* ========================================================================== */
00012 /* === analyze_worker ======================================================= */
00013 /* ========================================================================== */
00014 
00015 static Int analyze_worker /* returns KLU_OK or < 0 if error */
00016 (
00017     /* inputs, not modified */
00018     Int n,    /* A is n-by-n */
00019     Int Ap [ ],   /* size n+1, column pointers */
00020     Int Ai [ ],   /* size nz, row indices */
00021     Int nblocks,  /* # of blocks */
00022     Int Pbtf [ ], /* BTF row permutation */
00023     Int Qbtf [ ], /* BTF col permutation */
00024     Int R [ ],    /* size n+1, but only Rbtf [0..nblocks] is used */
00025     Int ordering, /* what ordering to use (0, 1, or 3 for this routine) */
00026 
00027     /* output only, not defined on input */
00028     Int P [ ],    /* size n */
00029     Int Q [ ],    /* size n */
00030     double Lnz [ ], /* size n, but only Lnz [0..nblocks-1] is used */
00031 
00032     /* workspace, not defined on input or output */
00033     Int Pblk [ ], /* size maxblock */
00034     Int Cp [ ],   /* size maxblock+1 */
00035     Int Ci [ ],   /* size MAX (nz+1, Cilen) */
00036     Int Cilen,    /* nz+1, or COLAMD_recommend(nz,n,n) for COLAMD */
00037     Int Pinv [ ], /* size maxblock */
00038 
00039     /* input/output */
00040     KLU_symbolic *Symbolic,
00041     KLU_common *Common
00042 )
00043 {
00044     double amd_Info [AMD_INFO], lnz, lnz1, flops, flops1 ;
00045     Int k1, k2, nk, k, block, oldcol, pend, newcol, result, pc, p, newrow,
00046   maxnz, nzoff, cstats [COLAMD_STATS], ok, err = KLU_INVALID ;
00047 
00048     /* ---------------------------------------------------------------------- */
00049     /* initializations */
00050     /* ---------------------------------------------------------------------- */
00051 
00052     /* compute the inverse of Pbtf */
00053 #ifndef NDEBUG
00054     for (k = 0 ; k < n ; k++)
00055     {
00056   P [k] = EMPTY ;
00057   Q [k] = EMPTY ;
00058   Pinv [k] = EMPTY ;
00059     }
00060 #endif
00061     for (k = 0 ; k < n ; k++)
00062     {
00063   ASSERT (Pbtf [k] >= 0 && Pbtf [k] < n) ;
00064   Pinv [Pbtf [k]] = k ;
00065     }
00066 #ifndef NDEBUG
00067     for (k = 0 ; k < n ; k++) ASSERT (Pinv [k] != EMPTY) ;
00068 #endif
00069     nzoff = 0 ;
00070     lnz = 0 ;
00071     maxnz = 0 ;
00072     flops = 0 ;
00073     Symbolic->symmetry = EMPTY ;  /* only computed by AMD */
00074 
00075     /* ---------------------------------------------------------------------- */
00076     /* order each block */
00077     /* ---------------------------------------------------------------------- */
00078 
00079     for (block = 0 ; block < nblocks ; block++)
00080     {
00081 
00082   /* ------------------------------------------------------------------ */
00083   /* the block is from rows/columns k1 to k2-1 */
00084   /* ------------------------------------------------------------------ */
00085 
00086   k1 = R [block] ;
00087   k2 = R [block+1] ;
00088   nk = k2 - k1 ;
00089   PRINTF (("BLOCK %d, k1 %d k2-1 %d nk %d\n", block, k1, k2-1, nk)) ;
00090 
00091   /* ------------------------------------------------------------------ */
00092   /* construct the kth block, C */
00093   /* ------------------------------------------------------------------ */
00094 
00095   Lnz [block] = EMPTY ;
00096   pc = 0 ;
00097   for (k = k1 ; k < k2 ; k++)
00098   {
00099       newcol = k-k1 ;
00100       Cp [newcol] = pc ;
00101       oldcol = Qbtf [k] ;
00102       pend = Ap [oldcol+1] ;
00103       for (p = Ap [oldcol] ; p < pend ; p++)
00104       {
00105     newrow = Pinv [Ai [p]] ;
00106     if (newrow < k1)
00107     {
00108         nzoff++ ;
00109     }
00110     else
00111     {
00112         /* (newrow,newcol) is an entry in the block */
00113         ASSERT (newrow < k2) ;
00114         newrow -= k1 ;
00115         Ci [pc++] = newrow ;
00116     }
00117       }
00118   }
00119   Cp [nk] = pc ;
00120   maxnz = MAX (maxnz, pc) ;
00121   ASSERT (KLU_valid (nk, Cp, Ci, NULL)) ;
00122 
00123   /* ------------------------------------------------------------------ */
00124   /* order the block C */
00125   /* ------------------------------------------------------------------ */
00126 
00127   if (nk <= 3)
00128   {
00129 
00130       /* -------------------------------------------------------------- */
00131       /* use natural ordering for tiny blocks (3-by-3 or less) */
00132       /* -------------------------------------------------------------- */
00133 
00134       for (k = 0 ; k < nk ; k++)
00135       {
00136     Pblk [k] = k ;
00137       }
00138       lnz1 = nk * (nk + 1) / 2 ;
00139       flops1 = nk * (nk - 1) / 2 + (nk-1)*nk*(2*nk-1) / 6 ;
00140       ok = TRUE ;
00141 
00142   }
00143   else if (ordering == 0)
00144   {
00145 
00146       /* -------------------------------------------------------------- */
00147       /* order the block with AMD (C+C') */
00148       /* -------------------------------------------------------------- */
00149 
00150       result = AMD_order (nk, Cp, Ci, Pblk, NULL, amd_Info) ;
00151       ok = (result >= AMD_OK) ;
00152       if (result == AMD_OUT_OF_MEMORY)
00153       {
00154     err = KLU_OUT_OF_MEMORY ;
00155       }
00156 
00157       /* account for memory usage in AMD */
00158       Common->mempeak = ( size_t) ( MAX (Common->mempeak,
00159     Common->memusage + amd_Info [AMD_MEMORY]) ) ;
00160 
00161       /* get the ordering statistics from AMD */
00162       lnz1 = (Int) (amd_Info [AMD_LNZ]) + nk ;
00163       flops1 = 2 * amd_Info [AMD_NMULTSUBS_LU] + amd_Info [AMD_NDIV] ;
00164       if (pc == maxnz)
00165       {
00166     /* get the symmetry of the biggest block */
00167     Symbolic->symmetry = amd_Info [AMD_SYMMETRY] ;
00168       }
00169 
00170   }
00171   else if (ordering == 1)
00172   {
00173 
00174       /* -------------------------------------------------------------- */
00175       /* order the block with COLAMD (C) */
00176       /* -------------------------------------------------------------- */
00177 
00178       /* order (and destroy) Ci, returning column permutation in Cp.
00179        * COLAMD "cannot" fail since the matrix has already been checked,
00180        * and Ci allocated. */
00181 
00182       ok = COLAMD (nk, nk, Cilen, Ci, Cp, NULL, cstats) ;
00183       lnz1 = EMPTY ;
00184       flops1 = EMPTY ;
00185 
00186       /* copy the permutation from Cp to Pblk */
00187       for (k = 0 ; k < nk ; k++)
00188       {
00189     Pblk [k] = Cp [k] ;
00190       }
00191 
00192   }
00193   else
00194   {
00195 
00196       /* -------------------------------------------------------------- */
00197       /* pass the block to the user-provided ordering function */
00198       /* -------------------------------------------------------------- */
00199 
00200       lnz1 = (Common->user_order) (nk, Cp, Ci, Pblk, Common) ;
00201       flops1 = EMPTY ;
00202       ok = (lnz1 != 0) ;
00203   }
00204 
00205   if (!ok)
00206   {
00207       return (err) ;  /* ordering method failed */
00208   }
00209 
00210   /* ------------------------------------------------------------------ */
00211   /* keep track of nnz(L) and flops statistics */
00212   /* ------------------------------------------------------------------ */
00213 
00214   Lnz [block] = lnz1 ;
00215   lnz = (lnz == EMPTY || lnz1 == EMPTY) ? EMPTY : (lnz + lnz1) ;
00216   flops = (flops == EMPTY || flops1 == EMPTY) ? EMPTY : (flops + flops1) ;
00217 
00218   /* ------------------------------------------------------------------ */
00219   /* combine the preordering with the BTF ordering */
00220   /* ------------------------------------------------------------------ */
00221 
00222   PRINTF (("Pblk, 1-based:\n")) ;
00223   for (k = 0 ; k < nk ; k++)
00224   {
00225       ASSERT (k + k1 < n) ;
00226       ASSERT (Pblk [k] + k1 < n) ;
00227       Q [k + k1] = Qbtf [Pblk [k] + k1] ;
00228   }
00229   for (k = 0 ; k < nk ; k++)
00230   {
00231       ASSERT (k + k1 < n) ;
00232       ASSERT (Pblk [k] + k1 < n) ;
00233       P [k + k1] = Pbtf [Pblk [k] + k1] ;
00234   }
00235     }
00236 
00237     PRINTF (("nzoff %d  Ap[n] %d\n", nzoff, Ap [n])) ;
00238     ASSERT (nzoff >= 0 && nzoff <= Ap [n]) ;
00239 
00240     /* return estimates of # of nonzeros in L including diagonal */
00241     Symbolic->lnz = lnz ;     /* EMPTY if COLAMD used */
00242     Symbolic->unz = lnz ;
00243     Symbolic->nzoff = nzoff ;
00244     Symbolic->est_flops = flops ;   /* EMPTY if COLAMD or user-ordering used */
00245     return (KLU_OK) ;
00246 }
00247 
00248 
00249 /* ========================================================================== */
00250 /* === order_and_analyze ==================================================== */
00251 /* ========================================================================== */
00252 
00253 /* Orders the matrix with or with BTF, then orders each block with AMD, COLAMD,
00254  * or the user ordering function.  Does not handle the natural or given
00255  * ordering cases. */
00256 
00257 static KLU_symbolic *order_and_analyze  /* returns NULL if error, or a valid
00258              KLU_symbolic object if successful */
00259 (
00260     /* inputs, not modified */
00261     Int n,    /* A is n-by-n */
00262     Int Ap [ ],   /* size n+1, column pointers */
00263     Int Ai [ ],   /* size nz, row indices */
00264     /* --------------------- */
00265     KLU_common *Common
00266 )
00267 {
00268     double work ;
00269     KLU_symbolic *Symbolic ;
00270     double *Lnz ;
00271     Int *Qbtf, *Cp, *Ci, *Pinv, *Pblk, *Pbtf, *P, *Q, *R ;
00272     Int nblocks, nz, block, maxblock, k1, k2, nk, do_btf, ordering, k, Cilen,
00273   *Work ;
00274 
00275     /* ---------------------------------------------------------------------- */
00276     /* allocate the Symbolic object, and check input matrix */
00277     /* ---------------------------------------------------------------------- */
00278 
00279     Symbolic = KLU_alloc_symbolic (n, Ap, Ai, Common) ;
00280     if (Symbolic == NULL)
00281     {
00282   return (NULL) ;
00283     }
00284     P = Symbolic->P ;
00285     Q = Symbolic->Q ;
00286     R = Symbolic->R ;
00287     Lnz = Symbolic->Lnz ;
00288     nz = Symbolic->nz ;
00289 
00290     ordering = Common->ordering ;
00291     if (ordering == 1)
00292     {
00293   /* COLAMD */
00294   Cilen = COLAMD_recommended (nz, n, n) ;
00295     }
00296     else if (ordering == 0 || (ordering == 3 && Common->user_order != NULL))
00297     {
00298   /* AMD or user ordering function */
00299   Cilen = nz+1 ;
00300     }
00301     else
00302     {
00303   /* invalid ordering */
00304   Common->status = KLU_INVALID ;
00305   KLU_free_symbolic (&Symbolic, Common) ;
00306   return (NULL) ;
00307     }
00308 
00309     /* AMD memory management routines */
00310     amesos_amd_malloc  = Common->malloc_memory ;
00311     amesos_amd_free    = Common->free_memory ;
00312     amesos_amd_calloc  = Common->calloc_memory ;
00313     amesos_amd_realloc = Common->realloc_memory ;
00314 
00315     /* ---------------------------------------------------------------------- */
00316     /* allocate workspace for BTF permutation */
00317     /* ---------------------------------------------------------------------- */
00318 
00319     Pbtf = KLU_malloc (n, sizeof (Int), Common) ;
00320     Qbtf = KLU_malloc (n, sizeof (Int), Common) ;
00321     if (Common->status < KLU_OK)
00322     {
00323   KLU_free (Pbtf, n, sizeof (Int), Common) ;
00324   KLU_free (Qbtf, n, sizeof (Int), Common) ;
00325   KLU_free_symbolic (&Symbolic, Common) ;
00326   return (NULL) ;
00327     }
00328 
00329     /* ---------------------------------------------------------------------- */
00330     /* get the common parameters for BTF and ordering method */
00331     /* ---------------------------------------------------------------------- */
00332 
00333     do_btf = Common->btf ;
00334     do_btf = (do_btf) ? TRUE : FALSE ;
00335     Symbolic->ordering = ordering ;
00336     Symbolic->do_btf = do_btf ;
00337     Symbolic->structural_rank = EMPTY ;
00338 
00339     /* ---------------------------------------------------------------------- */
00340     /* find the block triangular form (if requested) */
00341     /* ---------------------------------------------------------------------- */
00342 
00343     Common->work = 0 ;
00344 
00345     if (do_btf)
00346     {
00347   Work = KLU_malloc (5*n, sizeof (Int), Common) ;
00348   if (Common->status < KLU_OK)
00349   {
00350       /* out of memory */
00351       KLU_free (Pbtf, n, sizeof (Int), Common) ;
00352       KLU_free (Qbtf, n, sizeof (Int), Common) ;
00353       KLU_free_symbolic (&Symbolic, Common) ;
00354       return (NULL) ;
00355   }
00356 
00357   nblocks = BTF_order (n, Ap, Ai, Common->maxwork, &work, Pbtf, Qbtf, R,
00358     &(Symbolic->structural_rank), Work) ;
00359   Common->structural_rank = Symbolic->structural_rank ;
00360   Common->work += work ;
00361 
00362   KLU_free (Work, 5*n, sizeof (Int), Common) ;
00363 
00364   /* unflip Qbtf if the matrix does not have full structural rank */
00365   if (Symbolic->structural_rank < n)
00366   {
00367       for (k = 0 ; k < n ; k++)
00368       {
00369     Qbtf [k] = BTF_UNFLIP (Qbtf [k]) ;
00370       }
00371   }
00372 
00373   /* find the size of the largest block */
00374   maxblock = 1 ;
00375   for (block = 0 ; block < nblocks ; block++)
00376   {
00377       k1 = R [block] ;
00378       k2 = R [block+1] ;
00379       nk = k2 - k1 ;
00380       PRINTF (("block %d size %d\n", block, nk)) ;
00381       maxblock = MAX (maxblock, nk) ;
00382   }
00383     }
00384     else
00385     {
00386   /* BTF not requested */
00387   nblocks = 1 ;
00388   maxblock = n ;
00389   R [0] = 0 ;
00390   R [1] = n ;
00391   for (k = 0 ; k < n ; k++)
00392   {
00393       Pbtf [k] = k ;
00394       Qbtf [k] = k ;
00395   }
00396     }
00397 
00398     Symbolic->nblocks = nblocks ;
00399 
00400     PRINTF (("maxblock size %d\n", maxblock)) ;
00401     Symbolic->maxblock = maxblock ;
00402 
00403     /* ---------------------------------------------------------------------- */
00404     /* allocate more workspace, for analyze_worker */
00405     /* ---------------------------------------------------------------------- */
00406 
00407     Pblk = KLU_malloc (maxblock, sizeof (Int), Common) ;
00408     Cp   = KLU_malloc (maxblock + 1, sizeof (Int), Common) ;
00409     Ci   = KLU_malloc (MAX (Cilen, nz+1), sizeof (Int), Common) ;
00410     Pinv = KLU_malloc (n, sizeof (Int), Common) ;
00411 
00412     /* ---------------------------------------------------------------------- */
00413     /* order each block of the BTF ordering, and a fill-reducing ordering */
00414     /* ---------------------------------------------------------------------- */
00415 
00416     if (Common->status == KLU_OK)
00417     {
00418   PRINTF (("calling analyze_worker\n")) ;
00419   Common->status = analyze_worker (n, Ap, Ai, nblocks, Pbtf, Qbtf, R,
00420       ordering, P, Q, Lnz, Pblk, Cp, Ci, Cilen, Pinv, Symbolic, Common) ;
00421   PRINTF (("analyze_worker done\n")) ;
00422     }
00423 
00424     /* ---------------------------------------------------------------------- */
00425     /* free all workspace */
00426     /* ---------------------------------------------------------------------- */
00427 
00428     KLU_free (Pblk, maxblock, sizeof (Int), Common) ;
00429     KLU_free (Cp, maxblock+1, sizeof (Int), Common) ;
00430     KLU_free (Ci, MAX (Cilen, nz+1), sizeof (Int), Common) ;
00431     KLU_free (Pinv, n, sizeof (Int), Common) ;
00432     KLU_free (Pbtf, n, sizeof (Int), Common) ;
00433     KLU_free (Qbtf, n, sizeof (Int), Common) ;
00434 
00435     /* ---------------------------------------------------------------------- */
00436     /* return the symbolic object */
00437     /* ---------------------------------------------------------------------- */
00438 
00439     if (Common->status < KLU_OK)
00440     {
00441   KLU_free_symbolic (&Symbolic, Common) ;
00442     }
00443     return (Symbolic) ;
00444 }
00445 
00446 
00447 /* ========================================================================== */
00448 /* === KLU_analyze ========================================================== */
00449 /* ========================================================================== */
00450 
00451 KLU_symbolic *KLU_analyze /* returns NULL if error, or a valid
00452            KLU_symbolic object if successful */
00453 (
00454     /* inputs, not modified */
00455     Int n,    /* A is n-by-n */
00456     Int Ap [ ],   /* size n+1, column pointers */
00457     Int Ai [ ],   /* size nz, row indices */
00458     /* -------------------- */
00459     KLU_common *Common
00460 )
00461 {
00462 
00463     /* ---------------------------------------------------------------------- */
00464     /* get the control parameters for BTF and ordering method */
00465     /* ---------------------------------------------------------------------- */
00466 
00467     if (Common == NULL)
00468     {
00469   return (NULL) ;
00470     }
00471     Common->status = KLU_OK ;
00472     Common->structural_rank = EMPTY ;
00473 
00474     /* ---------------------------------------------------------------------- */
00475     /* order and analyze */
00476     /* ---------------------------------------------------------------------- */
00477 
00478     if (Common->ordering == 2)
00479     {
00480   /* natural ordering */
00481   return (KLU_analyze_given (n, Ap, Ai, NULL, NULL, Common)) ;
00482     }
00483     else
00484     {
00485   /* order with P and Q */
00486   return (order_and_analyze (n, Ap, Ai, Common)) ;
00487     }
00488 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines