From 0e2727113e8d48d0f161eac3f1b662c2a8b098c8 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Fri, 20 Nov 2009 13:51:58 +0000
Subject: [PATCH] Small changes in interface to make new solver possible.
---
AMDiS/src/OEMSolver.h | 32 +++++++++---------
AMDiS/src/Operator.cc | 34 +++----------------
AMDiS/src/Operator.h | 35 ++++++++++++++-----
AMDiS/src/SolverMatrix.h | 73 +++++++++++++++++++++++-----------------
4 files changed, 90 insertions(+), 84 deletions(-)
diff --git a/AMDiS/src/OEMSolver.h b/AMDiS/src/OEMSolver.h
index ad4fb7aa..a80dc29a 100644
--- a/AMDiS/src/OEMSolver.h
+++ b/AMDiS/src/OEMSolver.h
@@ -145,35 +145,35 @@ namespace AMDiS {
}
/// Solve a linear system for a vectorial problem.
- int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
- SystemVector& x,
- SystemVector& b)
+ virtual int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
+ SystemVector& x,
+ SystemVector& b)
{
- int ns= x.getSize(), // Number of systems.
- size= x.getUsedSize(); // Size of all DOFVectors
+ int ns = x.getSize(), // Number of systems.
+ size = x.getUsedSize(); // Size of all DOFVectors
// Copy vectors
mtl::dense_vector<value_type> xx(size), bb(size);
for (int i = 0, counter = 0; i < ns; i++) {
- DOFVector<double>::Iterator it_b(b.getDOFVector(i), USED_DOFS);
- DOFVector<double>::Iterator it_x(x.getDOFVector(i), USED_DOFS);
-
- for (it_b.reset(), it_x.reset(); !it_b.end(); ++it_b, ++it_x) {
- bb[counter] = *it_b;
- xx[counter] = *it_x;
- counter++;
- }
+ DOFVector<double>::Iterator it_b(b.getDOFVector(i), USED_DOFS);
+ DOFVector<double>::Iterator it_x(x.getDOFVector(i), USED_DOFS);
+
+ for (it_b.reset(), it_x.reset(); !it_b.end(); ++it_b, ++it_x) {
+ bb[counter] = *it_b;
+ xx[counter] = *it_x;
+ counter++;
+ }
}
-
+
// Solver on DOFVector for single system
int r = solveSystem(A.getMatrix(), xx, bb);
-
+
for (int i = 0, counter = 0; i < ns; i++) {
DOFVector<double>::Iterator it(x.getDOFVector(i), USED_DOFS);
for (it.reset(); !it.end(); ++it)
*it = xx[counter++];
}
-
+
return r;
}
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index a5757ecf..459c90d2 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -214,30 +214,6 @@ namespace AMDiS {
LALt[i][j] += factor * Lambda[i][k] * Lambda[j][l];
}
- void OperatorTerm::l1lt(const DimVec<WorldVector<double> >& Lambda,
- DimMat<double>& LALt,
- double factor) const
- {
- const int dimOfWorld = Global::getGeo(WORLD);
- int dim = LALt.getNumRows() - 1;
-
- for (int i = 0; i <= dim; i++) {
- double val = 0.0;
- for (int k = 0; k < dimOfWorld; k++)
- val += Lambda[i][k] * Lambda[i][k];
- val *= factor;
- LALt[i][i] += val;
- for (int j = i + 1; j <= dim; j++) {
- val = 0.0;
- for (int k = 0; k < dimOfWorld; k++)
- val += Lambda[i][k] * Lambda[j][k];
- val *= factor;
- LALt[i][j] += val;
- LALt[j][i] += val;
- }
- }
- }
-
Operator::Operator(Flag operatorType,
const FiniteElemSpace *row,
const FiniteElemSpace *col)
@@ -740,7 +716,7 @@ namespace AMDiS {
Vec2AndGradAtQP_FOT::Vec2AndGradAtQP_FOT(DOFVectorBase<double> *dv1,
DOFVectorBase<double> *dv2,
TertiaryAbstractFunction<double, double, double, WorldVector<double> > *f_, WorldVector<double> *b_)
- : FirstOrderTerm(8),
+ : FirstOrderTerm(f_->getDegree()),
vec1(dv1),
vec2(dv2),
f(f_),
@@ -829,7 +805,7 @@ namespace AMDiS {
Vec2AtQP_FOT::Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
WorldVector<double> *b_)
- : FirstOrderTerm(8), vec1(dv1), vec2(dv2), b(b_)
+ : FirstOrderTerm(0), vec1(dv1), vec2(dv2), b(b_)
{
auxFESpaces.push_back(dv1->getFESpace());
auxFESpaces.push_back(dv2->getFESpace());
@@ -837,7 +813,7 @@ namespace AMDiS {
Vec2AtQP_FOT::Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
int bIdx)
- : FirstOrderTerm(8), vec1(dv1), vec2(dv2)
+ : FirstOrderTerm(0), vec1(dv1), vec2(dv2)
{
bOne = bIdx;
auxFESpaces.push_back(dv1->getFESpace());
@@ -884,7 +860,7 @@ namespace AMDiS {
Vec3FctAtQP_FOT::Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,DOFVectorBase<double> *dv3,
BinaryAbstractFunction<double, double, double> *f_,
WorldVector<double> *bvec)
- : FirstOrderTerm(10),
+ : FirstOrderTerm(f_->getDegree()),
vec1(dv1),
vec2(dv2),
vec3(dv3),
@@ -899,7 +875,7 @@ namespace AMDiS {
Vec3FctAtQP_FOT::Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,DOFVectorBase<double> *dv3,
BinaryAbstractFunction<double, double, double> *f_,
int b)
- : FirstOrderTerm(10),
+ : FirstOrderTerm(f_->getDegree()),
vec1(dv1),
vec2(dv2),
vec3(dv3),
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index c274d643..00eb5f39 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -166,9 +166,28 @@ namespace AMDiS {
double factor);
/// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for A equal to the identity.
- void l1lt(const DimVec<WorldVector<double> >& Lambda,
- DimMat<double>& LALt,
- double factor) const;
+ inline void l1lt(const DimVec<WorldVector<double> >& Lambda,
+ DimMat<double>& LALt,
+ double factor) const
+ {
+ const int dim = LALt.getNumRows();
+
+ for (int i = 0; i < dim; i++) {
+ double val = 0.0;
+ for (int k = 0; k < dimOfWorld; k++)
+ val += Lambda[i][k] * Lambda[i][k];
+ val *= factor;
+ LALt[i][i] += val;
+ for (int j = i + 1; j < dim; j++) {
+ val = 0.0;
+ for (int k = 0; k < dimOfWorld; k++)
+ val += Lambda[i][k] * Lambda[j][k];
+ val *= factor;
+ LALt[i][j] += val;
+ LALt[j][i] += val;
+ }
+ }
+ }
/// Evaluation of \f$ \Lambda \cdot b\f$.
inline void lb(const DimVec<WorldVector<double> >& Lambda,
@@ -176,9 +195,9 @@ namespace AMDiS {
DimVec<double>& Lb,
double factor) const
{
- const int dim = Lambda.size() - 1;
+ const int dim = Lambda.getSize();
- for (int i = 0; i <= dim; i++) {
+ for (int i = 0; i < dim; i++) {
double val = 0.0;
for (int j = 0; j < dimOfWorld; j++)
@@ -193,7 +212,7 @@ namespace AMDiS {
DimVec<double>& Lb,
double factor) const
{
- const int dim = Lambda.size();
+ const int dim = Lambda.getSize();
for (int i = 0; i < dim; i++)
Lb[i] += Lambda[i][bOne] * factor;
@@ -207,9 +226,9 @@ namespace AMDiS {
DimVec<double>& Lb,
double factor) const
{
- const int dim = Lambda.size() - 1;
+ const int dim = Lambda.getSize();
- for (int i = 0; i <= dim; i++) {
+ for (int i = 0; i < dim; i++) {
double val = 0.0;
for (int j = 0; j < dimOfWorld; j++)
diff --git a/AMDiS/src/SolverMatrix.h b/AMDiS/src/SolverMatrix.h
index 43077aa7..f090a3ba 100644
--- a/AMDiS/src/SolverMatrix.h
+++ b/AMDiS/src/SolverMatrix.h
@@ -61,38 +61,49 @@ namespace AMDiS {
template <>
class SolverMatrix<Matrix<DOFMatrix*> >
{
- public :
- void setMatrix(const Matrix<DOFMatrix*>& A)
- {
- int ns= A.getSize();
- std::vector<int> block_starts(ns+1);
- block_starts[0]= 0;
- for (int i= 0; i < ns; ++i)
- block_starts[i+1]= block_starts[i] + A[i][i]->getFESpace()->getAdmin()->getUsedSize();
-
- matrix.change_dim(block_starts[ns], block_starts[ns]);
- set_to_zero(matrix);
-
- int nnz= 0;
- for (int rb= 0; rb < ns; ++rb)
- for (int cb= 0; cb < ns; ++cb)
- if (A[rb][cb])
- nnz+= A[rb][cb]->getBaseMatrix().nnz();
+ public :
+ void setMatrix(const Matrix<DOFMatrix*>& A)
+ {
+ int ns = A.getSize();
+ std::vector<int> block_starts(ns + 1);
+
+ block_starts[0] = 0;
+ for (int i= 0; i < ns; ++i)
+ block_starts[i+1]= block_starts[i] + A[i][i]->getFESpace()->getAdmin()->getUsedSize();
- DOFMatrix::inserter_type ins(matrix, int(1.2 * nnz / matrix.dim1()));
- for (int rb= 0; rb < ns; ++rb)
- for (int cb= 0; cb < ns; ++cb)
- if (A[rb][cb])
- ins[block_starts[rb]][block_starts[cb]] << A[rb][cb]->getBaseMatrix();
- }
-
- const DOFMatrix::base_matrix_type& getMatrix() const
- {
- return matrix;
- }
-
- private:
- DOFMatrix::base_matrix_type matrix;
+ matrix.change_dim(block_starts[ns], block_starts[ns]);
+ set_to_zero(matrix);
+
+ int nnz= 0;
+ for (int rb= 0; rb < ns; ++rb)
+ for (int cb= 0; cb < ns; ++cb)
+ if (A[rb][cb])
+ nnz+= A[rb][cb]->getBaseMatrix().nnz();
+
+ DOFMatrix::inserter_type ins(matrix, int(1.2 * nnz / matrix.dim1()));
+ for (int rb= 0; rb < ns; ++rb)
+ for (int cb= 0; cb < ns; ++cb)
+ if (A[rb][cb])
+ ins[block_starts[rb]][block_starts[cb]] << A[rb][cb]->getBaseMatrix();
+
+ originalMat = &A;
+ }
+
+ const DOFMatrix::base_matrix_type& getMatrix() const
+ {
+ return matrix;
+ }
+
+ const Matrix<DOFMatrix*>* getOriginalMat() const
+ {
+ return originalMat;
+ }
+
+ private:
+ DOFMatrix::base_matrix_type matrix;
+
+ const Matrix<DOFMatrix*>* originalMat;
+
};
--
GitLab