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

/*
 * PowerFlower class source.
 */
#include <vector>										// USES vector
#include "../PowerSystem/Interface/MatPowerInterface.h"	// USES MatPowerInterface
#include "../NonLinearSolver/NonLinearSolver.h"			// USES NonLinearSolver
#include "../NonLinearSolver/Preconditioner.h"			// USES Preconditioner
#include "PowerFlowProblem.h"							// USES PowerFlowProblem
#include "ForcingTerms.h"								// USES ForcingTerms
#include "PowerFlower.h"


POWERFLOWER_NAMESPACE_BEGIN

NONLINEARSOLVER_NAMESPACE_USING


struct SolveResult
{
	int solveCode;		// 0 if skipped, >0 if converged, <0 if diverged.
	double residualNorm;
	int nonlinearIterationCount;
	int linearIterationCount;
};


// Initialize static member objects
PowerFlower::ObjectListType PowerFlower::msObjectList;
LogEvents PowerFlower::msLogEvents;
const std::string PowerFlower::msPetscOptionFile = "-f";
const std::string PowerFlower::msPetscOptionPrecoMatrixType = "-pc_mat_type";
const std::string PowerFlower::msPetscOptionCAStartType = "-ca_start_type";
const std::string PowerFlower::msPetscOptionCAPrecoType = "-ca_preco_type";



/////////////////////////////// PUBLIC ///////////////////////////////////////


//============================= INIT/FINAL ===================================


#undef __FUNCT__
#define __FUNCT__ "PowerFlower::Initialize"
/* static */ bool PowerFlower::Initialize(int argc, char** ppArgv, char* pHelp)
{
	PetscErrorCode err;

	// Initialize PETSc if not already initialized or finalized
	if (IsInitialized())
	{
		err = PetscInfo(PETSC_NULL, "[Failed] PETSc is already initialized.\n"); CHKERRB(err);
		return false;
	}
	if (IsFinalized()) return false;

	// Initialize, get processor count (size), and get processor id (rank)
	err = PetscInitialize(&argc, &ppArgv, (char*)0, pHelp); CHKERRB(err);
	err = PetscInfo(PETSC_NULL, "PETSc Initialized.\n"); CHKERRB(err);
	PetscMPIInt mpiSize;	// Number of processors on which the code is running
	err = MPI_Comm_size(PETSC_COMM_WORLD, &mpiSize); CHKERRB(err);
	PetscMPIInt mpiRank;	// Rank of the processor that owns the PowerFlower object
	err = MPI_Comm_rank(PETSC_COMM_WORLD, &mpiRank); CHKERRB(err);

	// Check for sequential or parallel execution
	if (mpiSize > 1)
	{
		err = PetscInfo(PETSC_NULL, "[Failed] PowerFlower framework is not designed for parallel execution.\n"); CHKERRB(err);
		return false;
	}

	err = PetscInfo(PETSC_NULL, "Registering PowerFlower events.\n"); CHKERRB(err);
	if (!msLogEvents.Register()) return false;

	err = PetscInfo(PETSC_NULL, "PowerFlower framework initialized.\n"); CHKERRB(err);
	return true;
} // PowerFlower::Initialize


#undef __FUNCT__
#define __FUNCT__ "PowerFlower::IsInitialized"
/* static */ bool PowerFlower::IsInitialized(void)
{
	PetscBool petscInitialized;
	PetscInitialized(&petscInitialized);
	return petscInitialized;
	// ALTERTANIVELY USE: return PetscUtils::IsInitialized();
} // PowerFlower::IsInitialized


#undef __FUNCT__
#define __FUNCT__ "PowerFlower::Finalize"
/* static */ bool PowerFlower::Finalize(void)
{
	PetscErrorCode err;
	if (IsInitialized())
	{
		err = PetscInfo(PETSC_NULL, "Finalizing PowerFlower framework.\n"); CHKERRB(err);
	}

	// Clear all active objects.
	while (!msObjectList.empty())
	{
		msObjectList.back()->Clear();	// clear last object in list
		msObjectList.pop_back();		// delete last object from list
	}

	// Finalize PETSc
	if (IsInitialized())
	{
		err = PetscInfo(PETSC_NULL, "Finalizing PETSc.\n"); CHKERRB(err);
		err = PetscFinalize(); CHKERRB(err);
	}
	return true;
} // PowerFlower::Finalize


#undef __FUNCT__
#define __FUNCT__ "PowerFlower::IsFinalized"
/* static */ bool PowerFlower::IsFinalized(void)
{
	PetscBool petscFinalized;
	PetscFinalized(&petscFinalized);
	return petscFinalized;
} // PowerFlower::IsFinalized


//============================= LIFECYCLE ====================================


#undef __FUNCT__
#define __FUNCT__ "PowerFlower::PowerFlower"
PowerFlower::PowerFlower(void) : mpPowerSystem(NULL)
{
	if (IsInitialized())
	{
		PetscErrorCode err = PetscInfo(PETSC_NULL, "Creating PowerFlower object.\n"); CHKERRV(err);
	}
	msObjectList.push_back(this);
} // PowerFlower::PowerFlower


#undef __FUNCT__
#define __FUNCT__ "PowerFlower::~PowerFlower"
PowerFlower::~PowerFlower(void)
{
	Clear();
	msObjectList.remove(this);
	if (IsInitialized())
	{
		PetscErrorCode err = PetscInfo(PETSC_NULL, "Destroying PowerFlower object.\n"); CHKERRV(err);
	}
} // PowerFlower::~PowerFlower


//============================= OPERATIONS ===================================


#undef __FUNCT__
#define __FUNCT__ "PowerFlower::Clear"
void PowerFlower::Clear(void)
{
	if (IsInitialized())
	{
		PetscErrorCode err = PetscInfo(PETSC_NULL, "Clearing PowerFlower object data.\n"); CHKERRV(err);
	}
	if (mpPowerSystem != NULL) { delete mpPowerSystem; mpPowerSystem = NULL; }
} // PowerFlower::Clear


#undef __FUNCT__
#define __FUNCT__ "PowerFlower::ReadPowerSystem"
bool PowerFlower::ReadPowerSystem(std::string& rSourceName, PowerSystemExt::Interface& rInterface)
{
	if (mpPowerSystem == NULL) mpPowerSystem = new PowerSystemExt(&msLogEvents);
	if (!mpPowerSystem->Read(rSourceName, rInterface)) return false;
	if (!mpPowerSystem->BuildAdmittanceMatrix()) return false;
	if (!mpPowerSystem->BuildPowerVector()) return false;
	if (!mpPowerSystem->BuildVoltageVector()) return false;
	return true;
} // PowerFlower::ReadPowerSystem


#undef __FUNCT__
#define __FUNCT__ "PowerFlower::ReadPowerSystemFromMatPowerFile"
bool PowerFlower::ReadPowerSystemFromMatPowerFile(std::string& rSourceName)
{
	PowerSystem_Idema::MatPowerInterface interface;
	return ReadPowerSystem(rSourceName, interface);
} // PowerFlower::ReadPowerSystemFromMatPowerFile


#undef __FUNCT__
#define __FUNCT__ "PowerFlower::ReadPowerSystemFromMatPowerFile"
bool PowerFlower::ReadPowerSystemFromOptions(void)
{
	std::string fileName = "";
	PetscUtils::GetOptionsString(msPetscOptionFile, fileName, PETSC_MAX_PATH_LEN-1);
	if (fileName == "")
	{
		PetscErrorCode err = PetscInfo(PETSC_NULL, "[Failed] No input file specified. Use the command line option '-f filename'.\n"); CHKERRB(err);
		return false;
	}
	// Assume MatPower input file. No other input supported (yet).
	PowerSystem_Idema::MatPowerInterface interface;
	return ReadPowerSystem(fileName, interface);
} // PowerFlower::ReadPowerSystemFromMatPowerFile


#undef __FUNCT__
#define __FUNCT__ "PowerFlower::WritePowerSystemToString"
bool PowerFlower::WritePowerSystemToString(std::string& rString)
{
	if (mpPowerSystem == NULL) return false;
	rString = mpPowerSystem->ToString();
	return true;
} // PowerFlower::WritePowerSystemToString


#undef __FUNCT__
#define __FUNCT__ "PowerFlower::PowerFlow"
bool PowerFlower::PowerFlow(void) const
{
	msLogEvents.Begin(msLogEvents.GetPowerFlowEventId());
	if (!IsInitialized()) return false;

	PetscErrorCode err;
	if (mpPowerSystem == NULL)
	{
		err = PetscInfo(PETSC_NULL, "[Failed] No power system set.\n"); CHKERRB(err);
		msLogEvents.End(msLogEvents.GetPowerFlowEventId());
		return false;
	}

	// Set up nonlinear solver
	err = PetscInfo(PETSC_NULL, "Setting up nonlinear solver.\n"); CHKERRB(err);
	SolveResult result;
	NonLinearSolver solver;
	Preconditioner* pPreco = CreatePreconditionerFromOptions();	// :WARNING: needs to be called before solver is initialized, to catch options not supported by the nonlinear solver itself.
	if (!solver.Initialize())
	{
		err = PetscInfo(PETSC_NULL, "[Failed] Initializing nonlinear solver.\n"); CHKERRB(err);
		delete pPreco;
		msLogEvents.End(msLogEvents.GetPowerFlowEventId());
		return false;
	}

	// Set up power flow problem
	err = PetscInfo(PETSC_NULL, "Setting up power flow problem.\n"); CHKERRB(err);

	// :TODO: test if graph is fully connected?

	PowerFlowProblem problem(&msLogEvents);
	if (!problem.SetPowerSystem(*mpPowerSystem))
	{
		err = PetscInfo(PETSC_NULL, "[Failed] Setting up power flow problem.\n"); CHKERRB(err);
		delete pPreco;
		msLogEvents.End(msLogEvents.GetPowerFlowEventId());
		return false;
	}

	// :TODO: initial solution possibilities? Default: flat start.

	// Solve power flow problem
	err = PetscInfo(PETSC_NULL, "Solving power flow problem.\n"); CHKERRB(err);
	if (!SolvePowerFlowProblem(problem, NULL, &solver, pPreco, NULL, result))
	{
		err = PetscPrintf(PETSC_COMM_WORLD, " * [ERROR] Not solved.\n"); CHKERRB(err);
		delete pPreco;
		msLogEvents.End(msLogEvents.GetPowerFlowEventId());
		return false;
	}
	delete pPreco;

	// Print some data
	if (result.solveCode > 0)
	{
		err = PetscPrintf(PETSC_COMM_WORLD, " * converged:"); CHKERRB(err);
	}
	else
	{
		err = PetscPrintf(PETSC_COMM_WORLD, " * diverged :"); CHKERRB(err);
	}
	err = PetscPrintf(PETSC_COMM_WORLD, " %16.5e [%5.2f/%6.2f]\n", result.residualNorm, (double)result.nonlinearIterationCount, (double)result.linearIterationCount); CHKERRB(err);

	msLogEvents.End(msLogEvents.GetPowerFlowEventId());
	return true;
} // PowerFlower::PowerFlow


#undef __FUNCT__
#define __FUNCT__ "PowerFlower::ContingencyAnalysis"
bool PowerFlower::ContingencyAnalysis(void) const
{
	msLogEvents.Begin(msLogEvents.GetContingencyAnalysisEventId());
	if (!IsInitialized()) return false;

	PetscErrorCode err;
	if (mpPowerSystem == NULL)
	{
		err = PetscInfo(PETSC_NULL, "[Failed] No power system set.\n"); CHKERRB(err);
		msLogEvents.End(msLogEvents.GetContingencyAnalysisEventId());
		return false;
	}

	// Get contingency analysis options
	int initialSolutionType, preconditioningType;
	if (!GetContingencyAnalysisOptions(initialSolutionType, preconditioningType))
	{
		err = PetscInfo(PETSC_NULL, "[Failed] Setting contingency analysis options.\n"); CHKERRB(err);
		msLogEvents.End(msLogEvents.GetContingencyAnalysisEventId());
		return false;
	}
	if (initialSolutionType == 2)
	{
		err = PetscInfo(PETSC_NULL, "Setting contingency analysis to start from the solution of the base case.\n"); CHKERRB(err);
	}
	else
	{
		err = PetscInfo(PETSC_NULL, "Setting contingency analysis to use a flat start.\n"); CHKERRB(err);
	}

	// Get base case preconditioner from options
	std::auto_ptr<Preconditioner> pBasePreco(CreatePreconditionerFromOptions());	// :WARNING: CreatePreconditionerFromOptions() needs to be called before solver is initialized, to catch options not supported by the nonlinear solver itself.
	std::auto_ptr<PreconditionerFixed> pFixedPreco;		// Fixed preconditioner for all contingency cases (NULL pointer on creation)

	// Set up nonlinear solver
	err = PetscInfo(PETSC_NULL, "Setting up nonlinear solver.\n"); CHKERRB(err);
	SolveResult result;
	NonLinearSolver solver;
	if (!solver.Initialize())
	{
		err = PetscInfo(PETSC_NULL, "[Failed] Initializing nonlinear solver.\n"); CHKERRB(err);
		msLogEvents.End(msLogEvents.GetContingencyAnalysisEventId());
		return false;
	}

	/*** BASE CASE ***/

	err = PetscPrintf(PETSC_COMM_WORLD, "Base case power flow:\n"); CHKERRB(err);

	// Set up base case
	err = PetscInfo(PETSC_NULL, "Setting up base case.\n"); CHKERRB(err);

	// :TODO: test if graph is fully connected...

	// Remove inactive branches (because otherwise they would be activated when reactivating branches that were deactivated as contingency)
	err = PetscInfo(PETSC_NULL, "Removing inactive branches from the power flow problem.\n"); CHKERRB(err);
	PowerSystemExt::BranchIndexType removedCount = mpPowerSystem->RemoveInactiveBranches();
	err = PetscInfo1(PETSC_NULL, "Removed %d inactive branches.\n", removedCount); CHKERRB(err);

	PowerFlowProblem problem(&msLogEvents);
	if (!problem.SetPowerSystem(*mpPowerSystem))
	{
		err = PetscInfo(PETSC_NULL, "[Failed] Setting up power flow problem.\n"); CHKERRB(err);
		msLogEvents.End(msLogEvents.GetContingencyAnalysisEventId());
		return false;
	}

	// :TODO: Initial solution options for base case? Default is flat start.

	// Deal with contingency preconditioning that is not based on the solution of the base case
	if (preconditioningType == 0)			// Contingency cases use their own preco
	{
		err = PetscInfo(PETSC_NULL, "Setting contingency analysis to precondition each contingency seperately.\n"); CHKERRB(err);
	}
	else if (preconditioningType == 1)		// Contingency cases use preco based on base case in its intial solution
	{
		// :WARNING: Power flow problem objects needs to be done before setting fixedPreco
		err = PetscInfo(PETSC_NULL, "Setting contingency analysis to precondition all contingencies with the Jacobian of the base case in its initial solution.\n"); CHKERRB(err);
		pFixedPreco.reset(new PreconditionerFixed);
		pFixedPreco->SetMatrixCopy(problem.GetJacobianMatrix());
	}

	// Solve base case
	err = PetscInfo(PETSC_NULL, "Solving base case.\n"); CHKERRB(err);
	if (!SolvePowerFlowProblem(problem, NULL, &solver, pBasePreco.get(), NULL, result))
	{
		err = PetscPrintf(PETSC_COMM_WORLD, " * [ERROR] Not solved.\n"); CHKERRB(err);
		msLogEvents.End(msLogEvents.GetContingencyAnalysisEventId());
		return false;
	}

	// Print some data
	if (result.solveCode > 0)
	{
		err = PetscPrintf(PETSC_COMM_WORLD, " * converged:"); CHKERRB(err);
	}
	else
	{
		err = PetscPrintf(PETSC_COMM_WORLD, " * diverged :"); CHKERRB(err);
	}
	err = PetscPrintf(PETSC_COMM_WORLD, " %16.5e [%5.2f/%6.2f]\n", result.residualNorm, (double)result.nonlinearIterationCount, (double)result.linearIterationCount); CHKERRB(err);

	/*** CONTINGENCIES ***/

	// Set the initial solution based on the contingency analysis options
	std::auto_ptr<VoltageVector> pBaseSolution;			// Base solution (NULL pointer on creation)
	if (initialSolutionType == 2)
	{
		pBaseSolution.reset(new VoltageVector);
		*pBaseSolution = *problem.GetVoltageVector();	// Copy solution of the base case
	}

	// Deal with contingency preconditioning that is based on the solution of the base case
	if (preconditioningType == 2)		// Contingency cases use preco based on base case in its solution
	{
		err = PetscInfo(PETSC_NULL, "Setting contingency analysis to precondition all contingencies with the Jacobian of the base case in its solution.\n"); CHKERRB(err);
		pFixedPreco.reset(new PreconditionerFixed);
		pFixedPreco->SetMatrixCopy(problem.GetJacobianMatrix());
	}

	// Set preconditioner if it is fixed for all contingencies
	if (pFixedPreco.get() != NULL)
	{
		pBasePreco.reset();	// Done with base case preco, set to NULL so it can be passed to SolvePowerFlowProblem()
		if (!solver.SetPreconditioner(*pFixedPreco))		// Set fixed preco for all contingency cases
		{
			err = PetscInfo(PETSC_NULL, "[Failed] Setting preconditioner.\n"); CHKERRB(err);
			return false;
		}
	}



	// Set test focing terms
//	ForcingTermsTest testForcingTerms;



	/*** BRANCH CONTINGENCIES ***/

	err = PetscPrintf(PETSC_COMM_WORLD, "Branch contingencies:\n"); CHKERRB(err);

	std::vector<SolveResult> solveResults;

	// Get branch iterators for contingencies
	PowerSystemExt::ConstBranchIterator it = mpPowerSystem->GetBranchBegin();
	PowerSystemExt::ConstBranchIterator itEnd = mpPowerSystem->GetBranchEnd();
	PowerSystemExt::ConstBranchIterator itRollback, itNext;

	while (it != itEnd)
	{	// Branch contingency: disconnect all branches between the two busses that are
		// connected by the current branch, and solve the resulting power flow problem.
		err = PetscInfo(PETSC_NULL, "Setting up branch contingency case.\n"); CHKERRB(err);

		// Save the iterator of the current branch, so we can rollback the contingency disconnections later.
		itRollback = it;

		// Get bus indices of the current branch, such that bus1 < bus2.
		PowerSystemExt::BusIndexType bus1 = (it->GetFromBus() < it->GetToBus() ? it->GetFromBus() : it->GetToBus());
		PowerSystemExt::BusIndexType bus2 = (it->GetFromBus() < it->GetToBus() ? it->GetToBus() : it->GetFromBus());

		// Get the first iterator after all branches between the same busses as the current branch.
		itNext = mpPowerSystem->GetBranchNextUnique(it);

		// Disconnect all branches between between the two busses that are connected by the current branch.
		err = PetscInfo2(PETSC_NULL, "Disconnecting busses %d and %d.\n", bus1, bus2); CHKERRB(err);
		while (it != itNext)
		{
			mpPowerSystem->DeactivateBranch(it);
			++it;
		}

		// Test if disconnection did not cause two disconnected graphs with one of the two busses in each.
		bool isConnected = mpPowerSystem->IsConnected(bus1, bus2);
		if (isConnected)
		{
			err = PetscInfo(PETSC_NULL, "Solving branch contingency case.\n"); CHKERRB(err);
			if (!problem.UpdatePowerSystem(*mpPowerSystem))
			{
				err = PetscInfo(PETSC_NULL, "[Failed] Setting up power flow problem.\n"); CHKERRB(err);
				msLogEvents.End(msLogEvents.GetContingencyAnalysisEventId());
				return false;
			}

			// Solve contingency case
			if (!SolvePowerFlowProblem(problem, pBaseSolution.get(), &solver, pBasePreco.get(), NULL, result))
			{
				err = PetscInfo(PETSC_NULL, "[Failed] Solving power flow problem.\n"); CHKERRB(err);
				msLogEvents.End(msLogEvents.GetContingencyAnalysisEventId());
				return false;
			}
		}
		else
		{	// Graph no longer connected, can not solve this contingency case.
			err = PetscInfo2(PETSC_NULL, "Graph no longer connected. Skipping.\n", bus1, bus2); CHKERRB(err);

			// Set solve status to result info
			result.solveCode = 0;
			result.residualNorm = 0;
			result.nonlinearIterationCount = 0;
			result.linearIterationCount = 0;
		}

		// Store solve result
		solveResults.push_back(result);

		// Rollback the disconnections so the network is in original state again.
		it = itRollback;
		err = PetscInfo2(PETSC_NULL, "Reconnecting busses %d and %d.\n", bus1, bus2); CHKERRB(err);
		while (it != itNext)
		{
			mpPowerSystem->ActivateBranch(it);
			++it;
		}
	}

	// Print some data
	int skippedCount = 0;
	int convergedCount = 0;
	int convergedNonlinearIterationCount = 0;
	int convergedLinearIterationCount = 0;
	int divergedCount = 0;
	int divergedNonlinearIterationCount = 0;
	int divergedLinearIterationCount = 0;
	int contingencyCount = solveResults.size();

	std::vector<SolveResult>::const_iterator itResult = solveResults.begin();
	std::vector<SolveResult>::const_iterator itResultEnd = solveResults.end();
	while (itResult != itResultEnd)
	{
		if (itResult->solveCode > 0)
		{
			++convergedCount;
			convergedNonlinearIterationCount += itResult->nonlinearIterationCount;
			convergedLinearIterationCount += itResult->linearIterationCount;
		}
		else if (itResult->solveCode < 0)
		{
			++divergedCount;
			divergedNonlinearIterationCount += itResult->nonlinearIterationCount;
			divergedLinearIterationCount += itResult->linearIterationCount;
		}
		else
		{
			++skippedCount;
		}
		++itResult;
	}

	err = PetscPrintf(PETSC_COMM_WORLD, " * converged: %6d (%6.2f%) [%5.2f/%6.2f]\n", convergedCount, 100.0*convergedCount/contingencyCount, (convergedNonlinearIterationCount>0 ? (double)convergedNonlinearIterationCount/convergedCount : 0), (convergedNonlinearIterationCount>0 ? (double)convergedLinearIterationCount/convergedCount : 0)); CHKERRB(err);
	err = PetscPrintf(PETSC_COMM_WORLD, " * diverged : %6d (%6.2f%) [%5.2f/%6.2f]\n", divergedCount, 100.0*divergedCount/contingencyCount, (divergedNonlinearIterationCount>0 ? (double)divergedNonlinearIterationCount/divergedCount : 0), (divergedNonlinearIterationCount>0 ? (double)divergedLinearIterationCount/divergedCount : 0)); CHKERRB(err);
	err = PetscPrintf(PETSC_COMM_WORLD, " * skipped  : %6d (%6.2f%)\n", skippedCount, 100.0*skippedCount/contingencyCount); CHKERRB(err);

	msLogEvents.End(msLogEvents.GetContingencyAnalysisEventId());
	return true;
} // PowerFlower::ContingencyAnalysis


//============================= ACCESS     ===================================

/////////////////////////////// PRIVATE    ///////////////////////////////////


#undef __FUNCT__
#define __FUNCT__ "PowerFlower::GetContingencyAnalysisOptions"
bool PowerFlower::GetContingencyAnalysisOptions(int& rInitialSolutionType, int& rPrecondioningType) const
{
	PetscErrorCode err;
	std::string option;

	// Get the contingency analysis start type option string
	err = PetscUtils::GetOptionsString(msPetscOptionCAStartType, option, msPetscOptionMaxLength); CHKERRB(err);
	if (option == "")
	{
		rInitialSolutionType = 0;
	}
	else if (option == "flat")
	{
		rInitialSolutionType = 1;
	}
	else if (option == "basesolution")
	{
		rInitialSolutionType = 2;
	}
	else
	{
		rInitialSolutionType = 0;
		PetscError(PETSC_COMM_SELF,__LINE__,PETSC_FUNCTION_NAME,__FILE__,__SDIR__,PETSC_ERR_ARG_UNKNOWN_TYPE,PETSC_ERROR_INITIAL,"Unable to find requested contingency analysis start type %s",option.c_str());
		return false;
	}

	// Get the contingency analysis preconditioning type option string
	err = PetscUtils::GetOptionsString(msPetscOptionCAPrecoType, option, msPetscOptionMaxLength); CHKERRB(err);
	if (option == "")
	{
		rPrecondioningType = 0;
	}
	else if (option == "baseinitialjacobian")
	{
		rPrecondioningType = 1;
	}
	else if (option == "basesolutionjacobian")
	{
		rPrecondioningType = 2;
	}
	else
	{
		rPrecondioningType = 0;
		PetscError(PETSC_COMM_SELF,__LINE__,PETSC_FUNCTION_NAME,__FILE__,__SDIR__,PETSC_ERR_ARG_UNKNOWN_TYPE,PETSC_ERROR_INITIAL,"Unable to find requested contingency analysis preconditioning type %s",option.c_str());
		return false;
	}

	return true;
} // PowerFlower::GetContingencyAnalysisOptions


#undef __FUNCT__
#define __FUNCT__ "PowerFlower::CreatePreconditionerFromOptions"
Preconditioner* PowerFlower::CreatePreconditionerFromOptions(void) const
{
	// Get the preconditioner matrix type option string
	std::string option;
	PetscErrorCode err = PetscUtils::GetOptionsString(msPetscOptionPrecoMatrixType, option, msPetscOptionMaxLength);
	if (err != 0) return NULL;

	// Set up the preconditioner object
	Preconditioner* pPreco = NULL;
	if (option == "fdlf")
	{
		/* <begin> HAX THAT READS THE FDLF SPDM MATRIX FROM FILE */

		// Open file
		PetscViewer fd;
		std::string fdlfFile = "../data/fdlfspdm.petsc";
		err = PetscInfo1(PETSC_NULL, "Reading FDLF matrix from %s.\n", fdlfFile.c_str()); CHKERRB(err);
		err = PetscViewerBinaryOpen(PETSC_COMM_WORLD, fdlfFile.c_str(), FILE_MODE_READ, &fd); CHKERRB(err);

		// Create and read matrix
		Mat fdlf;
		err = MatCreate(PETSC_COMM_WORLD, &fdlf); CHKERRB(err);
		err = MatSetType(fdlf, MATSEQAIJ); CHKERRB(err);
		err = MatLoad(fdlf, fd); CHKERRB(err);

		// Close file
		err = PetscViewerDestroy(&fd); CHKERRB(err);

		// Create preco and copy FDLF matrix into the preco
		PreconditionerFixed* pPrecoFDLF = new PreconditionerFixed;
		if (!pPrecoFDLF->SetMatrixCopy(fdlf))
		{
			delete pPrecoFDLF;
			pPrecoFDLF = NULL;
		}
		err = MatDestroy(&fdlf); CHKERRB(err);
		pPreco = pPrecoFDLF;

		/* <end> HAX THAT READS THE FDLF SPDM MATRIX FROM FILE */

		// Remove the name-value pair from the options, so NonLinearSolver::CreatePreconditionerFromOptions
		// does not crash on the unknown value "fdlf", if it were to be called
		err = PetscOptionsClearValue(msPetscOptionPrecoMatrixType.c_str()); CHKERRB(err);
	}
	else
	{
		pPreco = NonLinearSolver::CreatePreconditionerFromOptions();
	}
	return pPreco;
} // PowerFlower::CreatePreconditionerFromOptions


#undef __FUNCT__
#define __FUNCT__ "PowerFlower::SolvePowerFlowProblem"
bool PowerFlower::SolvePowerFlowProblem(PowerFlowProblem& rProblem, VoltageVector* pInitialSolution, NonLinearSolver* pSolver, Preconditioner* pPreconditioner, ForcingTerms* pForcingTerms, SolveResult& rSolveResult) const
{
	msLogEvents.Begin(msLogEvents.GetPFSolveEventId());

	rSolveResult.solveCode = 0;		// Not solved (yet)
	rSolveResult.residualNorm = 0;
	rSolveResult.nonlinearIterationCount = 0;
	rSolveResult.linearIterationCount = 0;
	PetscErrorCode err;

	// Set up problem
	if (pInitialSolution && !rProblem.SetInitialSolution(*pInitialSolution))
	{	// Initial solution object provided, but failed setting it in the problem
		err = PetscInfo(PETSC_NULL, "[Failed] Setting initial solution.\n"); CHKERRB(err);
		msLogEvents.End(msLogEvents.GetPFSolveEventId());
		return false;
	}

	// Set up solver
	std::auto_ptr<NonLinearSolver> ownedSolver;
	if (!pSolver)
	{	// No solver provided, create solver object
		ownedSolver.reset(new NonLinearSolver);
		ownedSolver->Initialize();
		pSolver = ownedSolver.get();
	}
	if (pPreconditioner && !pSolver->SetPreconditioner(*pPreconditioner))
	{	// Preconditioner object provided, but failed setting it in the solver
		err = PetscInfo(PETSC_NULL, "[Failed] Setting preconditioner.\n"); CHKERRB(err);
		msLogEvents.End(msLogEvents.GetPFSolveEventId());
		return false;
	}
	if (pForcingTerms && !pSolver->SetForcingTerms(*pForcingTerms))
	{	// ForcingTerms object provided, but failed setting it in the solver
		err = PetscInfo(PETSC_NULL, "[Failed] Setting forcing terms.\n"); CHKERRB(err);
		msLogEvents.End(msLogEvents.GetPFSolveEventId());
		return false;
	}

	// Solve problem
	if (!pSolver->Solve(rProblem))
	{
		err = PetscInfo(PETSC_NULL, "[Failed] Solving power flow problem.\n"); CHKERRB(err);
		msLogEvents.End(msLogEvents.GetPFSolveEventId());
		return false;
	}

	// Set result info
	rSolveResult.solveCode = (pSolver->IsConverged() ? 1 : -1);
	rSolveResult.residualNorm = pSolver->GetFunctionNorm();
	rSolveResult.nonlinearIterationCount = pSolver->GetNonlinearIterationCount();
	rSolveResult.linearIterationCount = pSolver->GetLinearIterationCount();

	msLogEvents.End(msLogEvents.GetPFSolveEventId());
	return true;
} // PowerFlower::SolvePowerFlowProblem


POWERFLOWER_NAMESPACE_END
