Amesos Package Browser (Single Doxygen Collection) Development
amesos_paraklete_usolve_node.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === paraklete_usolve_node ================================================ */
00003 /* ========================================================================== */
00004 
00005 #include "amesos_paraklete_decl.h"
00006 
00007 /* Solve Ux=y 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_usolve_node
00016 (
00017     Int c,
00018     paraklete_numeric *LU,
00019     paraklete_symbolic *LUsymbolic,
00020     paraklete_common *Common
00021 )
00022 {
00023     double xj ;
00024     paraklete_node *LUnode ;
00025     cholmod_common *cm ;
00026     double *LUix, *Xchild, *Ux, *X ;
00027     Int *Child, *Childp, *Lost, *Lostp, *Cstart, *Uip, *Ulen,
00028   *Cn, *Qinv, *Ui, *Sched, *Cparent ;
00029     Int cp, nchild, child, cn, k, j, p, i, nfound, k1, k2, ulen,
00030   ci, nlost_in, cn_child, cn_nfound, myid, parent ;
00031     MPI (MPI_Status ms) ;
00032     MPI (MPI_Request req) ;
00033 
00034     /* ---------------------------------------------------------------------- */
00035     /* get local workspace */
00036     /* ---------------------------------------------------------------------- */
00037 
00038     cm = &(Common->cm) ;
00039     PR0 ((Common->file, "\n\n########################### Usolve NODE "ID"\n", c));
00040 
00041     /* ---------------------------------------------------------------------- */
00042     /* get the symbolic analysis of this node */
00043     /* ---------------------------------------------------------------------- */
00044 
00045     myid = Common->myid ;
00046     Childp = LUsymbolic->Childp ; /* Child [Childp [c]..Childp[c+1]-1] */
00047     Child  = LUsymbolic->Child ;  /* is list of children of node c */
00048     nchild = Childp [c+1] - Childp [c] ;    /* # of children of node c */
00049     Cn     = LUsymbolic->Cn ;   /* dimension of each node */
00050     Cstart = LUsymbolic->Cstart ;
00051     Sched = LUsymbolic->Sched ;
00052     Cparent = LUsymbolic->Cparent ;
00053 
00054     /* ---------------------------------------------------------------------- */
00055     /* get solution to Ly=Pb, and solution to Ux=y from parent */
00056     /* ---------------------------------------------------------------------- */
00057 
00058     LUnode = LU->LUnode [c] ;
00059     cn = LUnode->PK_NN ;
00060     nfound = LUnode->PK_NFOUND ;
00061     X = LUnode->X ;
00062     parent = Cparent [c] ;
00063 
00064     PR1 ((Common->file, "Usolve node "ID" cn "ID" nfound "ID"\n", c, cn, nfound)) ;
00065 
00066     if (parent != EMPTY && Sched [parent] != myid)
00067     {
00068   PR1 ((Common->file, "Recving usolve from "ID", size "ID"\n", Sched [parent],
00069     cn - nfound)) ;
00070   MPI (MPI_Irecv (X + nfound, cn-nfound, MPI_DOUBLE, Sched [parent],
00071     TAG0, MPI_COMM_WORLD, &req)) ;
00072   MPI (MPI_Wait (&req, &ms)) ;
00073     }
00074 
00075     /* ---------------------------------------------------------------------- */
00076     /* solve Ux=y at this node */
00077     /* ---------------------------------------------------------------------- */
00078 
00079     Uip = LUnode->up ;
00080     Ulen = LUnode->ulen ;
00081     LUix = LUnode->ix ;
00082 
00083     /* X1 = U2*X2 - X1 */
00084     for (j = cn-1 ; j >= nfound ; j--)
00085     {
00086   xj = X [j] ;
00087   GET_COLUMN (Uip, Ulen, LUix, j, Ui, Ux, ulen) ;
00088   for (p = 0 ; p < ulen ; p++)
00089   {
00090       X [Ui [p]] -= Ux [p] * xj ;
00091   }
00092     }
00093 
00094     /* solve U*X1 */
00095     for ( ; j >= 0 ; j--)
00096     {
00097   GET_COLUMN (Uip, Ulen, LUix, j, Ui, Ux, ulen) ;
00098   xj = X [j] / Ux [ulen-1] ;
00099   X [j] = xj ;
00100   for (p = 0 ; p < ulen-1 ; p++)
00101   {
00102       X [Ui [p]] -= Ux [p] * xj ;
00103   }
00104     }
00105 
00106     /* ---------------------------------------------------------------------- */
00107     /* get factors at this node */
00108     /* ---------------------------------------------------------------------- */
00109 
00110     nlost_in = 0 ;
00111     Lost  = LUnode->Lost ;
00112     Lostp = LUnode->Lostp ;
00113     nlost_in = LUnode->nlost_in ;
00114 
00115     /* The matrix factorized at this node is cn-by-cn. */
00116     cn = nlost_in + Cn [c] ;
00117     k1 = Cstart [c] ;
00118     k2 = Cstart [c+1] ;
00119 
00120     PR1 ((Common->file, "NODE "ID" nlost_in: "ID"\n", c, nlost_in)) ;
00121 
00122     Qinv = LUnode->Qinv ;
00123 
00124 #ifndef NDEBUG
00125     {
00126   Int npiv = LUnode->PK_NPIV ;
00127   Int *Qlocal = LUnode->Qlocal ;
00128   ASSERT (cn == LUnode->PK_NN) ;
00129   ASSERT (npiv == nlost_in + (k2 - k1)) ;
00130   for (k = 0 ; k < npiv ; k++)
00131   {
00132       i = Qlocal [k] ;
00133       ASSERT (i >= 0 && i < npiv) ;
00134       ASSERT (Qinv [i] == k) ;
00135   }
00136   DEBUG (CHOLMOD (print_perm) (Qlocal, npiv, npiv, "Qlocal at node c", cm));
00137   DEBUG (CHOLMOD (print_perm) (Qinv,   npiv, npiv, "Qinv at node c", cm)) ;
00138     }
00139 #endif
00140 
00141     /* ---------------------------------------------------------------------- */
00142     /* send solutions to each child */
00143     /* ---------------------------------------------------------------------- */
00144 
00145     for (cp = 0 ; cp < nchild ; cp++)
00146     {
00147   child = Child [Childp [c] + cp] ;
00148   cn_child = LU->LUnode [child]->PK_NN ;
00149   cn_nfound = LU->LUnode [child]->PK_NFOUND ;
00150   Xchild = LU->LUnode [child]->X + cn_nfound ;
00151 
00152   PR1 ((Common->file, "child "ID" cn "ID" nfound "ID"\n",
00153     child, cn_child, cn_nfound)) ;
00154 
00155   /* send the child its lost pivot cols */
00156   for (i = 0 ; i < Lost [cp] ; i++)
00157   {
00158       ci = i + Lostp [cp] ; /* ci is now "original" local index */
00159       k = Qinv [ci] ;   /* kth local pivot col */
00160       PR2 ((Common->file, "Xchild ["ID"] = %g from X ["ID"] (lost)\n",
00161         i, Xchild [i], k)) ;
00162       Xchild [i] = X [k] ;
00163   }
00164 
00165   /* send the child the rest of the solution from c */
00166   for ( ; i < Lost [cp] + (k2-k1) ; i++)
00167   {
00168       ci = i + (nlost_in - Lost [cp]) ;
00169       k = Qinv [ci] ;   /* kth local pivot col */
00170       PR2 ((Common->file, "Xchild ["ID"] = %g from X ["ID"] (cand)\n",
00171         i, Xchild [i], k)) ;
00172       Xchild [i] = X [k] ;
00173   }
00174 
00175   /* get the solutions of ancestors of c */
00176   for ( ; i < cn_child - cn_nfound ; i++)
00177   {
00178       k = i + (nlost_in - Lost [cp]) ;
00179       PR2 ((Common->file, "Xchild ["ID"] = %g from X ["ID"] (anc)\n",
00180         i, Xchild [i], k)) ;
00181       Xchild [i] = X [k] ;
00182   }
00183 
00184   PR1 ((Common->file, "Usolve: Sending from "ID" to child "ID"\n", c, child));
00185   if (Sched [child] != myid)
00186   {
00187       PR1 ((Common->file, "Sending to "ID", size "ID"\n", Sched [child], 
00188         cn_child - cn_nfound)) ;
00189       MPI (MPI_Isend (Xchild, cn_child - cn_nfound, MPI_DOUBLE,
00190       Sched [child], TAG0, MPI_COMM_WORLD, &req)) ;
00191       /* this request can be freed, because the work of this node is now
00192          done, and paraklete_usolve_node is followed by a barrier */
00193       MPI (MPI_Request_free (&req)) ;
00194   }
00195     }
00196 
00197     return (TRUE) ;
00198 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines