Commit 0e272711 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Small changes in interface to make new solver possible.

parent 06dbd114
......@@ -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;
}
......
......@@ -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),
......
......@@ -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++)
......
......@@ -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;
};
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment