ilu_mpi_pilu.c

00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2009) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #include "Euclid_dh.h"
00031 #include "Factor_dh.h"
00032 #include "Mat_dh.h"
00033 #include "ilu_dh.h"
00034 #include "Mem_dh.h"
00035 #include "Parser_dh.h"
00036 #include "Hash_dh.h"
00037 #include "getRow_dh.h"
00038 #include "SortedList_dh.h"
00039 #include "ExternalRows_dh.h"
00040 #include "SubdomainGraph_dh.h"
00041 
00042 
00043 static void iluk_symbolic_row_private (int localRow, int len, int *CVAL,
00044                        double *AVAL, ExternalRows_dh extRows,
00045                        SortedList_dh sList, Euclid_dh ctx,
00046                        bool debug);
00047 
00048 static void iluk_numeric_row_private (int new_row, ExternalRows_dh extRows,
00049                       SortedList_dh slist, Euclid_dh ctx,
00050                       bool debug);
00051 
00052 #undef __FUNC__
00053 #define __FUNC__ "iluk_mpi_pilu"
00054 void
00055 iluk_mpi_pilu (Euclid_dh ctx)
00056 {
00057   START_FUNC_DH int from = ctx->from, to = ctx->to;
00058   int i, m;
00059   int *n2o_row;         /* *o2n_col; */
00060   int *rp, *cval, *diag, *fill;
00061   int beg_row, beg_rowP, end_rowP;
00062   SubdomainGraph_dh sg = ctx->sg;
00063   int *CVAL, len, idx = 0, count;
00064   double *AVAL;
00065   REAL_DH *aval;
00066   Factor_dh F = ctx->F;
00067   SortedList_dh slist = ctx->slist;
00068   ExternalRows_dh extRows = ctx->extRows;
00069   bool bj, noValues, debug = false;
00070 
00071   /* for debugging */
00072   if (logFile != NULL && Parser_dhHasSwitch (parser_dh, "-debug_ilu"))
00073     debug = true;
00074   noValues = Parser_dhHasSwitch (parser_dh, "-noValues");
00075   bj = ctx->F->blockJacobi;
00076 
00077   m = F->m;
00078   rp = F->rp;
00079   cval = F->cval;
00080   fill = F->fill;
00081   diag = F->diag;
00082   aval = F->aval;
00083   /* work = ctx->work; */
00084 
00085   n2o_row = sg->n2o_row;
00086   /* o2n_col = sg->o2n_col; */
00087 
00088   if (from != 0)
00089     idx = rp[from];
00090 
00091   /* global numbers of first and last locally owned rows,
00092      with respect to A 
00093    */
00094   beg_row = sg->beg_row[myid_dh];
00095   /* end_row  = beg_row + sg->row_count[myid_dh]; */
00096 
00097   /* global number or first locally owned row, after reordering */
00098   beg_rowP = sg->beg_rowP[myid_dh];
00099   end_rowP = beg_rowP + sg->row_count[myid_dh];
00100 
00101 
00102   /* loop over rows to be factored (i references local rows) */
00103   for (i = from; i < to; ++i)
00104     {
00105 
00106       int row = n2o_row[i]; /* local row number */
00107       int globalRow = row + beg_row;    /* global row number */
00108 
00109       if (debug)
00110     {
00111       fprintf (logFile,
00112            "\nILU_pilu global: %i  old_Local: %i =========================================================\n",
00113            i + 1 + beg_rowP, row + 1);
00114     }
00115 
00116       EuclidGetRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
00117       CHECK_V_ERROR;
00118 
00119       if (debug)
00120     {
00121       int h;
00122       fprintf (logFile, "ILU_pilu  EuclidGetRow:\n");
00123       for (h = 0; h < len; ++h)
00124         fprintf (logFile, "    %i   %g\n", 1 + CVAL[h], AVAL[h]);
00125     }
00126 
00127 
00128       /* compute scaling value for row(i) */
00129       if (ctx->isScaled)
00130     {
00131       compute_scaling_private (i, len, AVAL, ctx);
00132       CHECK_V_ERROR;
00133     }
00134 
00135       SortedList_dhReset (slist, i);
00136       CHECK_V_ERROR;
00137 
00138       /* Compute symbolic factor for row(i);
00139          this also performs sparsification
00140        */
00141       iluk_symbolic_row_private (i, len, CVAL, AVAL,
00142                  extRows, slist, ctx, debug);
00143       CHECK_V_ERROR;
00144 
00145       /* enforce subdomain constraint */
00146       SortedList_dhEnforceConstraint (slist, sg);
00147 
00148       /* compute numeric factor for row */
00149       if (!noValues)
00150     {
00151       iluk_numeric_row_private (i, extRows, slist, ctx, debug);
00152       CHECK_V_ERROR;
00153     }
00154 
00155       EuclidRestoreRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
00156       CHECK_V_ERROR;
00157 
00158       /* Ensure adequate storage; reallocate, if necessary. */
00159       count = SortedList_dhReadCount (slist);
00160       CHECK_V_ERROR;
00161 
00162       /* Ensure adequate storage; reallocate, if necessary. */
00163       if (idx + count > F->alloc)
00164     {
00165       Factor_dhReallocate (F, idx, count);
00166       CHECK_V_ERROR;
00167       SET_INFO ("REALLOCATED from ilu_mpi_pilu");
00168       cval = F->cval;
00169       fill = F->fill;
00170       aval = F->aval;
00171     }
00172 
00173       /* Copy factor to permanent storage */
00174       if (bj)
00175     {           /* for debugging: blockJacobi case */
00176       int col;
00177       while (count--)
00178         {
00179           SRecord *sr = SortedList_dhGetSmallest (slist);
00180           CHECK_V_ERROR;
00181           col = sr->col;
00182           if (col >= beg_rowP && col < end_rowP)
00183         {
00184           cval[idx] = col;
00185           if (noValues)
00186             {
00187               aval[idx] = 0.0;
00188             }
00189           else
00190             {
00191               aval[idx] = sr->val;
00192             }
00193           fill[idx] = sr->level;
00194           ++idx;
00195         }
00196         }
00197     }
00198 
00199       if (debug)
00200     {
00201       fprintf (logFile, "ILU_pilu  ");
00202       while (count--)
00203         {
00204           SRecord *sr = SortedList_dhGetSmallest (slist);
00205           CHECK_V_ERROR;
00206           cval[idx] = sr->col;
00207           aval[idx] = sr->val;
00208           fill[idx] = sr->level;
00209           fprintf (logFile, "%i,%i,%g ; ", 1 + cval[idx], fill[idx],
00210                aval[idx]);
00211           ++idx;
00212         }
00213       fprintf (logFile, "\n");
00214     }
00215 
00216       else
00217     {
00218       while (count--)
00219         {
00220           SRecord *sr = SortedList_dhGetSmallest (slist);
00221           CHECK_V_ERROR;
00222           cval[idx] = sr->col;
00223           aval[idx] = sr->val;
00224           fill[idx] = sr->level;
00225           ++idx;
00226         }
00227     }
00228 
00229       /* add row-pointer to start of next row. */
00230       rp[i + 1] = idx;
00231 
00232       /* Insert pointer to diagonal */
00233       {
00234     int temp = rp[i];
00235     bool flag = true;
00236     while (temp < idx)
00237       {
00238         if (cval[temp] == i + beg_rowP)
00239           {
00240         diag[i] = temp;
00241         flag = false;
00242         break;
00243           }
00244         ++temp;
00245       }
00246     if (flag)
00247       {
00248         if (logFile != NULL)
00249           {
00250         int k;
00251         fprintf (logFile,
00252              "Failed to find diag in localRow %i (globalRow %i; ct= %i)\n   ",
00253              1 + i, i + 1 + beg_rowP, rp[i + 1] - rp[i]);
00254         for (k = rp[i]; k < rp[i + 1]; ++k)
00255           {
00256             fprintf (logFile, "%i ", cval[i] + 1);
00257           }
00258         fprintf (logFile, "\n\n");
00259           }
00260         sprintf (msgBuf_dh, "failed to find diagonal for localRow: %i",
00261              1 + i);
00262         SET_V_ERROR (msgBuf_dh);
00263       }
00264       }
00265 /*
00266     { int temp = rp[i]; 
00267       while (cval[temp] != i+beg_row) ++temp;
00268       diag[i] = temp;
00269     }
00270 */
00271 
00272       /* check for zero diagonal */
00273       if (!aval[diag[i]])
00274     {
00275       sprintf (msgBuf_dh, "zero diagonal in local row %i", i + 1);
00276       SET_V_ERROR (msgBuf_dh);
00277     }
00278 
00279     }
00280 
00281   /* adjust to local (zero) based, if block jacobi factorization */
00282   if (bj)
00283     {
00284       int nz = rp[m];
00285       for (i = 0; i < nz; ++i)
00286     cval[i] -= beg_rowP;
00287     }
00288 
00289 END_FUNC_DH}
00290 
00291 
00292 #undef __FUNC__
00293 #define __FUNC__ "iluk_symbolic_row_private"
00294 void
00295 iluk_symbolic_row_private (int localRow, int len, int *CVAL,
00296                double *AVAL, ExternalRows_dh extRows,
00297                SortedList_dh slist, Euclid_dh ctx, bool debug)
00298 {
00299   START_FUNC_DH int level = ctx->level, m = ctx->m;
00300   int beg_row = ctx->sg->beg_row[myid_dh];
00301   int beg_rowP = ctx->sg->beg_rowP[myid_dh];
00302   int *cval = ctx->F->cval, *diag = ctx->F->diag;
00303   int *rp = ctx->F->rp, *fill = ctx->F->fill;
00304   int j, node, col;
00305   int end_rowP = beg_rowP + m;
00306   int level_1, level_2;
00307   int *cvalPtr, *fillPtr;
00308   SRecord sr, *srPtr;
00309   REAL_DH scale, *avalPtr;
00310   double thresh = ctx->sparseTolA;
00311   bool wasInserted;
00312   int count = 0;
00313 
00314   scale = ctx->scale[localRow];
00315   ctx->stats[NZA_STATS] += (double) len;
00316 
00317   /* insert col indices in sorted linked list */
00318   sr.level = 0;
00319   for (j = 0; j < len; ++j)
00320     {
00321       sr.col = CVAL[j];
00322       sr.val = scale * AVAL[j];
00323 /*    if (fabs(sr.val) > thresh) { */
00324       wasInserted = SortedList_dhPermuteAndInsert (slist, &sr, thresh);
00325       CHECK_V_ERROR;
00326       if (wasInserted)
00327     ++count;
00328 /*    } */
00329       if (debug)
00330     {
00331       fprintf (logFile, "ILU_pilu   inserted from A: col= %i  val= %g\n",
00332            1 + CVAL[j], sr.val);
00333     }
00334     }
00335 
00336   /* ensure diagonal entry is inserted */
00337   sr.val = 0.0;
00338   sr.col = localRow + beg_rowP;
00339   srPtr = SortedList_dhFind (slist, &sr);
00340   CHECK_V_ERROR;
00341   if (srPtr == NULL)
00342     {
00343       SortedList_dhInsert (slist, &sr);
00344       CHECK_V_ERROR;
00345       ++count;
00346       if (debug)
00347     {
00348       fprintf (logFile, "ILU_pilu   inserted missing diagonal: %i\n",
00349            1 + localRow + beg_row);
00350     }
00351     }
00352   ctx->stats[NZA_USED_STATS] += (double) count;
00353 
00354   /* update row from previously factored rows */
00355   sr.val = 0.0;
00356   if (level > 0)
00357     {
00358       while (1)
00359     {
00360       srPtr = SortedList_dhGetSmallestLowerTri (slist);
00361       CHECK_V_ERROR;
00362       if (srPtr == NULL)
00363         break;
00364 
00365       node = srPtr->col;
00366 
00367       if (debug)
00368         {
00369           fprintf (logFile, "ILU_pilu   sf updating from row: %i\n",
00370                1 + srPtr->col);
00371         }
00372 
00373       level_1 = srPtr->level;
00374       if (level_1 < level)
00375         {
00376 
00377           /* case 1: locally owned row */
00378           if (node >= beg_rowP && node < end_rowP)
00379         {
00380           node -= beg_rowP;
00381           len = rp[node + 1] - diag[node] - 1;
00382           cvalPtr = cval + diag[node] + 1;
00383           fillPtr = fill + diag[node] + 1;
00384         }
00385 
00386           /* case 2: external row */
00387           else
00388         {
00389           len = 0;
00390           ExternalRows_dhGetRow (extRows, node, &len, &cvalPtr,
00391                      &fillPtr, &avalPtr);
00392           CHECK_V_ERROR;
00393           if (debug && len == 0)
00394             {
00395               fprintf (stderr,
00396                    "ILU_pilu  sf failed to get extern row: %i\n",
00397                    1 + node);
00398             }
00399         }
00400 
00401 
00402           /* merge in strict upper triangular portion of row */
00403           for (j = 0; j < len; ++j)
00404         {
00405           col = *cvalPtr++;
00406           level_2 = 1 + level_1 + *fillPtr++;
00407           if (level_2 <= level)
00408             {
00409               /* Insert new element, or update level if already inserted. */
00410               sr.col = col;
00411               sr.level = level_2;
00412               sr.val = 0.0;
00413               SortedList_dhInsertOrUpdate (slist, &sr);
00414               CHECK_V_ERROR;
00415             }
00416         }
00417         }
00418     }
00419     }
00420 END_FUNC_DH}
00421 
00422 
00423 #undef __FUNC__
00424 #define __FUNC__ "iluk_numeric_row_private"
00425 void
00426 iluk_numeric_row_private (int new_row, ExternalRows_dh extRows,
00427               SortedList_dh slist, Euclid_dh ctx, bool debug)
00428 {
00429   START_FUNC_DH int m = ctx->m;
00430   int beg_rowP = ctx->sg->beg_rowP[myid_dh];
00431   int end_rowP = beg_rowP + m;
00432   int len, row;
00433   int *rp = ctx->F->rp, *cval = ctx->F->cval, *diag = ctx->F->diag;
00434   REAL_DH *avalPtr, *aval = ctx->F->aval;
00435   int *cvalPtr;
00436   double multiplier, pc, pv;
00437   SRecord sr, *srPtr;
00438 
00439   /* note: non-zero entries from A were inserted in list during iluk_symbolic_row_private */
00440 
00441   SortedList_dhResetGetSmallest (slist);
00442   CHECK_V_ERROR;
00443   while (1)
00444     {
00445       srPtr = SortedList_dhGetSmallestLowerTri (slist);
00446       CHECK_V_ERROR;
00447       if (srPtr == NULL)
00448     break;
00449 
00450       /* update new_row's values from upper triangular portion of previously
00451          factored row
00452        */
00453       row = srPtr->col;
00454 
00455       if (row >= beg_rowP && row < end_rowP)
00456     {
00457       int local_row = row - beg_rowP;
00458 
00459       len = rp[local_row + 1] - diag[local_row];
00460       cvalPtr = cval + diag[local_row];
00461       avalPtr = aval + diag[local_row];
00462     }
00463       else
00464     {
00465       len = 0;
00466       ExternalRows_dhGetRow (extRows, row, &len, &cvalPtr,
00467                  NULL, &avalPtr);
00468       CHECK_V_ERROR;
00469       if (debug && len == 0)
00470         {
00471           fprintf (stderr, "ILU_pilu  failed to get extern row: %i\n",
00472                1 + row);
00473         }
00474 
00475     }
00476 
00477       if (len)
00478     {
00479       /* first, form and store pivot */
00480       sr.col = row;
00481       srPtr = SortedList_dhFind (slist, &sr);
00482       CHECK_V_ERROR;
00483       if (srPtr == NULL)
00484         {
00485           sprintf (msgBuf_dh,
00486                "find failed for sr.col = %i while factoring local row= %i \n",
00487                1 + sr.col, new_row + 1);
00488           SET_V_ERROR (msgBuf_dh);
00489         }
00490 
00491       pc = srPtr->val;
00492 
00493       if (pc != 0.0)
00494         {
00495           pv = *avalPtr++;
00496           --len;
00497           ++cvalPtr;
00498           multiplier = pc / pv;
00499           srPtr->val = multiplier;
00500 
00501           if (debug)
00502         {
00503           fprintf (logFile,
00504                "ILU_pilu   nf updating from row: %i; multiplier = %g\n",
00505                1 + srPtr->col, multiplier);
00506         }
00507 
00508           /* second, update from strict upper triangular portion of row */
00509           while (len--)
00510         {
00511           sr.col = *cvalPtr++;
00512           sr.val = *avalPtr++;
00513           srPtr = SortedList_dhFind (slist, &sr);
00514           CHECK_V_ERROR;
00515           if (srPtr != NULL)
00516             {
00517               srPtr->val -= (multiplier * sr.val);
00518             }
00519         }
00520         }
00521 
00522       else
00523         {
00524           if (debug)
00525         {
00526           fprintf (logFile,
00527                "ILU_pilu   NO UPDATE from row: %i; srPtr->val = 0.0\n",
00528                1 + srPtr->col);
00529         }
00530         }
00531 
00532     }
00533     }
00534 END_FUNC_DH}
 All Classes Files Functions Variables Enumerations Friends
Generated on Wed Apr 13 10:05:23 2011 for IFPACK by  doxygen 1.6.3