Amesos Package Browser (Single Doxygen Collection) Development
amesos_paraklete_lsolve_node.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === paraklete_lsolve_node ================================================ */
00003 /* ========================================================================== */
00004 
00005 #include "amesos_paraklete_decl.h"
00006 
00007 /* Solve Lx=b with node c of the separator tree.
00008  *
00009  * PARAKLETE version 0.3: parallel sparse LU factorization.  Nov 13, 2007
00010  * Copyright (C) 2007, Univ. of Florida.  Author: Timothy A. Davis
00011  * See License.txt for the Version 2.1 of the GNU Lesser General Public License
00012  * http://www.cise.ufl.edu/research/sparse
00013  */
00014 
00015 Int amesos_paraklete_lsolve_node
00016 (
00017     Int c,
00018     paraklete_numeric *LU,
00019     paraklete_symbolic *LUsymbolic,
00020     paraklete_common *Common
00021 )
00022 
00023 {
00024     double xj ;
00025     paraklete_node *LUnode ;
00026     cholmod_common *cm ;
00027     double *X, *B, *LUix, *Xchild, *Lx ;
00028     Int *Child, *Childp, *Lost, *Lostp, *Cstart, *Lip, *Llen,
00029   *Cn, *Pinv, *Li, *Cparent, *Sched ;
00030     Int cp, nchild, child, cn, k, j, p, i, nfound, k1, k2, llen,
00031   ci, nlost_in, cn_child, cn_nfound, parent, myid, pass ;
00032     MPI (MPI_Status ms) ;
00033     MPI (MPI_Request req) ;
00034 
00035     /* ---------------------------------------------------------------------- */
00036     /* get local workspace */
00037     /* ---------------------------------------------------------------------- */
00038 
00039     cm = &(Common->cm) ;
00040     PR0 ((Common->file, "\n\n########################## Lsolve NODE "ID"\n", c)) ;
00041 
00042     /* ---------------------------------------------------------------------- */
00043     /* get the symbolic analysis of this node */
00044     /* ---------------------------------------------------------------------- */
00045 
00046     myid = Common->myid ;
00047     Childp = LUsymbolic->Childp ; /* Child [Childp [c]..Childp[c+1]-1] */
00048     Child  = LUsymbolic->Child ;  /* is list of children of node c */
00049     nchild = Childp [c+1] - Childp [c] ;    /* # of children of node c */
00050     Cn     = LUsymbolic->Cn ;   /* dimension of each node */
00051     Cstart = LUsymbolic->Cstart ;
00052     Sched = LUsymbolic->Sched ;
00053     Cparent = LUsymbolic->Cparent ;
00054 
00055     /* ---------------------------------------------------------------------- */
00056     /* get factors at this node */
00057     /* ---------------------------------------------------------------------- */
00058 
00059     LUnode = LU->LUnode [c] ;
00060 
00061     Lip = LUnode->lp ;
00062     Llen = LUnode->llen ;
00063     LUix = LUnode->ix ;
00064     nfound = LUnode->PK_NFOUND ;
00065 
00066     nlost_in = 0 ;
00067     Lost  = LUnode->Lost ;
00068     Lostp = LUnode->Lostp ;
00069     nlost_in = LUnode->nlost_in ;
00070 
00071     /* The matrix factorized at this node is cn-by-cn. */
00072     cn = nlost_in + Cn [c] ;
00073     k1 = Cstart [c] ;
00074     k2 = Cstart [c+1] ;
00075 
00076     Pinv = LUnode->Pinv ;
00077 
00078 #ifndef NDEBUG
00079     {
00080   Int npiv = LUnode->PK_NPIV ;
00081   Int *Plocal = LUnode->Plocal ;
00082   ASSERT (cn == LUnode->PK_NN) ;
00083   ASSERT (npiv == nlost_in + (k2 - k1)) ;
00084   for (k = 0 ; k < npiv ; k++)
00085   {
00086       i = Plocal [k] ;
00087       ASSERT (i >= 0 && i < npiv) ;
00088       ASSERT (Pinv [i] == k) ;
00089   }
00090   DEBUG (CHOLMOD (print_perm) (Plocal, npiv, npiv, "Plocal at node c", cm));
00091   DEBUG (CHOLMOD (print_perm) (Pinv,   npiv, npiv, "Pinv at node c", cm)) ;
00092     }
00093 #endif
00094 
00095     /* ---------------------------------------------------------------------- */
00096     /* add up all the contributions to the right-hand-side */
00097     /* ---------------------------------------------------------------------- */
00098 
00099     X = LU->LUnode [c]->X ;
00100     PR1 ((Common->file, "Lsolve at Node "ID", cn "ID"\n", c, cn)) ;
00101     for (i = 0 ; i < cn ; i++)
00102     {
00103   X [i] = 0 ;
00104     }
00105 
00106     /* B [0..(k2-k1)-1] contains the right-hand-side corresponding to rows
00107      * Cstart [c] to Cstart [c+1]-1 of the global system (after permuted by
00108      * Cperm).
00109      */
00110     B = LUnode->B ;
00111 
00112     /* wait for B to arrive at this node */
00113     MPI (MPI_Wait (&(LU->LUnode [c]->req), &ms)) ;
00114 
00115     for (i = 0 ; i < k2-k1 ; i++)
00116     {
00117   ci = i + nlost_in ;
00118   k = Pinv [ci] ;
00119   PR2 ((Common->file, "orig B ["ID"] = %g goes to X ["ID"]\n", i, B [i], k)) ;
00120   ASSERT (k >= 0 && k < cn) ;
00121   X [k] = B [i] ;
00122     }
00123 
00124     /* scatter and add the contributions from each child */
00125     for (pass = 1 ; pass <= 2 ; pass++)
00126     {
00127   for (cp = 0 ; cp < nchild ; cp++)
00128   {
00129       child = Child [Childp [c] + cp] ;
00130       if (pass == 1)
00131       {
00132     /* do local contributions on first pass */
00133     if (Sched [child] != myid) continue ;
00134       }
00135       else
00136       {
00137     /* do non-local contributions on second pass */
00138     if (Sched [child] == myid) continue ;
00139     MPI (MPI_Wait (&(LUnode->Req [cp]), &ms)) ;
00140       }
00141 
00142       cn_child = LU->LUnode [child]->PK_NN ;
00143       cn_nfound = LU->LUnode [child]->PK_NFOUND ;
00144       Xchild = LU->LUnode [child]->X + cn_nfound ;
00145 
00146       PR1 ((Common->file, "child "ID" cn "ID" nfound "ID"\n",
00147         child, cn_child, cn_nfound)) ;
00148 
00149       /* get the contributions of child to its lost pivot rows */
00150       for (i = 0 ; i < Lost [cp] ; i++)
00151       {
00152     ci = i + Lostp [cp] ; /* ci is now "original" local index */
00153     k = Pinv [ci] ;   /* kth local pivot row */
00154     PR2 ((Common->file, "Xchild ["ID"] = %g goes to X ["ID"] (lost)\n",
00155       i, Xchild [i], k)) ;
00156     ASSERT (k >= 0 && k < cn) ;
00157     X [k] += Xchild [i] ;
00158       }
00159 
00160       /* get the contributions to candidate pivot rows of node c */
00161       for ( ; i < Lost [cp] + (k2-k1) ; i++)
00162       {
00163     ci = i + (nlost_in - Lost [cp]) ;
00164     k = Pinv [ci] ;   /* kth local pivot row */
00165     PR2 ((Common->file, "Xchild ["ID"] = %g goes to X ["ID"] (cand)\n",
00166       i, Xchild [i], k)) ;
00167     ASSERT (k >= 0 && k < cn) ;
00168     X [k] += Xchild [i] ;
00169       }
00170 
00171       /* get contributions to candidate pivot rows of ancestors of c */
00172       for ( ; i < cn_child - cn_nfound ; i++)
00173       {
00174     k = i + (nlost_in - Lost [cp]) ;
00175     PR2 ((Common->file, "Xchild ["ID"] = %g goes to X ["ID"] (anc)\n",
00176       i, Xchild [i], k)) ;
00177     ASSERT (k >= 0 && k < cn) ;
00178     X [k] += Xchild [i] ;
00179       }
00180   }
00181     }
00182 
00183     /* ---------------------------------------------------------------------- */
00184     /* solve Lx=b */
00185     /* ---------------------------------------------------------------------- */
00186 
00187     for (j = 0 ; j < nfound ; j++)
00188     {
00189   xj = X [j] ;
00190   GET_COLUMN (Lip, Llen, LUix, j, Li, Lx, llen) ;
00191   for (p = 1 ; p < llen ; p++)
00192   {
00193       X [Li [p]] -= Lx [p] * xj ;
00194   }
00195     }
00196 
00197     /* ---------------------------------------------------------------------- */
00198     /* send results to parent */
00199     /* ---------------------------------------------------------------------- */
00200 
00201     /* X [0..nfound-1] is part of the global solution, X [nfound..nn-1] is
00202      * the contribution to the parent and ancestors */
00203 
00204     parent = Cparent [c] ;
00205     if (parent != EMPTY && Sched [parent] != myid)
00206     {
00207   MPI (MPI_Isend (X, cn, MPI_DOUBLE, Sched [parent], TAG0, MPI_COMM_WORLD,
00208       &req)) ;
00209 
00210   /* this request can be freed, because this process is now done, and 
00211      paraklete_lsolve_node is followed by a barrier */
00212   MPI (MPI_Request_free (&req)) ;
00213     }
00214 
00215     return (TRUE) ;
00216 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines