Amesos Package Browser (Single Doxygen Collection) Development
amesos_paraklete_reanalyze.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === paraklete_reanalyze ================================================== */
00003 /* ========================================================================== */
00004 
00005 /* paraklete_analyze and paraklete_factorize have already been called.
00006  * This function constructs a new symbolic object that includes the original
00007  * fill-reducing nested dissection ordering (Cperm) and separator tree
00008  * from paraklete_analyze, combined with the partial-pivoting permutation
00009  * and lost-pivots from paraklete_factorize.
00010  *
00011  * All processes should do this, independently, since all processes have their
00012  * own complete copy of the LUsymbolic object.  Each process also has its own
00013  * copy of LU->P and LU->Q, the final row and column permutations from
00014  * paraklete_factorize.
00015  *
00016  * Everyone has the global P, Q, and LU->LUnode [*]->header.  The assignment
00017  * of rows/columns to the nodes can be finalized, to reflect the final
00018  * permutation.  All processes compute the new schedule.
00019  *
00020  * PARAKLETE version 0.3: parallel sparse LU factorization.  Nov 13, 2007
00021  * Copyright (C) 2007, Univ. of Florida.  Author: Timothy A. Davis
00022  * See License.txt for the Version 2.1 of the GNU Lesser General Public License
00023  * http://www.cise.ufl.edu/research/sparse
00024  */
00025 
00026 #include "amesos_paraklete_decl.h"
00027 
00028 paraklete_symbolic *amesos_paraklete_reanalyze
00029 (
00030     cholmod_sparse *A,      /* only root processor owns this */
00031     paraklete_numeric *LU,
00032     paraklete_symbolic *LUsymbolic,
00033     paraklete_common *Common
00034 )
00035 {
00036     paraklete_symbolic *LUsymbolic_new ;
00037     Int *Cparent, *Cstart, *Cn, *Cmember, *Pmember, *P, *Q, *Cnz, *Ap, *Ai ;
00038     Int i, j, k, c, n, ncomponents, jold, iold, parent, nparent, ci, cj, p ;
00039 
00040     /* TODO check inputs to make sure they are not NULL */
00041 
00042     /* ---------------------------------------------------------------------- */
00043     /* allocate the new symbolic object */
00044     /* ---------------------------------------------------------------------- */
00045 
00046     n = LUsymbolic->n ;
00047     ncomponents = LUsymbolic->ncomponents ;
00048     LUsymbolic_new = amesos_paraklete_alloc_symbolic (n, ncomponents, TRUE, Common) ;
00049 
00050     /* TODO check for out-of-memory here */
00051 
00052     /* ---------------------------------------------------------------------- */
00053     /* copy the parts of the separator tree that do not change with pivoting */
00054     /* ---------------------------------------------------------------------- */
00055 
00056     for (c = 0 ; c < ncomponents ; c++)
00057     {
00058   LUsymbolic_new->Child [c]  = LUsymbolic->Child [c] ;
00059     }
00060 
00061     for (c = 0 ; c <= ncomponents ; c++)
00062     {
00063   LUsymbolic_new->Childp [c] = LUsymbolic->Childp [c] ;
00064     }
00065 
00066     for (c = 0 ; c < ncomponents ; c++)
00067     {
00068   LUsymbolic_new->Sched [c] = LUsymbolic->Sched [c] ;
00069     }
00070 
00071     Cparent = LUsymbolic_new->Cparent ;
00072     for (c = 0 ; c < ncomponents ; c++)
00073     {
00074   Cparent [c] = LUsymbolic->Cparent [c] ;
00075     }
00076 
00077     /* ---------------------------------------------------------------------- */
00078     /* determine actual number of entries in LU factors at each node */
00079     /* ---------------------------------------------------------------------- */
00080 
00081     for (c = 0 ; c < ncomponents ; c++)
00082     {
00083   LUsymbolic_new->Clnz [c] = LU->LUnode [c]->lnz + LU->LUnode [c]->unz ;
00084     }
00085 
00086     /* ---------------------------------------------------------------------- */
00087     /* find the range of nodes in P*A*Q assigned to each node */
00088     /* ---------------------------------------------------------------------- */
00089 
00090     Cstart = LUsymbolic_new->Cstart ;
00091     Cstart [0] = 0 ;
00092     for (c = 0 ; c < ncomponents ; c++)
00093     {
00094   Cstart [c+1] = Cstart [c] + LU->LUnode [c]-> PK_NFOUND ;
00095     }
00096 
00097     /*
00098     if (Common->myid == 0)
00099   for (c = 0 ; c <= ncomponents ; c++)
00100       printf ("new Cstart "ID"\n", Cstart [c]) ;
00101     */
00102 
00103     /* ---------------------------------------------------------------------- */
00104     /* find size of C matrix at each node (all pivoting accounted for) */
00105     /* ---------------------------------------------------------------------- */
00106 
00107     Cn = LUsymbolic_new->Cn ;
00108     for (c = ncomponents - 1 ; c >= 0 ; c--)
00109     {
00110   parent = Cparent [c] ;
00111   nparent = (parent == EMPTY) ? 0 : Cn [parent] ;
00112   Cn [c] = (Cstart [c+1] - Cstart [c]) + nparent ;
00113   PR1 ((Common->file, "node "ID" new Cn: "ID"\n", c, Cn [c])) ;
00114     }
00115 
00116     /* ---------------------------------------------------------------------- */
00117     /* determine Cnz for each node */
00118     /* ---------------------------------------------------------------------- */
00119 
00120     Cmember = LUsymbolic_new->Cperm ; /* use Cperm as workspace */
00121     Pmember = LUsymbolic_new->Rperm ; /* use Rperm as workspace */
00122 
00123     Cnz = LUsymbolic_new->Cnz ;   /* new nonzero count for each node */
00124     Q = LU->Q ;       /* final column ordering, P*A*Q=L*U */
00125     P = LU->P ;       /* final row ordering */
00126 
00127     if (Common->myid == 0)
00128     {
00129   for (c = 0 ; c < ncomponents ; c++)
00130   {
00131       Cnz [c] = 0 ;
00132   }
00133 
00134   /* Cmember [k] = c if the kth diagonal entry of P*A*Q is in node c */
00135   for (c = 0 ; c < ncomponents ; c++)
00136   {
00137       for (k = Cstart [c] ; k < Cstart [c+1] ; k++)
00138       {
00139     Cmember [k] = c ;
00140       }
00141   }
00142 
00143   /*
00144   for (k = 0 ; k < n ; k++)
00145   {
00146       printf ("Cmember ["ID"] = "ID"\n", k, Cmember [k]) ;
00147   }
00148   for (k = 0 ; k < n ; k++)
00149   {
00150       printf ("P ["ID"] = "ID"\n", k, P [k]) ;
00151   }
00152   */
00153 
00154   /* Pmember [i] = c if A(i,i) becomes a pivot in node c */
00155   for (k = 0 ; k < n ; k++)
00156   {
00157       i = P [k] ;   /* kth row of P*A*Q is ith row of A */
00158       c = Cmember [k] ; /* kth row of P*A*Q is in node c */
00159       Pmember [i] = c ; /* A(i,i) becomes (P*A*Q)_kk, in node c */
00160   }
00161 
00162   /*
00163   for (i = 0 ; i < n ; i++)
00164   {
00165       printf ("Pmember ["ID"] = "ID"\n", i, Pmember [i]) ;
00166   }
00167   */
00168 
00169   /* count the entries in each node */
00170   Ap = A->p ;
00171   Ai = A->i ;
00172   for (j = 0 ; j < n ; j++)
00173   {
00174       jold = Q [j] ;  /* jth column of P*A*Q is column jold of A */
00175       cj = Cmember [j] ;  /* jth diagonal entry of P*A*Q in node cj */
00176       for (p = Ap [jold] ; p < Ap [jold+1] ; p++)
00177       {
00178     iold = Ai [p] ;   /* A(iold,jold) in original matrix */
00179     ci = Pmember [iold] ; /* A(iold,iold) in node ci */
00180     c = MIN (ci, cj) ;  /* determine node of (P*A*Q)_ij */
00181 
00182     /*
00183     printf ("A("ID","ID") PAQ("ID","ID") ci "ID" cj "ID" c "ID"\n",
00184       iold, jold, Ai [p], j, ci, cj, c) ;
00185     */
00186     Cnz [c]++ ;
00187       }
00188   }
00189 
00190   /*
00191   for (c = 0 ; c < ncomponents ; c++)
00192   {
00193       printf ("component "ID"\n", c) ;
00194       printf ("   new Cn "ID" Cnz "ID" Cstart "ID"\n",
00195       Cn [c], Cnz [c], Cstart [c]) ;
00196   }
00197   */
00198 
00199     }
00200 
00201     /* Cmember and Pmember workspace no longer needed */
00202 
00203     /* broadcast results to all processors */
00204     MPI (MPI_Bcast (Cnz, ncomponents, MPI_Int, TAG0, MPI_COMM_WORLD)) ;
00205 
00206     /* ---------------------------------------------------------------------- */
00207     /* copy the final P and Q */
00208     /* ---------------------------------------------------------------------- */
00209 
00210     /* column permutation Q */
00211     for (k = 0 ; k < n ; k++)
00212     {
00213   LUsymbolic_new->Cperm [k] = Q [k] ;
00214     }
00215 
00216     /* row permutation P and its inverse */
00217     for (k = 0 ; k < n ; k++)
00218     {
00219   i = P [k] ;
00220   LUsymbolic_new->RpermInv [i] = k ;
00221   LUsymbolic_new->Rperm [k] = i ;
00222     }
00223 
00224     return (LUsymbolic_new) ;
00225 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines