Commits (30)
......@@ -10,7 +10,7 @@ from .utilities.test_utilities import print_test_results
class Mpp:
def __init__(self, project_name='M++', executable='M++',
kernels=4, mute=False):
kernels=4, mute=False, fail=False):
self.project_name = project_name
self.executable = executable
......@@ -189,7 +189,6 @@ class Mpp:
with open(self.dm.mpi_tests_file_path, 'r') as file:
for index, line in enumerate(file.readlines()):
# remove /, \n
print(self.cwd)
absolute_path = line[1:-1]
self.cwd = "/" + absolute_path[:absolute_path.rfind('/')+1]
self.executable = absolute_path[absolute_path.rfind('/') + 1:]
......
class test_colors:
OK = '\033[92m' #GREEN
WARNING = '\033[93m' #YELLOW
FAIL = '\033[91m' #RED
RESET = '\033[0m' #RESET COLOR
def prRed(text):
print(test_colors.FAIL + text + test_colors.RESET)
def prGreen(text):
print(test_colors.OK + text + test_colors.RESET)
def print_test_results(duration, index, name, rc, total_num_of_tests):
if rc == 0:
message = f'{index + 1:>3}/{total_num_of_tests:<3}' \
f' Test {index + 1:>4}: {name} '
f' Test {index + 1:>3}: {name} '
message = message.ljust(64, '.')
message = message + ' Passed {:3.2f} sec'.format(duration)
print(message)
prGreen(message)
else:
message = f'{index + 1:>3}/{total_num_of_tests:<3}' \
f' Test {index + 1:>4}: {name} '
f' Test {index + 1:>3}: {name} '
message = message.ljust(64, '.')
message = message + ' Failed {:3.2f} sec'.format(duration)
print(message)
\ No newline at end of file
prRed(message)
\ No newline at end of file
......@@ -32,6 +32,8 @@ if __name__ == '__main__':
help='runs all registered tests in ctest')
parser.add_argument('--bench_tests', type=int, default=0,
help='runs all registered tests in mpp_bench')
parser.add_argument('-f', "--fail_state", type=int, default=1,
help='forwards the test fail to system')
args = parser.parse_args()
......@@ -42,13 +44,13 @@ if __name__ == '__main__':
project_name = line[line.rfind('=') + 1:line.rfind('\n')]
mpp = Mpp(project_name, executable=args.executable,
kernels=args.kernels, mute=bool(args.mute))
kernels=args.kernels, mute=bool(args.mute), fail=args.fail_state)
rc = 0
if args.mpi_tests:
mpp.run_mpi_tests(args.kernels)
if args.mpp_tests:
mpp.run_mpp_tests()
if args.bench_tests:
mpp.run_bench_tests(args.kernels)
rc = mpp.run_mpi_tests(args.kernels)
elif args.mpp_tests:
rc = mpp.run_mpp_tests()
elif args.bench_tests:
rc = mpp.run_bench_tests(args.kernels)
sys.exit(0 if rc == 0 else 1)
......@@ -6,7 +6,7 @@
class Distribution {
int verbose;
int verbose = 0;
int commSplit = 0;
......@@ -29,18 +29,16 @@ class Distribution {
int setProcsForCellsAndReturnMid(int rekDepth, int begin, int end);
public:
Distribution(int commSplit = 0, int verbose = 0) :
verbose(verbose),
commSplit(commSplit),
markedCells(PPM->Size(commSplit)) {
Distribution(int commSplit = 0) :
markedCells(PPM->Size(commSplit)),
commSplit(commSplit) {
config.get("DistributionVerbose", verbose);
config.get("Distribution", _distName);
}
Distribution(const std::string &distName, int commSplit = 0, int verbose = 0) :
verbose(verbose),
commSplit(commSplit),
Distribution(const std::string &distName, int commSplit = 0) :
markedCells(PPM->Size(commSplit)),
commSplit(commSplit),
_distName(distName) {
config.get("DistributionVerbose", verbose);
}
......
......@@ -6,7 +6,7 @@
class Overlap {
int verbose;
int verbose = 0;
int commSplit = 0;
......@@ -23,17 +23,15 @@ class Overlap {
void overlapLayer(Mesh &M, Vertices &bndVertices);
public:
Overlap(int commSplit = 0, int verbose = 0) :
verbose(verbose),
Overlap(int commSplit = 0) :
commSplit(commSplit) {
config.get("OverlapVerbose", verbose);
config.get("Overlap", _overlapName);
}
Overlap(const std::string &overlapName, int commSplit = 0, int verbose = 1) :
verbose(verbose),
commSplit(commSplit),
_overlapName(overlapName) {
Overlap(const std::string &overlapName, int commSplit = 0) :
_overlapName(overlapName),
commSplit(commSplit) {
config.get("OverlapVerbose", verbose);
}
......@@ -60,6 +58,7 @@ public:
void SlipFacesOverlap_cdd2d(Mesh &mesh);
};
// Todo remove as soon as spacetime is integrated
static void CellsOverlap_dG1(Mesh &M) {
std::map<Point, bool> clean_list;
......@@ -135,5 +134,4 @@ static void CellsOverlap_dG1(Mesh &M) {
}
}
#endif //OVERLAP_HPP
......@@ -712,16 +712,16 @@ void Mesh::PrintInfo() const {
PrintInfoEntry("Mesh Name", meshName, 1),
PrintInfoEntry("Distribution", distName, 1),
PrintInfoEntry("Overlap", olapName, 1),
PrintInfoEntry("Level", level),
PrintInfoEntry("pLevel", plevel),
PrintInfoEntry("Dimension", dim()),
PrintInfoEntry("Cells", CellCountGeometry()),
PrintInfoEntry("Overlap", OverlapCountGeometry()),
PrintInfoEntry("Vertices", VertexCountGeometry()),
PrintInfoEntry("Edges", EdgeCountGeometry()),
PrintInfoEntry("Faces", FaceCountGeometry()),
PrintInfoEntry("ProcSets", ProcSetsCountGeometryWithoutInfty()),
PrintInfoEntry("Mesh width [min, max]", MeshWidthStr()),
PrintInfoEntry("Level", level, 1),
PrintInfoEntry("pLevel", plevel, 1),
PrintInfoEntry("Dimension", dim(), 1),
PrintInfoEntry("Cells", CellCountGeometry(), 1),
PrintInfoEntry("Overlap", OverlapCountGeometry(), 1),
PrintInfoEntry("Vertices", VertexCountGeometry(), 1),
PrintInfoEntry("Edges", EdgeCountGeometry(), 1),
PrintInfoEntry("Faces", FaceCountGeometry(), 1),
PrintInfoEntry("ProcSets", ProcSetsCountGeometryWithoutInfty(), 1),
PrintInfoEntry("Mesh width [min, max]", MeshWidthStr(), 1),
PrintInfoEntry("Min cells on proc", PPM->Min(CellCount(), commSplit), 2),
PrintInfoEntry("Max cells on proc", PPM->Max(CellCount(), commSplit), 2),
PrintInfoEntry("Communicator split", commSplit, 2),
......
......@@ -131,7 +131,7 @@ Meshes::Meshes(const CoarseGeometry &CG,
std::map<string, std::list<int>> &bc_nodes,
std::map<string, int> &key_to_bc,
std::map<string, std::list<Point>> &surf,
std::map<string, int> &surf_id) : meshes(1), verbose(1) {
std::map<string, int> &surf_id) : meshes(1) {
initValues();
meshes[0] = std::make_unique<Mesh>(meshName, dist.Name(), olap.Name(),
plevel, 0, commSplit, scaling);
......@@ -144,7 +144,7 @@ Meshes::Meshes(const DataGeometry &OG,
std::map<string, std::list<int>> &bc_nodes,
std::map<string, int> &key_to_bc,
std::map<string, std::list<Point>> &surf,
std::map<string, int> &surf_id) : meshes(1), verbose(1) {
std::map<string, int> &surf_id) : meshes(1){
initValues();
meshes[0] = std::make_unique<Mesh>(meshName, dist.Name(), olap.Name(),
plevel, 0, commSplit, scaling);
......@@ -162,7 +162,6 @@ void Meshes::PrintInfo() const {
void Meshes::initValues() {
commSplit = 0;
scaling = 1.0;
verbose = 1;
fineLevel = 0;
plevel = fineLevel;
coarseLevel = plevel;
......
......@@ -12,15 +12,14 @@ class DataGeometry;
class CoarseGeometry;
class Meshes {
vector<std::unique_ptr<Mesh>> meshes{};
int verbose = 1;
double scaling;
int commSplit;
int verbose;
int plevel;
int fineLevel;
......
......@@ -12,14 +12,6 @@ private:
int _commSplit = 0;
int _distVerbose = 0;
int _olapVerbose = 0;
int _meshVerbose = 0;
int _meshesVerbose = 0;
std::string _distName = "RCB";
std::string _meshName = "Interval";
......@@ -62,26 +54,22 @@ public:
return *this;
}
MeshesCreator WithDistribute(const std::string &distName, int distVerbose = 0) {
_distVerbose = distVerbose;
MeshesCreator WithDistribute(const std::string &distName) {
_distName = distName;
return *this;
}
MeshesCreator WithoutDistribute(int distVerbose = 0) {
_distVerbose = distVerbose;
MeshesCreator WithoutDistribute() {
_distName = "NoDistribution";
return *this;
}
MeshesCreator WithOverlap(const std::string &olapName, int olapVerbose = 0) {
_olapVerbose = olapVerbose;
MeshesCreator WithOverlap(const std::string &olapName) {
_olapName = olapName;
return *this;
}
MeshesCreator WithoutOverlap(int olapVerbose = 0) {
_olapVerbose = olapVerbose;
MeshesCreator WithoutOverlap() {
_olapName = "NoOverlap";
return *this;
}
......@@ -94,8 +82,9 @@ public:
Meshes *Create() {
PPM->SplitCommunicators(_commSplit);
return new Meshes(coarseGeometry, _pLevel, _level, _commSplit,
Distribution(_distName, _commSplit, _distVerbose),
Overlap(_olapName, _commSplit, _olapVerbose), _meshName);
Distribution(_distName, _commSplit),
Overlap(_olapName, _commSplit),
_meshName);
}
};
......
......@@ -19,55 +19,64 @@ using DataSet = MeshPart<Point, DataContainer>;
*/
class DataMesh : public Mesh {
int chamber;
DataSet vertexData;
DataSet cellData;
int chamber;
DataSet vertexData;
DataSet cellData;
public:
/// Creates a new Mesh with a given chamber setting (e.g. only ventricle, fourchamber, ...).
DataMesh(double s = 1, int c = -1) : Mesh(s), chamber(c) {};
/// Creates a new Mesh with a given chamber setting (e.g. only ventricle, fourchamber, ...).
DataMesh(double s = 1, int c = -1) : Mesh(s), chamber(c) {};
/// returns the value corresponding to the chamber setting (e.g. only ventricle, fourchamber, ...).
int GetChamber() const override { return chamber; }
DataMesh(const std::string &meshName, const std::string &distName,
const std::string &olapName, int plevel, int level,
int commSplit, double scaling, int c = -1) :
Mesh(meshName, distName, olapName, plevel, level,
commSplit, scaling), chamber(c) {};
/// an identification if the mesh is a DataMesh
bool StoresData() const override { return true; }
/// returns the value corresponding to the chamber setting (e.g. only ventricle, fourchamber, ...).
int GetChamber() const override { return chamber; }
void Fill(const DataGeometry &geometry);
/// an identification if the mesh is a DataMesh
bool StoresData() const override { return true; }
void Fill(const string &geoName) override;
void Fill(const DataGeometry &geometry);
void Fill(const DataMesh &originMesh);
void Fill(const string &geoName) override;
void InsertData(Point p, DataContainer &data) override;
void Fill(const DataMesh &originMesh);
/// Inserts a new cell into the Mesh with DataContainers at each vertex
cell InsertCell(CELLTYPE tp, int flag, const vector<Point> &x, vector<DataContainer> &d) override;
void InsertData(Point p, DataContainer &data) override;
/// Removes an existing cell and the corresponding data
void RemoveCell(cell c) override;
/// Inserts a new cell into the Mesh with DataContainers at each vertex
cell InsertCell(CELLTYPE tp,
int flag,
const vector<Point> &x,
vector<DataContainer> &d) override;
std::unique_ptr<Mesh> Refine() override;
/// Removes an existing cell and the corresponding data
void RemoveCell(cell c) override;
/// Gets the DataContainer at the given Point
DataContainer &DataAt(const Point &z) const override {
return *(vertexData.Find(z)->second);
}
std::unique_ptr<Mesh> Refine() override;
/// Gets all DataContainers at the verteces of the cell c.
vector<DataContainer> DataIn(const cell &c) const override {
vector<DataContainer> data(c.Corners());
for (int i = 0; i < c.Corners(); i++)
data[i] = *(vertexData.Find(c.Corner(i))->second);
return data;
}
/// Gets the DataContainer at the given Point
DataContainer &DataAt(const Point &z) const override {
return *(vertexData.Find(z)->second);
}
/// Gets an averaged DataContainer at the center of the cell
DataContainer &DataAverage(const cell &c) const override {
return *(cellData.Find(c())->second);
};
/// Gets all DataContainers at the verteces of the cell c.
vector<DataContainer> DataIn(const cell &c) const override {
vector<DataContainer> data(c.Corners());
for (int i = 0; i < c.Corners(); i++)
data[i] = *(vertexData.Find(c.Corner(i))->second);
return data;
}
~DataMesh() override;
/// Gets an averaged DataContainer at the center of the cell
DataContainer &DataAverage(const cell &c) const override {
return *(cellData.Find(c())->second);
};
~DataMesh() override;
};
bool InChamber(int domain, int chamber);
......
......@@ -20,18 +20,18 @@ STDiscretizationT_DGDG_GLGL<T, sDim, tDim>::STDiscretizationT_DGDG_GLGL(
for (int degree = 0; degree <= 6; ++degree) {
auto shape = [degree](CELLTYPE c) { return createGaussLobattoShape<T, sDim, tDim>(c, degree);};
Quadrature *GL = new QintGaussLobattoT(degree + 1);
Quadrature *cell_quad = new QTensor({*GL, *GL});
const Quadrature &face_quad = GetQuadratureT<T, sDim, tDim>(FaceCellType(QUADRILATERAL), 2*maxSpace);
const Quadrature &cell_quad = GetGLQuadrature(QUADRILATERAL, degree + 1);
const Quadrature &face_quad = GetGLQuadratureT<T, sDim, tDim>(INTERVAL, maxSpace + 1);
IDiscretization::insert(QUADRILATERAL,
shape(QUADRILATERAL),
*cell_quad,
cell_quad,
face_quad,
this->cellIdx());
this->timeFaceQuad = GetQuadratureT<T, sDim, tDim>(INTERVAL, 2 * maxTime);
this->timeFaceQuad = GetGLQuadratureT<T, sDim, tDim>(INTERVAL, maxTime + 1);
const Quadrature &time_quad = GetGLQuadrature(INTERVAL, degree + 1);
IDiscretization::insert(INTERVAL,
shape(INTERVAL),
*GL,
time_quad,
this->timeFaceQuad, // use attribute for Time-Face-Quad
this->timeIdx());
/*this->insert(INTERVAL, createGaussLobattoShape<T, sDim, tDim>(INTERVAL, degree),
......@@ -49,6 +49,6 @@ void STDiscretizationT_DGDG_GLGL<T, sDim, tDim>::adaptQuadrature(int quadDeg) {
int maxTime = PPM->Max(max(quadDeg, this->get_max_time_deg()));
this->updateFaceQuadrature(2 * maxSpace, QUADRILATERAL, this->cellIdx());
this->timeFaceQuad = GetQuadratureT<T, sDim, tDim>(INTERVAL, 2 * maxTime);
this->timeFaceQuad = GetGLQuadratureT<T, sDim, tDim>(INTERVAL, maxTime + 1);
}
\ No newline at end of file
......@@ -109,6 +109,50 @@ const Quadrature &GetQuadrature(const CELLTYPE cellType,
exactUpTo2);
}
template<typename T, int sDim, int tDim>
const QuadratureT<T, sDim, tDim> &GetGLQuadratureT(const CELLTYPE cellType, int npcount) {
if (cellType == NONE) return quadEmpty;
if (cellType == POINT) return quadPoint;
int type = cellType;
std::string key= std::string("GL").append(std::to_string(type)).append("_")
.append(std::to_string(npcount));
if (Quads.empty()) Quads.reserve(1);
auto quad = Quads.find(key);
if (quad != Quads.end())
return *(quad->second);
switch (cellType) {
case INTERVAL:
Quads[key] = new QintGaussLobatto(npcount);
break;
case QUADRILATERAL: {
const QuadratureT<> &oneD = GetGLQuadratureT<T, sDim, tDim>(INTERVAL, npcount);
Quads[key] = new QTensor({oneD, oneD});
break;
}
case HEXAHEDRON: {
const QuadratureT<> &oneD = GetGLQuadratureT<T, sDim, tDim>(INTERVAL, npcount);
Quads[key] = new QTensor({oneD, oneD, oneD});
break;
}
default: Exit("No Quadrature implemented for given cell type!");
}
quad = Quads.find(key);
if (quad != Quads.end())
return *(quad->second);
Exit("Quadrature not implemented!");
}
const Quadrature &GetGLQuadrature(const CELLTYPE cellType, int npcount) {
return GetGLQuadratureT<double, SpaceDimension, TimeDimension>(cellType, npcount);
}
template<class T, int sDim, int tDim>
const QuadratureT<T, sDim, tDim> &GetQuadrature(const std::string &name) {
if (name == "Qempty")
......
......@@ -43,8 +43,9 @@ public:
this->w.resize(new_size);
if (quads.size() == 2) {
int n = 0;
for (int qX = 0; qX < quads[0].size(); ++qX) {
for (int qY = 0; qY < quads[1].size(); ++qY) {
for (int qY = 0; qY < quads[1].size(); ++qY) {
for (int qX = 0; qX < quads[0].size(); ++qX) {
this->z[n] = PointT<T, sDim, tDim>(quads[0].QPoint(qX)[0], quads[1].QPoint(qY)[0]);
this->w[n++] = quads[0].Weight(qX) * quads[1].Weight(qY);
}
......@@ -52,9 +53,9 @@ public:
this->name = "QTensor: [" + quads[0].Name() + ", " + quads[1].Name() + "]";
} else if (quads.size() == 3) {
int n = 0;
for (int qX = 0; qX < quads[0].size(); ++qX) {
for (int qZ = 0; qZ < quads[2].size(); ++qZ) {
for (int qY = 0; qY < quads[1].size(); ++qY) {
for (int qZ = 0; qZ < quads[2].size(); ++qZ) {
for (int qX = 0; qX < quads[0].size(); ++qX) {
this->z[n] = PointT<T, sDim, tDim>(quads[0].QPoint(qX)[0],
quads[1].QPoint(qY)[0],
quads[2].QPoint(qZ)[0]);
......@@ -142,6 +143,11 @@ const Quadrature &GetQuadrature(const CELLTYPE,
int exactUpTo1 = -1,
int exactUpTo2 = -1);
template<typename T, int sDim, int tDim>
const QuadratureT<T, sDim, tDim> &GetGLQuadratureT(const CELLTYPE, int npcount);
const Quadrature &GetGLQuadrature(const CELLTYPE, int npcount);
template<class T, int sDim, int tDim>
const QuadratureT<T, sDim, tDim> &GetQuadrature(const std::string &name);
......
......@@ -31,6 +31,7 @@ public:
return "GaussLobattoShape";
}
void NodalPoints(const Cell &c, vector<PointT<T, sDim, tDim>> &z) const override {
if (dim == 1) {
z = GaussLobattoNodalPointsInterval<T, sDim, tDim>(c, polynomialDegree + 1);
......
set(solvers_src
solver/preconditioner/PreconditionerBase.cpp
solver/preconditioner/Preconditioner.cpp
solver/preconditioner/CyclicPreconditioner.cpp
solver/preconditioner/LIB_PS_Solver.cpp
solver/linearsolver/LinearSolver.cpp
solver/Solver.cpp
......
#include "CyclicPreconditioner.hpp"
#include <memory>
#include "Preconditioner.hpp"
CyclicPreconditioner::CyclicPreconditioner(const std::string &prefix) {
vector<string> pc_names;
config.get(prefix + "CyclicPCNames", pc_names);
config.get(prefix + "CyclicPCIndices", cycle_indices);
for (const string &name: pc_names) {
preconditioners.push_back(std::shared_ptr<Preconditioner>(GetPC(name)));
}
}
CyclicPreconditioner::CyclicPreconditioner(
const vector<std::shared_ptr<Preconditioner>> &preconditioners,
const vector<int> &cycle_indices
) : preconditioners(preconditioners),
cycle_indices(cycle_indices) {
}
void CyclicPreconditioner::Construct(const Matrix &M) {
for (auto &PC : preconditioners) {
PC->Construct(M);
}
};
void CyclicPreconditioner::multiply(Vector &u, const Vector &b) const {
auto PC = preconditioners[cycle_indices[current_cycle_index]];
Date start;
PC->multiply(u, b);
Time t = Date() - start;
vout(1) << "[CPC] using " << PC->Name() << " t=" << t << endl;
current_cycle_index++;
current_cycle_index %= cycle_indices.size();
}
void CyclicPreconditioner::Destruct() {
for (auto &PC: preconditioners) {
PC->Destruct();
}
}
void CyclicPreconditioner::reset() {
current_cycle_index = 0;
}
string CyclicPreconditioner::Name() const {
std::stringstream ss;
ss << "CyclicPreconditioner[" << preconditioners[0]->Name();
for (int i = 1; i < preconditioners.size(); i++) {
ss << ", " << preconditioners[i]->Name();
}
ss << "](" << cycle_indices[0];
for (int i = 1; i < cycle_indices.size(); i++) {
ss << ", " << cycle_indices[i];
}
ss << ")";
return ss.str();
}
void CyclicPreconditioner::transpose() {
for (auto &PC: preconditioners) {
PC->transpose();
}
}
#ifndef CYCLICPRECONDITIONER_HPP
#define CYCLICPRECONDITIONER_HPP
#include "solver/preconditioner/PreconditionerBase.hpp"
#include <memory>
class CyclicPreconditioner : public Preconditioner {
private:
vector<std::shared_ptr<Preconditioner>> preconditioners;
vector<int> cycle_indices;
mutable int current_cycle_index = 0;
public:
void multiply(Vector &u, const Vector &b) const override;
explicit CyclicPreconditioner(const std::string &prefix = "");
CyclicPreconditioner(const vector<std::shared_ptr<Preconditioner>> &preconditioners,
const vector<int> &cycle_indices);
virtual void Construct(const Matrix &M) override;
virtual void Destruct() override;
void reset() override;
string Name() const override;
void transpose() override;
};
#endif //CYCLICPRECONDITIONER_HPP
......@@ -8,6 +8,8 @@
#include "LIB_PS_Solver.hpp"
#include "CyclicPreconditioner.hpp"
#include <set>
......@@ -1077,8 +1079,6 @@ public:
PC[k]->Construct(*A[k]);
vv[k - 1] = new Vector(G[k - 1]);
T[k - 1]->Project(A[k]->GetVector(), *vv[k - 1]);
// mout << OUT(A[k]->GetVector());
// mout << OUT(*vv[k-1]);
AA[k - 1] = new Matrix(*vv[k - 1]);
assemble.Initialize(*vv[k - 1]);
assemble.Jacobi(*vv[k - 1], *AA[k - 1]);
......
#ifndef _PRECONDITIONER_H_
#define _PRECONDITIONER_H_
#include "PreconditionerBase.hpp"
#include "TimeDate.hpp"
#include "Algebra.hpp"
#include "Assemble.hpp"
#include "Transfer.hpp"
#include "Sparse.hpp"
class Preconditioner : public Operator {
protected:
int verbose;
public:
virtual void Construct(const Matrix &) = 0;
virtual void Destruct() = 0;
Preconditioner() : verbose(0) {
config.get("PreconditionerVerbose", verbose);
}
virtual void multiply_transpose(Vector &u, const Vector &b) const {
multiply(u, b);
return;
Vector r(b);
Accumulate(r);
multiply(u, r);
MakeAdditive(u);
Accumulate(u);
}
virtual ~Preconditioner() {}
// transpose is used for the adjoint equation. For example by a multigrid preconditioner.
virtual void transpose() {}
// reset for preconditioners with state
virtual void reset() {}
virtual string Name() const = 0;
friend Logging &operator<<(Logging &s, const Preconditioner &PC) {
return s << PC.Name();
}
void SetVerbose(int v) { verbose = v; }
};
inline constAB<Operator, Vector>
operator*(const Preconditioner &PC, const Vector &v) {
return constAB<Operator, Vector>(PC, v);
}
inline constAB<Vector, Operator>
operator*(const Vector &v, const Preconditioner &PC) {
return constAB<Vector, Operator>(v, PC);
}
inline constAB<Operator, Vectors>
operator*(const Preconditioner &PC, const Vectors &v) {
return constAB<Operator, Vectors>(PC, v);
}
inline constAB<Vectors, Operator>
operator*(const Vectors &v, const Preconditioner &PC) {
return constAB<Vectors, Operator>(v, PC);
}
class Transfer;
......