ilu_mpi_bj.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 /* 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 
00041 int symbolic_row_private (int localRow, int beg_row, int end_row,
00042         int *list, int *marker, int *tmpFill,
00043         int len, int *CVAL, double *AVAL,
00044         int *o2n_col, Euclid_dh ctx);
00045 
00046 static int numeric_row_private (int localRow, int beg_row, int end_row,
00047         int len, int *CVAL, double *AVAL,
00048         REAL_DH * work, int *o2n_col, Euclid_dh ctx);
00049 
00050 
00051 /* all non-local column indices are discarded in symbolic_row_private() */
00052 #undef __FUNC__
00053 #define __FUNC__ "iluk_mpi_bj"
00054 void
00055 iluk_mpi_bj (Euclid_dh ctx)
00056 {
00057   START_FUNC_DH int *rp, *cval, *diag;
00058   int *CVAL;
00059   int i, j, len, count, col, idx = 0;
00060   int *list, *marker, *fill, *tmpFill;
00061   int temp, m, from = ctx->from, to = ctx->to;
00062   int *n2o_row, *o2n_col;
00063   int first_row, last_row;
00064   double *AVAL;
00065   REAL_DH *work, *aval;
00066   Factor_dh F = ctx->F;
00067   SubdomainGraph_dh sg = ctx->sg;
00068 
00069   if (ctx->F == NULL)
00070     {
00071       SET_V_ERROR ("ctx->F is NULL");
00072     }
00073   if (ctx->F->rp == NULL)
00074     {
00075       SET_V_ERROR ("ctx->F->rp is NULL");
00076     }
00077 
00078 /*  printf_dh("====================== starting iluk_mpi_bj; level= %i\n\n", ctx->level);
00079 */
00080 
00081   m = F->m;
00082   rp = F->rp;
00083   cval = F->cval;
00084   fill = F->fill;
00085   diag = F->diag;
00086   aval = F->aval;
00087   work = ctx->work;
00088 
00089   n2o_row = sg->n2o_row;
00090   o2n_col = sg->o2n_col;
00091 
00092   /* allocate and initialize working space */
00093   list = (int *) MALLOC_DH ((m + 1) * sizeof (int));
00094   CHECK_V_ERROR;
00095   marker = (int *) MALLOC_DH (m * sizeof (int));
00096   CHECK_V_ERROR;
00097   tmpFill = (int *) MALLOC_DH (m * sizeof (int));
00098   CHECK_V_ERROR;
00099   for (i = 0; i < m; ++i)
00100     {
00101       marker[i] = -1;
00102       work[i] = 0.0;
00103     }
00104 
00105   /*---------- main loop ----------*/
00106 
00107   /* global numbers of first and last locally owned rows,
00108      with respect to A 
00109    */
00110   first_row = sg->beg_row[myid_dh];
00111   last_row = first_row + sg->row_count[myid_dh];
00112   for (i = from; i < to; ++i)
00113     {
00114 
00115       int row = n2o_row[i]; /* local row number */
00116       int globalRow = row + first_row;  /* global row number */
00117 
00118       EuclidGetRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
00119       CHECK_V_ERROR;
00120 
00121       /* compute scaling value for row(i) */
00122       if (ctx->isScaled)
00123   {
00124     compute_scaling_private (i, len, AVAL, ctx);
00125     CHECK_V_ERROR;
00126   }
00127 
00128       /* Compute symbolic factor for row(i);
00129          this also performs sparsification
00130        */
00131       count = symbolic_row_private (i, first_row, last_row,
00132             list, marker, tmpFill,
00133             len, CVAL, AVAL, o2n_col, ctx);
00134       CHECK_V_ERROR;
00135 
00136       /* Ensure adequate storage; reallocate, if necessary. */
00137       if (idx + count > F->alloc)
00138   {
00139     Factor_dhReallocate (F, idx, count);
00140     CHECK_V_ERROR;
00141     SET_INFO ("REALLOCATED from lu_mpi_bj");
00142     cval = F->cval;
00143     fill = F->fill;
00144     aval = F->aval;
00145   }
00146 
00147       /* Copy factored symbolic row to permanent storage */
00148       col = list[m];
00149       while (count--)
00150   {
00151     cval[idx] = col;
00152     fill[idx] = tmpFill[col];
00153     ++idx;
00154     col = list[col];
00155   }
00156 
00157       /* add row-pointer to start of next row. */
00158       rp[i + 1] = idx;
00159 
00160       /* Insert pointer to diagonal */
00161       temp = rp[i];
00162       while (cval[temp] != i)
00163   ++temp;
00164       diag[i] = temp;
00165 
00166       /* compute numeric factor for current row */
00167       numeric_row_private (i, first_row, last_row,
00168          len, CVAL, AVAL, work, o2n_col, ctx);
00169       CHECK_V_ERROR EuclidRestoreRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
00170       CHECK_V_ERROR;
00171 
00172       /* Copy factored numeric row to permanent storage,
00173          and re-zero work vector
00174        */
00175       for (j = rp[i]; j < rp[i + 1]; ++j)
00176   {
00177     col = cval[j];
00178     aval[j] = work[col];
00179     work[col] = 0.0;
00180   }
00181 
00182       /* check for zero diagonal */
00183       if (!aval[diag[i]])
00184   {
00185     sprintf (msgBuf_dh, "zero diagonal in local row %i", i + 1);
00186     SET_V_ERROR (msgBuf_dh);
00187   }
00188     }
00189 
00190   FREE_DH (list);
00191   CHECK_V_ERROR;
00192   FREE_DH (tmpFill);
00193   CHECK_V_ERROR;
00194   FREE_DH (marker);
00195   CHECK_V_ERROR;
00196 
00197 END_FUNC_DH}
00198 
00199 
00200 
00201 /* Computes ILU(K) factor of a single row; returns fill 
00202    count for the row.  Explicitly inserts diag if not already 
00203    present.  On return, all column indices are local 
00204    (i.e, referenced to 0).
00205 */
00206 #undef __FUNC__
00207 #define __FUNC__ "symbolic_row_private"
00208 int
00209 symbolic_row_private (int localRow, int beg_row, int end_row,
00210           int *list, int *marker, int *tmpFill,
00211           int len, int *CVAL, double *AVAL,
00212           int *o2n_col, Euclid_dh ctx)
00213 {
00214   START_FUNC_DH int level = ctx->level, m = ctx->F->m;
00215   int *cval = ctx->F->cval, *diag = ctx->F->diag, *rp = ctx->F->rp;
00216   int *fill = ctx->F->fill;
00217   int count = 0;
00218   int j, node, tmp, col, head;
00219   int fill1, fill2;
00220   float val;
00221   double thresh = ctx->sparseTolA;
00222   REAL_DH scale;
00223 
00224   scale = ctx->scale[localRow];
00225   ctx->stats[NZA_STATS] += (double) len;
00226 
00227   /* Insert col indices in linked list, and values in work vector.
00228    * List[m] points to the first (smallest) col in the linked list.
00229    * Column values are adjusted from global to local numbering.
00230    */
00231   list[m] = m;
00232   for (j = 0; j < len; ++j)
00233     {
00234       tmp = m;
00235       col = *CVAL++;
00236       val = *AVAL++;
00237 
00238       /* throw out nonlocal columns */
00239       if (col >= beg_row && col < end_row)
00240   {
00241     col -= beg_row; /* adjust column to local zero-based */
00242     col = o2n_col[col]; /* permute column */
00243     if (fabs (scale * val) > thresh || col == localRow)
00244       {     /* sparsification */
00245         ++count;
00246         while (col > list[tmp])
00247     tmp = list[tmp];
00248         list[col] = list[tmp];
00249         list[tmp] = col;
00250         tmpFill[col] = 0;
00251         marker[col] = localRow;
00252       }
00253   }
00254     }
00255 
00256   /* insert diag if not already present */
00257   if (marker[localRow] != localRow)
00258     {
00259 /*     ctx->symbolicZeroDiags += 1; */
00260       tmp = m;
00261       while (localRow > list[tmp])
00262   tmp = list[tmp];
00263       list[localRow] = list[tmp];
00264       list[tmp] = localRow;
00265       tmpFill[localRow] = 0;
00266       marker[localRow] = localRow;
00267       ++count;
00268     }
00269   ctx->stats[NZA_USED_STATS] += (double) count;
00270 
00271   /* update row from previously factored rows */
00272   head = m;
00273   if (level > 0)
00274     {
00275       while (list[head] < localRow)
00276   {
00277     node = list[head];
00278     fill1 = tmpFill[node];
00279 
00280     if (fill1 < level)
00281       {
00282         for (j = diag[node] + 1; j < rp[node + 1]; ++j)
00283     {
00284       col = cval[j];
00285       fill2 = fill1 + fill[j] + 1;
00286 
00287       if (fill2 <= level)
00288         {
00289           /* if newly discovered fill entry, mark it as discovered;
00290            * if entry has level <= K, add it to the linked-list.
00291            */
00292           if (marker[col] < localRow)
00293       {
00294         tmp = head;
00295         marker[col] = localRow;
00296         tmpFill[col] = fill2;
00297         while (col > list[tmp])
00298           tmp = list[tmp];
00299         list[col] = list[tmp];
00300         list[tmp] = col;
00301         ++count;  /* increment fill count */
00302       }
00303 
00304           /* if previously-discovered fill, update the entry's level. */
00305           else
00306       {
00307         tmpFill[col] =
00308           (fill2 < tmpFill[col]) ? fill2 : tmpFill[col];
00309       }
00310         }
00311     }
00312       }
00313     head = list[head];  /* advance to next item in linked list */
00314   }
00315     }
00316 END_FUNC_VAL (count)}
00317 
00318 
00319 #undef __FUNC__
00320 #define __FUNC__ "numeric_row_private"
00321 int
00322 numeric_row_private (int localRow, int beg_row, int end_row,
00323          int len, int *CVAL, double *AVAL,
00324          REAL_DH * work, int *o2n_col, Euclid_dh ctx)
00325 {
00326   START_FUNC_DH double pc, pv, multiplier;
00327   int j, k, col, row;
00328   int *rp = ctx->F->rp, *cval = ctx->F->cval;
00329   int *diag = ctx->F->diag;
00330   double val;
00331   REAL_DH *aval = ctx->F->aval, scale;
00332 
00333   scale = ctx->scale[localRow];
00334 
00335   /* zero work vector */
00336   /* note: indices in col[] are already permuted, and are 
00337      local (zero-based)
00338    */
00339   for (j = rp[localRow]; j < rp[localRow + 1]; ++j)
00340     {
00341       col = cval[j];
00342       work[col] = 0.0;
00343     }
00344 
00345   /* init work vector with values from A */
00346   /* (note: some values may be na due to sparsification; this is O.K.) */
00347   for (j = 0; j < len; ++j)
00348     {
00349       col = *CVAL++;
00350       val = *AVAL++;
00351 
00352       if (col >= beg_row && col < end_row)
00353   {
00354     col -= beg_row; /* adjust column to local zero-based */
00355     col = o2n_col[col]; /* we permute the indices from A */
00356     work[col] = val * scale;
00357   }
00358     }
00359 
00360   for (j = rp[localRow]; j < diag[localRow]; ++j)
00361     {
00362       row = cval[j];
00363       pc = work[row];
00364 
00365       if (pc != 0.0)
00366   {
00367     pv = aval[diag[row]];
00368     multiplier = pc / pv;
00369     work[row] = multiplier;
00370 
00371     for (k = diag[row] + 1; k < rp[row + 1]; ++k)
00372       {
00373         col = cval[k];
00374         work[col] -= (multiplier * aval[k]);
00375       }
00376   }
00377     }
00378 
00379   /* check for zero or too small of a pivot */
00380 #if 0
00381   if (fabs (work[i]) <= pivotTol)
00382     {
00383       /* yuck! assume row scaling, and just stick in a value */
00384       aval[diag[i]] = pivotFix;
00385     }
00386 #endif
00387 
00388 END_FUNC_VAL (0)}
 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