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(0),
      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;
    }
  }


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

    // === Create vector of unique pointers to all coarse space maps. ===

    uniqueCoarseMap.clear();
    if (coarseSpaceMap.size()) {
      std::set<ParallelDofMapping*> tmp;
      for (map<int, ParallelDofMapping*>::iterator it = coarseSpaceMap.begin();
	   it != coarseSpaceMap.end(); ++it) {
	if (tmp.count(it->second) == 0) {	  
	  tmp.insert(it->second);
	  uniqueCoarseMap.push_back(it->second);
	}
      }
    }


    // === 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 < componentIthCoarseMap.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.
      vector<const FiniteElemSpace*> feSpaces = getComponentFeSpaces(seqMat);     
      interiorMap->setComputeMatIndex(!localMatrix);
      interiorMap->update(feSpaces);

      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]);
      MSG("REMOVE THIS!\n");
      MatSetOption(mat[i + 1][i + 1], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

      VecCreateMPI(mpiCommGlobal, nRankCoarseRows, nOverallCoarseRows, 
		   &vecSol[i + 1]);
      VecCreateMPI(mpiCommGlobal, nRankCoarseRows, nOverallCoarseRows, 
		   &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]);
	MSG("REMOVE THIS!\n");
	MatSetOption(mat[i][j], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
	  
	MatCreateAIJ(mpiCommGlobal,
		     nColsRankMat, nRowsRankMat,
		     nColsOverallMat, nRowsOverallMat,
		     0, nnz[j][i].dnnz, 0, nnz[j][i].onnz,
		     &mat[j][i]);
	MSG("REMOVE THIS!\n");
	MatSetOption(mat[j][i], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

      }
    }    
  }


  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);
    }
  }

}