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