Amesos Package Browser (Single Doxygen Collection) Development
amesos_paraklete_node.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === paraklete_factorize_node ============================================= */
00003 /* ========================================================================== */
00004 
00005 #include "amesos_paraklete_decl.h"
00006 
00007 /* Factorize node c of the separator tree:
00008  *  (1) obtain the Schur complements from the children of node c
00009  *  (2) sum up the Schur complements, and the input matrix at this node c
00010  *  (3) factorize node c
00011  *  (4) send the Schur complement of node c to its parent
00012  *
00013  * PARAKLETE version 0.3: parallel sparse LU factorization.  Nov 13, 2007
00014  * Copyright (C) 2007, Univ. of Florida.  Author: Timothy A. Davis
00015  * See License.txt for the Version 2.1 of the GNU Lesser General Public License
00016  * http://www.cise.ufl.edu/research/sparse
00017  */
00018 
00019 /* ========================================================================== */
00020 /* === paraklete_free_children ============================================== */
00021 /* ========================================================================== */
00022 
00023 /* Free the Schur complement of each child of node c.  Free the input matrix
00024  * A for this node c. */
00025 
00026 static void paraklete_free_children
00027 (
00028     Int c,
00029     paraklete_numeric *LU,
00030     paraklete_symbolic *LUsymbolic,
00031     paraklete_common *Common
00032 )
00033 {
00034     Int *Child, *Childp ;
00035     Int nchild, child, sn, cp ;
00036     cholmod_common *cm ;
00037 
00038     cm = &(Common->cm) ;
00039 
00040     /* free A for this node */
00041     CHOLMOD (free_sparse) (&(LU->LUnode [c]->A), cm) ;
00042 
00043     Child = LUsymbolic->Child ;
00044     Childp = LUsymbolic->Childp ;
00045     nchild = LU->LUnode [c]->nchild ;
00046     ASSERT (nchild == Childp [c+1] - Childp [c]) ;
00047 
00048     for (cp = 0 ; cp < nchild ; cp++)
00049     {
00050   /* free the Schur complement of the child */
00051   child = Child [Childp [c] + cp] ;
00052   sn = LU->LUnode [child]->PK_SN ;
00053 
00054   LU->LUnode [child]->sx = CHOLMOD (free) (
00055     LU->LUnode [child]->PK_SSIZE, sizeof (double),
00056     LU->LUnode [child]->sx, cm) ;
00057   LU->LUnode [child]->sp = CHOLMOD (free) (
00058     sn, sizeof (Int), LU->LUnode [child]->sp, cm) ;
00059   LU->LUnode [child]->slen = CHOLMOD (free) (
00060     sn, sizeof (Int), LU->LUnode [child]->slen, cm) ;
00061     }
00062 }
00063 
00064 
00065 /* ========================================================================== */
00066 /* === paraklete_send_to_parent ============================================= */
00067 /* ========================================================================== */
00068 
00069 /* Send the Schur complement or an error message to the parent */
00070 
00071 static void paraklete_send_to_parent
00072 (
00073     Int c,
00074     Int status,
00075     Int parent_id,
00076     paraklete_numeric *LU,
00077     paraklete_symbolic *LUsymbolic,
00078     paraklete_common *Common
00079 )
00080 {
00081     int ok ;
00082     cholmod_common *cm ;
00083     paraklete_node *LUnode ;
00084     MPI (MPI_Status ms) ;
00085 
00086     cm = &(Common->cm) ;
00087     LUnode = LU->LUnode [c] ;
00088 
00089     /* workspace W2 and C no longer needed */
00090     LUnode->W2 = CHOLMOD (free) (2*(LUnode->nlost_in), sizeof (Int),
00091       LUnode->W2, cm) ;
00092     CHOLMOD (free_sparse) (&(LUnode->C), cm) ;
00093 
00094     LUnode->PK_STATUS = status ;
00095 
00096     ok = TRUE ;
00097     if (parent_id != EMPTY && parent_id != Common->myid)
00098     {
00099   /* send header to process that owns the parent node (parent_id) */
00100   PR1 ((Common->file, "proc "ID" sends header parent proc "ID"\n",
00101     Common->myid, parent_id)) ;
00102   MPI (MPI_Send (LUnode->header, PK_HEADER, MPI_Int,
00103         parent_id, TAG0, MPI_COMM_WORLD)) ;
00104 
00105   /* wait to see if parent_id is OK */
00106   PR1 ((Common->file, "proc "ID" wait for parent "ID"\n",
00107     Common->myid, parent_id)) ;
00108   MPI (MPI_Recv (&ok, 1, MPI_INT, parent_id, TAG0, MPI_COMM_WORLD, &ms)) ;
00109 
00110   /* if status is not PK_OK, then parent_id will send back ok = FALSE */
00111   ASSERT (IMPLIES (status != PK_OK, !ok)) ;
00112 
00113   if (ok)
00114   {
00115       /* both parent_id and myid agree to send the Schur complement */
00116       /* TODO: send this as one or two messages */
00117       MPI (MPI_Rsend (LUnode->sp, LUnode->PK_SN, MPI_Int,
00118       parent_id, TAG0, MPI_COMM_WORLD)) ;
00119       MPI (MPI_Rsend (LUnode->slen, LUnode->PK_SN, MPI_Int,
00120       parent_id, TAG0, MPI_COMM_WORLD)) ;
00121       MPI (MPI_Rsend (LUnode->sx, LUnode->PK_SSIZE, MPI_DOUBLE,
00122       parent_id, TAG0, MPI_COMM_WORLD)) ;
00123       MPI (MPI_Rsend (LUnode->Pglobal, LUnode->PK_NPIV, MPI_Int,
00124       parent_id, TAG0, MPI_COMM_WORLD)) ;
00125       MPI (MPI_Rsend (LUnode->Qglobal, LUnode->PK_NPIV, MPI_Int,
00126       parent_id, TAG0, MPI_COMM_WORLD)) ;
00127   }
00128     }
00129 
00130     if (!ok || parent_id == EMPTY || parent_id != Common->myid)
00131     {
00132   /* free the Schur complement of node c if a failure occured, if node
00133    * c has no parent, or if the Schur complement was sent to a
00134    * different process (parent_id). */
00135   LUnode->sx = CHOLMOD (free) (LUnode->PK_SSIZE,
00136     sizeof (double), LUnode->sx, cm) ;
00137   LUnode->sp = CHOLMOD (free) (LUnode->PK_SN,
00138     sizeof (Int), LUnode->sp, cm) ;
00139   LUnode->slen = CHOLMOD (free) (LUnode->PK_SN,
00140     sizeof (Int), LUnode->slen, cm) ;
00141     }
00142 
00143     if (!ok)
00144     {
00145   /* free the Schur complements of the children if a failure occured */
00146   paraklete_free_children (c, LU, LUsymbolic, Common) ;
00147     }
00148 }
00149 
00150 
00151 /* ========================================================================== */
00152 /* === paraklete_factorize_node ============================================= */
00153 /* ========================================================================== */
00154 
00155 Int amesos_paraklete_factorize_node
00156 (
00157     Int c,
00158     paraklete_numeric *LU,
00159     paraklete_symbolic *LUsymbolic,
00160     paraklete_common *Common
00161 )
00162 {
00163     paraklete_node *LUnode ;
00164     cholmod_common *cm ;
00165     cholmod_sparse *C ;
00166     double *Cx, *W, *Sx, *Ax, *S ;
00167     Int *Child, *Childp, *Lost, *Lostp, *Plost_in, *Qlost_in, *Cp, *Ci, *Flag,
00168   *Sp, *Si, *Slen, *Ap, *Ai, *Cn, *Plocal, *Qlocal, *Pglobal, *Qglobal,
00169   *Sched, *Pinv ;
00170     Int cp, cnz, clnz, nchild, k1, k2, child, cn, npiv, k, j, p, i, nfound,
00171   mark, cj, ci, nlost_in, len, sn, myid, ssize, parent, nmessages,
00172   status, parent_id ;
00173     int ok ;
00174     size_t s ;
00175 
00176     MPI (MPI_Status ms) ;
00177     MPI (MPI_Request req [5]) ;
00178 
00179     /* ---------------------------------------------------------------------- */
00180     /* get local workspace */
00181     /* ---------------------------------------------------------------------- */
00182 
00183     cm = &(Common->cm) ;
00184     PR0 ((Common->file, "\n\n######################## FACTORIZE NODE "ID"\n", c));
00185 
00186     /* ---------------------------------------------------------------------- */
00187     /* get the symbolic analysis of this node */
00188     /* ---------------------------------------------------------------------- */
00189 
00190     cnz    = LUsymbolic->Cnz [c] ;  /* # entries in A for this node c */
00191     Childp = LUsymbolic->Childp ; /* Child [Childp [c]..Childp[c+1]-1] */
00192     Child  = LUsymbolic->Child ;  /* is list of children of node c */
00193     clnz   = LUsymbolic->Clnz [c] ; /* est. # entries in L for node c */
00194     nchild = Childp [c+1] - Childp [c] ;    /* # of children of node c */
00195     k1     = LUsymbolic->Cstart [c] ; /* global index of first pivot cand. */
00196     k2     = LUsymbolic->Cstart [c+1] ; /* global index of last piv plus one */
00197     Cn     = LUsymbolic->Cn ;   /* dimension of each node */
00198     Sched  = LUsymbolic->Sched ;
00199     parent = LUsymbolic->Cparent [c] ;
00200     myid = Common->myid ;
00201     parent_id = (parent == EMPTY) ? EMPTY : Sched [parent] ;
00202 
00203     PR0 ((Common->file, "proc "ID" at node "ID", clnz: "ID"\n", myid, c, clnz)) ;
00204 
00205     /* ---------------------------------------------------------------------- */
00206     /* get the arrowhead of the input matrix to factorize at this node */
00207     /* ---------------------------------------------------------------------- */
00208 
00209     LUnode = LU->LUnode [c] ;
00210     ASSERT (nchild == LUnode->nchild) ;
00211     Ap = LUnode->A->p ;     /* A matrix of dimension Cn [c] */
00212     Ai = LUnode->A->i ;
00213     Ax = LUnode->A->x ;
00214     DEBUG (CHOLMOD (print_sparse) (LUnode->A, "Arrowhead", cm)) ;
00215 
00216 #ifndef NDEBUG
00217     PR1 ((Common->file, "proc "ID" node "ID" nchild "ID"\n", myid, c, nchild)) ;
00218     {
00219   Int cc ;
00220   for (cc = 0 ; cc < LUsymbolic->ncomponents ; cc++)
00221   {
00222       PR2 ((Common->file,
00223     "proc "ID": node "ID" sched "ID" cparent "ID" childp ["ID":"ID"]\n",
00224     myid, cc, Sched [cc], LUsymbolic->Cparent [cc], Childp [cc],
00225     Childp [cc+1]-1)) ;
00226       for (cp = 0 ; cp < Childp [cc+1] - Childp [cc] ; cp++)
00227       {
00228     PR2 ((Common->file, "node "ID": child node "ID"\n",
00229           cc, Child [Childp[c] + cp])) ;
00230       }
00231   }
00232     }
00233 #endif
00234 
00235     /* ---------------------------------------------------------------------- */
00236     /* post a non-blocking receive for the header information from each child */
00237     /* ---------------------------------------------------------------------- */
00238 
00239     PR1 ((Common->file, "proc "ID" posting recv at node "ID", nchild "ID" status "ID"\n",
00240       myid, c, nchild, cm->status)) ;
00241     nmessages = 0 ;
00242     for (cp = 0 ; cp < nchild ; cp++)
00243     {
00244   child = Child [Childp [c] + cp] ;
00245   PR2 ((Common->file, "proc "ID" child "ID" owned by "ID"\n",
00246         myid, child, Sched [child])) ;
00247   if (Sched [child] != myid)
00248   {
00249       PR2 ((Common->file, "parent proc "ID" awaits header from "ID"\n",
00250         myid, Sched [child])) ;
00251       MPI (MPI_Irecv (LU->LUnode [child]->header, PK_HEADER, MPI_Int,
00252         Sched [child], TAG0, MPI_COMM_WORLD, &(LUnode->Req [cp]))) ;
00253       nmessages++ ;
00254   }
00255   else
00256   {
00257       MPI (LUnode->Req [cp] = MPI_REQUEST_NULL) ;
00258   }
00259     }
00260 
00261 #ifdef NMPI
00262     /* all nodes in the tree are scheduled to process zero, if no MPI */
00263     ASSERT (nmessages == 0) ;
00264 #endif
00265 
00266     /* ---------------------------------------------------------------------- */
00267     /* get header and Schur complement from each child */
00268     /* ---------------------------------------------------------------------- */
00269 
00270     ok = TRUE ;
00271     for (k = 0 ; k < nmessages ; k++)
00272     {
00273 
00274   /* ------------------------------------------------------------------ */
00275   /* get header from a child who is ready to send its Schur complement*/
00276   /* ------------------------------------------------------------------ */
00277 
00278   cp = 0 ;
00279   MPI (MPI_Waitany (nchild, LUnode->Req, &cp, &ms)) ;
00280 
00281   child = Child [Childp [c] + cp] ;
00282   ASSERT (Sched [child] != myid) ;
00283 
00284   status = LU->LUnode [child]->PK_STATUS ;
00285   sn = LU->LUnode [child]->PK_SN ;
00286   nfound = LU->LUnode [child]->PK_NFOUND ;
00287   npiv = LU->LUnode [child]->PK_NPIV ;
00288   ssize = LU->LUnode [child]->PK_SSIZE ;
00289   cn = LU->LUnode [child]->PK_NN ;
00290 
00291   ok = ok && (status == PK_OK) ;
00292 
00293   /* ------------------------------------------------------------------ */
00294   /* allocate space for Schur complement of the child */
00295   /* ------------------------------------------------------------------ */
00296 
00297   if (ok)
00298   {
00299       /* allocate space for next message: the Schur complement */
00300       LU->LUnode [child]->sp = CHOLMOD (malloc) (sn,
00301         sizeof (Int), cm) ;
00302       LU->LUnode [child]->slen = CHOLMOD (malloc) (sn,
00303         sizeof (Int), cm) ;
00304       LU->LUnode [child]->Pglobal = CHOLMOD (malloc) (npiv,
00305         sizeof (Int), cm) ;
00306       LU->LUnode [child]->Qglobal = CHOLMOD (malloc) (npiv,
00307         sizeof (Int), cm) ;
00308       LU->LUnode [child]->sx = CHOLMOD (malloc) (ssize,
00309         sizeof (double), cm) ;
00310 
00311       /* allocate space for foward/backsolve */
00312       LU->LUnode [child]->X = CHOLMOD (malloc) (cn,
00313         sizeof (double), cm) ;
00314   }
00315 
00316         if (cm->status != CHOLMOD_OK)
00317         {
00318             /* TODO return NULL, and broadcast error to all processes */
00319             PARAKLETE_ERROR (PK_OUT_OF_MEMORY, "out of memory") ;
00320         }
00321 
00322   /* check if we ran out of memory */
00323   ok = ok && (cm->status == CHOLMOD_OK) ;
00324 
00325   /* ------------------------------------------------------------------ */
00326   /* post receives for the Schur complement from the child */
00327   /* ------------------------------------------------------------------ */
00328 
00329   if (ok)
00330   {
00331       /* both parent and child agree to send the Schur complement */
00332       /* TODO: fold this information into fewer messages */
00333       MPI (MPI_Irecv (LU->LUnode [child]->sp, sn, MPI_Int,
00334       Sched [child], TAG0, MPI_COMM_WORLD, &(req [0]))) ;
00335       MPI (MPI_Irecv (LU->LUnode [child]->slen, sn, MPI_Int,
00336       Sched [child], TAG0, MPI_COMM_WORLD, &(req [1]))) ;
00337       MPI (MPI_Irecv (LU->LUnode [child]->sx, ssize, MPI_DOUBLE,
00338       Sched [child], TAG0, MPI_COMM_WORLD, &(req [2]))) ;
00339       MPI (MPI_Irecv (LU->LUnode [child]->Pglobal, npiv, MPI_Int,
00340       Sched [child], TAG0, MPI_COMM_WORLD, &(req [3]))) ;
00341       MPI (MPI_Irecv (LU->LUnode [child]->Qglobal, npiv, MPI_Int,
00342       Sched [child], TAG0, MPI_COMM_WORLD, &(req [4]))) ;
00343   }
00344 
00345   /* ------------------------------------------------------------------ */
00346   /* tell child that we are ready to receive its Schur complement */
00347   /* ------------------------------------------------------------------ */
00348 
00349   PR1 ((Common->file, "parent proc "ID" replies to proc "ID"\n",
00350         myid, Sched [child])) ;
00351   MPI (MPI_Send (&ok, 1, MPI_INT, Sched [child], TAG0, MPI_COMM_WORLD)) ;
00352 
00353   /* ------------------------------------------------------------------ */
00354   /* wait for the Schur complement to be received */
00355   /* ------------------------------------------------------------------ */
00356 
00357   if (ok)
00358   {
00359       MPI (MPI_Waitall (5, req, MPI_STATUSES_IGNORE)) ;
00360   }
00361     }
00362 
00363     PR1 ((Common->file, "proc "ID" node "ID" arrowhead again\n", myid, c)) ;
00364     DEBUG (CHOLMOD (print_sparse) (LUnode->A, "Arrowhead again", cm)) ;
00365 
00366     /* ---------------------------------------------------------------------- */
00367     /* report failure to parent of c, if a failure occured */
00368     /* ---------------------------------------------------------------------- */
00369 
00370     if (!ok)
00371     {
00372   PR0 ((Common->file, "proc "ID", node "ID", report failure to parent: "ID"\n",
00373     myid, c, parent_id)) ;
00374   /* 
00375         printf ("proc "ID" out of memory at node "ID"\n", myid, c) ;
00376         */
00377   paraklete_send_to_parent (c, PK_OUT_OF_MEMORY, parent_id,
00378     LU, LUsymbolic, Common) ;
00379   return (FALSE) ;
00380     }
00381 
00382     DEBUG (CHOLMOD (print_sparse) (LUnode->A, "Arrowhead yet again", cm)) ;
00383 
00384     /* ---------------------------------------------------------------------- */
00385     /* get lost pivots and Schur complement from each child */
00386     /* ---------------------------------------------------------------------- */
00387 
00388     nlost_in = 0 ;
00389 
00390     Lost  = LUnode->Lost ;
00391     Lostp = LUnode->Lostp ;
00392     LUnode->nchild = nchild ;
00393 
00394     s = cnz ;
00395     for (cp = 0 ; cp < nchild ; cp++)
00396     {
00397   /* find the failed pivot rows/cols of this child and add them
00398    * to the pivot candidate sets of this node s */
00399   child = Child [Childp [c] + cp] ;
00400   Lostp [cp] = nlost_in ;
00401   Lost [cp] = LU->LUnode [child]->PK_NLOST ;
00402   PR1 ((Common->file, "child "ID" lost "ID" \n", child, Lost [cp])) ;
00403   nlost_in += Lost [cp] ;
00404   s = CHOLMOD (add_size_t) (s, LU->LUnode [child]->PK_SNZ, &ok) ;
00405         if (!ok)
00406         {
00407             /* TODO broadcast the error to all processes */
00408             PARAKLETE_ERROR (PK_TOO_LARGE, "problem too large") ;
00409         }
00410     }
00411     cnz = s ;
00412     Lostp [nchild] = nlost_in ;
00413     LUnode->nlost_in = nlost_in ;
00414 
00415     /* The matrix to factorize at this node is cn-by-cn.  Up to npiv pivots
00416      * can be found in this node */
00417     cn = nlost_in + Cn [c] ;
00418     npiv = nlost_in + (k2 - k1) ;
00419     LUnode->PK_NN = cn ;
00420 
00421     LUnode->W2 = CHOLMOD (malloc) (nlost_in, 2*sizeof (Int), cm) ;
00422     Plost_in = LUnode->W2 ;   /* size nlost_in */
00423     Qlost_in = LUnode->W2 + nlost_in ;  /* size nlost_in */
00424 
00425     /* get workspace */
00426     CHOLMOD (allocate_work) (cn, 3*cn, cn, cm) ;
00427 
00428     DEBUG (CHOLMOD (print_sparse) (LUnode->A, "Arrowhead once again", cm)) ;
00429 
00430     /* ---------------------------------------------------------------------- */
00431     /* C = original entries in this node, plus Schur complements of children */
00432     /* ---------------------------------------------------------------------- */
00433 
00434     /* TODO: assemble the Schur complements but do not compute C, just to
00435      * determine cnz.  Then do the CHOLMOD (allocate_sparse), below: */
00436 
00437     C = CHOLMOD (allocate_sparse) (cn, cn, cnz, FALSE, TRUE, 0, TRUE, cm) ;
00438     LUnode->C = C ;
00439 
00440     if (cm->status != CHOLMOD_OK)
00441     {
00442   /* out of memory; tell the parent that we failed */
00443         PARAKLETE_ERROR (PK_OUT_OF_MEMORY, "out of memory") ;
00444   paraklete_send_to_parent (c, PK_OUT_OF_MEMORY, parent_id,
00445     LU, LUsymbolic, Common) ;
00446   return (FALSE) ;
00447     }
00448 
00449     Cp = C->p ;
00450     Ci = C->i ;
00451     Cx = C->x ;
00452 
00453     /* ---------------------------------------------------------------------- */
00454     /* concatenate and remap the lost pivot columns of each child */
00455     /* ---------------------------------------------------------------------- */
00456 
00457     k = 0 ;
00458     cnz = 0 ;
00459     for (cp = 0 ; cp < nchild ; cp++)
00460     {
00461   /* get the Schur complement of the child */
00462   child = Child [Childp [c] + cp] ;
00463   sn = LU->LUnode [child]->PK_SN ;
00464 
00465   PR1 ((Common->file, "\nconcatenate lost child "ID", Lost [cp] "ID"\n",
00466         child, Lost [cp])) ;
00467   S    = LU->LUnode [child]->sx ;
00468   Sp   = LU->LUnode [child]->sp ;
00469   Slen = LU->LUnode [child]->slen ;
00470   nfound = LU->LUnode [child]->PK_NFOUND ;
00471 
00472   PR1 ((Common->file, "Cn [c] is "ID"\n", Cn [c])) ;
00473   PR1 ((Common->file, "child is "ID"\n", child)) ;
00474   PR1 ((Common->file, "Lost[cp] is "ID"\n", Lost [cp])) ;
00475   PR1 ((Common->file, "sn "ID"  Lost[cp]+Cn[cn] "ID"\n",
00476       sn, Lost [cp] + Cn [c]));
00477 
00478   ASSERT (sn == Lost [cp] + Cn [c]) ;
00479 
00480   for (j = 0 ; j < Lost [cp] ; j++)
00481   {
00482       /* column j of the Schur complement of the child becomes column
00483        * k of C */
00484       PR2 ((Common->file, "lost child col "ID" becomes col "ID" of C\n",
00485     j, k)) ;
00486       Cp [k] = cnz ;
00487       GET_COLUMN (Sp, Slen, S, j, Si, Sx, len) ;
00488       for (p = 0 ; p < len ; p++)
00489       {
00490     i = Si [p] ;
00491     ci = i + ((i < Lost [cp]) ? Lostp [cp] : (nlost_in - Lost[cp]));
00492     Ci [cnz] = ci ;
00493     Cx [cnz] = Sx [p] ;
00494     PR3 ((Common->file,
00495       "  Lost child col: row "ID" newrow "ID" value %g\n",
00496       i, ci, Sx [p])) ;
00497     cnz++ ;
00498       }
00499       /* get the lost pivot row/column from the child */
00500       Plost_in [k] = LU->LUnode [child]->Pglobal [nfound+j] ;
00501       Qlost_in [k] = LU->LUnode [child]->Qglobal [nfound+j] ;
00502       k++ ;
00503   }
00504     }
00505     ASSERT (k == nlost_in) ;
00506 
00507     /* ---------------------------------------------------------------------- */
00508     /* assemble original entries and Schur complements of each child */
00509     /* ---------------------------------------------------------------------- */
00510 
00511     Flag = cm->Flag ;
00512     W = cm->Xwork ;
00513     /* mark = CHOLMOD (clear_flag) (cm) ; */
00514     CHOLMOD_CLEAR_FLAG (cm) ;
00515     mark = cm->mark ;
00516 
00517     ASSERT (cm->xworksize >= (size_t) cn) ;
00518     DEBUG (for (cj = 0 ; cj < cn ; cj++) ASSERT (W [cj] == 0)) ;
00519 
00520     for (j = 0 ; j < Cn [c] ; j++)
00521     {
00522 
00523   /* Column j is nominal column index, if there were no failed pivots.
00524    * The column is shifted over to accomodate incoming lost pivot columns
00525    * from the children, to obtain the new column cj. */
00526 
00527   cj = j + nlost_in ;
00528   Cp [cj] = cnz ;
00529   PR2 ((Common->file, "\nAdd column "ID" of Schur\n", cj)) ;
00530 
00531   /* ------------------------------------------------------------------ */
00532   /* add up all the contributions to column cj of C */
00533   /* ------------------------------------------------------------------ */
00534 
00535   /* scatter the original entries of column j.  A is the input arrowhead
00536    * for this node, with row/col indices already mapped to the local
00537    * row/col index space (in range 0 to Cn [c]). */
00538   for (p = Ap [j] ; p < Ap [j+1] ; p++)
00539   {
00540       i = Ai [p] ;
00541       ci = i + nlost_in ;
00542       PR3 ((Common->file, "scatter original ("ID","ID") %g to ("ID","ID")\n",
00543       i, j, Ax [p], ci, cj)) ;
00544       Flag [ci] = mark ;
00545       Ci [cnz++] = ci ;
00546       W [ci] = Ax [p] ;
00547   }
00548 
00549   /* scatter and add the contributions from each child */
00550   for (cp = 0 ; cp < nchild ; cp++)
00551   {
00552       /* get column j+Lost[cp] of the Schur complement of the child */
00553       child = Child [Childp [c] + cp] ;
00554       S    = LU->LUnode [child]->sx ;
00555       Sp   = LU->LUnode [child]->sp ;
00556       Slen = LU->LUnode [child]->slen ;
00557       GET_COLUMN (Sp, Slen, S, j + Lost [cp], Si, Sx, len) ;
00558       for (p = 0 ; p < len ; p++)
00559       {
00560     i = Si [p] ;
00561     ci = i + ((i < Lost [cp]) ? Lostp [cp] : (nlost_in - Lost[cp]));
00562     if (Flag [ci] < mark)
00563     {
00564         Flag [ci] = mark ;
00565         Ci [cnz++] = ci ;
00566     }
00567     W [ci] += Sx [p] ;
00568       }
00569   }
00570 
00571   /* gather the results */
00572   PR2 ((Common->file, "gather "ID" to "ID"-1\n", Cp [cj], cnz)) ;
00573   for (p = Cp [cj] ; p < cnz ; p++)
00574   {
00575       ci = Ci [p] ;
00576       Cx [p] = W [ci] ;
00577       W [ci] = 0 ;
00578       PR3 ((Common->file, ""ID": %g\n", ci, Cx [p])) ;
00579   }
00580   DEBUG (for (cj = 0 ; cj < cn ; cj++) ASSERT (W [cj] == 0)) ;
00581 
00582   /* clear Flag array */
00583   /* mark = CHOLMOD (clear_flag) (cm) ; */
00584   CHOLMOD_CLEAR_FLAG (cm)  ;
00585   mark = cm->mark ;
00586     }
00587 
00588     Cp [cn] = cnz ;
00589 
00590     /* shrink C to be just large enough */
00591     CHOLMOD (reallocate_sparse) (cnz, C, cm) ;
00592     ASSERT (cm->status == CHOLMOD_OK) ;
00593 
00594     /* ---------------------------------------------------------------------- */
00595     /* A for this node, and Schur complements of children, no longer needed */
00596     /* ---------------------------------------------------------------------- */
00597 
00598     paraklete_free_children (c, LU, LUsymbolic, Common) ;
00599 
00600     /* ---------------------------------------------------------------------- */
00601     /* allocate the LU factors for this node */
00602     /* ---------------------------------------------------------------------- */
00603 
00604     /* LUnode->lusize = 1.2 * (2 * ((3*clnz)/2 + 2)) ; */
00605     /* clnz is already checked for integer overflow */
00606     s = clnz ;
00607     s = CHOLMOD (mult_size_t) (s, 3, &ok) / 2 ;       /* (3*clnz)/2 */
00608     s = CHOLMOD (add_size_t)  (s, 2, &ok) ;           /* add 2 */
00609     s = CHOLMOD (mult_size_t) (s, 2, &ok) ;           /* times 2 */
00610     s = CHOLMOD (add_size_t)  (s, s/5, &ok) ;         /* add s/5 */
00611     LUnode->lusize = s ;
00612 
00613     LUnode->PK_NPIV = npiv ;
00614 
00615     LUnode->llen = CHOLMOD (malloc) (npiv, sizeof (Int), cm) ;
00616     LUnode->lp   = CHOLMOD (malloc) (npiv, sizeof (Int), cm) ;
00617 
00618     LUnode->ulen = CHOLMOD (malloc) (cn, sizeof (Int), cm) ;
00619     LUnode->up   = CHOLMOD (malloc) (cn, sizeof (Int), cm) ;
00620 
00621     LUnode->Plocal  = CHOLMOD (malloc) (npiv, sizeof (Int), cm) ;
00622     LUnode->Pglobal = CHOLMOD (malloc) (npiv, sizeof (Int), cm) ;
00623     LUnode->Qlocal  = CHOLMOD (malloc) (npiv, sizeof (Int), cm) ;
00624     LUnode->Qglobal = CHOLMOD (malloc) (npiv, sizeof (Int), cm) ;
00625 
00626     LUnode->Pinv = CHOLMOD (malloc) (npiv, sizeof (Int), cm) ;
00627     LUnode->Qinv = CHOLMOD (malloc) (npiv, sizeof (Int), cm) ;
00628 
00629     /* this block of memory may grow in size */
00630     LUnode->ix = CHOLMOD (malloc) (LUnode->lusize, sizeof (double), cm) ;
00631     PR1 ((Common->file, "allocated LUnode->ix: size "ID"\n", LUnode->lusize)) ;
00632 
00633     if (cm->status != CHOLMOD_OK)
00634     {
00635   /* out of memory */
00636         PARAKLETE_ERROR (PK_OUT_OF_MEMORY, "out of memory") ;
00637   paraklete_send_to_parent (c, PK_OUT_OF_MEMORY, parent_id,
00638     LU, LUsymbolic, Common) ;
00639   return (FALSE) ;
00640     }
00641 
00642     /* ---------------------------------------------------------------------- */
00643     /* factorize P*C*Q into L*U+S */
00644     /* ---------------------------------------------------------------------- */
00645 
00646     ASSERT (CHOLMOD (print_sparse) (C, "C = A + sum (Schur)", cm)) ;
00647     ok = amesos_paraklete_kernel (C, LUnode, Common) ;
00648 
00649     if (!ok)
00650     {
00651   /* out of memory */
00652         PARAKLETE_ERROR (PK_OUT_OF_MEMORY, "out of memory") ;
00653   paraklete_send_to_parent (c, PK_OUT_OF_MEMORY, parent_id,
00654     LU, LUsymbolic, Common) ;
00655   return (FALSE) ;
00656     }
00657 
00658     PR0 ((Common->file, "============= NODE "ID" NPIV "ID" FOUND "ID" LOST "ID"\n", c,
00659       LUnode->PK_NPIV, LUnode->PK_NFOUND,
00660       LUnode->PK_NPIV - LUnode->PK_NFOUND)) ;
00661 
00662     /* TODO ensure the kernel does this */
00663     cm->mark = EMPTY ;
00664     /* CHOLMOD (clear_flag) (cm) ; */
00665     CHOLMOD_CLEAR_FLAG (cm) 
00666 
00667     Plocal = LUnode->Plocal ;
00668     Pinv = LUnode->Pinv ;
00669     for (k = 0 ; k < npiv ; k++)
00670     {
00671   i = Plocal [k] ;
00672   ASSERT (i >= 0 && i < npiv) ;
00673   ASSERT (Pinv [i] == k) ;
00674     }
00675 
00676     /* ---------------------------------------------------------------------- */
00677     /* determine the global pivot ordering (not including LUsymbolic->Cperm) */
00678     /* ---------------------------------------------------------------------- */
00679 
00680     Plocal = LUnode->Plocal ;
00681     Qlocal = LUnode->Qlocal ;
00682 
00683     Pglobal = LUnode->Pglobal ;
00684     Qglobal = LUnode->Qglobal ;
00685 
00686     DEBUG (CHOLMOD (print_perm) (Plocal, npiv, npiv, "Node P, local", cm)) ;
00687     DEBUG (CHOLMOD (print_perm) (Qlocal, npiv, npiv, "Node Q, local", cm)) ;
00688 
00689     nfound = LUnode->PK_NFOUND ;
00690 
00691     for (k = 0 ; k < npiv ; k++)
00692     {
00693   i = Plocal [k] ;
00694   if (i < nlost_in)
00695   {
00696       /* row i was the ith incoming lost row, from the children */
00697       Pglobal [k] = Plost_in [i] ;
00698   }
00699   else
00700   {
00701       /* row i was an original candidate for this node.  It was shifted
00702        * by nlost_in to obtain the local index.  Then k1 needs to be added
00703        * to obtain the global index. */
00704       Pglobal [k] = i - nlost_in + k1 ;
00705   }
00706     }
00707     for (k = 0 ; k < npiv ; k++)
00708     {
00709   j = Qlocal [k] ;
00710   if (j < nlost_in)
00711   {
00712       /* col j was the jth incoming lost col, from the children */
00713       Qglobal [k] = Qlost_in [j] ;
00714   }
00715   else
00716   {
00717       /* col j was an original candidate for this node.  It was shifted
00718        * by nlost_in to obtain the local index.  Then k1 needs to be added
00719        * to obtain the global index. */
00720       Qglobal [k] = j - nlost_in + k1 ;
00721   }
00722     }
00723 
00724     LUnode->PK_NLOST = npiv - nfound ;
00725 
00726     DEBUG (CHOLMOD (print_perm) (Pglobal, npiv, LU->n, "Node P, global", cm)) ;
00727     DEBUG (CHOLMOD (print_perm) (Qglobal, npiv, LU->n, "Node Q, global", cm)) ;
00728 
00729     /* ---------------------------------------------------------------------- */
00730     /* allocate space for subsequent forward/backsolve */
00731     /* ---------------------------------------------------------------------- */
00732 
00733     LU->LUnode [c]->X = CHOLMOD (malloc) (cn, sizeof (double), cm) ;
00734     ASSERT (k2-k1 == LUnode->nk) ;
00735     LU->LUnode [c]->B = CHOLMOD (malloc) (k2-k1, sizeof (double), cm) ;
00736 
00737     if (cm->status != CHOLMOD_OK)
00738     {
00739   /* out of memory */
00740         PARAKLETE_ERROR (PK_OUT_OF_MEMORY, "out of memory") ;
00741   paraklete_send_to_parent (c, PK_OUT_OF_MEMORY, parent_id,
00742     LU, LUsymbolic, Common) ;
00743   return (FALSE) ;
00744     }
00745 
00746     /* ---------------------------------------------------------------------- */
00747     /* send Schur complement to the parent */
00748     /* ---------------------------------------------------------------------- */
00749 
00750     PR0 ((Common->file, "proc "ID" done at node "ID"\n", myid, c)) ;
00751 
00752     paraklete_send_to_parent (c, PK_OK, parent_id, LU, LUsymbolic, Common) ;
00753 
00754     PR0 ((Common->file, "proc "ID" done send to parent at node "ID"\n", myid, c)) ;
00755 
00756     return (TRUE) ;
00757 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines