Added support for element-by-element preconditioner using option pc = mat

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1457 e10b68d5-8a6e-419e-a041-bce267b0401d

src/LinAlg/PETScMatrix.C

@@ -15,6 +15,8 @@
 #ifdef HAS_PETSC
 #include "SAM.h"
 #include "petscversion.h"
+#include "petscis.h"
+#include "petscsys.h"
 
 #if PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 2
 #include "petscpcmg.h"
@@ -186,6 +188,9 @@ PETScMatrix::PETScMatrix(const LinSolParams& spar) : solParams(spar)
     KSPSetNullSpace(ksp,nsp);
   }
   LinAlgInit::increfs();
+
+  elmIS = 0;
+  ISsize = 0;
 }
 
 
@@ -204,6 +209,9 @@ PETScMatrix::PETScMatrix (const PETScMatrix& B) : solParams(B.solParams)
     KSPSetNullSpace(ksp,nsp);
   }
   LinAlgInit::increfs();
+
+  elmIS = 0;
+  ISsize = 0;
 }
 
 
@@ -219,6 +227,10 @@ PETScMatrix::~PETScMatrix ()
   // Deallocation of matrix object.
   MatDestroy(PETSCMANGLE(A));   
   LinAlgInit::decrefs();
+
+  for (int i = 0;i < ISsize;i++)
+    ISDestroy(elmIS[i]);
+  delete elmIS;
 }
 
 
@@ -445,6 +457,10 @@ void PETScMatrix::initAssembly (const SAM& sam)
   // Get number of equations in linear system
   const PetscInt neq = sam.getNoEquations();
 
+  // RUNAR
+  if (!strncasecmp(solParams.getPreconditioner(),"mat",3))
+    std::cout << "EBE = " << this->makeElementIS(sam) << std::endl;
+
   // Set correct number of rows and columns for matrix.
   MatSetSizes(A,neq,neq,PETSC_DECIDE,PETSC_DECIDE);
   MPI_Barrier(PETSC_COMM_WORLD);
@@ -577,6 +593,8 @@ bool PETScMatrix::multiply (const SystemVector& B, SystemVector& C)
 
 bool PETScMatrix::solve (SystemVector& B, bool newLHS)
 {
+  Mat AebeI;
+
   const PETScVector* Bptr = dynamic_cast<PETScVector*>(&B);
   if (!Bptr)
     return false;
@@ -586,7 +604,12 @@ bool PETScMatrix::solve (SystemVector& B, bool newLHS)
   VecCopy(Bptr->getVector(),x);
   
   // Has lefthand side changed?
-  if (newLHS)
+  if (!strncasecmp(solParams.getPreconditioner(),"mat",3)) {
+    if (!this->makeEBEpreconditioner(A,&AebeI))
+      return false;
+    KSPSetOperators(ksp,A,AebeI,SAME_NONZERO_PATTERN);
+  }
+  else if (newLHS)
     KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);
   else
     KSPSetOperators(ksp,A,A,SAME_PRECONDITIONER);
@@ -600,6 +623,9 @@ bool PETScMatrix::solve (SystemVector& B, bool newLHS)
   PetscPrintf(PETSC_COMM_WORLD,"\n Iterations for %s = %D\n",solParams.getMethod(),its);
   VecDestroy(PETSCMANGLE(x));
 
+  if (ISsize > 0)
+    MatDestroy(AebeI);
+
   return true;
 }
 
@@ -731,4 +757,87 @@ real PETScMatrix::Linfnorm () const
   return norm;
 }
 
+
+bool PETScMatrix::makeElementIS(const SAM& sam)
+{
+  std::vector<PetscInt> meen;
+
+  ISsize = sam.getNoElms();
+
+  elmIS = new IS[ISsize];
+  for (int e = 1;e <= ISsize;e++) {
+    if (!sam.getElmEqns(meen,e))
+      return false;
+    std::sort(meen.begin(),meen.end());
+
+    for (size_t i = 0;i < meen.size();i++)
+      meen[i]--;
+
+    int ndof = 0;
+    for (size_t i = 0;i < meen.size();i++)
+      if (meen[i] >= 0)
+	ndof++;
+
+    PetscInt l2g[ndof];
+    ndof = 0;
+    for (size_t i = 0;i < meen.size();i++)
+      if (meen[i] >= 0)
+	l2g[ndof++] = meen[i];
+
+    ISCreateGeneral(PETSC_COMM_SELF,ndof,l2g,&(elmIS[e-1]));
+  }
+   
+  return true;
+}
+
+
+bool PETScMatrix::makeEBEpreconditioner(const Mat A, Mat* AeI)
+{
+  PetscInt         nedof;
+  const PetscInt*  indx;
+  Vector           vals;
+  std::vector<int> locidx;
+
+  if (!elmIS)
+    return false;
+
+  Mat* subMats;
+  MatGetSubMatrices(A,ISsize,elmIS,elmIS,MAT_INITIAL_MATRIX,&subMats);
+
+  MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,AeI);
+  MatZeroEntries(*AeI);
+
+  Matrix elmMat;
+  for (PetscInt e = 0;e < ISsize;e++) {
+    ISGetSize(elmIS[e],&nedof);
+
+    locidx.resize(nedof);
+    for (int i = 0;i < nedof;i++)
+      locidx[i] = i;
+
+    vals.resize(nedof*nedof,true);
+    MatGetValues(subMats[e],nedof,&(locidx.front()),nedof,&(locidx.front()),&(vals.front()));
+
+    int dof = 0;
+    elmMat.resize(nedof,nedof);
+    for (int i = 1;i <= nedof;i++)
+      for (int j = 1;j <= nedof;j++)
+	elmMat(i,j) = vals[dof++];
+
+    utl::invert(elmMat);
+    elmMat.transpose();
+
+    ISGetIndices(elmIS[e],&indx);
+    MatSetValues(*AeI,nedof,indx,nedof,indx,elmMat.ptr(),ADD_VALUES);
+    ISRestoreIndices(elmIS[e],&indx);
+  }
+
+   MatDestroyMatrices(ISsize,&subMats);
+
+   MatAssemblyBegin(*AeI,MAT_FINAL_ASSEMBLY);
+   MatAssemblyEnd(*AeI,MAT_FINAL_ASSEMBLY);
+
+   return true;
+}
+
 #endif // HAS_PETSC

src/LinAlg/PETScMatrix.h

@@ -227,6 +227,14 @@ private:
   KSP                 ksp;       //!< Linear solver
   MatNullSpace        nsp;       //!< Null-space of linear operator
   const LinSolParams& solParams; //!< Linear solver parameters
+  IS*                 elmIS;     //!< Element index sets
+  PetscInt            ISsize;    //!< Number of index sets/elements
+
+  //! \brief Constructs index set needed for element-by-element preconditioner
+  bool makeElementIS(const SAM& sam);
+
+  //! \brief Constructs the EBE preconditioner of the given matrix
+  bool makeEBEpreconditioner(const Mat A, Mat* AeI);
 
 #else // dummy implementation when PETSc is not included
   virtual SystemMatrix* copy() const { return 0; }