//
// 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 "parallel/PetscSolver.h"
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"

namespace AMDiS {

  using namespace std;

  void PetscSolver::printSolutionInfo(AdaptInfo *adaptInfo,
				      bool iterationCounter, 
				      bool residual)
  {
    FUNCNAME("PetscSolver::printSolutionInfo()");

    if (iterationCounter) {
      int iterations = 0;
      KSPGetIterationNumber(solver, &iterations);
      MSG("  Number of iterations: %d\n", iterations);
      adaptInfo->setSolverIterations(iterations);
    }

    if (residual) {
      double norm = 0.0;
      MatMult(petscMatrix, petscSolVec, petscTmpVec);
      VecAXPY(petscTmpVec, -1.0, petscRhsVec);
      VecNorm(petscTmpVec, NORM_2, &norm);
      MSG("  Residual norm: %e\n", norm);
    }
  }


  void PetscSolver::copyVec(Vec& originVec, Vec& destVec, 
			    vector<int>& originIndex, vector<int>& destIndex)
  {
    FUNCNAME("PetscSolver::copyVec()");

    IS originIs, destIs;
    ISCreateGeneral(PETSC_COMM_WORLD, 
		    originIndex.size(), 
		    &(originIndex[0]),
		    PETSC_USE_POINTER,
		    &originIs);

    ISCreateGeneral(PETSC_COMM_WORLD, 
		    destIndex.size(), 
		    &(destIndex[0]),
		    PETSC_USE_POINTER,
		    &destIs);

    VecScatter scatter;
    VecScatterCreate(originVec, originIs, destVec, destIs, &scatter);
    VecScatterBegin(scatter, originVec, destVec,
		    INSERT_VALUES, SCATTER_FORWARD);
    VecScatterEnd(scatter, originVec, destVec,
		  INSERT_VALUES, SCATTER_FORWARD);

    ISDestroy(&originIs);
    ISDestroy(&destIs);    
    VecScatterDestroy(&scatter);
  }


  vector<const FiniteElemSpace*> PetscSolver::getFeSpaces(Matrix<DOFMatrix*> *mat)
  {
    FUNCNAME("PetscSolver::getFeSpaces()");

    int nComponents = mat->getNumRows();
    vector<const FiniteElemSpace*> result(nComponents);

    for (int i = 0; i < nComponents; i++) 
      for (int j = 0; j < nComponents; j++)
	if ((*mat)[i][j]) {
	  result[i] = (*mat)[i][j]->getRowFeSpace();
	  break;
	}

#if (DEBUG != 0)
    // === In debug mode, we test if all FE spaces of the matrix are also ===
    // === considered by the mesh distributor.                            ===
    checkFeSpaces(result);
#endif

    return result;
  }


  vector<const FiniteElemSpace*> PetscSolver::getFeSpaces(SystemVector *vec)
  {
    FUNCNAME("PetscSolver::getFeSpaces()");

    int nComponents = vec->getSize();
    vector<const FiniteElemSpace*> result(nComponents);

    for (int i = 0; i < nComponents; i++) 
      result[i] = vec->getFeSpace(i);

#if (DEBUG != 0)
    // === In debug mode, we test if all FE spaces of the matrix are also ===
    // === considered by the mesh distributor.                            ===
    checkFeSpaces(result);
#endif

    return result;
  }


  void PetscSolver::checkFeSpaces(vector<const FiniteElemSpace*>& feSpaces)
  {
    FUNCNAME("PetscSolver::checkFeSpaces()");

    for (unsigned int i = 0; i < feSpaces.size(); i++) {
      bool found = false;
      for (unsigned int j = 0; j < meshDistributor->getFeSpaces().size(); j++)
	if (feSpaces[i] == meshDistributor->getFeSpace(j)) {
	  found = true;
	  break;
	}
      TEST_EXIT(found)("FE spaces not found in mesh distributor!\n");
    }
  }

}