Commit bc5b945f authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Serveral bug fixes for parallelization.

parent 87712cf8
......@@ -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"
......
......@@ -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) {
......
......@@ -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
......
......@@ -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]];
......
......@@ -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);
......
......@@ -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();
......
......@@ -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;
......@@ -586,6 +586,9 @@ namespace AMDiS {
delete [] sendBuffers[i];
MatDestroy(petscMatrix);
VecDestroy(petscRhsVec);
VecDestroy(petscSolVec);
VecDestroy(petscTmpVec);
}
......@@ -611,16 +614,9 @@ 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;
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);
......@@ -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.
......
......@@ -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,
......@@ -1492,11 +1500,39 @@ namespace AMDiS {
}
}
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,12 +1549,26 @@ namespace AMDiS {
DegreeOfFreedom dof;
int dimOfWorld = Global::getGeo(WORLD);
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;
nnzValues.resize(nComponents);
rhsValues.resize(nComponents);
solValues.resize(nComponents);
for (int i = 0; i < nComponents; i++)
nnzValues[i].resize(nComponents);
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(filename.c_str());
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++)
......@@ -1564,26 +1614,157 @@ namespace AMDiS {
DegreeOfFreedom rowHere = coordToDof[rowCoords];
DegreeOfFreedom colHere = coordToDof[colCoords];
double valueHere = (dofmatrix->getBaseMatrix())[rowHere][colHere];
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;
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];
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];
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);
}
}
......@@ -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.