diff --git a/AMDiS/src/MacroReader.cc b/AMDiS/src/MacroReader.cc
index a513b1278e01a746be1f707f88216ed4f2a2e332..b3f245aba4d2888c5e1156c35715ad4b4fcf44bf 100644
--- a/AMDiS/src/MacroReader.cc
+++ b/AMDiS/src/MacroReader.cc
@@ -475,7 +475,7 @@ namespace AMDiS {
FUNCNAME("MacroInfo::readAMDiSMacro");
FILE *file;
int dim;
- int dow, nv, ne, i, j, k;
+ int dow, nv, ne, j, k;
double dbl;
char name[128], line[256];
int line_no, n_keys, i_key, sort_key[N_KEYS], nv_key, ne_key;
@@ -600,15 +600,13 @@ namespace AMDiS {
break;
case 4:
fscanf(file, "%*s %*s");
- for (i = 0; i < nv; i++)
- {
- for (j = 0; j <Global::getGeo(WORLD) ; j++)
- {
- TEST_EXIT(fscanf(file, "%lf", &dbl) == 1)
- ("error while reading coordinates, check file %s\n", name);
- coords[i][j] = dbl;
- }
+ for (int i = 0; i < nv; i++) {
+ for (j = 0; j <Global::getGeo(WORLD) ; j++) {
+ TEST_EXIT(fscanf(file, "%lf", &dbl) == 1)
+ ("error while reading coordinates, check file %s\n", name);
+ coords[i][j] = dbl;
}
+ }
key_def[4] = true;
break;
case 5:
@@ -617,7 +615,7 @@ namespace AMDiS {
/* global number of vertices on a single element */
/****************************************************************************/
- for (i = 0; i < ne; i++) {
+ for (int i = 0; i < ne; i++) {
TEST_EXIT(read_indices(file, *ind))
("cannot read vertex indices of element %d in file %s\n", i, name);
@@ -632,8 +630,7 @@ namespace AMDiS {
/****************************************************************************/
/* MEL boundary pointers */
/****************************************************************************/
-
- for (i = 0; i < ne; i++) {
+ for (int i = 0; i < ne; i++) {
// boundary information of ith element
TEST_EXIT(read_indices(file, *ind))
@@ -658,7 +655,7 @@ namespace AMDiS {
/* else init them by a call of fill_mel_neighbour() */
/****************************************************************************/
neigh_set = true;
- for (i = 0; i < ne; i++) {
+ for (int i = 0; i < ne; i++) {
// neighbour information about ith element
if (read_indices(file, *ind))
@@ -682,7 +679,7 @@ namespace AMDiS {
if (dim == 2 || dim == 1)
ERROR("there is no element type in 2d and 2d; ignoring data for elType\n");
- for (i = 0; i < ne; i++) {
+ for (int i = 0; i < ne; i++) {
TEST_EXIT(fscanf(file, "%d", &j) == 1)
("cannot read elType of element %d in file %s\n",
i, name);
@@ -700,11 +697,11 @@ namespace AMDiS {
int numFaces = mesh->getGeo(FACE);
int numEdgesAtBoundary = 0;
- for(k = 1; k < dim; k++) {
+ for (k = 1; k < dim; k++) {
numEdgesAtBoundary += k;
}
- for (i = 0; i < ne; i++) {
+ for (int i = 0; i < ne; i++) {
TEST_EXIT(read_indices(file, *ind))
("cannot read boundary projector of element %d in file %s\n", i, name);
@@ -736,7 +733,7 @@ namespace AMDiS {
/* MEL regions */
/****************************************************************************/
- for (i = 0; i < ne; i++) {
+ for (int i = 0; i < ne; i++) {
TEST_EXIT(fscanf(file, "%d", &j) == 1)
("cannot read region of element %d in file %s\n", i, name);
if(j >= 0) {
@@ -751,7 +748,7 @@ namespace AMDiS {
break;
case 11:
fscanf(file, "%*s %*s");
- for (i = 0; i < ne; i++) {
+ for (int i = 0; i < ne; i++) {
TEST_EXIT(read_indices(file, *ind))
("cannot read surface regions of element %d in file %s\n", i, name);
@@ -807,10 +804,8 @@ namespace AMDiS {
void MacroInfo::dirichletBoundary()
{
- int i, k;
-
- for (i = 0; i < static_cast<int>( mel.size()); i++) {
- for (k = 0; k < mesh->getGeo(NEIGH); k++) {
+ for (int i = 0; i < static_cast<int>( mel.size()); i++) {
+ for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
if (mel[i]->neighbour[k])
mel[i]->boundary[k] = INTERIOR;
else
diff --git a/AMDiS/src/OEMSolver.h b/AMDiS/src/OEMSolver.h
index bc70a618dd02213835776673a0af9c8cd113a8e9..c938b21616c37bab41402c169dcb4ecd718a08bf 100644
--- a/AMDiS/src/OEMSolver.h
+++ b/AMDiS/src/OEMSolver.h
@@ -135,7 +135,9 @@ namespace AMDiS {
/** \brief
* Returns solvers \ref name.
*/
- inline const ::std::string& getName() { return name; };
+ inline const ::std::string& getName() {
+ return name;
+ };
/** \brief
* Returns \ref tolerance
@@ -158,10 +160,6 @@ namespace AMDiS {
return residual;
};
- // inline Preconditioner<VectorType> *getLeftPrecon() { return leftPrecon; };
-
- // inline Preconditioner<VectorType> *getRightPrecon() { return rightPrecon; };
-
/** \} */
// ===== setting-methods ======================================================
@@ -173,32 +171,28 @@ namespace AMDiS {
/** \brief
* Sets \ref tolerance
*/
- inline void setTolerance(double tol)
- {
+ inline void setTolerance(double tol) {
tolerance = tol;
};
/** \brief
* Sets \ref relative
*/
- inline void setRelative(bool rel)
- {
+ inline void setRelative(bool rel) {
relative = rel;
};
/** \brief
* Sets \ref max_iter
*/
- inline void setMaxIterations(int i)
- {
+ inline void setMaxIterations(int i) {
max_iter = i;
};
/** \brief
* Sets \ref info
*/
- inline void setInfo(int i)
- {
+ inline void setInfo(int i) {
info = i;
};
@@ -211,14 +205,6 @@ namespace AMDiS {
return vectorCreator;
};
- // inline void setLeftPrecon(Preconditioner<VectorType> *precon) {
- // leftPrecon = precon;
- // };
-
- // inline void setRightPrecon(Preconditioner<VectorType> *precon) {
- // rightPrecon = precon;
- // };
-
/** \} */
protected:
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 7b588164f31677d0614a6d085380060f176e20a5..b965e1c9c76fe0007fb83fff6bb018361ab53b4b 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -192,19 +192,19 @@ namespace AMDiS {
GET_PARAMETER(0, name_ + "->mesh", &meshName);
TEST_EXIT(meshName != "")("no mesh name spezified\n");
int dim = 0;
+ ::std::cout << name_ + "->dim" << ::std::endl;
GET_PARAMETER(0, name_ + "->dim", "%d", &dim);
TEST_EXIT(dim)("no problem dimension spezified!\n");
- int i, refSet;
- for(i = 0; i < numComponents_; i++) {
+ for (int i = 0; i < numComponents_; i++) {
sprintf(number, "%d", i);
- refSet = -1;
- GET_PARAMETER(0, name_ + "->refinement set[" + number + "]","%d", &refSet);
- if(refSet < 0) {
+ int refSet = -1;
+ GET_PARAMETER(0, name_ + "->refinement set[" + number + "]", "%d", &refSet);
+ if (refSet < 0) {
WARNING("refinement set for component %d not set\n", i);
refSet = 0;
}
- if(meshForRefinementSet[refSet] == NULL) {
+ if (meshForRefinementSet[refSet] == NULL) {
Mesh *newMesh = NEW Mesh(meshName, dim);
meshForRefinementSet[refSet] = newMesh;
meshes_.push_back(newMesh);
@@ -237,13 +237,13 @@ namespace AMDiS {
char number[3];
::std::map< ::std::pair<Mesh*, int>, FiniteElemSpace*> feSpaceMap;
- int i, dim = -1;
+ int dim = -1;
GET_PARAMETER(0, name_ + "->dim", "%d", &dim);
TEST_EXIT(dim != -1)("no problem dimension spezified!\n");
componentSpaces_.resize(numComponents_, NULL);
- for(i = 0; i < numComponents_; i++) {
+ for (int i = 0; i < numComponents_; i++) {
sprintf(number, "%d", i);
GET_PARAMETER(0, name_ + "->polynomial degree[" + number + "]","%d", °ree);
@@ -263,7 +263,7 @@ namespace AMDiS {
}
// create dof admin for vertex dofs if neccessary
- for (i = 0; i < static_cast<int>(meshes_.size()); i++) {
+ for (int i = 0; i < static_cast<int>(meshes_.size()); i++) {
if (meshes_[i]->getNumberOfDOFs(VERTEX) == 0) {
DimVec<int> ln_dof(meshes_[i]->getDim(), DEFAULT_VALUE, 0);
ln_dof[VERTEX]= 1;
@@ -892,13 +892,10 @@ namespace AMDiS {
PeriodicBC *periodic = new PeriodicBC(type, feSpace);
- // TEST_EXIT(mesh->getPeriodicBCMap()[type] == NULL)
- // ("periodic condition already set in mesh\n");
- // mesh->getPeriodicBCMap()[type] = periodic;
-
- if(systemMatrix_ && (*systemMatrix_)[row][col])
+ if (systemMatrix_ && (*systemMatrix_)[row][col])
(*systemMatrix_)[row][col]->getBoundaryManager()->addBoundaryCondition(periodic);
- if(rhs_)
+
+ if (rhs_)
rhs_->getDOFVector(row)->getBoundaryManager()->
addBoundaryCondition(periodic);
}
diff --git a/AMDiS/src/UmfPackSolver.h b/AMDiS/src/UmfPackSolver.h
index 19282bc187e32dbbcb2f920e4ff6967d6a292189..7b1bd476c4dc1f90feb047649135a022663419e7 100644
--- a/AMDiS/src/UmfPackSolver.h
+++ b/AMDiS/src/UmfPackSolver.h
@@ -102,6 +102,20 @@ namespace AMDiS {
* These vectors are justed to calculate the final residual of the solution.
*/
VectorType *r, *p;
+
+ /** \brief
+ * Stores the result of umfpack_di_symbolic (the symbolic analysis which is needed
+ * for the numeric factorization).
+ */
+ void *symbolic_;
+
+ /** \brief
+ * If the symbolic analysis should be done only once (for example, if the matrices
+ * to solve have all the same pattern because of no adaptivity), this variable
+ * is one, otherwise 0 and the symbolic analysis will be performed at every call
+ * of solveSystem().
+ */
+ int storeSymbolic_;
};
}
diff --git a/AMDiS/src/UmfPackSolver.hh b/AMDiS/src/UmfPackSolver.hh
index 9fa5313bf94bd16c8bcfc8f4be4384f422252262..f1d7f0cdadf629b2a7c232691b41f3e6eb8a11f0 100644
--- a/AMDiS/src/UmfPackSolver.hh
+++ b/AMDiS/src/UmfPackSolver.hh
@@ -11,11 +11,20 @@ namespace AMDiS {
template<typename VectorType>
UmfPackSolver<VectorType>::UmfPackSolver(::std::string name)
- : OEMSolver<VectorType>(name)
- {}
+ : OEMSolver<VectorType>(name),
+ symbolic_(NULL),
+ storeSymbolic_(0)
+ {
+ GET_PARAMETER(0, name + "->store symbolic", "%d", &storeSymbolic_);
+ }
template<typename VectorType>
- UmfPackSolver<VectorType>::~UmfPackSolver() {}
+ UmfPackSolver<VectorType>::~UmfPackSolver()
+ {
+ if (storeSymbolic_ && (symbolic_ != NULL)) {
+ umfpack_di_free_symbolic(&symbolic_);
+ }
+ }
template<typename VectorType>
int UmfPackSolver<VectorType>::solveSystem(MatVecMultiplier<VectorType> *matVec,
@@ -105,7 +114,7 @@ namespace AMDiS {
}
- void *Symbolic, *Numeric;
+ void *numeric;
double Control[UMFPACK_CONTROL];
double Info[UMFPACK_INFO];
int status = 0;
@@ -116,24 +125,32 @@ namespace AMDiS {
umfpack_di_defaults(Control);
// Run UMFPack.
- status = umfpack_di_symbolic(newMatrixSize, newMatrixSize, Ap, Ai, Ax, &Symbolic, Control, Info);
- if (!status == UMFPACK_OK) {
- ERROR_EXIT("UMFPACK Error in function umfpack_di_symbolic");
- }
-
- status = umfpack_di_numeric(Ap, Ai, Ax, Symbolic, &Numeric, Control, Info);
- if (!status == UMFPACK_OK) {
- ERROR_EXIT("UMFPACK Error in function umfpack_di_numberic");
+ if (!(storeSymbolic_) ||
+ (storeSymbolic_ && (symbolic_ == NULL))) {
+
+ status = umfpack_di_symbolic(newMatrixSize, newMatrixSize, Ap, Ai, Ax, &symbolic_, Control, Info);
+ if (!status == UMFPACK_OK) {
+ ERROR_EXIT("UMFPACK Error in function umfpack_di_symbolic");
+ }
+ ::std::cout << "symbolic durchgefuehrt" << ::std::endl;
+ } else {
+ ::std::cout << "symbolic ausgelassen" << ::std::endl;
}
- status = umfpack_di_solve(UMFPACK_A, Ap, Ai, Ax, xvec, bvec, Numeric, Control, Info);
+ status = umfpack_di_numeric(Ap, Ai, Ax, symbolic_, &numeric, Control, Info);
+ if (!status == UMFPACK_OK) {
+ ERROR_EXIT("UMFPACK Error in function umfpack_di_numeric");
+ }
+
+ status = umfpack_di_solve(UMFPACK_A, Ap, Ai, Ax, xvec, bvec, numeric, Control, Info);
if (!status == UMFPACK_OK) {
ERROR_EXIT("UMFPACK Error in function umfpack_di_solve");
}
-
-
- umfpack_di_free_symbolic(&Symbolic);
- umfpack_di_free_numeric(&Numeric);
+
+ if (!storeSymbolic_) {
+ umfpack_di_free_symbolic(&symbolic_);
+ }
+ umfpack_di_free_numeric(&numeric);
// Copy and resort solution.