ParallelCoarseSpaceMatVec.cc

//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


#include "AMDiS.h"
#include "parallel/ParallelCoarseSpaceMatVec.h"
#include "parallel/MatrixNnzStructure.h"

namespace AMDiS {

  using namespace std;

  ParallelCoarseSpaceMatVec::ParallelCoarseSpaceMatVec()
    : interiorMap(NULL),
      lastMeshNnz(-1),
      alwaysCreateNnzStructure(false),
      meshDistributor(NULL),
      subdomainLevel(0),
      rStartInterior(0),
      nGlobalOverallInterior(0)
  {
    Parameters::get("parallel->always create nnz structure", 
		    alwaysCreateNnzStructure);
  }


  void ParallelCoarseSpaceMatVec::setCoarseSpaceDofMapping(ParallelDofMapping *coarseDofs, 
							   int component)
  {
    FUNCNAME("ParallelCoarseSpaceMatVec::setCoarseSpaceDofMapping()");

    TEST_EXIT_DBG(coarseDofs)("Should not happen!\n");

    if (component == -1) {
      // === Set coarse space for all components. ===

      coarseSpaceMap.clear();

      int nComponents = coarseDofs->getNumberOfComponents();
      for (int i = 0; i < nComponents; i++)
	coarseSpaceMap[i] = coarseDofs;
    } else {
      // === Set coarse space for just one component. ===

      coarseSpaceMap[component] = coarseDofs;
    }

    if (find(uniqueCoarseMap.begin(), uniqueCoarseMap.end(), coarseDofs) == 
	uniqueCoarseMap.end())
      uniqueCoarseMap.push_back(coarseDofs);
  }


  void ParallelCoarseSpaceMatVec::prepare()
  {
    FUNCNAME("ParallelCoarseSpaceMatVec:prepare()");

    TEST_EXIT(uniqueCoarseMap.size() <= 2)
      ("Not yet implemented for more than two coarse spaces!\n");

    // === Create pointers to PETSc matrix and vector objects. ===

    int nCoarseMap = uniqueCoarseMap.size();
    mat.resize(nCoarseMap + 1);
    for (int i = 0; i < nCoarseMap + 1; i++)
      mat[i].resize(nCoarseMap + 1);

    vecSol.resize(nCoarseMap + 1);
    vecRhs.resize(nCoarseMap + 1);
   
    // === Create map from component number to its coarse space map. ===

    componentIthCoarseMap.resize(coarseSpaceMap.size());
    for (unsigned int i = 0; i < coarseSpaceMap.size(); i++) {
      bool found = false;
      for (int j = 0; j < nCoarseMap; j++) {
	if (coarseSpaceMap[i] == uniqueCoarseMap[j]) {
	  componentIthCoarseMap[i] = j;
	  found = true;
	  break;
	}
      }
      
      TEST_EXIT_DBG(found)("Should not happen!\n");
    }
  }


  void ParallelCoarseSpaceMatVec::createMatVec(Matrix<DOFMatrix*>& seqMat)
  {
    FUNCNAME("ParallelCoarseSpaceMatVec::createMatVec()");

    // === Prepare coarse space information and generate the correct number ===
    // === of empty PETSc matrix and vector objects.                        ===

    prepare();

    
    // === Update subdomain data (mostly required for multi-level methods) ===

    updateSubdomainData();


    // === If required, recompute non zero structure of the matrix. ===

    bool localMatrix = (coarseSpaceMap.size() && subdomainLevel == 0);

    if (checkMeshChange()) {
      // Mesh has been changed, recompute interior DOF mapping.
      interiorMap->setComputeMatIndex(!localMatrix);
      interiorMap->update();

      int nMat = uniqueCoarseMap.size() + 1;
      nnz.resize(nMat);
      for (int i = 0; i < nMat; i++) {
	nnz[i].resize(nMat);
	for (int j = 0; j < nMat; j++)
	  nnz[i][j].clear();
      }

      nnz[0][0].create(seqMat, mpiCommGlobal, *interiorMap,
		       (coarseSpaceMap.size() == 0 ? &(meshDistributor->getPeriodicMap()) : NULL),
		       meshDistributor->getElementObjectDb(),
		       localMatrix);

      for (int i = 0; i < nMat; i++) {
	for (int j = 0; j < nMat; j++) {
	  if (i == 0 && j == 0)
	    continue;

	  ParallelDofMapping &rowMap = 
	    (i == 0 ? *interiorMap : *(uniqueCoarseMap[i - 1]));
	  ParallelDofMapping &colMap =
	    (j == 0 ? *interiorMap : *(uniqueCoarseMap[j - 1]));

	  nnz[i][j].create(seqMat, mpiCommGlobal, rowMap, colMap, NULL,
			   meshDistributor->getElementObjectDb());
	}
      }
    }


    // === Create PETSc matrices and PETSc vectors with the computed  ===
    // === nnz data structure.                                        ===
    
    int nRankInteriorRows = interiorMap->getRankDofs();
    int nOverallInteriorRows = interiorMap->getOverallDofs();
    
    if (localMatrix) {
      MatCreateSeqAIJ(mpiCommLocal, nRankInteriorRows, nRankInteriorRows,
		      0, nnz[0][0].dnnz,
		      &mat[0][0]);
      MatSetOption(mat[0][0], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
    } else {
      MatCreateAIJ(mpiCommGlobal, nRankInteriorRows, nRankInteriorRows, 
		   nOverallInteriorRows, nOverallInteriorRows,
		   0, nnz[0][0].dnnz, 0, nnz[0][0].onnz,		   
		   &mat[0][0]);
      MatSetOption(mat[0][0], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
    }

    VecCreateMPI(mpiCommGlobal, nRankInteriorRows, nOverallInteriorRows, &vecSol[0]);
    VecCreateMPI(mpiCommGlobal, nRankInteriorRows, nOverallInteriorRows, &vecRhs[0]);

    int nCoarseMap = uniqueCoarseMap.size();
    for (int i = 0; i < nCoarseMap; i++) {
      ParallelDofMapping* cMap = uniqueCoarseMap[i];
      
      int nRankCoarseRows = cMap->getRankDofs();
      int nOverallCoarseRows = cMap->getOverallDofs();
      
      MatCreateAIJ(mpiCommGlobal,
		   nRankCoarseRows, nRankCoarseRows,
		   nOverallCoarseRows, nOverallCoarseRows,
		   0, nnz[i + 1][i + 1].dnnz, 0, nnz[i + 1][i + 1].onnz,
		   &mat[i + 1][i + 1]);
      cMap->createVec(vecSol[i + 1]);
      cMap->createVec(vecRhs[i + 1]);
    }

    for (int i = 0; i < nCoarseMap + 1; i++) {
      for (int j = 0; j < nCoarseMap + 1; j++) {
	if (i == j)
	  continue;

	int nRowsRankMat = (i == 0 ? nRankInteriorRows : uniqueCoarseMap[i - 1]->getRankDofs());
	int nRowsOverallMat = (i == 0 ? nOverallInteriorRows : uniqueCoarseMap[i - 1]->getOverallDofs());

	int nColsRankMat = (j == 0 ? nRankInteriorRows : uniqueCoarseMap[j - 1]->getRankDofs());
	int nColsOverallMat = (j == 0 ? nOverallInteriorRows : uniqueCoarseMap[j - 1]->getOverallDofs());

	MatCreateAIJ(mpiCommGlobal,
		     nRowsRankMat, nColsRankMat,
		     nRowsOverallMat, nColsOverallMat,
		     0, nnz[i][j].dnnz, 0, nnz[i][j].onnz,
		     &mat[i][j]);	  
	MatCreateAIJ(mpiCommGlobal,
		     nColsRankMat, nRowsRankMat,
		     nColsOverallMat, nRowsOverallMat,
		     0, nnz[j][i].dnnz, 0, nnz[j][i].onnz,
		     &mat[j][i]);
      }
    }    
  }


  void ParallelCoarseSpaceMatVec::matDestroy()
  {
    FUNCNAME("ParallelCoarseSpaceMatVec::matDestroy()");

    int nMatrix = mat.size();
    for (int i = 0; i < nMatrix; i++)
      for (int j = 0; j < nMatrix; j++)
	MatDestroy(&mat[i][j]);
  }


  void ParallelCoarseSpaceMatVec::vecDestroy()
  {
    FUNCNAME("ParallelCoarseSpaceMatVec::vecDestroy()");

    int nVec = vecSol.size();
    for (int i = 0; i < nVec; i++) {
      VecDestroy(&vecSol[i]);
      VecDestroy(&vecRhs[i]);
    }
  }


  void ParallelCoarseSpaceMatVec::matAssembly()
  {
    FUNCNAME("ParallelCoarseSpaceMatVec::matAssembly()");

    int nMatrix = mat.size();
    for (int i = 0; i < nMatrix; i++) {
      for (int j = 0; j < nMatrix; j++) {
	MatAssemblyBegin(mat[i][j], MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(mat[i][j], MAT_FINAL_ASSEMBLY);  
      }
    }
  }


  void ParallelCoarseSpaceMatVec::vecRhsAssembly()
  {
    FUNCNAME("ParallelCoarseSpaceMatVec::vecRhsAssembly()");

    int nVec = vecRhs.size();
    for (int i = 0; i < nVec; i++) {
      VecAssemblyBegin(vecRhs[i]);
      VecAssemblyEnd(vecRhs[i]);
    }
  }


  void ParallelCoarseSpaceMatVec::vecSolAssembly()
  {
    FUNCNAME("ParallelCoarseSpaceMatVec::vecSolAssembly()");

    int nVec = vecRhs.size();
    for (int i = 0; i < nVec; i++) {
      VecAssemblyBegin(vecSol[i]);
      VecAssemblyEnd(vecSol[i]);
    }
  }


  bool ParallelCoarseSpaceMatVec::checkMeshChange()
  {
    FUNCNAME("ParallelCoarseSpaceMatVec::checkMeshChange()");

    int recvAllValues = 0;
    int sendValue = 
      static_cast<int>(meshDistributor->getLastMeshChangeIndex() != lastMeshNnz);
    mpiCommGlobal.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);

    if (recvAllValues != 0 || alwaysCreateNnzStructure) {
      lastMeshNnz = meshDistributor->getLastMeshChangeIndex();
      return true;
    }

    return false;
  }


  void ParallelCoarseSpaceMatVec::updateSubdomainData()
  {
    FUNCNAME("ParallelCoarseSpaceMatVec::updateSubdomainData()");

    if (mpiCommLocal.Get_size() == 1) {
      rStartInterior = 0;
      nGlobalOverallInterior = interiorMap->getOverallDofs();
    } else {
      int groupRowsInterior = 0;
      if (mpiCommLocal.Get_rank() == 0)
	groupRowsInterior = interiorMap->getOverallDofs();

      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
			   rStartInterior, nGlobalOverallInterior);

      int tmp = 0;
      if (mpiCommLocal.Get_rank() == 0)
	tmp = rStartInterior;

      mpiCommLocal.Allreduce(&tmp, &rStartInterior, 1, MPI_INT, MPI_SUM);
    }
  }

}