/*
$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];
}
}
*/
syntax highlighted by Code2HTML, v. 0.9.1