From bc5b945f8e70fe454ba8f1efe815755adec2ec39 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Wed, 2 Sep 2009 13:49:58 +0000
Subject: [PATCH] Serveral bug fixes for parallelization.
---
AMDiS/libtool | 18 +--
AMDiS/src/DOFVector.hh | 17 +--
AMDiS/src/Element.cc | 14 +-
AMDiS/src/ElementDofIterator.h | 2 +-
AMDiS/src/Mesh.cc | 30 ++++
AMDiS/src/Mesh.h | 7 +
AMDiS/src/ParallelDomainBase.cc | 70 +++++----
AMDiS/src/ProblemVec.cc | 263 +++++++++++++++++++++++++++-----
AMDiS/src/ProblemVec.h | 30 ++--
AMDiS/src/SubAssembler.cc | 9 +-
AMDiS/src/ZeroOrderAssembler.cc | 30 ++--
11 files changed, 349 insertions(+), 141 deletions(-)
diff --git a/AMDiS/libtool b/AMDiS/libtool
index 41ace2e6..9fbd8125 100755
--- a/AMDiS/libtool
+++ b/AMDiS/libtool
@@ -82,13 +82,13 @@ AR="ar"
AR_FLAGS="cru"
# A C compiler.
-LTCC="gcc"
+LTCC="/usr/lib/openmpi/1.2.7-gcc//bin/mpicc"
# LTCC compiler flags.
LTCFLAGS="-g -O2"
# A language-specific compiler.
-CC="gcc"
+CC="/usr/lib/openmpi/1.2.7-gcc//bin/mpicc"
# Is the compiler the GNU C compiler?
with_gcc=yes
@@ -171,7 +171,7 @@ dlopen_self=unknown
dlopen_self_static=unknown
# Compiler flag to prevent dynamic linking.
-link_static_flag="-static"
+link_static_flag=""
# Compiler flag to turn off builtin functions.
no_builtin_flag=" -fno-builtin"
@@ -6798,13 +6798,13 @@ AR="ar"
AR_FLAGS="cru"
# A C compiler.
-LTCC="gcc"
+LTCC="/usr/lib/openmpi/1.2.7-gcc//bin/mpicc"
# LTCC compiler flags.
LTCFLAGS="-g -O2"
# A language-specific compiler.
-CC="g++"
+CC="/usr/lib/openmpi/1.2.7-gcc//bin/mpiCC"
# Is the compiler the GNU C compiler?
with_gcc=yes
@@ -6887,7 +6887,7 @@ dlopen_self=unknown
dlopen_self_static=unknown
# Compiler flag to prevent dynamic linking.
-link_static_flag="-static"
+link_static_flag=""
# Compiler flag to turn off builtin functions.
no_builtin_flag=" -fno-builtin"
@@ -6954,11 +6954,11 @@ predeps=""
# Dependencies to place after the objects being linked to create a
# shared library.
-postdeps="-lstdc++ -lm -lgcc_s -lc -lgcc_s"
+postdeps="-lmpi_cxx -lmpi -lopen-rte -lopen-pal -ldl -lnsl -lutil -ldl -lstdc++ -lm -lgcc_s -lpthread -lc -lgcc_s"
# The library search path used internally by the compiler when linking
# a shared library.
-compiler_lib_search_path="-L/usr/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-linux/4.1.2/../../.."
+compiler_lib_search_path="-L/usr/lib/openmpi/1.2.7-gcc/lib -L/usr/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-linux/4.1.2/../../.."
# Method to check whether dependent libraries are shared objects.
deplibs_check_method="pass_all"
@@ -7103,7 +7103,7 @@ AR="ar"
AR_FLAGS="cru"
# A C compiler.
-LTCC="gcc"
+LTCC="/usr/lib/openmpi/1.2.7-gcc//bin/mpicc"
# LTCC compiler flags.
LTCFLAGS="-g -O2"
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index c761c559..198715e6 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -432,15 +432,13 @@ namespace AMDiS {
{
const BasisFunction *basFct = traverseVector->getFESpace()->getBasisFcts();
const DOFAdmin* admin = traverseVector->getFESpace()->getAdmin();
- const DegreeOfFreedom *dof = basFct->getLocalIndices(const_cast<Element *>(elinfo->getElement()),
- admin, NULL);
- const T *inter_val = const_cast<BasisFunction*>(basFct)->interpol(elinfo,
- 0,
- NULL,
- traverseVector->interFct,
- NULL);
- int number = basFct->getNumber();
- for (int i = 0; i < number; i++)
+ const DegreeOfFreedom *dof =
+ basFct->getLocalIndices(const_cast<Element*>(elinfo->getElement()), admin, NULL);
+ const T *inter_val =
+ const_cast<BasisFunction*>(basFct)->interpol(elinfo, 0, NULL,
+ traverseVector->interFct, NULL);
+ int nBasFcts = basFct->getNumber();
+ for (int i = 0; i < nBasFcts; i++)
(*traverseVector)[dof[i]] = inter_val[i];
return 0;
@@ -959,7 +957,6 @@ namespace AMDiS {
static T* localVec = NULL;
static int localVecSize = 0;
const DOFAdmin* admin = feSpace->getAdmin();
-
T *result;
if (d) {
diff --git a/AMDiS/src/Element.cc b/AMDiS/src/Element.cc
index 6934a84c..97c66eb7 100644
--- a/AMDiS/src/Element.cc
+++ b/AMDiS/src/Element.cc
@@ -258,13 +258,14 @@ namespace AMDiS {
void Element::newDOFFct2(const DOFAdmin* admin)
{
- int i, j, k, n0, nd, nd0;
+ int i, j, k, n0, nd0;
DegreeOfFreedom *ldof;
int vertices = mesh->getGeo(VERTEX);
int edges = mesh->getGeo(EDGE);
int faces = mesh->getGeo(FACE);
- if (nd = admin->getNumberOfDOFs(VERTEX)) {
+ int nd = admin->getNumberOfDOFs(VERTEX);
+ if (nd) {
nd0 = admin->getNumberOfPreDOFs(VERTEX);
n0 = admin->getMesh()->getNode(VERTEX);
for (i = 0; i < vertices; i++) {
@@ -273,7 +274,8 @@ namespace AMDiS {
}
if (mesh->getDim() > 1) {
- if (nd = admin->getNumberOfDOFs(EDGE)) {
+ nd = admin->getNumberOfDOFs(EDGE);
+ if (nd) {
nd0 = admin->getNumberOfPreDOFs(EDGE);
n0 = admin->getMesh()->getNode(EDGE);
for (i = 0; i < edges; i++) {
@@ -283,7 +285,8 @@ namespace AMDiS {
}
if (mesh->getDim() == 3) {
- if (nd = admin->getNumberOfDOFs(FACE)) {
+ nd = admin->getNumberOfDOFs(FACE);
+ if (nd) {
nd0 = admin->getNumberOfPreDOFs(FACE);
n0 = admin->getMesh()->getNode(FACE);
for (i = 0; i < faces; i++) {
@@ -292,7 +295,8 @@ namespace AMDiS {
}
}
- if (nd = admin->getNumberOfDOFs(CENTER)) {
+ nd = admin->getNumberOfDOFs(CENTER);
+ if (nd) {
nd0 = admin->getNumberOfPreDOFs(CENTER);
n0 = admin->getMesh()->getNode(CENTER);
// only one center
diff --git a/AMDiS/src/ElementDofIterator.h b/AMDiS/src/ElementDofIterator.h
index a64cee19..718632a2 100644
--- a/AMDiS/src/ElementDofIterator.h
+++ b/AMDiS/src/ElementDofIterator.h
@@ -62,7 +62,7 @@ namespace AMDiS {
bool nextStrict();
/// Returns the dof index of the current dof.
- inline const DegreeOfFreedom getDof()
+ inline DegreeOfFreedom getDof()
{
if (inOrder)
return dofs[node0 + elementPos][n0 + orderPosition[dofPos]];
diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc
index 3e75bcad..07b11ec9 100644
--- a/AMDiS/src/Mesh.cc
+++ b/AMDiS/src/Mesh.cc
@@ -24,6 +24,8 @@
#include "Projection.h"
#include "ElInfoStack.h"
+#include "mpi.h"
+
namespace AMDiS {
#define TIME_USED(f,s) ((double)((s)-(f))/(double)CLOCKS_PER_SEC)
@@ -1194,6 +1196,34 @@ namespace AMDiS {
macroFileInfo = NULL;
}
+ void Mesh::computeMatrixFillin(const FiniteElemSpace *feSpace,
+ std::vector<int> &nnzInRow, int &overall, int &average)
+ {
+ std::map<DegreeOfFreedom, int> dofCounter;
+
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(this, -1, Mesh::CALL_LEAF_EL);
+ ElementDofIterator elDofIter(feSpace);
+ while (elInfo) {
+ elDofIter.reset(elInfo->getElement());
+ do {
+ DegreeOfFreedom dof = elDofIter.getDof();
+ if (dofCounter.count(dof) == 0) {
+ dofCounter[dof] = 1;
+ } else {
+ dofCounter[dof]++;
+ }
+ } while (elDofIter.next());
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ overall = 0;
+ for (std::map<DegreeOfFreedom, int>::iterator it = dofCounter.begin();
+ it != dofCounter.end(); ++it) {
+ overall += it->second * 15;
+ }
+ }
+
int Mesh::calcMemoryUsage()
{
int result = sizeof(Mesh);
diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h
index 68dc9399..015e5aed 100644
--- a/AMDiS/src/Mesh.h
+++ b/AMDiS/src/Mesh.h
@@ -564,6 +564,13 @@ namespace AMDiS {
///
void clearMacroFileInfo();
+ /** \brief
+ * Traverse this mesh to compute the number of non zero elements the assembled
+ * matrix will have in each row.
+ */
+ void computeMatrixFillin(const FiniteElemSpace *feSpace,
+ std::vector<int> &nnzInRow, int &overall, int &average);
+
///
int calcMemoryUsage();
diff --git a/AMDiS/src/ParallelDomainBase.cc b/AMDiS/src/ParallelDomainBase.cc
index c77f2bb7..08839103 100644
--- a/AMDiS/src/ParallelDomainBase.cc
+++ b/AMDiS/src/ParallelDomainBase.cc
@@ -144,23 +144,8 @@ namespace AMDiS {
DbgTestCommonDofs(true);
#endif
-
- // === Create petsc matrix. ===
-
nRankRows = nRankDOFs * nComponents;
nOverallRows = nOverallDOFs * nComponents;
-
- VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
- VecSetSizes(petscRhsVec, nRankRows, nOverallRows);
- VecSetType(petscRhsVec, VECMPI);
-
- VecCreate(PETSC_COMM_WORLD, &petscSolVec);
- VecSetSizes(petscSolVec, nRankRows, nOverallRows);
- VecSetType(petscSolVec, VECMPI);
-
- VecCreate(PETSC_COMM_WORLD, &petscTmpVec);
- VecSetSizes(petscTmpVec, nRankRows, nOverallRows);
- VecSetType(petscTmpVec, VECMPI);
}
@@ -244,21 +229,12 @@ namespace AMDiS {
for (icursor_type icursor = begin<nz>(cursor),
icend = end<nz>(cursor); icursor != icend; ++icursor)
if (value(*icursor) != 0.0) {
- if (mpiRank == 0 && r == 0)
- std::cout << "C = " << col(*icursor)
- << " = " << mapLocalGlobalDOFs[col(*icursor)]
- << " = " << mapLocalGlobalDOFs[col(*icursor)] * dispMult + dispAddCol
- << " -> " << value(*icursor)
- << std::endl;
-
cols.push_back(mapLocalGlobalDOFs[col(*icursor)] * dispMult + dispAddCol);
values.push_back(value(*icursor));
}
MatSetValues(petscMatrix, 1, &r, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES);
}
-
- std::cout << "==========" << std::endl;
}
@@ -281,6 +257,18 @@ namespace AMDiS {
MatSetSizes(petscMatrix, nRankRows, nRankRows, nOverallRows, nOverallRows);
MatSetType(petscMatrix, MATAIJ);
+ VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
+ VecSetSizes(petscRhsVec, nRankRows, nOverallRows);
+ VecSetType(petscRhsVec, VECMPI);
+
+ VecCreate(PETSC_COMM_WORLD, &petscSolVec);
+ VecSetSizes(petscSolVec, nRankRows, nOverallRows);
+ VecSetType(petscSolVec, VECMPI);
+
+ VecCreate(PETSC_COMM_WORLD, &petscTmpVec);
+ VecSetSizes(petscTmpVec, nRankRows, nOverallRows);
+ VecSetType(petscTmpVec, VECMPI);
+
setDofMatrix(mat);
MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
@@ -299,6 +287,18 @@ namespace AMDiS {
clock_t first = clock();
+ VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
+ VecSetSizes(petscRhsVec, nRankRows, nOverallRows);
+ VecSetType(petscRhsVec, VECMPI);
+
+ VecCreate(PETSC_COMM_WORLD, &petscSolVec);
+ VecSetSizes(petscSolVec, nRankRows, nOverallRows);
+ VecSetType(petscSolVec, VECMPI);
+
+ VecCreate(PETSC_COMM_WORLD, &petscTmpVec);
+ VecSetSizes(petscTmpVec, nRankRows, nOverallRows);
+ VecSetType(petscTmpVec, VECMPI);
+
using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
namespace traits= mtl::traits;
typedef DOFMatrix::base_matrix_type Matrix;
@@ -481,7 +481,7 @@ namespace AMDiS {
}
i++;
- }
+ }
MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows,
0, d_nnz, 0, o_nnz, &petscMatrix);
@@ -499,7 +499,7 @@ namespace AMDiS {
for (int i = 0; i < nComponents; i++)
for (int j = 0; j < nComponents; j++)
- if ((*mat)[i][j])
+ if ((*mat)[i][j])
setDofMatrix((*mat)[i][j], nComponents, i, j);
MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
@@ -586,6 +586,9 @@ namespace AMDiS {
delete [] sendBuffers[i];
MatDestroy(petscMatrix);
+ VecDestroy(petscRhsVec);
+ VecDestroy(petscSolVec);
+ VecDestroy(petscTmpVec);
}
@@ -611,17 +614,10 @@ namespace AMDiS {
for (int i = 0; i < nComponents; i++) {
DOFVector<double> *dofvec = vec.getDOFVector(i);
- for (int j = 0; j < nRankDOFs; j++) {
- if (i == 0 && mapLocalToDofIndex[j] == 0)
- std::cout << "ZUGRIFF AUF: " << j * nComponents + i << std::endl;
-
- (*dofvec)[mapLocalToDofIndex[j]] = vecPointer[j * nComponents + i];
- }
+ for (int j = 0; j < nRankDOFs; j++)
+ (*dofvec)[mapLocalToDofIndex[j]] = vecPointer[j * nComponents + i];
}
- if (mpiRank == 0)
- std::cout << "RESULT = " << (*(vec.getDOFVector(0)))[0] << std::endl;
-
VecRestoreArray(petscSolVec, &vecPointer);
synchVectors(vec);
@@ -637,6 +633,9 @@ namespace AMDiS {
MSG(" Residual norm: %e\n", norm);
MatDestroy(petscMatrix);
+ VecDestroy(petscRhsVec);
+ VecDestroy(petscSolVec);
+ VecDestroy(petscTmpVec);
KSPDestroy(solver);
}
@@ -933,7 +932,6 @@ namespace AMDiS {
int i = 0;
for (DofContainer::iterator dofIt = rankAllDofs.begin();
dofIt != rankAllDofs.end(); ++dofIt) {
-
rankDofsNewLocalIndex[*dofIt] = i;
// First, we set all dofs in ranks partition to be owend by the rank. Later,
// the dofs in ranks partition that are owned by other rank are set to false.
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 7fe5b59a..a4fbd4c6 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -32,7 +32,7 @@ namespace AMDiS {
Flag adoptFlag)
{
FUNCNAME("ProblemVec::initialize()");
-
+
// === create meshes ===
if (meshes.size() != 0) {
WARNING("meshes already created\n");
@@ -646,6 +646,14 @@ namespace AMDiS {
int nnz = 0;
for (int i = 0; i < nComponents; i++) {
+
+#if 0
+ std::vector<int> nnzInRow;
+ int overallNnz, averageNnz;
+ componentSpaces[i]->getMesh()->computeMatrixFillin(componentSpaces[i],
+ nnzInRow, overallNnz, averageNnz);
+#endif
+
MSG("%d DOFs for %s\n",
componentSpaces[i]->getAdmin()->getUsedSize(),
componentSpaces[i]->getName().c_str());
@@ -1384,7 +1392,7 @@ namespace AMDiS {
{
const BasisFunction *basisFcts = componentSpaces[0]->getBasisFcts();
WorldVector<double> coords;
- ElementDofIterator elDofIter(componentSpaces[0]);
+ ElementDofIterator elDofIter(componentSpaces[0], true);
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(componentMeshes[0], -1,
@@ -1456,7 +1464,7 @@ namespace AMDiS {
for (DofToCoord::iterator it = dofToCoords.begin(); it != dofToCoords.end(); ++it) {
out << it->first << std::endl;
for (int j = 0; j < dimOfWorld; j++)
- out << (it->second)[j] << std::endl;
+ out << (it->second)[j] << std::endl;
}
// Write the matrices.
@@ -1491,12 +1499,40 @@ namespace AMDiS {
TEST_EXIT(testNnz == nnz)("NNZ does not fit together!\n");
}
}
+
+
+ for (int i = 0; i < nComponents; i++) {
+ out << rhs->getDOFVector(i)->getUsedSize() << std::endl;
+ DOFIterator<double> it(rhs->getDOFVector(i), USED_DOFS);
+ int c = 0;
+ for (it.reset(); !it.end(); ++it) {
+ out << c << " " << *it << std::endl;
+ c++;
+ }
+
+ TEST_EXIT(c == rhs->getDOFVector(i)->getUsedSize())
+ ("RHS NNZ does not fit together!\n");
+ }
+
+
+ for (int i = 0; i < nComponents; i++) {
+ out << solution->getDOFVector(i)->getUsedSize() << std::endl;
+ DOFIterator<double> it(solution->getDOFVector(i), USED_DOFS);
+ int c = 0;
+ for (it.reset(); !it.end(); ++it) {
+ out << c << " " << *it << std::endl;
+ c++;
+ }
+
+ TEST_EXIT(c == rhs->getDOFVector(i)->getUsedSize())
+ ("RHS NNZ does not fit together!\n");
+ }
out.close();
}
-
- void ProblemVec::readAndCompareDbgMatrix(std::string filename)
+
+ void ProblemVec::readAndCompareDbgMatrix(std::vector<std::string> filenames)
{
using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
namespace traits = mtl::traits;
@@ -1513,42 +1549,92 @@ namespace AMDiS {
DegreeOfFreedom dof;
int dimOfWorld = Global::getGeo(WORLD);
- // Open file and read the number of DOFs the file contains.
- std::ifstream in(filename.c_str());
- int nReadDofs;
- in >> nReadDofs;
+ std::vector<std::vector<std::map<std::pair<int, int>, double> > > nnzValues;
+ std::vector<std::map<int, double> > rhsValues;
+ std::vector<std::map<int, double> > solValues;
- // Read all DOF indices and their world coordinates.
- for (int i = 0; i < nReadDofs; i++) {
- in >> dof;
- for (int j = 0; j < dimOfWorld; j++)
- in >> coords[j];
- dofToCoord[dof] = coords;
- }
+ nnzValues.resize(nComponents);
+ rhsValues.resize(nComponents);
+ solValues.resize(nComponents);
+ for (int i = 0; i < nComponents; i++)
+ nnzValues[i].resize(nComponents);
- // Now we traverse all matrices and check them.
- for (int i = 0; i < nComponents; i++) {
- for (int j = 0; j < nComponents; j++) {
- DOFMatrix *dofmatrix = (*systemMatrix)[i][j];
-
- if (!dofmatrix)
- continue;
+ for (std::vector<std::string>::iterator fileIt = filenames.begin();
+ fileIt != filenames.end(); ++fileIt) {
+ // Open file and read the number of DOFs the file contains.
+ std::ifstream in(fileIt->c_str());
+ int nReadDofs;
+ in >> nReadDofs;
+
+ // Read all DOF indices and their world coordinates.
+ dofToCoord.clear();
+ for (int i = 0; i < nReadDofs; i++) {
+ in >> dof;
+ for (int j = 0; j < dimOfWorld; j++)
+ in >> coords[j];
+ dofToCoord[dof] = coords;
+ }
+
+ // Now we traverse all matrices and check them.
+ for (int i = 0; i < nComponents; i++) {
+ for (int j = 0; j < nComponents; j++) {
+ DOFMatrix *dofmatrix = (*systemMatrix)[i][j];
+
+ if (!dofmatrix)
+ continue;
+
+ int readNnz;
+ in >> readNnz;
+
+ // Read each entry in file and check it with the corresponding matrix value
+ // in this problem.
+ for (int k = 0; k < readNnz; k++) {
+ DegreeOfFreedom row, col;
+ double value;
+
+ in >> row;
+ in >> col;
+ in >> value;
+
+ WorldVector<double> rowCoords = dofToCoord[row];
+ WorldVector<double> colCoords = dofToCoord[col];
+
+ if (coordToDof.find(rowCoords) == coordToDof.end()) {
+ std::cout << "CANNOT ROW FIND" << std::endl;
+ rowCoords.print();
+ exit(0);
+ }
+
+ if (coordToDof.find(colCoords) == coordToDof.end()) {
+ std::cout << "CANNOT COL FIND" << std::endl;
+ rowCoords.print();
+ exit(0);
+ }
+
+ DegreeOfFreedom rowHere = coordToDof[rowCoords];
+ DegreeOfFreedom colHere = coordToDof[colCoords];
+
+ std::pair<int, int> rowcol = make_pair(rowHere, colHere);
+ if (nnzValues[i][j].count(rowcol) == 0)
+ nnzValues[i][j][rowcol] = value;
+ else
+ nnzValues[i][j][rowcol] += value;
+ }
+ }
+ }
+
+ for (int i = 0; i < nComponents; i++) {
int readNnz;
in >> readNnz;
-
- // Read each entry in file and check it with the corresponding matrix value
- // in this problem.
for (int k = 0; k < readNnz; k++) {
- DegreeOfFreedom row, col;
+ int row;
double value;
-
+
in >> row;
- in >> col;
in >> value;
WorldVector<double> rowCoords = dofToCoord[row];
- WorldVector<double> colCoords = dofToCoord[col];
if (coordToDof.find(rowCoords) == coordToDof.end()) {
std::cout << "CANNOT ROW FIND" << std::endl;
@@ -1556,34 +1642,129 @@ namespace AMDiS {
exit(0);
}
- if (coordToDof.find(colCoords) == coordToDof.end()) {
- std::cout << "CANNOT COL FIND" << std::endl;
+ DegreeOfFreedom rowHere = coordToDof[rowCoords];
+
+ if (rhsValues[i].count(rowHere) == 0)
+ rhsValues[i][rowHere] = value;
+ else
+ rhsValues[i][rowHere] += value;
+ }
+ }
+
+ for (int i = 0; i < nComponents; i++) {
+ int readNnz;
+ in >> readNnz;
+ for (int k = 0; k < readNnz; k++) {
+ int row;
+ double value;
+
+ in >> row;
+ in >> value;
+
+ WorldVector<double> rowCoords = dofToCoord[row];
+
+ if (coordToDof.find(rowCoords) == coordToDof.end()) {
+ std::cout << "CANNOT ROW FIND" << std::endl;
rowCoords.print();
exit(0);
}
DegreeOfFreedom rowHere = coordToDof[rowCoords];
- DegreeOfFreedom colHere = coordToDof[colCoords];
- double valueHere = (dofmatrix->getBaseMatrix())[rowHere][colHere];
- double diff = fabs(valueHere - value);
+ if (solValues[i].count(rowHere) == 0) {
+ solValues[i][rowHere] = value;
+ } else {
+ double diff = fabs(solValues[i][rowHere] - value);
+ if (diff > 1e-8) {
+ std::cout << "DIFFERENT values in solution vector!" << std::endl;
+ exit(0);
+ }
+ }
+ }
+ }
+
+ in.close();
+
+ std::cout << "File read!" << std::endl;
+ }
+
+ std::cout.precision(15);
+
+ for (int i = 0; i < nComponents; i++) {
+ for (int j = 0; j < nComponents; j++) {
+ DOFMatrix *dofmatrix = (*systemMatrix)[i][j];
+
+ if (!dofmatrix)
+ continue;
+
+ for (std::map<std::pair<int, int>, double>::iterator nnzIt =
+ nnzValues[i][j].begin(); nnzIt != nnzValues[i][j].end(); ++nnzIt) {
+
+ int row = nnzIt->first.first;
+ int col = nnzIt->first.second;
+ double value = nnzIt->second;
+
+ double valueHere = (dofmatrix->getBaseMatrix())[row][col];
+ double diff = fabs(valueHere - value);
+
// And so we check the difference between the value read from file and the
// corresponding value in the matrix.
- if (diff > 1e-10) {
+ if (diff > 1e-8) {
std::cout << "WRONG value in matrix[" << i << "][" << j << "]!" << std::endl
- << " DOF in other matrix [" << row << "][" << col << "] = " << value << std::endl
- << " DOF in this matrix [" << rowHere << "][" << colHere << "] = " << valueHere << std::endl;
-
-
+ << " DOF in other matrix[" << row << "][" << col << "] = " << value << std::endl
+ << " DOF in this matrix[" << row << "][" << col << "] = " << valueHere << std::endl;
+
+
exit(0);
}
+ }
+ }
+ }
+ for (int i = 0; i < nComponents; i++) {
+ for (std::map<int, double>::iterator rhsIt = rhsValues[i].begin();
+ rhsIt != rhsValues[i].end(); ++rhsIt) {
+
+ int row = rhsIt->first;
+ double value = rhsIt->second;
+
+ double valueHere = (*(rhs->getDOFVector(i)))[row];
+ double diff = fabs(valueHere - value);
+
+ if (diff > 1e-8) {
+ std::cout << "WRONG value in rhs[" << i << "]!" << std::endl
+ << " DOF in other rhs[" << row << "] = " << value << std::endl
+ << " DOF in this rhs[" << row << "] = " << valueHere << std::endl;
+
+ exit(0);
+ }
+ }
+ }
+
+ for (int i = 0; i < nComponents; i++) {
+ for (std::map<int, double>::iterator solIt = solValues[i].begin();
+ solIt != solValues[i].end(); ++solIt) {
+
+ int row = solIt->first;
+ double value = solIt->second;
+
+ double valueHere = (*(solution->getDOFVector(i)))[row];
+ double diff = fabs(valueHere - value);
+
+ if (diff > 1e-8) {
+ std::cout << "WRONG value in sol[" << i << "]!" << std::endl
+ << " DOF in other sol[" << row << "] = " << value << std::endl
+ << " DOF in this sol[" << row << "] = " << valueHere << std::endl;
+
+ exit(0);
}
}
}
+ std::cout << "FINISHED COMPARING!" << std::endl;
- in.close();
+ exit(0);
}
+
}
diff --git a/AMDiS/src/ProblemVec.h b/AMDiS/src/ProblemVec.h
index f4081696..7a4de800 100644
--- a/AMDiS/src/ProblemVec.h
+++ b/AMDiS/src/ProblemVec.h
@@ -473,6 +473,21 @@ namespace AMDiS {
return fileWriters;
}
+ /** \brief
+ * This function writes the system matrix and the system vector to a file
+ * for debugging. The entries of both, the matrix and the vector, are not
+ * indicated by the dof indices, but by the world coordinates of the dofs.
+ * This makes it possible to compare dof matrices, where the dof indices
+ * are different, e.g., when using domain decomposition parallelization.
+ */
+ void writeDbgMatrix(std::string filename);
+
+ /** \brief
+ * Reads a file, that was created using the function \ref writeDbgMatrix, and
+ * compares the date with the system matrix of this problem.
+ */
+ void readAndCompareDbgMatrix(std::vector<std::string> filenames);
+
protected:
/** \brief
* If the exact solution is known, the problem can compute the exact
@@ -494,21 +509,6 @@ namespace AMDiS {
*/
void createCoordToDofMap(CoordToDof &dofMap);
- /** \brief
- * This function writes the system matrix and the system vector to a file
- * for debugging. The entries of both, the matrix and the vector, are not
- * indicated by the dof indices, but by the world coordinates of the dofs.
- * This makes it possible to compare dof matrices, where the dof indices
- * are different, e.g., when using domain decomposition parallelization.
- */
- void writeDbgMatrix(std::string filename);
-
- /** \brief
- * Reads a file, that was created using the function \ref writeDbgMatrix, and
- * compares the date with the system matrix of this problem.
- */
- void readAndCompareDbgMatrix(std::string filename);
-
protected:
/// Name of this problem.
diff --git a/AMDiS/src/SubAssembler.cc b/AMDiS/src/SubAssembler.cc
index ba8d3706..e616e470 100644
--- a/AMDiS/src/SubAssembler.cc
+++ b/AMDiS/src/SubAssembler.cc
@@ -179,11 +179,10 @@ namespace AMDiS {
if (opt && !quad && sameFESpaces) {
const BasisFunction *psi = owner->rowFESpace->getBasisFcts();
const BasisFunction *phi = owner->colFESpace->getBasisFcts();
- if (vec->getFESpace()->getBasisFcts() == psi) {
+ if (vec->getFESpace()->getBasisFcts() == psi)
psiFast = updateFastQuadrature(psiFast, psi, INIT_PHI);
- } else if(vec->getFESpace()->getBasisFcts() == phi) {
- phiFast = updateFastQuadrature(phiFast, phi, INIT_PHI);
- }
+ else if(vec->getFESpace()->getBasisFcts() == phi)
+ phiFast = updateFastQuadrature(phiFast, phi, INIT_PHI);
}
// calculate new values
@@ -200,7 +199,7 @@ namespace AMDiS {
} else {
vec->getVecAtQPs(elInfo, localQuad, NULL, values);
}
-
+
valuesAtQPs[vec]->valid = true;
return values;
diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc
index a8d6f083..c7975367 100644
--- a/AMDiS/src/ZeroOrderAssembler.cc
+++ b/AMDiS/src/ZeroOrderAssembler.cc
@@ -49,10 +49,8 @@ namespace AMDiS {
sort(assTerms.begin(), assTerms.end());
- if ((opTerms == assTerms) &&
- ((*subAssemblers)[i]->getQuadrature() == quad)) {
+ if (opTerms == assTerms && (*subAssemblers)[i]->getQuadrature() == quad)
return dynamic_cast<ZeroOrderAssembler*>((*subAssemblers)[i]);
- }
}
}
@@ -69,11 +67,10 @@ namespace AMDiS {
if (!optimized) {
newAssembler = new StandardZOA(op, assembler, quad);
} else {
- if (pwConst) {
- newAssembler = new PrecalcZOA(op, assembler, quad);
- } else {
- newAssembler = new FastQuadZOA(op, assembler, quad);
- }
+ if (pwConst)
+ newAssembler = new PrecalcZOA(op, assembler, quad);
+ else
+ newAssembler = new FastQuadZOA(op, assembler, quad);
}
subAssemblers->push_back(newAssembler);
@@ -95,9 +92,8 @@ namespace AMDiS {
int myRank = omp_get_thread_num();
std::vector<OperatorTerm*>::iterator termIt;
- for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
+ for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt)
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
- }
if (symmetric) {
TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
@@ -106,9 +102,8 @@ namespace AMDiS {
c[iq] *= elInfo->getDet();
// calculate phi at QPs only once!
- for (int i = 0; i < nCol; i++) {
+ for (int i = 0; i < nCol; i++)
phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
- }
for (int i = 0; i < nRow; i++) {
double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
@@ -125,15 +120,13 @@ namespace AMDiS {
c[iq] *= elInfo->getDet();
// calculate phi at QPs only once!
- for (int i = 0; i < nCol; i++) {
+ for (int i = 0; i < nCol; i++)
phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
- }
for (int i = 0; i < nRow; i++) {
double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
- for (int j = 0; j < nCol; j++) {
+ for (int j = 0; j < nCol; j++)
mat[i][j] += quadrature->getWeight(iq) * c[iq] * psival * phival[j];
- }
}
}
}
@@ -143,12 +136,11 @@ namespace AMDiS {
{
int nPoints = quadrature->getNumPoints();
std::vector<double> c(nPoints, 0.0);
-
int myRank = omp_get_thread_num();
+
std::vector<OperatorTerm*>::iterator termIt;
- for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
+ for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt)
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
- }
for (int iq = 0; iq < nPoints; iq++) {
c[iq] *= elInfo->getDet();
--
GitLab