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

#include "precondsg.h"

#include "utils.h"
#include "numerics.h"

#include "sysmat.h"
#include "sysmatsp.h"
#include "family.h"
#include "dirichlet.h"

#include "cmdpars.h"
extern CmdPars Cmd;


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


void DummyPreconditioner:: initialize(SystemMatrix* /*A*/, Vector<Num>& /*x*/, 
				      Vector<Num>& /*b*/)
{ 
    if (infoLinSystem) cout << " -- Preconditioner: none\n";
    // FORALL(x,i) if (dirichletBCs->isSet(i)) x[i] = dirichletBCs->value(i);
}

void DummyPreconditioner:: invert(Vector<Num>& lhs, SystemMatrix* /*A*/, 
				  Vector<Num>& rhs)
{
    int i;
    FORALL(lhs,i) lhs[i] = rhs[i]; 		// dummy procedure
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void Jacobi:: initialize(SystemMatrix* /*A*/, Vector<Num>& /*x*/, Vector<Num>& /*b*/)
{
    if (infoLinSystem) cout << " -- Jacobi Preconditioner\n";
}

void Jacobi:: invert(Vector<Num>& lhs, SystemMatrix* /*A*/, Vector<Num>& rhs)
{
    AP->DiagDiv(lhs,rhs);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void SGS:: initialize(SystemMatrix* /*A*/, Vector<Num>& /*x*/, Vector<Num>& /*b*/)
{
    if (infoLinSystem)  cout << " -- SGS Preconditioner\n";
}

void SGS:: invert(Vector<Num>& lhs, SystemMatrix* /*A*/, Vector<Num>& rhs)
{
    AP->Fm1(lhs,rhs);
    AP->DiagMult(lhs,lhs);
    AP->FmT(lhs,lhs);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

void SSOR:: initParameters()
{
    if (!Cmd.get("Omega",&omega)) omega = AP->omegaOpt(singleGrid);
}

void SSOR:: initialize(SystemMatrix* /*A*/, Vector<Num>& /*x*/, Vector<Num>& /*b*/)
{
    if (infoLinSystem) 
         cout << " -- SSOR Preconditioner (omega = " << omega <<")\n";
}

void SSOR:: invert(Vector<Num>& lhs, SystemMatrix* /*A*/, Vector<Num>& rhs)
{
    // Real factor = omega*(2.0-omega);

    AP->Fm1(lhs,rhs,omega);
    AP->DiagMult(lhs,lhs);              	// AP->DiagMult(lhs,lhs,factor);
    AP->FmT(lhs,lhs,omega);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


TrSSOR:: TrSSOR() : aux(1)   { }


void TrSSOR:: initParameters()
{
    if (!Cmd.get("Omega",&omega)) omega = AP->omegaOpt(singleGrid);
}

void TrSSOR:: initialize(SystemMatrix* A, Vector<Num>& x, Vector<Num>& b)
{
    int i;

    if (infoLinSystem) 
         cout << " -- TR-SSOR Preconditioner (omega = " << omega <<")\n";

    aux.resize(A->Dim());

    FORALL(aux,i) aux[i] = x[i];
    A->FT (x,aux,omega);
    A->Fm1(b,b,omega);
}

void TrSSOR:: close(SystemMatrix* A, Vector<Num>& x, Vector<Num>& b)
{
    int i;
    FORALL(aux,i) aux[i] = b[i];
    A->F(b,aux,omega);
    A->FmT(x,x,omega);
    aux.resize(1);
}


void TrSSOR:: invert(Vector<Num>& lhs, SystemMatrix* A, Vector<Num>& rhs)
{
    A->DiagMult(lhs,rhs);
}

void TrSSOR:: AMult(Vector<Num>& lhs, SystemMatrix* A, Vector<Num>& rhs)
{
    int i;
    A->FmT(aux,rhs,omega);

    A->DiagMult(lhs,aux,omega-2.0);
    FORALL(lhs,i) lhs[i] += rhs[i];

    A->Fm1(lhs,lhs,omega);
    FORALL(lhs,i) lhs[i] += aux[i]; 

    if (!equal(omega,1.0)) FORALL(lhs,i) lhs[i] /= omega;
}

void TrSSOR:: ATMult(Vector<Num>& lhs, SystemMatrix* A, Vector<Num>& rhs)
{
    int i;
    A->Fm1(aux,rhs,omega);

    A->DiagMult(lhs,aux);
    FORALL(lhs,i) lhs[i] = rhs[i] - lhs[i];

    A->FmT(lhs,lhs,omega);
    FORALL(lhs,i) lhs[i] += aux[i]; 
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


HB:: HB() 
{ 
    A0      = 0;
    AM1     = 0;
    AHBDiag = 0;
    maxLevel = -1;
}

HB:: ~HB()
{
    delete A0;
    delete AHBDiag;
    delete AM1;
}
//-------------------------------------------------------------------------


void HB:: update(SystemMatrix* A, FamilyTree* familyTree0,
		 DirichletBCs* dirichletBCs0)
{
    int i;

    dirichletBCs = dirichletBCs0;
    familyTree = familyTree0;
    ++maxLevel;

    delete AM1;
    AM1 = A;

    if (maxLevel == 0) 
    {
	if (A->DirectSolution())  { A0 = A;  AM1 = 0; }
	else
	{
	    A0  = 0;
	    AHBDiag = new Vector<Num>(1,A->Dim());
	    FORALL(*AHBDiag,i) (*AHBDiag)[i] = A->Diag(i);
	}
    }
    else if (maxLevel == 1 && A0)
    {
	AHBDiag = new Vector<Num>(A0->Dim()+1, A->Dim());
	FORALL(*AHBDiag,i) (*AHBDiag)[i] = A->Diag(i);
    }
    else
    {
	int prevDim = AHBDiag->high();
	AHBDiag->extendAndCopy(A->Dim());
	for (i=prevDim+1; i<=A->Dim(); ++i) (*AHBDiag)[i] = A->Diag(i);
    }
}
//-------------------------------------------------------------------------

void HB:: initialize(SystemMatrix* /*A*/, Vector<Num>& /*x*/, Vector<Num>& /*b*/)
{
    if (infoLinSystem) cout << " -- HB Preconditioner\n";
}
//-------------------------------------------------------------------------


void HB:: invert(Vector<Num>& lhs, SystemMatrix* /*A*/, Vector<Num>& rhs)
{
    int i, level;
    int max = familyTree->MaxLevel();
    int min = familyTree->MinLevel();

    FORALL(lhs,i) lhs[i] = rhs[i];

    for (level=max; level>=min; --level) 
    {
	familyTree->rhsToHB(lhs,level);
	FORALL(lhs,i) { if (dirichletBCs->isSet(i)) lhs[i] = 0.0; }
    }

    if (A0)  A0->FBSubst(lhs,lhs);
    FORALL(*AHBDiag,i) lhs[i] /= (*AHBDiag)[i];   
		
    for (level=min; level<=max; ++level) 
    {
	familyTree->solToNB(lhs,level);
	FORALL(lhs,i) { if (dirichletBCs->isSet(i)) lhs[i] = 0.0; }
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
/*
   ILU::  ILU() : AILU(0) 
   { 
   weight = 0.0;	Cmd.get("ILUWeight", &weight);
   }
   
   ILU:: ~ILU() { delete AILU; }
   
   void ILU:: initialize(SystemMatrix* A, Vector<Num>& x, Vector<Num>& b)
   {
   if (infoLinSystem) cout << " -- ILU Precond. (w=" << weight << ")\n";
   }
  
   void ILU:: update(SystemMatrix* APNew, FamilyTree* familyTree0,
   DirichletBCs* dirichletBCs0)  
   { 
   static Real accTime = 0.0;
   Timer timer, accTimer;
   
   dirichletBCs = dirichletBCs0;
   
   if (APNew != AP) 
   {
   delete AP;
   AP = APNew;
   }
   if (APNew->type() != SystemMatrix::Sparse)  
   { 
   cout << "\n*** No SparseMatrix handed to ILU Preconditioner\n";
   cout.flush(); abort();
   }
   // cannot be factorized
   delete AILU;
   AILU = new SparseMatrix(APNew->SpaceDim(), APNew);   // copy matrix to AILU
   
   AILU->ILUDecompose(weight);  
   
   // initParameters();
   
   
   if (timeLinSystem)    { cout << "\n\tILU-Factorization: "; timer.cpu(); }
   if (accTimeLinSystem) 
   {
   accTime += accTimer.cpu(false);
   cout << Form("\tAccumulated time:  %1.2f sec.\n", accTime); 
   }
   }
   
   void ILU:: invert(Vector<Num>& lhs, SystemMatrix* A, Vector<Num>& rhs)
   {
   AILU->ILUFBSubst(lhs,rhs);		// forward-backward-substitution
   }
   */
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


ILU:: ILU()
{ 
    dropTol = 1e-3;	Cmd.get("DropTolerance", &dropTol);
    
    if (equal(dropTol, 0.0))
    {
        cout << "\n*** ILU:: DropTolerance must be > 0 \n";
        cout.flush(); abort();
    }
}
//-------------------------------------------------------------------------

ILU:: ~ILU() { }

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

void ILU:: initialize(SystemMatrix* /*A*/, Vector<Num>& /*x*/, Vector<Num>& /*b*/)
{
    if (infoLinSystem) 
      cout << Form(" -- ILU-Preconditioner (dropTol=%1.3g)\n", dropTol);

    AP->ILUDecompose(dropTol);  
}
//-------------------------------------------------------------------------


void ILU:: invert(Vector<Num>& lhs, SystemMatrix* /*A*/, Vector<Num>& rhs)
{
    AP->ILUFBSubst(lhs,rhs);		// forward-backward-substitution
}