amesos_paraklete_btf.c

Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === Paraklete/paraklete_btf.c ============================================ */
00003 /* ========================================================================== */
00004 
00005 /* User-callable functions that combine the functions of Paraklete,
00006  * KLU, and BTF.  See pb.c for examples on how to use these functions.
00007  *
00008  * PARAKLETE version 0.3: parallel sparse LU factorization.  Nov 13, 2007
00009  * Copyright (C) 2007, Univ. of Florida.  Author: Timothy A. Davis
00010  */
00011 
00012 #include "amesos_paraklete_decl.h"
00013 
00014 /* ========================================================================== */
00015 /* === paraklete_btf_bcast_symbolic ========================================= */
00016 /* ========================================================================== */
00017 
00018 #ifndef NMPI
00019 static Int paraklete_btf_bcast_symbolic
00020 (
00021     paraklete_btf_symbolic **LU_btf_symbolic_handle,
00022     paraklete_common *Common
00023 )
00024 {
00025     paraklete_btf_symbolic *LU_btf_symbolic = NULL ;
00026     Int n, nblocks, header [4] ;
00027     int ok, all_ok ;
00028     cholmod_common *cm ;
00029     KLU_common *km ;
00030 
00031     cm = &(Common->cm) ;
00032     km = &(Common->km) ;
00033 
00034     n = EMPTY ;
00035     nblocks = EMPTY ;
00036 
00037     /* broadcast number of diagonal blocks in the BTF form, or -1 if failure */
00038     if (Common->myid == 0)
00039     {
00040   LU_btf_symbolic = *LU_btf_symbolic_handle ;
00041   if (LU_btf_symbolic != NULL)
00042   {
00043       n = LU_btf_symbolic->n ;
00044       nblocks = LU_btf_symbolic->nblocks ;
00045   }
00046     }
00047     else
00048     {
00049   /* other processes do not yet have the BTF symbolic object */
00050   *LU_btf_symbolic_handle = NULL ;
00051     }
00052 
00053     /* broadcast the size of the object, or -1 if a failure occurred */
00054     header [0] = n ;
00055     header [1] = nblocks ;
00056     MPI (MPI_Bcast (&header, 2, MPI_Int, TAG0, MPI_COMM_WORLD)) ;
00057     n = header [0] ;
00058     nblocks = header [1] ;
00059     if (n == EMPTY)
00060     {
00061   return (FALSE) ;
00062     }
00063 
00064     if (Common->myid != 0)
00065     {
00066   /* allocate the copy of this analysis on the slave */
00067   LU_btf_symbolic = amesos_paraklete_btf_alloc_symbolic (n, nblocks, Common) ;
00068   *LU_btf_symbolic_handle = LU_btf_symbolic ;
00069         /* TODO return if failed, and broadcast error to all processes */
00070         /* PARAKLETE_ERROR already called */
00071     }
00072 
00073     /* all processes find out if any one process fails to allocate memory */
00074     ok = (cm->status == CHOLMOD_OK) && (LU_btf_symbolic != NULL) ;
00075     all_ok = ok ;
00076     MPI (MPI_Allreduce (&ok, &all_ok, 1, MPI_INT, MPI_LAND, MPI_COMM_WORLD)) ;
00077     if (!all_ok)
00078     {
00079   /* out of memory; inform all processes */
00080   PR0 ((Common->file, "proc "ID" all fail in analyze\n", Common->myid)) ;
00081   amesos_paraklete_btf_free_symbolic (&LU_btf_symbolic, Common) ;
00082   *LU_btf_symbolic_handle = NULL ;
00083   return (FALSE) ;
00084     }
00085 
00086     /* broadcast the contents of the LU_btf_symbolic object (excluding
00087      * the LUsymbolic analysis of each diagonal block) */
00088     MPI (MPI_Bcast (LU_btf_symbolic->Mem_n, 3*n+1, MPI_Int, TAG0,
00089         MPI_COMM_WORLD)) ;
00090     return (TRUE) ;
00091 }
00092 #endif
00093 
00094 
00095 /* ========================================================================== */
00096 /* === paraklete_btf_alloc_symbolic ========================================= */
00097 /* ========================================================================== */
00098 
00099 paraklete_btf_symbolic *amesos_paraklete_btf_alloc_symbolic
00100 (
00101     Int n,
00102     Int nblocks,
00103     paraklete_common *Common
00104 )
00105 {
00106     cholmod_common *cm ;
00107     KLU_common *km ;
00108     paraklete_btf_symbolic *LU_btf_symbolic ;
00109     Int *p ;
00110     int ok = TRUE ;
00111     size_t n3 ;
00112 
00113     cm = &(Common->cm) ;
00114     km = &(Common->km) ;
00115 
00116     /* n3 = 3*n+1 */
00117     n3 = CHOLMOD (mult_size_t) (n, 3, &ok) ;
00118     n3 = CHOLMOD (add_size_t) (n3, 1, &ok) ;
00119     if (!ok) PARAKLETE_ERROR (PK_TOO_LARGE, "problem too large") ;
00120 
00121     LU_btf_symbolic = CHOLMOD (malloc) (1, sizeof (paraklete_btf_symbolic), cm) ;
00122     if (cm->status != CHOLMOD_OK)
00123     {
00124         PARAKLETE_ERROR (PK_OUT_OF_MEMORY, "out of memory") ;
00125   return (NULL) ;
00126     }
00127 
00128     LU_btf_symbolic->n = n ;
00129     LU_btf_symbolic->nblocks = nblocks ;
00130     LU_btf_symbolic->cnz = 0 ;
00131     LU_btf_symbolic->fnz = 0 ;
00132     p = LU_btf_symbolic->Mem_n = CHOLMOD (calloc) (n3, sizeof (Int), cm) ;
00133     LU_btf_symbolic->Qbtf = p ;         /* size n */
00134     LU_btf_symbolic->Pbinv = p + n ;        /* size n */
00135     LU_btf_symbolic->Rbtf = p + 2*n ;       /* size n+1 */
00136     /* only nblocks of the LUsymbolic array is used: */
00137     LU_btf_symbolic->LUsymbolic = CHOLMOD (calloc) (n, sizeof (void *), cm) ;
00138 
00139     if (cm->status != CHOLMOD_OK)
00140     {
00141         /* TODO free object and return NULL if malloc fails */
00142         PARAKLETE_ERROR (PK_OUT_OF_MEMORY, "out of memory") ;
00143     }
00144     return (LU_btf_symbolic) ;
00145 }
00146 
00147 
00148 /* ========================================================================== */
00149 /* === paraklete_btf_analyze ================================================ */
00150 /* ========================================================================== */
00151 
00152 paraklete_btf_symbolic *amesos_paraklete_btf_analyze
00153 (
00154     cholmod_sparse *A,        /* matrix to analyze */ 
00155     paraklete_common *Common
00156 )
00157 {
00158     double work = 0 ;
00159     Int *Pbtf = NULL, *Qbtf = NULL, *Rbtf = NULL, *Work = NULL, *Ap = NULL,
00160   *Ai = NULL, *Cp = NULL, *Ci = NULL, *Pbinv = NULL ;
00161     cholmod_sparse *C = NULL ;
00162     cholmod_common *cm ;
00163     KLU_common *km ;
00164     paraklete_btf_symbolic *LU_btf_symbolic = NULL ;
00165     void **LUsymbolic ;
00166     Int nblocks, n = 0, nmatch, block, k, k1, k2, inew, jold, anz = 0, fnz,
00167   cnz, p, nk ;
00168 
00169     /* ---------------------------------------------------------------------- */
00170     /* root processor finds the BTF ordering */
00171     /* ---------------------------------------------------------------------- */
00172 
00173     cm = &(Common->cm) ;
00174     km = &(Common->km) ;
00175 
00176     if (Common->myid == 0)
00177     {
00178   Ap = A->p ;
00179   Ai = A->i ;
00180   n = A->nrow ;
00181   anz = Ap [n] ;
00182 
00183   LU_btf_symbolic = amesos_paraklete_btf_alloc_symbolic (n, n, Common) ;
00184         if (LU_btf_symbolic == NULL)
00185         {
00186             /* TODO return NULL if malloc fails */
00187             /* PARAKLETE_ERROR already called */ ;
00188         }
00189 
00190   Qbtf = LU_btf_symbolic->Qbtf ;        /* size n */
00191   Rbtf = LU_btf_symbolic->Rbtf ;        /* size n+1 */
00192   Pbinv = LU_btf_symbolic->Pbinv ;      /* size n */
00193   LUsymbolic = LU_btf_symbolic->LUsymbolic ;  /* size n */
00194 
00195   /* ------------------------------------------------------------------ */
00196   /* find the BTF ordering */
00197   /* ------------------------------------------------------------------ */
00198 
00199   Work = CHOLMOD (malloc) (n, 6*sizeof (Int), cm) ;
00200   Pbtf = Work + 5*n ;
00201         if (cm->status != CHOLMOD_OK)
00202         {
00203             /* TODO return if malloc fails */
00204             PARAKLETE_ERROR (PK_OUT_OF_MEMORY, "out of memory") ;
00205         }
00206 
00207   nblocks = BTF_order (n, Ap, Ai, 0, &work, Pbtf, Qbtf, Rbtf, &nmatch,
00208       Work) ;
00209 
00210   /* Pbinv = inverse of Pbtf */
00211   for (k = 0 ; k < n ; k++)
00212   {
00213       Pbinv [Pbtf [k]] = k ;
00214   }
00215 
00216   CHOLMOD (free) (n, 6*sizeof (Int), Work, cm) ;
00217   LU_btf_symbolic->nblocks = nblocks ;
00218     }
00219 
00220     /* ---------------------------------------------------------------------- */
00221     /* broadcast the BTF information */
00222     /* ---------------------------------------------------------------------- */
00223 
00224     MPI (paraklete_btf_bcast_symbolic (&LU_btf_symbolic, Common)) ;
00225     /* TODO return if broadcast fails */
00226     Rbtf = LU_btf_symbolic->Rbtf ;
00227     LUsymbolic = LU_btf_symbolic->LUsymbolic ;
00228     nblocks = LU_btf_symbolic->nblocks ;
00229 
00230     /* ---------------------------------------------------------------------- */
00231     /* symbolic analysis of diagonal blocks of A(p,q) */
00232     /* ---------------------------------------------------------------------- */
00233 
00234     if (Common->myid == 0)
00235     {
00236   /* C = empty n-by-n sparse matrix */
00237   fnz = 0 ;
00238   C = CHOLMOD (allocate_sparse) (n, n, anz, FALSE, TRUE, 0, CHOLMOD_PATTERN,
00239       cm) ;
00240         if (cm->status != CHOLMOD_OK)
00241         {
00242             /* TODO return if malloc fails */
00243             PARAKLETE_ERROR (PK_OUT_OF_MEMORY, "out of memory") ;
00244         }
00245   Cp = C->p ;
00246   Ci = C->i ;
00247   LU_btf_symbolic->cnz = 0 ;
00248     }
00249 
00250     for (block = 0 ; block < nblocks ; block++)
00251     {
00252   k1 = Rbtf [block] ;
00253   k2 = Rbtf [block+1] ;
00254   nk = k2 - k1 ;
00255 
00256   if (Common->myid == 0)
00257   {
00258       cnz = 0 ;
00259       for (k = k1 ; k < k2 ; k++)
00260       {
00261     Cp [k - k1] = cnz ;       /* note change of index */
00262     jold = Qbtf [k] ;
00263     for (p = Ap [jold] ; p < Ap [jold+1] ; p++)
00264     {
00265         inew = Pbinv [Ai [p]] ;
00266         if (inew < k1)
00267         {
00268       /* this entry goes in the off-diagonal part */
00269       fnz++ ;
00270         }
00271         else
00272         {
00273       /* this entry goes in the diagonal block */
00274       Ci [cnz++] = inew - k1 ;  /* note the change of index */
00275         }
00276     }
00277       }
00278       LU_btf_symbolic->cnz = MAX (LU_btf_symbolic->cnz, cnz) ;
00279       Cp [nk] = cnz ;
00280       C->nrow = nk ;
00281       C->ncol = nk ;
00282   }
00283 
00284   if (nk == 1)
00285   {
00286       /* singleton; nothing to do */
00287       LUsymbolic [block] = NULL ;
00288   }
00289         else if (nk < 1000)
00290         {
00291             if (Common->myid == 0)
00292             {
00293                 /* use KLU on process 0 */
00294                 LUsymbolic [block] = (void *) KLU_analyze (nk, Cp, Ci, km) ;
00295                 if (km->status != KLU_OK)
00296                 {
00297                     /* TODO return if KLU_analyze failed; broadcast the error */
00298                     PARAKLETE_ERROR (PK_UNKNOWN, "KLU analyze failed") ;
00299                 }
00300             }
00301             else
00302             {
00303                 /* small block; nothing to do for other processors*/
00304                 LUsymbolic [block] = NULL ;
00305             }
00306         }
00307   else
00308   {
00309             /* parallel factorization of a large block */
00310       /* process 0 does the work and broadcasts it to all the others */
00311       LUsymbolic [block] = (void *) amesos_paraklete_analyze (C, Common) ;
00312             if (LUsymbolic [block] == NULL)
00313             {
00314           /* TODO return if analyze fails and broadcast the error */
00315                 /* note that PARAKLETE_ERROR has already been called */
00316             }
00317   }
00318     }
00319 
00320     /* ---------------------------------------------------------------------- */
00321     /* free workspace sparse matrix C */
00322     /* ---------------------------------------------------------------------- */
00323 
00324     if (Common->myid == 0)
00325     {
00326         /* process 0 finalizes the object and frees its workspace */
00327   LU_btf_symbolic->fnz = fnz ;
00328   C->nrow = n ;
00329   C->ncol = n ;
00330   CHOLMOD (free_sparse) (&C, cm) ;
00331     }
00332 
00333     /* ---------------------------------------------------------------------- */
00334     /* return the result */
00335     /* ---------------------------------------------------------------------- */
00336 
00337     return (LU_btf_symbolic) ;
00338 }
00339 
00340 
00341 /* ========================================================================== */
00342 /* === paraklete_btf_factorize ============================================== */
00343 /* ========================================================================== */
00344 
00345 paraklete_btf_numeric *amesos_paraklete_btf_factorize
00346 (
00347     cholmod_sparse *A,        /* matrix to analyze */ 
00348     paraklete_btf_symbolic *LU_btf_symbolic,  /* symbolic analysis */
00349     paraklete_common *Common
00350 )
00351 {
00352     cholmod_sparse *F = NULL, *C = NULL ;
00353     cholmod_common *cm ;
00354     KLU_common *km ;
00355     double *Ax = NULL, *Cx = NULL, *Fx = NULL, *Singleton = NULL ;
00356     Int *Qbtf = NULL, *Rbtf = NULL, *Ap = NULL, *Ai = NULL, *Cp = NULL,
00357   *Ci = NULL, *Pbinv = NULL, *Fp = NULL, *Fi = NULL ;
00358     Int fnz = 0, cnz = 0, nblocks, n, k, k1, k2, block, jold, inew, p, nk ;
00359     int ok = TRUE, all_ok ;
00360     void **LUsymbolic = NULL;
00361     paraklete_btf_numeric *LU_btf_numeric = NULL ;
00362     void **LUnumeric = NULL ;
00363 
00364     /* ---------------------------------------------------------------------- */
00365     /* get inputs and allocate result */
00366     /* ---------------------------------------------------------------------- */
00367 
00368     cm = &(Common->cm) ;
00369     km = &(Common->km) ;
00370 
00371     LU_btf_numeric = CHOLMOD (malloc) (1, sizeof (paraklete_btf_numeric), cm) ;
00372     if (cm->status != CHOLMOD_OK)
00373     {
00374         PARAKLETE_ERROR (PK_OUT_OF_MEMORY, "out of memory") ;
00375   return (NULL) ;
00376     }
00377 
00378     if (Common->myid == 0)
00379     {
00380   fnz = LU_btf_symbolic->fnz ;
00381   cnz = LU_btf_symbolic->cnz ;
00382   n = A->nrow ;
00383   F = CHOLMOD (allocate_sparse) (n, n, fnz, FALSE, TRUE, 0,
00384       CHOLMOD_REAL, cm) ;
00385   C = CHOLMOD (allocate_sparse) (n, n, cnz, FALSE, TRUE, 0,
00386       CHOLMOD_REAL, cm) ;
00387         if (cm->status != CHOLMOD_OK)
00388         {
00389       /* TODO free LU_btf_numeric and return NULL;
00390              * broadcast error to all processes */
00391             PARAKLETE_ERROR (PK_OUT_OF_MEMORY, "out of memory") ;
00392         }
00393     }
00394 
00395     n = LU_btf_symbolic->n ;
00396     nblocks = LU_btf_symbolic->nblocks ;
00397     Rbtf = LU_btf_symbolic->Rbtf ;
00398     Qbtf = LU_btf_symbolic->Qbtf ;
00399     Pbinv = LU_btf_symbolic->Pbinv ;
00400     LUsymbolic = LU_btf_symbolic->LUsymbolic ;
00401 
00402     if (Common->myid == 0)
00403     {
00404   fnz = 0 ;
00405   Fp = F->p ;
00406   Fi = F->i ;
00407   Fx = F->x ;
00408   Ap = A->p ;
00409   Ai = A->i ;
00410   Ax = A->x ;
00411   Cp = C->p ;
00412   Ci = C->i ;
00413   Cx = C->x ;
00414   Singleton = CHOLMOD (calloc) (nblocks, sizeof (double), cm) ;
00415         if (cm->status != CHOLMOD_OK)
00416         {
00417       /* TODO free LU_btf_numeric and return NULL;
00418              * broadcast error to all processes */
00419             PARAKLETE_ERROR (PK_OUT_OF_MEMORY, "out of memory") ;
00420         }
00421     }
00422 
00423     /* all processes do this */
00424     LU_btf_numeric->Singleton = Singleton ;
00425     LU_btf_numeric->LUnumeric = LUnumeric =
00426   CHOLMOD (calloc) (nblocks, sizeof (void *), cm) ;
00427     LU_btf_numeric->F = F ;
00428     LU_btf_numeric->nblocks = nblocks ;
00429 
00430     if (cm->status != CHOLMOD_OK)
00431     {
00432         /* TODO free LU_btf_numeric and return NULL;
00433          * broadcast error to all processes */
00434         PARAKLETE_ERROR (PK_OUT_OF_MEMORY, "out of memory") ;
00435     }
00436 
00437     /* ---------------------------------------------------------------------- */
00438     /* factorize each diagonal block, and construct F */
00439     /* ---------------------------------------------------------------------- */
00440 
00441     for (block = 0 ; block < nblocks ; block++)
00442     {
00443   k1 = Rbtf [block] ;
00444   k2 = Rbtf [block+1] ;
00445   nk = k2 - k1 ;
00446 
00447   if (Common->myid == 0)
00448   {
00449       cnz = 0 ;
00450       for (k = k1 ; k < k2 ; k++)
00451       {
00452     Fp [k] = fnz ;
00453     Cp [k - k1] = cnz ;       /* note change of index */
00454     jold = Qbtf [k] ;
00455     for (p = Ap [jold] ; p < Ap [jold+1] ; p++)
00456     {
00457         inew = Pbinv [Ai [p]] ;
00458         if (inew < k1)
00459         {
00460       /* this entry goes in the off-diagonal part */
00461       Fi [fnz] = inew ;
00462       Fx [fnz] = Ax [p] ;
00463       fnz++ ;
00464         }
00465         else
00466         {
00467       /* this entry goes in the diagonal block */
00468       Ci [cnz] = inew - k1 ;  /* note the change of index */
00469       Cx [cnz] = Ax [p] ;
00470       cnz++ ;
00471         }
00472     }
00473       }
00474             Cp [nk] = cnz ;
00475             C->nrow = nk ;
00476             C->ncol = nk ;
00477   }
00478 
00479   if (nk == 1)
00480   {
00481       /* singleton */
00482       if (Common->myid == 0)
00483       {
00484     Singleton [block] = Cx [0] ;
00485       }
00486             LUnumeric [block] = NULL ;
00487   }
00488         else if (nk < 1000)
00489         {
00490             if (Common->myid == 0)
00491             {
00492                 /* use KLU on process 0 */
00493                 LUnumeric [block] = (void *) KLU_factor (Cp, Ci, Cx,
00494                     (KLU_symbolic *) LUsymbolic [block], km) ;
00495                 if (km->status != KLU_OK)
00496                 {
00497                     /* TODO return if KLU failed; broadcast the error */
00498                     PARAKLETE_ERROR (PK_UNKNOWN, "KLU factorize failed") ;
00499                     ok = FALSE ;
00500                 }
00501             }
00502             else
00503             {
00504                 /* small block; nothing to do for other processors */
00505                 LUnumeric [block] = NULL ;
00506             }
00507         }
00508   else
00509   {
00510             /* parallel factorization of a large block */
00511       LUnumeric [block] = (void *) amesos_paraklete_factorize (C,
00512                     (paraklete_symbolic *) LUsymbolic [block], Common) ;
00513             if (LUnumeric [block] == NULL)
00514             {
00515                 /* PARAKLETE failed */
00516                 ok = FALSE ;
00517             }
00518   }
00519     }
00520 
00521     if (Common->myid == 0)
00522     {
00523         /* process 0 frees its workspace */
00524   Fp [n] = fnz ;
00525   C->nrow = n ;
00526   C->ncol = n ;
00527   CHOLMOD (free_sparse) (&C, cm) ;
00528     }
00529 
00530     /* all processes find out if any one process fails */
00531     all_ok = ok ;
00532     MPI (MPI_Allreduce (&ok, &all_ok, 1, MPI_INT, MPI_LAND, MPI_COMM_WORLD)) ;
00533     if (!all_ok)
00534     {
00535   /* one of the blocks failed; all processes free result */
00536   amesos_paraklete_btf_free_numeric (&LU_btf_numeric, Common) ;
00537     }
00538     return (LU_btf_numeric) ;
00539 }
00540 
00541 
00542 /* ========================================================================== */
00543 /* === paraklete_btf_solve ================================================== */
00544 /* ========================================================================== */
00545 
00546 Int amesos_paraklete_btf_solve                 /* TRUE if OK, FALSE otherwise */
00547 (
00548     paraklete_btf_numeric *LU_btf_numeric,  /* numeric factorization */
00549     paraklete_btf_symbolic *LU_btf_symbolic,  /* symbolic analysis */
00550     double *B,        /* right-hand-side; soln on output */
00551     paraklete_common *Common
00552 )
00553 {
00554     double wk ;
00555     cholmod_common *cm ;
00556     KLU_common *km ;
00557     cholmod_sparse *F ;
00558     void **LUsymbolic ;
00559     void **LUnumeric ;
00560     double *Fx = NULL, *W = NULL, *Singleton ;
00561     Int *Qbtf, *Pbinv, *Rbtf, *Fp = NULL, *Fi = NULL ;
00562     Int iold, n, jnew, k, k1, k2, nk, p, nblocks, block ;
00563 
00564     /* ---------------------------------------------------------------------- */
00565     /* get inputs and allocate workspace */
00566     /* ---------------------------------------------------------------------- */
00567 
00568     cm = &(Common->cm) ;
00569     km = &(Common->km) ;
00570 
00571     nblocks = LU_btf_symbolic->nblocks ;
00572     n = LU_btf_symbolic->n ;
00573     Rbtf = LU_btf_symbolic->Rbtf ;
00574     Qbtf = LU_btf_symbolic->Qbtf ;
00575     Pbinv = LU_btf_symbolic->Pbinv ;
00576     LUsymbolic = LU_btf_symbolic->LUsymbolic ;
00577     LUnumeric = LU_btf_numeric->LUnumeric ;
00578     Singleton = LU_btf_numeric->Singleton ;
00579     F = LU_btf_numeric->F ;
00580     if (Common->myid == 0)
00581     {
00582         /* only the master process has the off-diagonal blocks F */
00583         Fp = F->p ;
00584         Fi = F->i ;
00585         Fx = F->x ;
00586     }
00587 
00588     /* ---------------------------------------------------------------------- */
00589     /* permute right-hand side; W = Pbtf*B */
00590     /* ---------------------------------------------------------------------- */
00591 
00592     if (Common->myid == 0)
00593     {
00594         /* only the master process does this */
00595         W = CHOLMOD (malloc) (n, sizeof (double), cm) ;
00596         if (cm->status != CHOLMOD_OK)
00597         {
00598             /* TODO return NULL, and broadcast error to all processes */
00599             PARAKLETE_ERROR (PK_OUT_OF_MEMORY, "out of memory") ;
00600         }
00601         for (iold = 0 ; iold < n ; iold++)
00602         {
00603             W [Pbinv [iold]] = B [iold] ;
00604         }
00605     }
00606 
00607     /* ---------------------------------------------------------------------- */
00608     /* block backsolve */
00609     /* ---------------------------------------------------------------------- */
00610 
00611     for (block = nblocks-1 ; block >= 0 ; block--)
00612     {
00613   k1 = Rbtf [block] ;
00614   k2 = Rbtf [block+1] ;
00615   nk = k2 - k1 ;
00616   if (nk == 1)
00617   {
00618       /* singleton */
00619             if (Common->myid == 0)
00620             {
00621                 W [k1] /= Singleton [block] ;
00622             }
00623   }
00624         else if (nk < 1000)
00625         {
00626             if (Common->myid == 0)
00627             {
00628                 /* use KLU on process 0 */
00629                 KLU_solve (
00630                     (KLU_symbolic *) LUsymbolic [block],
00631                     (KLU_numeric  *) LUnumeric  [block], nk, 1, W+k1, km) ;
00632 
00633                 if (km->status != KLU_OK)
00634                 {
00635                     /* TODO return NULL, and broadcast error to all processes */
00636                     PARAKLETE_ERROR (PK_UNKNOWN, "KLU solve failed") ;
00637                 }
00638 
00639             }
00640             else
00641             {
00642                 /* small block; nothing to do for other processors*/
00643                 LUnumeric [block] = NULL ;
00644             }
00645         }
00646   else
00647   {
00648       /* solve the diagonal block */
00649       amesos_paraklete_solve (
00650                 (paraklete_numeric *)  LUnumeric  [block], 
00651                 (paraklete_symbolic *) LUsymbolic [block], W+k1, Common) ;
00652             /* TODO check for error here */
00653   }
00654   /* backsolve for off-diagonal entries */
00655         if (Common->myid == 0)
00656         {
00657             for (k = k1 ; k < k2 ; k++)
00658             {
00659                 wk = W [k] ;
00660                 for (p = Fp [k] ; p < Fp [k+1] ; p++)
00661                 {
00662                     W [Fi [p]] -= Fx [p] * wk ;
00663                 }
00664             }
00665         }
00666     }
00667 
00668     /* ---------------------------------------------------------------------- */
00669     /* permute solution; B = Qbtf'*W */
00670     /* ---------------------------------------------------------------------- */
00671 
00672     if (Common->myid == 0)
00673     {
00674         for (jnew = 0 ; jnew < n ; jnew++)
00675         {
00676             B [Qbtf [jnew]] = W [jnew] ;
00677         }
00678         CHOLMOD (free) (n, sizeof (double), W, cm) ;
00679     }
00680 
00681     return (TRUE) ;
00682 }
00683 
00684 
00685 /* ========================================================================== */
00686 /* === paraklete_btf_free_symbolic ========================================== */
00687 /* ========================================================================== */
00688 
00689 void amesos_paraklete_btf_free_symbolic
00690 (
00691     paraklete_btf_symbolic **LU_btf_symbolic_handle,
00692     paraklete_common *Common
00693 )
00694 {
00695     paraklete_btf_symbolic *LU_btf_symbolic ;
00696     void **LUsymbolic ;
00697     Int n, block, nblocks, k1, k2, nk, *Rbtf ;
00698     cholmod_common *cm ;
00699     KLU_common *km ;
00700 
00701     if (LU_btf_symbolic_handle == NULL)
00702     {
00703   return ;
00704     }
00705     LU_btf_symbolic = *LU_btf_symbolic_handle ; 
00706     if (LU_btf_symbolic == NULL)
00707     {
00708   *LU_btf_symbolic_handle = NULL ;
00709   return ;
00710     }
00711     cm = &(Common->cm) ;
00712     km = &(Common->km) ;
00713 
00714     LUsymbolic = LU_btf_symbolic->LUsymbolic ;
00715     Rbtf = LU_btf_symbolic->Rbtf ;
00716     n = LU_btf_symbolic->n ;
00717     nblocks = LU_btf_symbolic->nblocks ;
00718     for (block = 0 ; block < nblocks ; block++)
00719     {
00720   k2 = Rbtf [block+1] ;
00721   k1 = Rbtf [block] ;
00722   nk = k2 - k1 ;
00723         if (nk < 1000)
00724         {
00725             KLU_symbolic *S ;
00726             S = (KLU_symbolic *) LUsymbolic [block] ;
00727             KLU_free_symbolic (&S, &(Common->km)) ;
00728         }
00729         else
00730         {
00731             paraklete_symbolic *S ;
00732             S = (paraklete_symbolic *) LUsymbolic [block] ;
00733             amesos_paraklete_free_symbolic (&S, Common) ;
00734         }
00735         LUsymbolic [block] = NULL ;
00736     }
00737     CHOLMOD (free) (n, sizeof (void *), LU_btf_symbolic->LUsymbolic, cm) ;
00738     CHOLMOD (free) (3*n+1, sizeof (Int), LU_btf_symbolic->Mem_n, cm) ;
00739     *LU_btf_symbolic_handle = CHOLMOD (free) (1, sizeof (paraklete_btf_symbolic),
00740   LU_btf_symbolic, cm) ;
00741 }
00742 
00743 
00744 /* ========================================================================== */
00745 /* === paraklete_btf_free_numeric =========================================== */
00746 /* ========================================================================== */
00747 
00748 void amesos_paraklete_btf_free_numeric
00749 (
00750     paraklete_btf_numeric **LU_btf_numeric_handle,
00751     paraklete_common *Common
00752 )
00753 {
00754     paraklete_btf_numeric *LU_btf_numeric ;
00755     void **LUnumeric ;
00756     Int block, nblocks ;
00757     cholmod_common *cm ;
00758 
00759     if (LU_btf_numeric_handle == NULL)
00760     {
00761   return ;
00762     }
00763     LU_btf_numeric = *LU_btf_numeric_handle ; 
00764     if (LU_btf_numeric == NULL)
00765     {
00766   *LU_btf_numeric_handle = NULL ;
00767   return ;
00768     }
00769     cm = &(Common->cm) ;
00770 
00771     LUnumeric = LU_btf_numeric->LUnumeric ;
00772     nblocks = LU_btf_numeric->nblocks ;
00773 
00774     for (block = 0 ; block < nblocks ; block++)
00775     {
00776         paraklete_numeric *N ;
00777         N = (paraklete_numeric *) LUnumeric [block] ;
00778         if (N != NULL)
00779         {
00780             if (N->magic == PARAKLETE_MAGIC)
00781             {
00782                 /* this block was factorized with Paraklete */
00783                 amesos_paraklete_free_numeric (&N, Common) ;
00784             }
00785             else
00786             {
00787                 /* this block was factorized with KLU */
00788                 KLU_numeric *K ;
00789                 K = (KLU_numeric *) LUnumeric [block] ;
00790                 KLU_free_numeric (&K, &(Common->km)) ;
00791             }
00792         }
00793         LUnumeric [block] = NULL ;
00794     }
00795 
00796     CHOLMOD (free) (nblocks, sizeof (double), LU_btf_numeric->Singleton, cm) ;
00797     CHOLMOD (free) (nblocks, sizeof (void *), LUnumeric, cm) ;
00798     CHOLMOD (free_sparse) (&(LU_btf_numeric->F), cm) ;
00799 
00800     *LU_btf_numeric_handle = CHOLMOD (free) (1, sizeof (paraklete_btf_numeric),
00801   LU_btf_numeric, cm) ;
00802 }

Generated on Sat Mar 14 12:22:51 2009 for Amesos Package Browser (Single Doxygen Collection) by  doxygen 1.3.9.1