ilu_seq.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 /* to do: re-integrate fix-smalll-pivots */
00031 
00032 #include "ilu_dh.h"
00033 #include "Mem_dh.h"
00034 #include "Parser_dh.h"
00035 #include "Euclid_dh.h"
00036 #include "getRow_dh.h"
00037 #include "Factor_dh.h"
00038 #include "SubdomainGraph_dh.h"
00039 
00040 static bool check_constraint_private (Euclid_dh ctx, int b, int j);
00041 
00042 static int symbolic_row_private (int localRow,
00043                  int *list, int *marker, int *tmpFill,
00044                  int len, int *CVAL, double *AVAL,
00045                  int *o2n_col, Euclid_dh ctx, bool debug);
00046 
00047 static int numeric_row_private (int localRow,
00048                 int len, int *CVAL, double *AVAL,
00049                 REAL_DH * work, int *o2n_col, Euclid_dh ctx,
00050                 bool debug);
00051 
00052 
00053 #undef __FUNC__
00054 #define __FUNC__ "compute_scaling_private"
00055 void
00056 compute_scaling_private (int row, int len, double *AVAL, Euclid_dh ctx)
00057 {
00058   START_FUNC_DH double tmp = 0.0;
00059   int j;
00060 
00061   for (j = 0; j < len; ++j)
00062     tmp = MAX (tmp, fabs (AVAL[j]));
00063   if (tmp)
00064     {
00065       ctx->scale[row] = 1.0 / tmp;
00066     }
00067 END_FUNC_DH}
00068 
00069 #if 0
00070 
00071 /* not used ? */
00072 #undef __FUNC__
00073 #define __FUNC__ "fixPivot_private"
00074 double
00075 fixPivot_private (int row, int len, float *vals)
00076 {
00077   START_FUNC_DH int i;
00078   float max = 0.0;
00079   bool debug = false;
00080 
00081   for (i = 0; i < len; ++i)
00082     {
00083       float tmp = fabs (vals[i]);
00084       max = MAX (max, tmp);
00085     }
00086 END_FUNC_VAL (max * ctxPrivate->pivotFix)}
00087 
00088 #endif
00089 
00090 
00091 
00092 
00093 #undef __FUNC__
00094 #define __FUNC__ "iluk_seq"
00095 void
00096 iluk_seq (Euclid_dh ctx)
00097 {
00098   START_FUNC_DH int *rp, *cval, *diag;
00099   int *CVAL;
00100   int i, j, len, count, col, idx = 0;
00101   int *list, *marker, *fill, *tmpFill;
00102   int temp, m, from = ctx->from, to = ctx->to;
00103   int *n2o_row, *o2n_col, beg_row, beg_rowP;
00104   double *AVAL;
00105   REAL_DH *work, *aval;
00106   Factor_dh F = ctx->F;
00107   SubdomainGraph_dh sg = ctx->sg;
00108   bool debug = false;
00109 
00110   if (logFile != NULL && Parser_dhHasSwitch (parser_dh, "-debug_ilu"))
00111     debug = true;
00112 
00113   m = F->m;
00114   rp = F->rp;
00115   cval = F->cval;
00116   fill = F->fill;
00117   diag = F->diag;
00118   aval = F->aval;
00119   work = ctx->work;
00120   count = rp[from];
00121 
00122   if (sg == NULL)
00123     {
00124       SET_V_ERROR ("subdomain graph is NULL");
00125     }
00126 
00127   n2o_row = ctx->sg->n2o_row;
00128   o2n_col = ctx->sg->o2n_col;
00129   beg_row = ctx->sg->beg_row[myid_dh];
00130   beg_rowP = ctx->sg->beg_rowP[myid_dh];
00131 
00132   /* allocate and initialize working space */
00133   list = (int *) MALLOC_DH ((m + 1) * sizeof (int));
00134   CHECK_V_ERROR;
00135   marker = (int *) MALLOC_DH (m * sizeof (int));
00136   CHECK_V_ERROR;
00137   tmpFill = (int *) MALLOC_DH (m * sizeof (int));
00138   CHECK_V_ERROR;
00139   for (i = 0; i < m; ++i)
00140     marker[i] = -1;
00141 
00142   /* working space for values */
00143   for (i = 0; i < m; ++i)
00144     work[i] = 0.0;
00145 
00146 /*    printf_dh("====================== starting iluk_seq; level= %i\n\n", ctx->level); 
00147 */
00148 
00149 
00150   /*---------- main loop ----------*/
00151 
00152   for (i = from; i < to; ++i)
00153     {
00154       int row = n2o_row[i]; /* local row number */
00155       int globalRow = row + beg_row;    /* global row number */
00156 
00157 /*fprintf(logFile, "--------------------------------- localRow= %i\n", 1+i);
00158 */
00159 
00160       if (debug)
00161     {
00162       fprintf (logFile,
00163            "ILU_seq ================================= starting local row: %i, (global= %i) level= %i\n",
00164            i + 1, i + 1 + sg->beg_rowP[myid_dh], ctx->level);
00165     }
00166 
00167       EuclidGetRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
00168       CHECK_V_ERROR;
00169 
00170       /* compute scaling value for row(i) */
00171       if (ctx->isScaled)
00172     {
00173       compute_scaling_private (i, len, AVAL, ctx);
00174       CHECK_V_ERROR;
00175     }
00176 
00177       /* Compute symbolic factor for row(i);
00178          this also performs sparsification
00179        */
00180       count = symbolic_row_private (i, list, marker, tmpFill,
00181                     len, CVAL, AVAL, o2n_col, ctx, debug);
00182       CHECK_V_ERROR;
00183 
00184       /* Ensure adequate storage; reallocate, if necessary. */
00185       if (idx + count > F->alloc)
00186     {
00187       Factor_dhReallocate (F, idx, count);
00188       CHECK_V_ERROR;
00189       SET_INFO ("REALLOCATED from ilu_seq");
00190       cval = F->cval;
00191       fill = F->fill;
00192       aval = F->aval;
00193     }
00194 
00195       /* Copy factored symbolic row to permanent storage */
00196       col = list[m];
00197       while (count--)
00198     {
00199       cval[idx] = col;
00200       fill[idx] = tmpFill[col];
00201       ++idx;
00202 /*fprintf(logFile, "  col= %i\n", 1+col);
00203 */
00204       col = list[col];
00205     }
00206 
00207       /* add row-pointer to start of next row. */
00208       rp[i + 1] = idx;
00209 
00210       /* Insert pointer to diagonal */
00211       temp = rp[i];
00212       while (cval[temp] != i)
00213     ++temp;
00214       diag[i] = temp;
00215 
00216 /*fprintf(logFile, "  diag[i]= %i\n", diag);
00217 */
00218 
00219       /* compute numeric factor for current row */
00220       numeric_row_private (i, len, CVAL, AVAL, work, o2n_col, ctx, debug);
00221       CHECK_V_ERROR EuclidRestoreRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
00222       CHECK_V_ERROR;
00223 
00224       /* Copy factored numeric row to permanent storage,
00225          and re-zero work vector
00226        */
00227       if (debug)
00228     {
00229       fprintf (logFile, "ILU_seq:  ");
00230       for (j = rp[i]; j < rp[i + 1]; ++j)
00231         {
00232           col = cval[j];
00233           aval[j] = work[col];
00234           work[col] = 0.0;
00235           fprintf (logFile, "%i,%i,%g ; ", 1 + cval[j], fill[j], aval[j]);
00236           fflush (logFile);
00237         }
00238       fprintf (logFile, "\n");
00239     }
00240       else
00241     {
00242       for (j = rp[i]; j < rp[i + 1]; ++j)
00243         {
00244           col = cval[j];
00245           aval[j] = work[col];
00246           work[col] = 0.0;
00247         }
00248     }
00249 
00250       /* check for zero diagonal */
00251       if (!aval[diag[i]])
00252     {
00253       sprintf (msgBuf_dh, "zero diagonal in local row %i", i + 1);
00254       SET_V_ERROR (msgBuf_dh);
00255     }
00256     }
00257 
00258   FREE_DH (list);
00259   CHECK_V_ERROR;
00260   FREE_DH (tmpFill);
00261   CHECK_V_ERROR;
00262   FREE_DH (marker);
00263   CHECK_V_ERROR;
00264 
00265   /* adjust column indices back to global */
00266   if (beg_rowP)
00267     {
00268       int start = rp[from];
00269       int stop = rp[to];
00270       for (i = start; i < stop; ++i)
00271     cval[i] += beg_rowP;
00272     }
00273 
00274   /* for debugging: this is so the Print methods will work, even if
00275      F hasn't been fully factored
00276    */
00277   for (i = to + 1; i < m; ++i)
00278     rp[i] = 0;
00279 
00280 END_FUNC_DH}
00281 
00282 
00283 #undef __FUNC__
00284 #define __FUNC__ "iluk_seq_block"
00285 void
00286 iluk_seq_block (Euclid_dh ctx)
00287 {
00288   START_FUNC_DH int *rp, *cval, *diag;
00289   int *CVAL;
00290   int h, i, j, len, count, col, idx = 0;
00291   int *list, *marker, *fill, *tmpFill;
00292   int temp, m;
00293   int *n2o_row, *o2n_col, *beg_rowP, *n2o_sub, blocks;
00294   int *row_count, *dummy = NULL, dummy2[1];
00295   double *AVAL;
00296   REAL_DH *work, *aval;
00297   Factor_dh F = ctx->F;
00298   SubdomainGraph_dh sg = ctx->sg;
00299   bool bj = false, constrained = false;
00300   int discard = 0;
00301   int gr = -1;          /* globalRow */
00302   bool debug = false;
00303 
00304   if (logFile != NULL && Parser_dhHasSwitch (parser_dh, "-debug_ilu"))
00305     debug = true;
00306 
00307 /*fprintf(stderr, "====================== starting iluk_seq_block; level= %i\n\n", ctx->level); 
00308 */
00309 
00310   if (!strcmp (ctx->algo_par, "bj"))
00311     bj = true;
00312   constrained = !Parser_dhHasSwitch (parser_dh, "-unconstrained");
00313 
00314   m = F->m;
00315   rp = F->rp;
00316   cval = F->cval;
00317   fill = F->fill;
00318   diag = F->diag;
00319   aval = F->aval;
00320   work = ctx->work;
00321 
00322   if (sg != NULL)
00323     {
00324       n2o_row = sg->n2o_row;
00325       o2n_col = sg->o2n_col;
00326       row_count = sg->row_count;
00327       /* beg_row   = sg->beg_row ; */
00328       beg_rowP = sg->beg_rowP;
00329       n2o_sub = sg->n2o_sub;
00330       blocks = sg->blocks;
00331     }
00332 
00333   else
00334     {
00335       dummy = (int *) MALLOC_DH (m * sizeof (int));
00336       CHECK_V_ERROR;
00337       for (i = 0; i < m; ++i)
00338     dummy[i] = i;
00339       n2o_row = dummy;
00340       o2n_col = dummy;
00341       dummy2[0] = m;
00342       row_count = dummy2;
00343       /* beg_row   = 0; */
00344       beg_rowP = dummy;
00345       n2o_sub = dummy;
00346       blocks = 1;
00347     }
00348 
00349   /* allocate and initialize working space */
00350   list = (int *) MALLOC_DH ((m + 1) * sizeof (int));
00351   CHECK_V_ERROR;
00352   marker = (int *) MALLOC_DH (m * sizeof (int));
00353   CHECK_V_ERROR;
00354   tmpFill = (int *) MALLOC_DH (m * sizeof (int));
00355   CHECK_V_ERROR;
00356   for (i = 0; i < m; ++i)
00357     marker[i] = -1;
00358 
00359   /* working space for values */
00360   for (i = 0; i < m; ++i)
00361     work[i] = 0.0;
00362 
00363   /*---------- main loop ----------*/
00364 
00365   for (h = 0; h < blocks; ++h)
00366     {
00367       /* 1st and last row in current block, with respect to A */
00368       int curBlock = n2o_sub[h];
00369       int first_row = beg_rowP[curBlock];
00370       int end_row = first_row + row_count[curBlock];
00371 
00372       if (debug)
00373     {
00374       fprintf (logFile, "\n\nILU_seq BLOCK: %i @@@@@@@@@@@@@@@ \n",
00375            curBlock);
00376     }
00377 
00378       for (i = first_row; i < end_row; ++i)
00379     {
00380       int row = n2o_row[i];
00381       ++gr;
00382 
00383       if (debug)
00384         {
00385           fprintf (logFile,
00386                "ILU_seq  global: %i  local: %i =================================\n",
00387                1 + gr, 1 + i - first_row);
00388         }
00389 
00390 /*prinft("first_row= %i  end_row= %i\n", first_row, end_row);
00391 */
00392 
00393       EuclidGetRow (ctx->A, row, &len, &CVAL, &AVAL);
00394       CHECK_V_ERROR;
00395 
00396       /* compute scaling value for row(i) */
00397       if (ctx->isScaled)
00398         {
00399           compute_scaling_private (i, len, AVAL, ctx);
00400           CHECK_V_ERROR;
00401         }
00402 
00403       /* Compute symbolic factor for row(i);
00404          this also performs sparsification
00405        */
00406       count = symbolic_row_private (i, list, marker, tmpFill,
00407                     len, CVAL, AVAL, o2n_col, ctx, debug);
00408       CHECK_V_ERROR;
00409 
00410       /* Ensure adequate storage; reallocate, if necessary. */
00411       if (idx + count > F->alloc)
00412         {
00413           Factor_dhReallocate (F, idx, count);
00414           CHECK_V_ERROR;
00415           SET_INFO ("REALLOCATED from ilu_seq");
00416           cval = F->cval;
00417           fill = F->fill;
00418           aval = F->aval;
00419         }
00420 
00421       /* Copy factored symbolic row to permanent storage */
00422       col = list[m];
00423       while (count--)
00424         {
00425 
00426           /* constrained pilu */
00427           if (constrained && !bj)
00428         {
00429           if (col >= first_row && col < end_row)
00430             {
00431               cval[idx] = col;
00432               fill[idx] = tmpFill[col];
00433               ++idx;
00434             }
00435           else
00436             {
00437               if (check_constraint_private (ctx, curBlock, col))
00438             {
00439               cval[idx] = col;
00440               fill[idx] = tmpFill[col];
00441               ++idx;
00442             }
00443               else
00444             {
00445               ++discard;
00446             }
00447             }
00448           col = list[col];
00449         }
00450 
00451           /* block jacobi case */
00452           else if (bj)
00453         {
00454           if (col >= first_row && col < end_row)
00455             {
00456               cval[idx] = col;
00457               fill[idx] = tmpFill[col];
00458               ++idx;
00459             }
00460           else
00461             {
00462               ++discard;
00463             }
00464           col = list[col];
00465         }
00466 
00467           /* general case */
00468           else
00469         {
00470           cval[idx] = col;
00471           fill[idx] = tmpFill[col];
00472           ++idx;
00473           col = list[col];
00474         }
00475         }
00476 
00477       /* add row-pointer to start of next row. */
00478       rp[i + 1] = idx;
00479 
00480       /* Insert pointer to diagonal */
00481       temp = rp[i];
00482       while (cval[temp] != i)
00483         ++temp;
00484       diag[i] = temp;
00485 
00486       /* compute numeric factor for current row */
00487       numeric_row_private (i, len, CVAL, AVAL, work, o2n_col, ctx, debug);
00488       CHECK_V_ERROR EuclidRestoreRow (ctx->A, row, &len, &CVAL, &AVAL);
00489       CHECK_V_ERROR;
00490 
00491       /* Copy factored numeric row to permanent storage,
00492          and re-zero work vector
00493        */
00494       if (debug)
00495         {
00496           fprintf (logFile, "ILU_seq: ");
00497           for (j = rp[i]; j < rp[i + 1]; ++j)
00498         {
00499           col = cval[j];
00500           aval[j] = work[col];
00501           work[col] = 0.0;
00502           fprintf (logFile, "%i,%i,%g ; ", 1 + cval[j], fill[j],
00503                aval[j]);
00504         }
00505           fprintf (logFile, "\n");
00506         }
00507 
00508       /* normal operation */
00509       else
00510         {
00511           for (j = rp[i]; j < rp[i + 1]; ++j)
00512         {
00513           col = cval[j];
00514           aval[j] = work[col];
00515           work[col] = 0.0;
00516         }
00517         }
00518 
00519       /* check for zero diagonal */
00520       if (!aval[diag[i]])
00521         {
00522           sprintf (msgBuf_dh, "zero diagonal in local row %i", i + 1);
00523           SET_V_ERROR (msgBuf_dh);
00524         }
00525     }
00526     }
00527 
00528 /*  printf("bj= %i  constrained= %i  discarded= %i\n", bj, constrained, discard); */
00529 
00530   if (dummy != NULL)
00531     {
00532       FREE_DH (dummy);
00533       CHECK_V_ERROR;
00534     }
00535   FREE_DH (list);
00536   CHECK_V_ERROR;
00537   FREE_DH (tmpFill);
00538   CHECK_V_ERROR;
00539   FREE_DH (marker);
00540   CHECK_V_ERROR;
00541 
00542 END_FUNC_DH}
00543 
00544 
00545 
00546 /* Computes ILU(K) factor of a single row; returns fill 
00547    count for the row.  Explicitly inserts diag if not already 
00548    present.  On return, all column indices are local 
00549    (i.e, referenced to 0).
00550 */
00551 #undef __FUNC__
00552 #define __FUNC__ "symbolic_row_private"
00553 int
00554 symbolic_row_private (int localRow,
00555               int *list, int *marker, int *tmpFill,
00556               int len, int *CVAL, double *AVAL,
00557               int *o2n_col, Euclid_dh ctx, bool debug)
00558 {
00559   START_FUNC_DH int level = ctx->level, m = ctx->F->m;
00560   int *cval = ctx->F->cval, *diag = ctx->F->diag, *rp = ctx->F->rp;
00561   int *fill = ctx->F->fill;
00562   int count = 0;
00563   int j, node, tmp, col, head;
00564   int fill1, fill2, beg_row;
00565   double val;
00566   double thresh = ctx->sparseTolA;
00567   REAL_DH scale;
00568 
00569   scale = ctx->scale[localRow];
00570   ctx->stats[NZA_STATS] += (double) len;
00571   beg_row = ctx->sg->beg_row[myid_dh];
00572 
00573   /* Insert col indices in linked list, and values in work vector.
00574    * List[m] points to the first (smallest) col in the linked list.
00575    * Column values are adjusted from global to local numbering.
00576    */
00577   list[m] = m;
00578   for (j = 0; j < len; ++j)
00579     {
00580       tmp = m;
00581       col = *CVAL++;
00582       col -= beg_row;       /* adjust to zero based */
00583       col = o2n_col[col];   /* permute the column */
00584       val = *AVAL++;
00585       val *= scale;     /* scale the value */
00586 
00587       if (fabs (val) > thresh || col == localRow)
00588     {           /* sparsification */
00589       ++count;
00590       while (col > list[tmp])
00591         tmp = list[tmp];
00592       list[col] = list[tmp];
00593       list[tmp] = col;
00594       tmpFill[col] = 0;
00595       marker[col] = localRow;
00596     }
00597     }
00598 
00599   /* insert diag if not already present */
00600   if (marker[localRow] != localRow)
00601     {
00602       tmp = m;
00603       while (localRow > list[tmp])
00604     tmp = list[tmp];
00605       list[localRow] = list[tmp];
00606       list[tmp] = localRow;
00607       tmpFill[localRow] = 0;
00608       marker[localRow] = localRow;
00609       ++count;
00610     }
00611   ctx->stats[NZA_USED_STATS] += (double) count;
00612 
00613   /* update row from previously factored rows */
00614   head = m;
00615   if (level > 0)
00616     {
00617       while (list[head] < localRow)
00618     {
00619       node = list[head];
00620       fill1 = tmpFill[node];
00621 
00622       if (debug)
00623         {
00624           fprintf (logFile, "ILU_seq   sf updating from row: %i\n",
00625                1 + node);
00626         }
00627 
00628       if (fill1 < level)
00629         {
00630           for (j = diag[node] + 1; j < rp[node + 1]; ++j)
00631         {
00632           col = cval[j];
00633           fill2 = fill1 + fill[j] + 1;
00634 
00635           if (fill2 <= level)
00636             {
00637               /* if newly discovered fill entry, mark it as discovered;
00638                * if entry has level <= K, add it to the linked-list.
00639                */
00640               if (marker[col] < localRow)
00641             {
00642               tmp = head;
00643               marker[col] = localRow;
00644               tmpFill[col] = fill2;
00645               while (col > list[tmp])
00646                 tmp = list[tmp];
00647               list[col] = list[tmp];
00648               list[tmp] = col;
00649               ++count;  /* increment fill count */
00650             }
00651 
00652               /* if previously-discovered fill, update the entry's level. */
00653               else
00654             {
00655               tmpFill[col] =
00656                 (fill2 < tmpFill[col]) ? fill2 : tmpFill[col];
00657             }
00658             }
00659         }
00660         }           /* fill1 < level  */
00661       head = list[head];    /* advance to next item in linked list */
00662     }
00663     }
00664 END_FUNC_VAL (count)}
00665 
00666 
00667 #undef __FUNC__
00668 #define __FUNC__ "numeric_row_private"
00669 int
00670 numeric_row_private (int localRow,
00671              int len, int *CVAL, double *AVAL,
00672              REAL_DH * work, int *o2n_col, Euclid_dh ctx, bool debug)
00673 {
00674   START_FUNC_DH double pc, pv, multiplier;
00675   int j, k, col, row;
00676   int *rp = ctx->F->rp, *cval = ctx->F->cval;
00677   int *diag = ctx->F->diag;
00678   int beg_row;
00679   double val;
00680   REAL_DH *aval = ctx->F->aval, scale;
00681 
00682   scale = ctx->scale[localRow];
00683   beg_row = ctx->sg->beg_row[myid_dh];
00684 
00685   /* zero work vector */
00686   /* note: indices in col[] are already permuted. */
00687   for (j = rp[localRow]; j < rp[localRow + 1]; ++j)
00688     {
00689       col = cval[j];
00690       work[col] = 0.0;
00691     }
00692 
00693   /* init work vector with values from A */
00694   /* (note: some values may be na due to sparsification; this is O.K.) */
00695   for (j = 0; j < len; ++j)
00696     {
00697       col = *CVAL++;
00698       col -= beg_row;
00699       val = *AVAL++;
00700       col = o2n_col[col];   /* note: we permute the indices from A */
00701       work[col] = val * scale;
00702     }
00703 
00704 
00705 
00706 /*fprintf(stderr, "local row= %i\n", 1+localRow);
00707 */
00708 
00709 
00710   for (j = rp[localRow]; j < diag[localRow]; ++j)
00711     {
00712       row = cval[j];        /* previously factored row */
00713       pc = work[row];
00714 
00715 
00716       pv = aval[diag[row]]; /* diagonal of previously factored row */
00717 
00718 /*
00719 if (pc == 0.0 || pv == 0.0) {
00720 fprintf(stderr, "pv= %g; pc= %g\n", pv, pc);
00721 }
00722 */
00723 
00724       if (pc != 0.0 && pv != 0.0)
00725     {
00726       multiplier = pc / pv;
00727       work[row] = multiplier;
00728 
00729       if (debug)
00730         {
00731           fprintf (logFile,
00732                "ILU_seq   nf updating from row: %i; multiplier= %g\n",
00733                1 + row, multiplier);
00734         }
00735 
00736       for (k = diag[row] + 1; k < rp[row + 1]; ++k)
00737         {
00738           col = cval[k];
00739           work[col] -= (multiplier * aval[k]);
00740         }
00741     }
00742       else
00743     {
00744       if (debug)
00745         {
00746           fprintf (logFile,
00747                "ILU_seq   nf NO UPDATE from row %i; pc = %g; pv = %g\n",
00748                1 + row, pc, pv);
00749         }
00750     }
00751     }
00752 
00753   /* check for zero or too small of a pivot */
00754 #if 0
00755   if (fabs (work[i]) <= pivotTol)
00756     {
00757       /* yuck! assume row scaling, and just stick in a value */
00758       aval[diag[i]] = pivotFix;
00759     }
00760 #endif
00761 
00762 END_FUNC_VAL (0)}
00763 
00764 
00765 /*-----------------------------------------------------------------------*
00766  * ILUT starts here
00767  *-----------------------------------------------------------------------*/
00768 int ilut_row_private (int localRow, int *list, int *o2n_col, int *marker,
00769               int len, int *CVAL, double *AVAL,
00770               REAL_DH * work, Euclid_dh ctx, bool debug);
00771 
00772 #undef __FUNC__
00773 #define __FUNC__ "ilut_seq"
00774 void
00775 ilut_seq (Euclid_dh ctx)
00776 {
00777   START_FUNC_DH int *rp, *cval, *diag, *CVAL;
00778   int i, len, count, col, idx = 0;
00779   int *list, *marker;
00780   int temp, m, from, to;
00781   int *n2o_row, *o2n_col, beg_row, beg_rowP;
00782   double *AVAL, droptol;
00783   REAL_DH *work, *aval, val;
00784   Factor_dh F = ctx->F;
00785   SubdomainGraph_dh sg = ctx->sg;
00786   bool debug = false;
00787 
00788   if (logFile != NULL && Parser_dhHasSwitch (parser_dh, "-debug_ilu"))
00789     debug = true;
00790 
00791   m = F->m;
00792   rp = F->rp;
00793   cval = F->cval;
00794   diag = F->diag;
00795   aval = F->aval;
00796   work = ctx->work;
00797   from = ctx->from;
00798   to = ctx->to;
00799   count = rp[from];
00800   droptol = ctx->droptol;
00801 
00802   if (sg == NULL)
00803     {
00804       SET_V_ERROR ("subdomain graph is NULL");
00805     }
00806 
00807   n2o_row = ctx->sg->n2o_row;
00808   o2n_col = ctx->sg->o2n_col;
00809   beg_row = ctx->sg->beg_row[myid_dh];
00810   beg_rowP = ctx->sg->beg_rowP[myid_dh];
00811 
00812 
00813   /* allocate and initialize working space */
00814   list = (int *) MALLOC_DH ((m + 1) * sizeof (int));
00815   CHECK_V_ERROR;
00816   marker = (int *) MALLOC_DH (m * sizeof (int));
00817   CHECK_V_ERROR;
00818   for (i = 0; i < m; ++i)
00819     marker[i] = -1;
00820   rp[0] = 0;
00821 
00822   /* working space for values */
00823   for (i = 0; i < m; ++i)
00824     work[i] = 0.0;
00825 
00826   /* ----- main loop start ----- */
00827   for (i = from; i < to; ++i)
00828     {
00829       int row = n2o_row[i]; /* local row number */
00830       int globalRow = row + beg_row;    /* global row number */
00831       EuclidGetRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
00832       CHECK_V_ERROR;
00833 
00834       /* compute scaling value for row(i) */
00835       compute_scaling_private (i, len, AVAL, ctx);
00836       CHECK_V_ERROR;
00837 
00838       /* compute factor for row i */
00839       count = ilut_row_private (i, list, o2n_col, marker,
00840                 len, CVAL, AVAL, work, ctx, debug);
00841       CHECK_V_ERROR;
00842 
00843       EuclidRestoreRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
00844       CHECK_V_ERROR;
00845 
00846       /* Ensure adequate storage; reallocate, if necessary. */
00847       if (idx + count > F->alloc)
00848     {
00849       Factor_dhReallocate (F, idx, count);
00850       CHECK_V_ERROR;
00851       SET_INFO ("REALLOCATED from ilu_seq");
00852       cval = F->cval;
00853       aval = F->aval;
00854     }
00855 
00856       /* Copy factored row to permanent storage,
00857          apply 2nd drop test,
00858          and re-zero work vector
00859        */
00860       col = list[m];
00861       while (count--)
00862     {
00863       val = work[col];
00864       if (col == i || fabs (val) > droptol)
00865         {
00866           cval[idx] = col;
00867           aval[idx++] = val;
00868           work[col] = 0.0;
00869         }
00870       col = list[col];
00871     }
00872 
00873       /* add row-pointer to start of next row. */
00874       rp[i + 1] = idx;
00875 
00876       /* Insert pointer to diagonal */
00877       temp = rp[i];
00878       while (cval[temp] != i)
00879     ++temp;
00880       diag[i] = temp;
00881 
00882       /* check for zero diagonal */
00883       if (!aval[diag[i]])
00884     {
00885       sprintf (msgBuf_dh, "zero diagonal in local row %i", i + 1);
00886       SET_V_ERROR (msgBuf_dh);
00887     }
00888     }               /* --------- main loop end --------- */
00889 
00890   /* adjust column indices back to global */
00891   if (beg_rowP)
00892     {
00893       int start = rp[from];
00894       int stop = rp[to];
00895       for (i = start; i < stop; ++i)
00896     cval[i] += beg_rowP;
00897     }
00898 
00899   FREE_DH (list);
00900   FREE_DH (marker);
00901 END_FUNC_DH}
00902 
00903 
00904 #undef __FUNC__
00905 #define __FUNC__ "ilut_row_private"
00906 int
00907 ilut_row_private (int localRow, int *list, int *o2n_col, int *marker,
00908           int len, int *CVAL, double *AVAL,
00909           REAL_DH * work, Euclid_dh ctx, bool debug)
00910 {
00911   START_FUNC_DH Factor_dh F = ctx->F;
00912   int j, col, m = ctx->m, *rp = F->rp, *cval = F->cval;
00913   int tmp, *diag = F->diag;
00914   int head;
00915   int count = 0, beg_row;
00916   double val;
00917   double mult, *aval = F->aval;
00918   double scale, pv, pc;
00919   double droptol = ctx->droptol;
00920   double thresh = ctx->sparseTolA;
00921 
00922   scale = ctx->scale[localRow];
00923   ctx->stats[NZA_STATS] += (double) len;
00924   beg_row = ctx->sg->beg_row[myid_dh];
00925 
00926 
00927   /* Insert col indices in linked list, and values in work vector.
00928    * List[m] points to the first (smallest) col in the linked list.
00929    * Column values are adjusted from global to local numbering.
00930    */
00931   list[m] = m;
00932   for (j = 0; j < len; ++j)
00933     {
00934       tmp = m;
00935       col = *CVAL++;
00936       col -= beg_row;       /* adjust to zero based */
00937       col = o2n_col[col];   /* permute the column */
00938       val = *AVAL++;
00939       val *= scale;     /* scale the value */
00940 
00941       if (fabs (val) > thresh || col == localRow)
00942     {           /* sparsification */
00943       ++count;
00944       while (col > list[tmp])
00945         tmp = list[tmp];
00946       list[col] = list[tmp];
00947       list[tmp] = col;
00948       work[col] = val;
00949       marker[col] = localRow;
00950     }
00951     }
00952 
00953   /* insert diag if not already present */
00954   if (marker[localRow] != localRow)
00955     {
00956       tmp = m;
00957       while (localRow > list[tmp])
00958     tmp = list[tmp];
00959       list[localRow] = list[tmp];
00960       list[tmp] = localRow;
00961       marker[localRow] = localRow;
00962       ++count;
00963     }
00964 
00965   /* update current row from previously factored rows */
00966   head = m;
00967   while (list[head] < localRow)
00968     {
00969       int row = list[head];
00970 
00971       /* get the multiplier, and apply 1st drop tolerance test */
00972       pc = work[row];
00973       if (pc != 0.0)
00974     {
00975       pv = aval[diag[row]]; /* diagonal (pivot) of previously factored row */
00976       mult = pc / pv;
00977 
00978       /* update localRow from previously factored "row" */
00979       if (fabs (mult) > droptol)
00980         {
00981           work[row] = mult;
00982 
00983           for (j = diag[row] + 1; j < rp[row + 1]; ++j)
00984         {
00985           col = cval[j];
00986           work[col] -= (mult * aval[j]);
00987 
00988           /* if col isn't already present in the linked-list, insert it.  */
00989           if (marker[col] < localRow)
00990             {
00991               marker[col] = localRow;   /* mark the column as known fill */
00992               tmp = head;   /* insert in list [this and next 3 lines] */
00993               while (col > list[tmp])
00994             tmp = list[tmp];
00995               list[col] = list[tmp];
00996               list[tmp] = col;
00997               ++count;  /* increment fill count */
00998             }
00999         }
01000         }
01001     }
01002       head = list[head];    /* advance to next item in linked list */
01003     }
01004 
01005 END_FUNC_VAL (count)}
01006 
01007 
01008 #undef __FUNC__
01009 #define __FUNC__ "check_constraint_private"
01010 bool
01011 check_constraint_private (Euclid_dh ctx, int p1, int j)
01012 {
01013   START_FUNC_DH bool retval = false;
01014   int i, p2;
01015   int *nabors, count;
01016   SubdomainGraph_dh sg = ctx->sg;
01017 
01018   if (sg == NULL)
01019     {
01020       SET_ERROR (-1, "ctx->sg == NULL");
01021     }
01022 
01023   p2 = SubdomainGraph_dhFindOwner (ctx->sg, j, true);
01024 
01025 
01026   nabors = sg->adj + sg->ptrs[p1];
01027   count = sg->ptrs[p1 + 1] - sg->ptrs[p1];
01028 
01029 /*
01030 printf("p1= %i, p2= %i;  p1's nabors: ", p1, p2);
01031 for (i=0; i<count; ++i) printf("%i ", nabors[i]);
01032 printf("\n");
01033 */
01034 
01035   for (i = 0; i < count; ++i)
01036     {
01037 /* printf("  @@@ next nabor= %i\n", nabors[i]);
01038 */
01039       if (nabors[i] == p2)
01040     {
01041       retval = true;
01042       break;
01043     }
01044     }
01045 
01046 END_FUNC_VAL (retval)}
 All Classes Files Functions Variables Enumerations Friends
Generated on Wed Apr 13 10:05:23 2011 for IFPACK by  doxygen 1.6.3