ilu_mpi_pilu.c

Go to the documentation of this file.
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 Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:35 2011 for Ifpack Package Browser (Single Doxygen Collection) by  doxygen 1.6.3