ilu_mpi_pilu.c

Go to the documentation of this file.

/*@HEADER
// ***********************************************************************
//
//       Ifpack: Object-Oriented Algebraic Preconditioner Package
//                 Copyright (2002) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
//@HEADER
*/

#include "Euclid_dh.h"
#include "Factor_dh.h"
#include "Mat_dh.h"
#include "ilu_dh.h"
#include "Mem_dh.h"
#include "Parser_dh.h"
#include "Hash_dh.h"
#include "getRow_dh.h"
#include "SortedList_dh.h"
#include "ExternalRows_dh.h"
#include "SubdomainGraph_dh.h"


static void iluk_symbolic_row_private (int localRow, int len, int *CVAL,
               double *AVAL, ExternalRows_dh extRows,
               SortedList_dh sList, Euclid_dh ctx,
               bool debug);

static void iluk_numeric_row_private (int new_row, ExternalRows_dh extRows,
              SortedList_dh slist, Euclid_dh ctx,
              bool debug);

#undef __FUNC__
#define __FUNC__ "iluk_mpi_pilu"
void
iluk_mpi_pilu (Euclid_dh ctx)
{
  START_FUNC_DH int from = ctx->from, to = ctx->to;
  int i, m;
  int *n2o_row;     /* *o2n_col; */
  int *rp, *cval, *diag, *fill;
  int beg_row, beg_rowP, end_rowP;
  SubdomainGraph_dh sg = ctx->sg;
  int *CVAL, len, idx = 0, count;
  double *AVAL;
  REAL_DH *aval;
  Factor_dh F = ctx->F;
  SortedList_dh slist = ctx->slist;
  ExternalRows_dh extRows = ctx->extRows;
  bool bj, noValues, debug = false;

  /* for debugging */
  if (logFile != NULL && Parser_dhHasSwitch (parser_dh, "-debug_ilu"))
    debug = true;
  noValues = Parser_dhHasSwitch (parser_dh, "-noValues");
  bj = ctx->F->blockJacobi;

  m = F->m;
  rp = F->rp;
  cval = F->cval;
  fill = F->fill;
  diag = F->diag;
  aval = F->aval;
  /* work = ctx->work; */

  n2o_row = sg->n2o_row;
  /* o2n_col = sg->o2n_col; */

  if (from != 0)
    idx = rp[from];

  /* global numbers of first and last locally owned rows,
     with respect to A 
   */
  beg_row = sg->beg_row[myid_dh];
  /* end_row  = beg_row + sg->row_count[myid_dh]; */

  /* global number or first locally owned row, after reordering */
  beg_rowP = sg->beg_rowP[myid_dh];
  end_rowP = beg_rowP + sg->row_count[myid_dh];


  /* loop over rows to be factored (i references local rows) */
  for (i = from; i < to; ++i)
    {

      int row = n2o_row[i]; /* local row number */
      int globalRow = row + beg_row;  /* global row number */

      if (debug)
  {
    fprintf (logFile,
       "\nILU_pilu global: %i  old_Local: %i =========================================================\n",
       i + 1 + beg_rowP, row + 1);
  }

      EuclidGetRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
      CHECK_V_ERROR;

      if (debug)
  {
    int h;
    fprintf (logFile, "ILU_pilu  EuclidGetRow:\n");
    for (h = 0; h < len; ++h)
      fprintf (logFile, "    %i   %g\n", 1 + CVAL[h], AVAL[h]);
  }


      /* compute scaling value for row(i) */
      if (ctx->isScaled)
  {
    compute_scaling_private (i, len, AVAL, ctx);
    CHECK_V_ERROR;
  }

      SortedList_dhReset (slist, i);
      CHECK_V_ERROR;

      /* Compute symbolic factor for row(i);
         this also performs sparsification
       */
      iluk_symbolic_row_private (i, len, CVAL, AVAL,
         extRows, slist, ctx, debug);
      CHECK_V_ERROR;

      /* enforce subdomain constraint */
      SortedList_dhEnforceConstraint (slist, sg);

      /* compute numeric factor for row */
      if (!noValues)
  {
    iluk_numeric_row_private (i, extRows, slist, ctx, debug);
    CHECK_V_ERROR;
  }

      EuclidRestoreRow (ctx->A, globalRow, &len, &CVAL, &AVAL);
      CHECK_V_ERROR;

      /* Ensure adequate storage; reallocate, if necessary. */
      count = SortedList_dhReadCount (slist);
      CHECK_V_ERROR;

      /* Ensure adequate storage; reallocate, if necessary. */
      if (idx + count > F->alloc)
  {
    Factor_dhReallocate (F, idx, count);
    CHECK_V_ERROR;
    SET_INFO ("REALLOCATED from ilu_mpi_pilu");
    cval = F->cval;
    fill = F->fill;
    aval = F->aval;
  }

      /* Copy factor to permanent storage */
      if (bj)
  {     /* for debugging: blockJacobi case */
    int col;
    while (count--)
      {
        SRecord *sr = SortedList_dhGetSmallest (slist);
        CHECK_V_ERROR;
        col = sr->col;
        if (col >= beg_rowP && col < end_rowP)
    {
      cval[idx] = col;
      if (noValues)
        {
          aval[idx] = 0.0;
        }
      else
        {
          aval[idx] = sr->val;
        }
      fill[idx] = sr->level;
      ++idx;
    }
      }
  }

      if (debug)
  {
    fprintf (logFile, "ILU_pilu  ");
    while (count--)
      {
        SRecord *sr = SortedList_dhGetSmallest (slist);
        CHECK_V_ERROR;
        cval[idx] = sr->col;
        aval[idx] = sr->val;
        fill[idx] = sr->level;
        fprintf (logFile, "%i,%i,%g ; ", 1 + cval[idx], fill[idx],
           aval[idx]);
        ++idx;
      }
    fprintf (logFile, "\n");
  }

      else
  {
    while (count--)
      {
        SRecord *sr = SortedList_dhGetSmallest (slist);
        CHECK_V_ERROR;
        cval[idx] = sr->col;
        aval[idx] = sr->val;
        fill[idx] = sr->level;
        ++idx;
      }
  }

      /* add row-pointer to start of next row. */
      rp[i + 1] = idx;

      /* Insert pointer to diagonal */
      {
  int temp = rp[i];
  bool flag = true;
  while (temp < idx)
    {
      if (cval[temp] == i + beg_rowP)
        {
    diag[i] = temp;
    flag = false;
    break;
        }
      ++temp;
    }
  if (flag)
    {
      if (logFile != NULL)
        {
    int k;
    fprintf (logFile,
       "Failed to find diag in localRow %i (globalRow %i; ct= %i)\n   ",
       1 + i, i + 1 + beg_rowP, rp[i + 1] - rp[i]);
    for (k = rp[i]; k < rp[i + 1]; ++k)
      {
        fprintf (logFile, "%i ", cval[i] + 1);
      }
    fprintf (logFile, "\n\n");
        }
      sprintf (msgBuf_dh, "failed to find diagonal for localRow: %i",
         1 + i);
      SET_V_ERROR (msgBuf_dh);
    }
      }
/*
    { int temp = rp[i]; 
      while (cval[temp] != i+beg_row) ++temp;
      diag[i] = temp;
    }
*/

      /* check for zero diagonal */
      if (!aval[diag[i]])
  {
    sprintf (msgBuf_dh, "zero diagonal in local row %i", i + 1);
    SET_V_ERROR (msgBuf_dh);
  }

    }

  /* adjust to local (zero) based, if block jacobi factorization */
  if (bj)
    {
      int nz = rp[m];
      for (i = 0; i < nz; ++i)
  cval[i] -= beg_rowP;
    }

END_FUNC_DH}


#undef __FUNC__
#define __FUNC__ "iluk_symbolic_row_private"
void
iluk_symbolic_row_private (int localRow, int len, int *CVAL,
         double *AVAL, ExternalRows_dh extRows,
         SortedList_dh slist, Euclid_dh ctx, bool debug)
{
  START_FUNC_DH int level = ctx->level, m = ctx->m;
  int beg_row = ctx->sg->beg_row[myid_dh];
  int beg_rowP = ctx->sg->beg_rowP[myid_dh];
  int *cval = ctx->F->cval, *diag = ctx->F->diag;
  int *rp = ctx->F->rp, *fill = ctx->F->fill;
  int j, node, col;
  int end_rowP = beg_rowP + m;
  int level_1, level_2;
  int *cvalPtr, *fillPtr;
  SRecord sr, *srPtr;
  REAL_DH scale, *avalPtr;
  double thresh = ctx->sparseTolA;
  bool wasInserted;
  int count = 0;

  scale = ctx->scale[localRow];
  ctx->stats[NZA_STATS] += (double) len;

  /* insert col indices in sorted linked list */
  sr.level = 0;
  for (j = 0; j < len; ++j)
    {
      sr.col = CVAL[j];
      sr.val = scale * AVAL[j];
/*    if (fabs(sr.val) > thresh) { */
      wasInserted = SortedList_dhPermuteAndInsert (slist, &sr, thresh);
      CHECK_V_ERROR;
      if (wasInserted)
  ++count;
/*    } */
      if (debug)
  {
    fprintf (logFile, "ILU_pilu   inserted from A: col= %i  val= %g\n",
       1 + CVAL[j], sr.val);
  }
    }

  /* ensure diagonal entry is inserted */
  sr.val = 0.0;
  sr.col = localRow + beg_rowP;
  srPtr = SortedList_dhFind (slist, &sr);
  CHECK_V_ERROR;
  if (srPtr == NULL)
    {
      SortedList_dhInsert (slist, &sr);
      CHECK_V_ERROR;
      ++count;
      if (debug)
  {
    fprintf (logFile, "ILU_pilu   inserted missing diagonal: %i\n",
       1 + localRow + beg_row);
  }
    }
  ctx->stats[NZA_USED_STATS] += (double) count;

  /* update row from previously factored rows */
  sr.val = 0.0;
  if (level > 0)
    {
      while (1)
  {
    srPtr = SortedList_dhGetSmallestLowerTri (slist);
    CHECK_V_ERROR;
    if (srPtr == NULL)
      break;

    node = srPtr->col;

    if (debug)
      {
        fprintf (logFile, "ILU_pilu   sf updating from row: %i\n",
           1 + srPtr->col);
      }

    level_1 = srPtr->level;
    if (level_1 < level)
      {

        /* case 1: locally owned row */
        if (node >= beg_rowP && node < end_rowP)
    {
      node -= beg_rowP;
      len = rp[node + 1] - diag[node] - 1;
      cvalPtr = cval + diag[node] + 1;
      fillPtr = fill + diag[node] + 1;
    }

        /* case 2: external row */
        else
    {
      len = 0;
      ExternalRows_dhGetRow (extRows, node, &len, &cvalPtr,
           &fillPtr, &avalPtr);
      CHECK_V_ERROR;
      if (debug && len == 0)
        {
          fprintf (stderr,
             "ILU_pilu  sf failed to get extern row: %i\n",
             1 + node);
        }
    }


        /* merge in strict upper triangular portion of row */
        for (j = 0; j < len; ++j)
    {
      col = *cvalPtr++;
      level_2 = 1 + level_1 + *fillPtr++;
      if (level_2 <= level)
        {
          /* Insert new element, or update level if already inserted. */
          sr.col = col;
          sr.level = level_2;
          sr.val = 0.0;
          SortedList_dhInsertOrUpdate (slist, &sr);
          CHECK_V_ERROR;
        }
    }
      }
  }
    }
END_FUNC_DH}


#undef __FUNC__
#define __FUNC__ "iluk_numeric_row_private"
void
iluk_numeric_row_private (int new_row, ExternalRows_dh extRows,
        SortedList_dh slist, Euclid_dh ctx, bool debug)
{
  START_FUNC_DH int m = ctx->m;
  int beg_rowP = ctx->sg->beg_rowP[myid_dh];
  int end_rowP = beg_rowP + m;
  int len, row;
  int *rp = ctx->F->rp, *cval = ctx->F->cval, *diag = ctx->F->diag;
  REAL_DH *avalPtr, *aval = ctx->F->aval;
  int *cvalPtr;
  double multiplier, pc, pv;
  SRecord sr, *srPtr;

  /* note: non-zero entries from A were inserted in list during iluk_symbolic_row_private */

  SortedList_dhResetGetSmallest (slist);
  CHECK_V_ERROR;
  while (1)
    {
      srPtr = SortedList_dhGetSmallestLowerTri (slist);
      CHECK_V_ERROR;
      if (srPtr == NULL)
  break;

      /* update new_row's values from upper triangular portion of previously
         factored row
       */
      row = srPtr->col;

      if (row >= beg_rowP && row < end_rowP)
  {
    int local_row = row - beg_rowP;

    len = rp[local_row + 1] - diag[local_row];
    cvalPtr = cval + diag[local_row];
    avalPtr = aval + diag[local_row];
  }
      else
  {
    len = 0;
    ExternalRows_dhGetRow (extRows, row, &len, &cvalPtr,
         NULL, &avalPtr);
    CHECK_V_ERROR;
    if (debug && len == 0)
      {
        fprintf (stderr, "ILU_pilu  failed to get extern row: %i\n",
           1 + row);
      }

  }

      if (len)
  {
    /* first, form and store pivot */
    sr.col = row;
    srPtr = SortedList_dhFind (slist, &sr);
    CHECK_V_ERROR;
    if (srPtr == NULL)
      {
        sprintf (msgBuf_dh,
           "find failed for sr.col = %i while factoring local row= %i \n",
           1 + sr.col, new_row + 1);
        SET_V_ERROR (msgBuf_dh);
      }

    pc = srPtr->val;

    if (pc != 0.0)
      {
        pv = *avalPtr++;
        --len;
        ++cvalPtr;
        multiplier = pc / pv;
        srPtr->val = multiplier;

        if (debug)
    {
      fprintf (logFile,
         "ILU_pilu   nf updating from row: %i; multiplier = %g\n",
         1 + srPtr->col, multiplier);
    }

        /* second, update from strict upper triangular portion of row */
        while (len--)
    {
      sr.col = *cvalPtr++;
      sr.val = *avalPtr++;
      srPtr = SortedList_dhFind (slist, &sr);
      CHECK_V_ERROR;
      if (srPtr != NULL)
        {
          srPtr->val -= (multiplier * sr.val);
        }
    }
      }

    else
      {
        if (debug)
    {
      fprintf (logFile,
         "ILU_pilu   NO UPDATE from row: %i; srPtr->val = 0.0\n",
         1 + srPtr->col);
    }
      }

  }
    }
END_FUNC_DH}