Amesos Package Browser (Single Doxygen Collection) Development
amesos_paraklete_solve.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === paraklete_solve ====================================================== */
00003 /* ========================================================================== */
00004 
00005 #include "amesos_paraklete_decl.h"
00006 
00007 /* paraklete_solve (LU, LUsymbolic, B, Common) solves Ly=Pb, where P is
00008  * the initial fill-reducing ordering (P=Cperm).
00009  *
00010  * The solution is left in LU->LUnode [c]->X, for each node c.  Next, it solves
00011  * UQ'x=y, where the solution is written into the input array B,
00012  * and where Q is the final column ordering.
00013  *
00014  * This solve is performed using the original LUsymbolic object, not the one
00015  * returned from paraklete_reanalyze.
00016  *
00017  * PARAKLETE version 0.3: parallel sparse LU factorization.  Nov 13, 2007
00018  * Copyright (C) 2007, Univ. of Florida.  Author: Timothy A. Davis
00019  * See License.txt for the Version 2.1 of the GNU Lesser General Public License
00020  * http://www.cise.ufl.edu/research/sparse
00021  */
00022 
00023 Int amesos_paraklete_solve
00024 (
00025     /* inputs, not modified */
00026     paraklete_numeric *LU,
00027     paraklete_symbolic *LUsymbolic,
00028     double *B,
00029     paraklete_common *Common
00030 )
00031 {
00032     cholmod_common *cm ;
00033     double *Blocal, *X, *W, *X2 ;
00034     Int *Cperm, *Rperm, *Cstart, *Q, *Sched, *Child, *Childp ;
00035     Int i, n, ncomponents, k, c, k1, k2, nfound, myid, cp, nchild, child ;
00036     MPI (MPI_Status ms) ;
00037 
00038     /* ---------------------------------------------------------------------- */
00039     /* get inputs */
00040     /* ---------------------------------------------------------------------- */
00041 
00042     cm = &(Common->cm) ;
00043     n = LU->n ;
00044     W = LU->W ;     /* only non-NULL on the root process 0 */
00045     myid = Common->myid ;
00046 
00047     ncomponents = LUsymbolic->ncomponents ;
00048     Cperm  = LUsymbolic->Cperm ;
00049     Rperm  = LUsymbolic->Rperm ;
00050     Rperm  = (Rperm == NULL) ? Cperm : Rperm ;
00051     Cstart = LUsymbolic->Cstart ;
00052     Sched  = LUsymbolic->Sched ;
00053 
00054     /* ---------------------------------------------------------------------- */
00055     /* W = Rperm*B, on process 0 only */
00056     /* ---------------------------------------------------------------------- */
00057 
00058     if (myid == 0)
00059     {
00060   for (k = 0 ; k < n ; k++)
00061   {
00062       W [k] = B [Rperm [k]] ;
00063   }
00064     }
00065 
00066     /* ---------------------------------------------------------------------- */
00067     /* distribute the permuted B to the nodes */
00068     /* ---------------------------------------------------------------------- */
00069 
00070     for (c = 0 ; c < ncomponents ; c++)
00071     {
00072   k1 = Cstart [c] ;
00073   k2 = Cstart [c+1] ;
00074 
00075   Blocal = LU->LUnode [c]->B ;
00076 
00077   /* send Blocal to node c */
00078   MPI (LU->LUnode [c]->req = MPI_REQUEST_NULL) ;
00079   if (myid == 0)
00080   {
00081       if (Sched [c] == myid)
00082       {
00083     for (i = 0 ; i < k2-k1 ; i++)
00084     {
00085         Blocal [i] = W [k1 + i] ;
00086     }
00087       }
00088       else
00089       {
00090     MPI (MPI_Isend (W + k1, k2-k1, MPI_DOUBLE, Sched [c],
00091       /* TAG: */ c, MPI_COMM_WORLD, &(LU->LUnode [c]->req))) ;
00092       }
00093   }
00094   else
00095   {
00096       if (Sched [c] == myid)
00097       {
00098     MPI (MPI_Irecv (Blocal, k2-k1, MPI_DOUBLE, 0,
00099       /* TAG: */ c, MPI_COMM_WORLD, &(LU->LUnode [c]->req))) ;
00100       }
00101   }
00102     }
00103 
00104     /* ---------------------------------------------------------------------- */
00105     /* forward solve: post a receive for X from each non-local child */
00106     /* ---------------------------------------------------------------------- */
00107 
00108     Childp = LUsymbolic->Childp ; /* Child [Childp [c]..Childp[c+1]-1] */
00109     Child  = LUsymbolic->Child ;  /* is list of children of node c */
00110 
00111     for (c = 0 ; c < ncomponents ; c++)
00112     {
00113   if (Sched [c] == myid)
00114   {
00115       nchild = Childp [c+1] - Childp [c] ;
00116       for (cp = 0 ; cp < nchild ; cp++)
00117       {
00118     child = Child [Childp [c] + cp] ;
00119     if (Sched [child] != myid)
00120     {
00121         MPI (MPI_Irecv (LU->LUnode [child]->X,
00122       LU->LUnode [child]->PK_NN, MPI_DOUBLE, Sched [child],
00123       TAG0, MPI_COMM_WORLD, &(LU->LUnode [c]->Req [cp]))) ;
00124     }
00125     else
00126     {
00127         MPI (LU->LUnode [c]->Req [cp] = MPI_REQUEST_NULL) ;
00128     }
00129       }
00130   }
00131     }
00132 
00133     /* ---------------------------------------------------------------------- */
00134     /* forward solve: Ly=b */
00135     /* ---------------------------------------------------------------------- */
00136 
00137     for (c = 0 ; c < ncomponents ; c++)
00138     {
00139   if (Sched [c] == myid)
00140   {
00141       amesos_paraklete_lsolve_node (c, LU, LUsymbolic, Common) ;
00142   }
00143     }
00144 
00145     MPI (MPI_Barrier (MPI_COMM_WORLD)) ;
00146 
00147     /* ---------------------------------------------------------------------- */
00148     /* backsolve: Ux=y */
00149     /* ---------------------------------------------------------------------- */
00150 
00151     for (c = ncomponents-1 ; c >= 0 ; c--)
00152     {
00153   if (Sched [c] == myid)
00154   {
00155       amesos_paraklete_usolve_node (c, LU, LUsymbolic, Common) ;
00156   }
00157     }
00158 
00159     MPI (MPI_Barrier (MPI_COMM_WORLD)) ;
00160 
00161     /* ---------------------------------------------------------------------- */
00162     /* get the permuted solution from each node */
00163     /* ---------------------------------------------------------------------- */
00164 
00165     Q = LU->Q ;
00166     k = 0 ;
00167     for (c = 0 ; c < ncomponents ; c++)
00168     {
00169   X = LU->LUnode [c]->X ;
00170   nfound = LU->LUnode [c]->PK_NFOUND ;
00171   if (myid == 0)
00172   {
00173       PR1 ((Common->file, "get soln, node c="ID", nfound "ID"\n", c, nfound));
00174       /* get X from node c */
00175       if (Sched [c] != myid)
00176       {
00177     PR1 ((Common->file, "recv node "ID" from "ID"\n", c, Sched [c])) ;
00178     MPI (MPI_Recv (W, nfound, MPI_DOUBLE, Sched [c],
00179           TAG0, MPI_COMM_WORLD, &ms)) ;
00180     X2 = W ;
00181       }
00182       else
00183       {
00184     PR1 ((Common->file, "I own it already\n")) ;
00185     X2 = X ;
00186       }
00187       PR1 ((Common->file, "got it from Sched [c] = "ID"\n", Sched [c])) ;
00188       for (i = 0 ; i < nfound ; i++)
00189       {
00190     B [Q [k]] = X2 [i] ;
00191     PR2 ((Common->file, "X ["ID"] is global B ["ID"] %g\n",
00192         i, k, X2 [i])) ;
00193     k++ ;
00194       }
00195   }
00196   else
00197   {
00198       if (Sched [c] == myid)
00199       {
00200     PR1 ((Common->file,
00201         "send soln, node c = "ID", myid "ID" nfound "ID"\n",
00202         c, myid, nfound)) ;
00203     MPI (MPI_Send (X, nfound, MPI_DOUBLE, 0, TAG0, MPI_COMM_WORLD));
00204       }
00205   }
00206     }
00207 
00208     return (TRUE) ;
00209 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines