/*
 * Copyright (C) 2008-2012 Delft University of Technology - All rights reserved.
 *
 * Author: Reijer Idema
 */

/*
 * Preconditioner class source.
 */
#include "Preconditioner.h"


NONLINEARSOLVER_NAMESPACE_BEGIN

//============================= ABC ==========================================

#undef __FUNCT__
#define __FUNCT__ "Preconditioner::Preconditioner"
Preconditioner::Preconditioner(void) : mName(""), mpProblem(NULL)
{
} // Preconditioner::Preconditioner

#undef __FUNCT__
#define __FUNCT__ "Preconditioner::~Preconditioner"
/* virtual */ Preconditioner::~Preconditioner(void)
{ // Always provide a virtual destructor on base classes, even if it is empty.
} // Preconditioner::~Preconditioner

#undef __FUNCT__
#define __FUNCT__ "Preconditioner::Initialize"
/* virtual */ PetscErrorCode Preconditioner::Initialize(NonLinearProblem& rProblem)
{
	mpProblem = &rProblem;
	return 0;
} // Preconditioner::Initialize

/* virtual */ std::string Preconditioner::GetName(void) const { return mName; }

//============================= JACOBIAN ===============================

#undef __FUNCT__
#define __FUNCT__ "PreconditionerJacobian::PreconditionerJacobian"
PreconditionerJacobian::PreconditionerJacobian(void)
{ // Base constructor is automatically called
	mName = "jacobian";
} // PreconditionerJacobian::PreconditionerJacobian

#undef __FUNCT__
#define __FUNCT__ "PreconditionerJacobian::GetMatrix"
/* virtual */ Mat PreconditionerJacobian::GetMatrix(void)
{
	if (mpProblem == NULL) return NULL;
	return mpProblem->GetJacobianMatrix();
} // PreconditionerJacobian::GetMatrix

#undef __FUNCT__
#define __FUNCT__ "PreconditionerJacobian::Calculate"
/* virtual */ PetscErrorCode PreconditionerJacobian::Calculate(SNES snes, PetscUtils::ConstVec pX, PetscUtils::ConstMat pJ, Mat* ppM, MatStructure* pFlag)
{
	// Check if the SNES preconditioner matrix points to the correct data
	// If not it can probably be fixed with a simple: *ppM = mpMatrix;
	if (*ppM != pJ) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_GENERAL, "SNES preconditioner matrix no longer points to the Jacobian");
	*pFlag = (mpProblem->IsJacobianMatrixStructureSame() ? SAME_NONZERO_PATTERN : DIFFERENT_NONZERO_PATTERN);
	return 0;
} // PreconditionerJacobian::Calculate

//============================= INITIAL JACOBIAN =============================

#undef __FUNCT__
#define __FUNCT__ "PreconditionerInitialJacobian::PreconditionerInitialJacobian"
PreconditionerInitialJacobian::PreconditionerInitialJacobian(void) : mpOwnedMatrix(PETSC_NULL)
{ // Base constructor is automatically called
	mName = "initial jacobian";
} // PreconditionerInitialJacobian::PreconditionerInitialJacobian

#undef __FUNCT__
#define __FUNCT__ "PreconditionerInitialJacobian::~PreconditionerInitialJacobian"
/* virtual */ PreconditionerInitialJacobian::~PreconditionerInitialJacobian(void)
{ // Base destructor is automatically called
	if (mpOwnedMatrix != PETSC_NULL)
	{
		PetscErrorCode err = MatDestroy(&mpOwnedMatrix); mpOwnedMatrix = PETSC_NULL; CHKERRV(err);
	}
} // PreconditionerInitialJacobian::~PreconditionerInitialJacobian

#undef __FUNCT__
#define __FUNCT__ "PreconditionerInitialJacobian::GetMatrix"
/* virtual */ Mat PreconditionerInitialJacobian::GetMatrix(void)
{
	return mpOwnedMatrix;
} // PreconditionerInitialJacobian::GetMatrix

#undef __FUNCT__
#define __FUNCT__ "PreconditionerInitialJacobian::Initialize"
/* virtual */ PetscErrorCode PreconditionerInitialJacobian::Initialize(NonLinearProblem& rProblem)
{
	PetscErrorCode err;
	if (mpOwnedMatrix != PETSC_NULL)
	{
		err = MatDestroy(&mpOwnedMatrix); mpOwnedMatrix = PETSC_NULL; CHKERRQ(err);
	}

	Mat pInitialJacobian = rProblem.GetJacobianMatrix();
	if (pInitialJacobian == PETSC_NULL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_GENERAL, "Jacobian matrix is NULL");

	// Copy initial jacobian matrix
	err = MatDuplicate(pInitialJacobian, MAT_COPY_VALUES, &mpOwnedMatrix); CHKERRQ(err);
	// :WARNING: Not assembling crashes LU factorization within SNESSolve(), although ILU seems to work fine without.
	err = MatAssemblyBegin(mpOwnedMatrix, MAT_FINAL_ASSEMBLY); CHKERRQ(err);
	err = MatAssemblyEnd(mpOwnedMatrix, MAT_FINAL_ASSEMBLY); CHKERRQ(err);

	return Preconditioner::Initialize(rProblem);
} // PreconditionerInitialJacobian::Initialize

#undef __FUNCT__
#define __FUNCT__ "PreconditionerInitialJacobian::Calculate"
/* virtual */ PetscErrorCode PreconditionerInitialJacobian::Calculate(SNES snes, PetscUtils::ConstVec pX, PetscUtils::ConstMat pJ, Mat* ppM, MatStructure* pFlag)
{
	// Check if the SNES preconditioner matrix points to the correct data
	// If not it can probably be fixed with a simple: *ppM = mpMatrix;
	if (*ppM != mpOwnedMatrix) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_GENERAL, "SNES preconditioner matrix no longer points to the correct data");
	*pFlag = SAME_PRECONDITIONER;
	return 0;
} // PreconditionerInitialJacobian::Calculate

//============================= FIXED ========================================

#undef __FUNCT__
#define __FUNCT__ "PreconditionerFixed::PreconditionerFixed"
PreconditionerFixed::PreconditionerFixed(void) :
	mpMatrix(PETSC_NULL), mpOwnedMatrix(PETSC_NULL)
{ // Base constructor is automatically called
	mName = "fixed";
} // PreconditionerFixed::PreconditionerFixed

#undef __FUNCT__
#define __FUNCT__ "PreconditionerFixed::~PreconditionerFixed"
/* virtual */ PreconditionerFixed::~PreconditionerFixed(void)
{ // Base destructor is automatically called
	if (mpOwnedMatrix != PETSC_NULL)
	{
		PetscErrorCode err = MatDestroy(&mpOwnedMatrix); mpOwnedMatrix = PETSC_NULL; CHKERRV(err);
	}
} // PreconditionerFixed::~PreconditionerFixed

#undef __FUNCT__
#define __FUNCT__ "PreconditionerFixed::SetMatrix"
/* virtual */ bool PreconditionerFixed::SetMatrix(Mat pMat)
{
	// Destroy current matrix if it exists and is owned
	if (mpOwnedMatrix != PETSC_NULL)
	{
		PetscErrorCode err = MatDestroy(&mpOwnedMatrix); mpOwnedMatrix = PETSC_NULL; CHKERRB(err);
	}
	// Set preconditioner
	mpMatrix = pMat;
	return true;
} // PreconditionerFixed::SetMatrix

#undef __FUNCT__
#define __FUNCT__ "PreconditionerFixed::SetMatrixCopy"
/* virtual */ bool PreconditionerFixed::SetMatrixCopy(Mat pMat)
{
	PetscErrorCode err;
	// Destroy current matrix if it exists and is owned
	if (mpOwnedMatrix != PETSC_NULL)
	{
		err = MatDestroy(&mpOwnedMatrix); mpOwnedMatrix = PETSC_NULL; CHKERRB(err);
	}
	// Copy provided matrix
	err = MatDuplicate(pMat, MAT_COPY_VALUES, &mpOwnedMatrix); CHKERRB(err);
	// :WARNING: Not assembling crashes LU factorization within SNESSolve(), although ILU seems to work fine without.
	err = MatAssemblyBegin(mpOwnedMatrix, MAT_FINAL_ASSEMBLY); CHKERRB(err);
	err = MatAssemblyEnd(mpOwnedMatrix, MAT_FINAL_ASSEMBLY); CHKERRB(err);
	// Set copy as preconditioner
	mpMatrix = mpOwnedMatrix;
	return true;
} // PreconditionerFixed::SetMatrixCopy

#undef __FUNCT__
#define __FUNCT__ "PreconditionerFixed::GetMatrix"
/* virtual */ Mat PreconditionerFixed::GetMatrix(void)
{
	return mpMatrix;
} // PreconditionerFixed::GetMatrix

#undef __FUNCT__
#define __FUNCT__ "PreconditionerFixed::Calculate"
/* virtual */ PetscErrorCode PreconditionerFixed::Calculate(SNES snes, PetscUtils::ConstVec pX, PetscUtils::ConstMat pJ, Mat* ppM, MatStructure* pFlag)
{
	// Check if the SNES preconditioner matrix points to the correct data
	// If not it can probably be fixed with a simple: *ppM = mpMatrix;
	if (*ppM != mpMatrix) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_GENERAL, "SNES preconditioner matrix no longer points to the correct data");
	*pFlag = SAME_PRECONDITIONER;
	return 0;
} // PreconditionerFixed::Calculate

//============================================================================

NONLINEARSOLVER_NAMESPACE_END
