Skip to content

GitLab

Commit 0b5440f3 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

pfc preconditioner updated

parent 9f5bf49d
Hide whitespace changes
Inline Side-by-side
extensions/preconditioner/BlockPreconditioner.h
View file @ 0b5440f3
......@@ -57,6 +57,12 @@ namespace AMDiS {
TEST_EXIT(false)("Noch nicht implementiert!\n");
}
mtl::irange& getRows(size_t i)
{
return rows[i];
}
protected:
template<typename SolverType, typename RunnerType>
......
extensions/preconditioner/MTLPreconPfc.h
View file @ 0b5440f3
......@@ -19,39 +19,11 @@
#include "AMDiS.h"
#include "BlockPreconditioner.h"
// #include "TimeTracer.h"
using namespace std;
using namespace AMDiS;
template< typename Base >
struct PfcSchurMatrix
{
PfcSchurMatrix() : super(NULL) {};
PfcSchurMatrix(Base* super) : super(super)
{
m = num_rows(super->getM());
y1.change_dim(m);
y2.change_dim(m);
y3.change_dim(m);
}
template <typename VectorIn, typename VectorOut, typename Assign>
void mult(const VectorIn& b, VectorOut& x, Assign) const;
template <typename VectorIn>
mtl::vector::mat_cvec_multiplier<PfcSchurMatrix<Base>, VectorIn> operator*(const VectorIn& v) const
{ return mtl::vector::mat_cvec_multiplier<PfcSchurMatrix<Base>, VectorIn>(*this, v); }
int m;
private:
Base* super;
mutable MTLVector y1;
mutable MTLVector y2;
mutable MTLVector y3;
};
using namespace AMDiS::MTLTypes;
struct MTLPreconPfc : BlockPreconditioner
{
......@@ -70,26 +42,12 @@ struct MTLPreconPfc : BlockPreconditioner
MTLPreconPfc(std::string name)
: BlockPreconditioner(),
name(name),
tau(NULL),
solverM(NULL), solverM2(NULL), solverMpL(NULL), runnerM(NULL), runnerM2(NULL), runnerMpL(NULL),
Pid(NULL), Pdiag(NULL),
maxIterS(30), restartS(100)
{
createSubSolver(name + "->subsolver M", solverM, runnerM, "cg");
createSubSolver(name + "->subsolver M2", solverM2, runnerM2, "cg");
createSubSolver(name + "->subsolver MpL", solverMpL, runnerMpL, "cg");
Parameters::get(name + "->subsolver S->solver->max iteration", maxIterS);
Parameters::get(name + "->subsolver S->solver->restart", restartS);
}
tau(NULL)
{ }
virtual ~MTLPreconPfc()
{
delete solverM;
delete solverMpL;
if (Pid) { delete Pid; Pid = NULL; }
if (Pdiag) { delete Pdiag; Pdiag = NULL; }
exit();
}
void setData(double* tau_)
......@@ -100,13 +58,11 @@ struct MTLPreconPfc : BlockPreconditioner
virtual void init(const typename super::BlockMatrix& A, const MTLTypes::MTLMatrix& fullMatrix);
virtual void exit()
{
runnerM->exit();
runnerM2->exit();
runnerMpL->exit();
if (Pid) { delete Pid; Pid = NULL; }
if (Pdiag) { delete Pdiag; Pdiag = NULL; }
{
if (PId) { delete PId; PId = NULL; }
if (PDiagM) { delete PDiagM; PDiagM = NULL; }
if (PDiagMpL) { delete PDiagMpL; PDiagMpL = NULL; }
if (PDiagMpL2) { delete PDiagMpL2; PDiagMpL2 = NULL; }
}
virtual void solve(const MTLVector& b, MTLVector& x) const;
......@@ -114,79 +70,58 @@ struct MTLPreconPfc : BlockPreconditioner
template <typename VectorX, typename VectorB>
void solveM(VectorX& x, const VectorB& b, int nIter = -1) const
{
runnerM->solve(getM(), x, b);
// TimeTracer t("solveM");
itl::basic_iteration<double> iter(b, 5, 0, 1.e-6);
x = 0.0;
itl::cg(getM(), x, b, *PDiagM, *PId, iter);
// std::cout << "................ residual(M) = " << solverM->getResidual() << "\n";
}
template <typename VectorX, typename VectorB>
void solveM2(VectorX& x, const VectorB& b) const
void solveMpL(VectorX& x, const VectorB& b) const
{
runnerM2->solve(getM(), x, b);
// TimeTracer t("solveMpL");
itl::basic_iteration<double> iter(b, 20, 0, 1.e-6);
x = 0.0;
itl::cg(MpL, x, b, *PDiagMpL, *PId, iter);
// std::cout << "............. residual(MpK) = " << solverMpL->getResidual() << "\n";
}
template <typename VectorX, typename VectorB>
void solveMpK(VectorX& x, const VectorB& b) const
void solveMpL2(VectorX& x, const VectorB& b) const
{
runnerMpL->solve(MpL, x, b);
}
template <typename SchurMatrix, typename VectorX, typename VectorB>
void solveS(const SchurMatrix& S, VectorX& x, const VectorB& b) const
{
itl::basic_iteration<double> iter(b, maxIterS, 0, 0);
// TimeTracer t("solveMpL2");
itl::basic_iteration<double> iter(b, 10, 0, 1.e-6);
x = 0.0;
itl::fgmres_householder(S, x, b, *Pdiag, iter, restartS); // S*x = b
itl::cg(MpL2, x, b, *PDiagMpL2, *PId, iter);
// std::cout << "............. residual(MpK) = " << solverMpL->getResidual() << "\n";
}
const MTLMatrix& getM() const { return A->getSubMatrix(2,2); }
const MTLMatrix& getL() const { return A->getSubMatrix(2,1); }
double getTau() const { return *tau; }
protected:
MTLMatrix MpL;
PfcSchurMatrix<MTLPreconPfc> S;
MTLMatrix MpL2;
mutable MTLVector y0;
mutable MTLVector y1;
mutable MTLVector tmp;
LinearSolver* solverM;
LinearSolver* solverM2;
LinearSolver* solverMpL;
MTL4Runner<MTLTypes::MTLMatrix, MTLTypes::MTLVector>* runnerM;
MTL4Runner<MTLTypes::MTLMatrix, MTLTypes::MTLVector>* runnerM2;
MTL4Runner<MTLTypes::MTLMatrix, MTLTypes::MTLVector>* runnerMpL;
itl::pc::identity<MTLTypes::MTLMatrix, double> *Pid;
itl::pc::diagonal<MTLTypes::MTLMatrix, double> *Pdiag;
itl::pc::identity<MTLTypes::MTLMatrix, double> *PId;
itl::pc::diagonal<MTLTypes::MTLMatrix, double> *PDiagM;
itl::pc::diagonal<MTLTypes::MTLMatrix, double> *PDiagMpL;
itl::pc::diagonal<MTLTypes::MTLMatrix, double> *PDiagMpL2;
std::string name;
double* tau;
int maxIterS;
int restartS;
};
// some necessary functions that must be available for matrix-types
inline std::size_t size(const PfcSchurMatrix<MTLPreconPfc>& A) { return A.m * A.m; }
inline std::size_t num_rows(const PfcSchurMatrix<MTLPreconPfc>& A) { return A.m; }
inline std::size_t num_cols(const PfcSchurMatrix<MTLPreconPfc>& A) { return A.m; }
namespace mtl {
template <>
struct Collection<PfcSchurMatrix<MTLPreconPfc> >
{
typedef double value_type;
typedef int size_type;
};
namespace ashape {
template <> struct ashape_aux<PfcSchurMatrix<MTLPreconPfc> >
{ typedef nonscal type; };
}
}
#include "MTLPreconPfc.hh"
#endif // MTL_PFC_PRECONDITIONER_H
......
extensions/preconditioner/MTLPreconPfc.hh
View file @ 0b5440f3
......@@ -18,33 +18,53 @@
void MTLPreconPfc::init(const BlockMatrix& A_, const MTLMatrix& fullMatrix_)
{
assert(tau != NULL);
// TimeTracer t0("MTLPreconPfc::init");
super::init(A_, fullMatrix_);
// helper-matrix MpL = M + sqrt(tau)*L
MpL.change_dim(num_rows(getM()), num_cols(getM()));
MpL = getM() + sqrt(*tau)*getL();
typedef typename mtl::Collection<MTLMatrix>::size_type size_type;
typedef typename mtl::Collection<MTLMatrix>::value_type value_type;
size_type n = num_rows(getM());
double delta = sqrt(getTau());
double eps = sqrt(0.5);
// helper-matrix MpL = M + delta*L
MpL.change_dim(n, n);
MpL = getM() + delta*getL();
// helper-matrix MpL = M + eps*sqrt(delta)*L
MpL2.change_dim(n, n);
MpL2 = getM() + ( eps*sqrt(delta) ) * getL();
// schur-komplement-Matrix S = M - 2sqrt(tau)*L + sqrt(tau)*L*M^(-1)*L*M^(-1)*L
S = PfcSchurMatrix<MTLPreconPfc>(this);
PId = new itl::pc::identity<MTLTypes::MTLMatrix, double>(getM());
PDiagM = new itl::pc::diagonal<MTLTypes::MTLMatrix, double>(getM());
PDiagMpL = new itl::pc::diagonal<MTLTypes::MTLMatrix, double>(MpL);
PDiagMpL2 = new itl::pc::diagonal<MTLTypes::MTLMatrix, double>(MpL2);
// {
// mtl::io::matrix_market_ostream outM("M.mtx");
// mtl::io::matrix_market_ostream outL("L.mtx");
// mtl::io::matrix_market_ostream outMpL("MpL.mtx");
// mtl::io::matrix_market_ostream outMpL2("MpL2.mtx");
// outM << getM();
// outL << getL();
// outMpL << MpL;
// outMpL2 << MpL2;
// }
// temporary variables
y0.change_dim(num_rows(getM()));
y1.change_dim(num_rows(getM()));
tmp.change_dim(num_rows(getM()));
// init sub-solvers
Pid = new itl::pc::identity<MTLTypes::MTLMatrix, double>(getM());
Pdiag = new itl::pc::diagonal<MTLTypes::MTLMatrix, double>(MpL);
runnerM->init(A_, getM());
runnerM2->init(A_, getM());
runnerMpL->init(A_, MpL);
tmp.change_dim(num_rows(getM()));
}
void MTLPreconPfc::solve(const MTLVector& b, MTLVector& x) const
{ FUNCNAME("MTLPreconPfc::solve()");
// TimeTracer t0("PreconPfc::solve");
x.change_dim(num_rows(b));
const MTLVector b0(b[rows[0]]);
......@@ -57,41 +77,20 @@ void MTLPreconPfc::solve(const MTLVector& b, MTLVector& x) const
double delta = sqrt(getTau());
solveM(y0, b0);
solveM(y0, b0); // M*y0 = b0
y1 = getL() * y0; // y1 = K*y0
tmp = b1 - getTau()*y1; // tmp := b1 - tau*y1
solveMpK(y1, tmp);
x0 = y0 + (1.0/delta)*y1;// x0 = y0 + (1/delta)*y1
solveMpL(y1, tmp); // (M + delta*K) * y1 = tmp
x0 = y0 + (1.0/delta)*y1; // x0 = y0 + (1/delta)*y1
tmp = getM() * y1; // tmp := M*y1
solveS(S, x1, tmp);
solveMpL2(y1, tmp); // (M+eps*sqrt(delta)K) * y1 = tmp
tmp = getM() * y1; // tmp := M*y1
solveMpL2(x1, tmp); // (M+eps*sqrt(delta)K) * x1 = tmp
x0-= (1.0/delta)*x1; // x0 = x0 - (1/delta)*x1 = y0 + (1/delta)*(y1 - x1)
x0-= (1.0/delta)*x1; // x0 = x0 - (1/delta)*x1 = y0 + (1/delta)*(y1 - x1)
tmp = b2 - getL() * x1; // tmp := b2 - K*x1
solveM2(x2, tmp);
}
template<>
template <typename VectorIn, typename VectorOut, typename Assign>
void PfcSchurMatrix<MTLPreconPfc>::mult(const VectorIn& b, VectorOut& x, Assign) const
{ FUNCNAME("PfcSchurMatrix<MTLPreconPfc>::mult()");
int m = num_rows(super->getM());
assert(int(size(b)) == m);
assert(size(b) == size(x));
double delta = sqrt(super->getTau());
y1 = super->getL() * b;
super->solveM(y2, y1);
y3 = super->getL() * y2;
y2 = super->getM() * b;
y2-= (2.0*delta) * y1;
y2+= delta * y3;
Assign::apply(x, y2);
solveM(x2, tmp);
}
extensions/preconditioner/PetscPreconPfc.cc
View file @ 0b5440f3
......@@ -15,41 +15,16 @@
*
******************************************************************************/
#include "PetscPreconPfc.h"
// #include "TimeTracer.h"
namespace AMDiS {
using namespace std;
int pcPfcMatMultSchur(Mat mat, Vec b, Vec x)
{
void *ctx;
MatShellGetContext(mat, &ctx);
PfcData* data = static_cast<PfcData*>(ctx);
Vec y1, y2, y3;
VecDuplicate(b, &y1);
VecDuplicate(b, &y2);
VecDuplicate(b, &y3);
MatMult(data->matK, b, y1); // y1 := K*b
KSPSolve(data->kspM, y1, y2); // M*y2 = y1
MatMult(data->matK, y2, y3); // y3 := K*y2
// KSPSolve(data->kspM, y3, y2); // M*y2 = y3
// MatMult(data->matK, y2, y3); // y3 := K*y2
VecScale(y3, data->delta); // y3 = delta*y3
MatMultAdd(data->matM, b, y3, x); // x = M*b + y3
VecAXPY(x, -2.0*(data->delta), y1); // x := M*b - 2*delta*K*b + delta*B*b
VecDestroy(&y1);
VecDestroy(&y2);
VecDestroy(&y3);
return 0;
}
/// solve Pfc Preconditioner
PetscErrorCode pcPfcShell(PC pc, Vec b, Vec x) // solve Px=b
{ FUNCNAME("pcPfcShell()");
void *ctx;
PCShellGetContext(pc, &ctx);
PfcData* data = static_cast<PfcData*>(ctx);
......@@ -75,64 +50,62 @@ namespace AMDiS {
VecAYPX(tmp, -(data->tau), b2); // tmp := b2 - tau*tmp
KSPSolve(data->kspMpK, tmp, y2); // (M+sqrt(tau)K)*y2 = tmp
MatMult(data->matM, y2, tmp); // tmp := M*y2
KSPSolve(data->kspMpK2, tmp, x2); // MpK2*x2 = tmp
MatMult(data->matM, x2, tmp); // tmp := M*x2
KSPSolve(data->kspMpK2, tmp, x2); // MpK2*x2 = tmp
KSPSolve(data->kspSchur, tmp, x2); // S*x2 = tmp
VecCopy(x2, x1); // x1 := x2
VecAXPBYPCZ(x1, 1.0, 1.0/(data->delta), -1.0/(data->delta), y1, y2); // x1 = 1*y1 + factor*y2 - factor*x1
MatMult(data->matK, x2, tmp); // tmp := K*x2
VecAYPX(tmp, -1.0, b3); // tmp := b3 - tmp
KSPSolve(data->kspM2, tmp, x3); // M*x3 = tmp
KSPSolve(data->kspM, tmp, x3); // M*x3 = tmp
VecDestroy(&y1);
VecDestroy(&y2);
VecDestroy(&tmp);
return 0;
}
void PetscPreconPfc::init(const BlockMatrix& A, const Mat& matrix)
{
assert(tau != NULL);
dynamic_cast<PetscSolver<PetscPreconPfc>*>(&oem)->setNestedVectors(true);
PetscRunner::init(A, matrix);
// init shell preconditioner
PCSetType(getPc(), PCSHELL);
PCShellSetApply(getPc(), pcPfcShell);
PCShellSetContext(getPc(), &data);
double delta = sqrt(*tau);
double eps = sqrt(0.5);
MatNestGetSubMat(matrix, 2, 2, &data.matM);
MatNestGetSubMat(matrix, 2, 1, &data.matK);
MatDuplicate(data.matM, MAT_COPY_VALUES, &MpK);
MatAXPY(MpK, sqrt(*tau), data.matK, SAME_NONZERO_PATTERN);
MatAXPY(MpK, delta, data.matK, SAME_NONZERO_PATTERN);
// schur-komplete-Matrix
PetscInt numRows, numCols;
MatGetSize(data.matM, &numRows, &numCols);
MatCreateShell(PETSC_COMM_SELF, numRows, numCols,
PETSC_DETERMINE, PETSC_DETERMINE,
&data, &matSchur);
MatShellSetOperation(matSchur, MATOP_MULT, (void(*)(void))pcPfcMatMultSchur);
KSPCreate(PETSC_COMM_SELF, &data.kspSchur);
KSPSetOperators(data.kspSchur, matSchur, MpK, SAME_NONZERO_PATTERN);
setSolver(data.kspSchur, "schur_", KSPFGMRES, PCJACOBI, 1e-6, 1e-8, 4);
KSPSetFromOptions(data.kspSchur);
MatDuplicate(data.matM, MAT_COPY_VALUES, &MpK2);
MatAXPY(MpK2, eps*sqrt(delta), data.matK, SAME_NONZERO_PATTERN);
// init sub-solvers
createSubSolver(data.kspM, data.matM, "mass_");
createSubSolver(data.kspM2, data.matM, "mass2_");
createSubSolver(data.kspMpK, MpK, "MpK_");
createSubSolver(data.kspMpK2, MpK2, "MpK2_");
data.delta = sqrt(*tau);
data.tau = (*tau);
data.delta = delta;
data.tau = *tau;
setSolver(data.kspM, "mass_", KSPCG, PCJACOBI, 0.0, 1e-14, 5);
setSolver(data.kspM2, "mass2_", KSPCG, PCJACOBI, 0.0, 1e-14, 15);
setSolver(data.kspMpK, "MpK_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 5);
setSolver(data.kspM, "mass_", KSPCG, PCJACOBI, 0.0, 1e-8, 5);
setSolver(data.kspMpK, "MpK_", KSPCG, PCJACOBI, 0.0, 1e-6, 20);
setSolver(data.kspMpK2, "MpK2_", KSPCG, PCJACOBI, 0.0, 1e-6, 10);
}
......@@ -141,12 +114,11 @@ namespace AMDiS {
FUNCNAME("PetscPreconPfc::exit()");
MatDestroy(&MpK);
MatDestroy(&matSchur);
MatDestroy(&MpK2);
KSPDestroy(&data.kspM);
KSPDestroy(&data.kspM2);
KSPDestroy(&data.kspMpK);
KSPDestroy(&data.kspSchur);
KSPDestroy(&data.kspMpK2);
PetscRunner::exit();
}
......
extensions/preconditioner/PetscPreconPfc.h
View file @ 0b5440f3
......@@ -25,7 +25,7 @@ namespace AMDiS {
using namespace std;
struct PfcData {
KSP kspM, kspM2, kspMpK, kspSchur;
KSP kspM, kspMpK, kspMpK2;
Mat matM, matK;
double delta, tau;
};
......@@ -48,7 +48,7 @@ namespace AMDiS {
private:
Mat MpK;
Mat matSchur;
Mat MpK2;
PfcData data;
double *tau;
......
extensions/preconditioner/PetscSolverPfc.cc
View file @ 0b5440f3
......@@ -22,33 +22,6 @@
namespace AMDiS { namespace Parallel {
using namespace std;
PetscErrorCode pcPfcMatMultSchur(Mat mat, Vec b, Vec x)
{
void *ctx;
MatShellGetContext(mat, &ctx);
PfcData* data = static_cast<PfcData*>(ctx);
Vec y1, y2, y3;
VecDuplicate(b, &y1);
VecDuplicate(b, &y2);
VecDuplicate(b, &y3);
MatMult(data->matK, b, y1); // y1 := K*b
KSPSolve(data->kspM, y1, y2); // M*y2 = y1
MatMult(data->matK, y2, y3); // y3 := K*y2
// KSPSolve(data->kspM, y3, y2); // M*y2 = y3
// MatMult(data->matK, y2, y3); // y3 := K*y2
VecScale(y3, data->delta); // y3 = delta*y3
MatMultAdd(data->matM, b, y3, x); // x = M*b + y3
VecAXPY(x, -2.0*(data->delta), y1); // x := M*b - 2*delta*K*b + delta*B*b
VecDestroy(&y1);
VecDestroy(&y2);
VecDestroy(&y3);
PetscFunctionReturn(0);
}
/// solve Pfc Preconditioner
PetscErrorCode pcPfcShell(PC pc, Vec b, Vec x) // solve Px=b
......@@ -79,21 +52,25 @@ namespace AMDiS { namespace Parallel {
VecAYPX(tmp, -(data->tau), b2); // tmp := b2 - tau*tmp
KSPSolve(data->kspMpK, tmp, y2); // (M+sqrt(tau)K)*y2 = tmp
MatMult(data->matM, y2, tmp); // tmp := M*y2
KSPSolve(data->kspMpK2, tmp, x2); // MpK2*x2 = tmp
MatMult(data->matM, x2, tmp); // tmp := M*x2
KSPSolve(data->kspMpK2, tmp, x2); // MpK2*x2 = tmp
KSPSolve(data->kspSchur, tmp, x2); // S*x2 = tmp
VecCopy(x2, x1); // x1 := x2
VecAXPBYPCZ(x1, 1.0, 1.0/(data->delta), -1.0/(data->delta), y1, y2); // x1 = 1*y1 + factor*y2 - factor*x1
MatMult(data->matK, x2, tmp); // tmp := K*x2
VecAYPX(tmp, -1.0, b3); // tmp := b3 - tmp
KSPSolve(data->kspM2, tmp, x3); // M*x3 = tmp
KSPSolve(data->kspM, tmp, x3); // M*x3 = tmp
VecDestroy(&y1);
VecDestroy(&y2);
VecDestroy(&tmp);
PetscFunctionReturn(0);
return 0;
}
void PetscSolverPfc::initSolver(KSP &ksp)
......@@ -124,6 +101,9 @@ namespace AMDiS { namespace Parallel {
PCShellSetApply(getPc(), pcPfcShell);
PCShellSetContext(getPc(), &data);
double delta = sqrt(*tau);
double eps = sqrt(0.5);
const FiniteElemSpace *feSpace = componentSpaces[0];
// create mass-matrix
......@@ -138,14 +118,10 @@ namespace AMDiS { namespace Parallel {
data.matM = solverM->getMatInterior();
data.kspM = solverM->getSolver();
solverM2 = createSubSolver(0, "M2_");
solverM2->fillPetscMatrix(&matrixM);
data.kspM2 = solverM2->getSolver();
// create laplace-matrix
// create MpK-matrix
DOFMatrix matrixMpK(feSpace, feSpace);
Operator laplaceOp2(feSpace, feSpace);
Simple_SOT sot2(sqrt(*tau));
Simple_SOT sot2(delta);
laplaceOp2.addTerm(&zot);
laplaceOp2.addTerm(&sot2);
matrixMpK.assembleOperator(laplaceOp2);
......@@ -155,6 +131,20 @@ namespace AMDiS { namespace Parallel {
matMpK = solverMpK->getMatInterior();
data.kspMpK = solverMpK->getSolver();
// create MpK2-matrix
DOFMatrix matrixMpK2(feSpace, feSpace);
Operator laplaceOp3(feSpace, feSpace);
Simple_SOT sot3(eps*sqrt(delta));
laplaceOp3.addTerm(&zot);
laplaceOp3.addTerm(&sot3);
matrixMpK2.assembleOperator(laplaceOp3);
solverMpK2 = createSubSolver(0, "MpK2_");
solverMpK2->fillPetscMatrix(&matrixMpK2);
matMpK2 = solverMpK2->