/*
 $Id: linsystem.cc,v 1.3 1996/10/11 09:09:57 roitzsch Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "linsystem.h"

#include "sysmat.h"
#include "connect.h"
#include "dirichlet.h"

#include "cmdpars.h"
extern CmdPars Cmd;

//-------------------------------------------------------------------------

LinSystem:: LinSystem(int spaceDim0, symmetryType sym, 
		      Preconditioner* 	    preCond0,
		      const DirichletBCs*   dirichletBCs0)

	: symmetry(sym), spaceDim(spaceDim0), 
	  A(0), b(1), bSave(1), preCond(preCond0), dirichletBCs(dirichletBCs0),
	  statistic(0), status(True)
{ 

    extPrecFactor = 1.0;  Cmd.get("extPrecFactor",&extPrecFactor);
    maxOrthoGMRes = 15;   Cmd.get("maxOrthoGMRes",&maxOrthoGMRes);


    infoLinSystem = 0;  	Cmd.get("infoLinSystem", &infoLinSystem);
    timeLinSystem = 0;  	Cmd.get("timeLinSystem", &timeLinSystem);
    accTimeLinSystem = 0;  	Cmd.get("accTimeLinSystem", &accTimeLinSystem);

    level0DirectSolverLimit = 25000; 
    directSolverLimit 	    = 0; 

    Cmd.get("level0DirectSolver",&level0DirectSolverLimit);
    Cmd.get("directSolver",      &directSolverLimit);

    solver = cgOmin;
    if      (Cmd.isSet("linSolver", "PCG"))      solver = cgOmin;
    else if (Cmd.isSet("linSolver", "CGOmin"))   solver = cgOmin;
    else if (Cmd.isSet("linSolver", "CGOdir"))   solver = cgOdir;
    else if (Cmd.isSet("linSolver", "CROmin"))   solver = crOmin;
    else if (Cmd.isSet("linSolver", "CRODir"))   solver = crOdir;
    else if (Cmd.isSet("linSolver", "LSQCG"))    solver = lsqCG;
    else if (Cmd.isSet("linSolver", "DdCG"))     solver = ddCG;
    else if (Cmd.isSet("linSolver", "StdBiCG"))  solver = biCG;
    else if (Cmd.isSet("linSolver", "BiCGStab")) solver = biCGStab;
    else if (Cmd.isSet("linSolver", "CGS"))  	 solver = cgs;
    else if (Cmd.isSet("linSolver", "GMRes"))    solver = gmRes;
    else if (Cmd.isSet("linSolver", "Relax"))    solver = relax;
    else if (Cmd.isSet("linSolver", "nonLinRelax")) solver = nonLinRelax;
    else if (Cmd.isSet("linSolver", "Test1"))    solver = testSolver1;
    else if (Cmd.isSet("linSolver", "Test2"))    solver = testSolver2;


    convTest0 = residual;
    if      (Cmd.isSet("convergenceTest","ccgDd"))    convTest0 = ccgDd;
    else if (Cmd.isSet("convergenceTest","ccgDB"))    convTest0 = ccgDB;
    else if (Cmd.isSet("convergenceTest","ci"))       convTest0 = ci;
    else if (Cmd.isSet("convergenceTest","residual")) convTest0 = residual;
    else if (Cmd.isSet("convergenceTest","decay")) convTest0 = decayOfResidual;
    else if (Cmd.isSet("convergenceTest","vectorIter")) 
    						   convTest0 = vectorIteration;
}
//-------------------------------------------------------------------------

LinSystem:: ~LinSystem() { }

//-------------------------------------------------------------------------

void LinSystem:: notImplemented(const char* /*s*/)
{
    cout << "\n*** class LinSystem -- " << " not implemented\n";
    cout.flush(); abort();
}
//-------------------------------------------------------------------------

void LinSystem:: formNewEigenSystem(Num /*lambda*/)
{
    notImplemented("formNewEigenSystem");
}
//-------------------------------------------------------------------------

int LinSystem:: Dim() const { return A->Dim(); }

//-------------------------------------------------------------------------

void LinSystem:: update(SystemMatrix* ANew)
{
    A = ANew;				// old matrix handled by preconditioner
    b.resize(A->Dim());			// right hand side vector
}
//-------------------------------------------------------------------------

void LinSystem:: reset()
{
    int i;
    A->reset(); 
    FORALL(b,i) b[i] = 0.0;
}
//-------------------------------------------------------------------------

void LinSystem:: setRhs(Vector<Num>& rhs)
{
    int i;
    if (b.h != rhs.h) b.resize(rhs.h);
    FORALL(b,i) b[i] = rhs[i];
}
//-------------------------------------------------------------------------

// --			res = b - A*x

void LinSystem:: Residual(Vector<Num>& res, Vector<Num>& x) const
{
    int i;
    A->Mult(res,x);
    FORALL(b,i) res[i] = b[i] - res[i];
}
//-------------------------------------------------------------------------

Real LinSystem:: normRhs() const 
{
    int i;
    static Real tiny = machMin(Real(0.0));

    Num sum = 0.0;
    FORALL(b,i)  if (!dirichletBCs->isSet(i)) sum += conj(b[i])*b[i];

    if (sum <= tiny) 
    {
	Vector<Num> r(b.high()), x(b.high());
	FORALL(x,i)  
	{
	    if (dirichletBCs->isSet(i)) x[i] = dirichletBCs->value(i);
	    else x[i] = 0.0;
	}
	Residual(r,x);
	sum = 0.0;
	FORALL(r,i) sum += conj(r[i])*r[i];

	cout << Form("\n\t* norm of rhs == 0; set to |r0| = %1.4g\n", 
		       sqrt(Abs(sum)));

	if (sum <= tiny) 
	{
	    sum = 0.0;
	    FORALL(b,i) sum += conj(b[i])*b[i];
	    cout << Form("\n\t* norm of r0 == 0; set to |b| = %1.4g \n",
			   sqrt(Abs(sum)));
	}
    }
    
    return sqrt(Abs(sum));
}
//-------------------------------------------------------------------------


void LinSystem:: store(const Matrix<Num>& AElem, const Vector<Num>& bElem, 
		       const Vector<int>& globalNodes)
{
    int i, row;

    FORALL(globalNodes,i) if ((row=globalNodes[i])) b[row] += bElem[i];

    A->store(AElem, globalNodes);
}
//-------------------------------------------------------------------------


void LinSystem:: setDirichletBCs(DirichletBCs& dirichletBCs0)
{ 
    int i;

    bSave.resize(1,b.high());
    FORALL(bSave,i) bSave[i] = b[i];   // save original rhs vector

    A->setDirichletBCs(dirichletBCs0, b, bSave, Fct0, E0); 
}
//-------------------------------------------------------------------------

void LinSystem:: memSpaceInfo()
{
    A->memSpace(1);
}
//-------------------------------------------------------------------------

void LinSystem:: print()
{
    cout << "\nLinear System\n";
    A->print();
    cout << "\n right hand side:\n"; 
    b.print();
    cout << "\n";
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

/*

SparseLinSystem:: SparseLinSystem (int spaceDim, symmetryType sym,
				   Preconditioner* 	 preCond,
				   const DirichletBCs*   dirichletBCs)

	: LinSystem(spaceDim, sym, preCond, dirichletBCs)
{  } 
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

SparseEigenSystem:: SparseEigenSystem (int spaceDim, symmetryType sym,
				       const Interface*     interface,
				       const DirichletBCs*  dirichletBCs)

    : SparseLinSystem(spaceDim, sym, interface, dirichletBCs, solver, convTest)
{  
    AA = BB = 0;
}
//-------------------------------------------------------------------------

SparseEigenSystem:: ~SparseEigenSystem()
{
    delete AA;
    delete BB;
}
//-------------------------------------------------------------------------

void SparseEigenSystem:: reset()
{
    AA->reset();
    BB->reset();
}
//-------------------------------------------------------------------------


void SparseEigenSystem:: update()
{
    int i;

    refLevel  = interface->refLevel;
    noOfNodes = interface->noOfNodes[refLevel];

    if (Cmd.isSet("PreMat","AA")) 
    {
	preCond->dispose(AA);   
	delete BB; 
	delete A;
    }
    else if (Cmd.isSet("PreMat","BB"))
    { 
	preCond->dispose(BB);   
	delete AA; 
	delete A;
    }
    else 	
    { 
	preCond->dispose(A); 
	delete AA; 
	delete BB;
    }

    A = new SparseMatrix(*interface->cPattern, symmetry);
    A->reset();
    AA = new SparseMatrix(A);
    BB = new SparseMatrix(A);
}
//-------------------------------------------------------------------------


void SparseEigenSystem:: updateAssembledSystem()
{
    if (Cmd.isSet("PreMat","AA")) 	preCond->update(AA, interface);
    else if (Cmd.isSet("PreMat","BB")) 	preCond->update(BB, interface);
    else 				preCond->update(A,  interface);
}
//-------------------------------------------------------------------------


void SparseEigenSystem:: formNewEigenSystem(Num lambda)
{
    int i, j, k;

    FORALL(A->D,i)  A->D[i] = AA->D[i] - lambda*BB->D[i];

    FORALL(A->D,i) 
	    for (k=A->end[i-1]+1; k<=A->end[i]; ++k)
	    {
		(*A->L)[k] = (*AA->L)[k] - lambda * (*BB->L)[k];

		if (symmetry == asym) 
		(*A->U)[k] = (*AA->U)[k] - lambda * (*BB->U)[k];
	    }
}
*/