Liebe Gitlab-Nutzer, lieber Gitlab-Nutzer, es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Ein Anmelden über dieses erzeugt ein neues Konto. Das alte Konto ist über den Reiter "Standard" erreichbar. Die Administratoren

Dear Gitlab user, it is now possible to log in to our service using the ZIH login/LDAP. Logging in via this will create a new account. The old account can be accessed via the "Standard" tab. The administrators

Commit 2aff8c3b authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

some corrections in the solvers

parent 0b5440f3
......@@ -255,6 +255,9 @@ int i ;
void *Symbolic, *Numeric ;
double stats [2];
umfpack_tic (stats);
(void) umfpack_di_symbolic (n, n, Ap, Ai, Ax, &Symbolic, null, null) ;
(void) umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, null, null) ;
......@@ -265,6 +268,7 @@ umfpack_di_free_symbolic (&Symbolic) ;
umfpack_di_free_numeric (&Numeric) ;
umfpack_toc (stats);
return (0) ;
......@@ -276,18 +280,30 @@ return (0) ;
OUTPUT_VARIABLE CHOLMOD_OUT)
if(NOT CHOLMOD_TEST)
find_library(CHOLMOD_LIB cholmod)
find_library(COLAMD_LIB colamd)
if(CHOLMOD_LIB)
list(APPEND AMDIS_LIBRARIES ${CHOLMOD_LIB})
else()
message(FATAL_ERROR "your umfpack seems to need cholmod, but cmake could not find it")
endif()
find_library(COLAMD_LIB colamd)
if(COLAMD_LIB)
list(APPEND AMDIS_LIBRARIES ${COLAMD_LIB})
else()
message(FATAL_ERROR "your umfpack seems to need colamd, but cmake could not find it")
endif()
try_compile(CHOLMOD_TEST2 ${_CHOLMOD_TEST_DIR}/build ${_CHOLMOD_TEST_DIR} cholmodTest
OUTPUT_VARIABLE CHOLMOD_OUT)
if(NOT CHOLMOD_TEST2)
find_library(SUITESPARSECONFIG_LIB libsuitesparseconfig)
if(SUITESPARSECONFIG_LIB)
list(APPEND AMDIS_LIBRARIES ${SUITESPARSECONFIG_LIB})
else()
message(FATAL_ERROR "your umfpack seems to need suitesparseconfig, but cmake could not find it")
endif()
endif()
endif()
endif(AMDIS_NEED_UMFPACK)
......
......@@ -53,8 +53,11 @@ elseif(AMDIS_NEED_SEQ_PETSC)
endif()
if(AMDIS_NEED_HYPRE)
find_package(HYPRE REQUIRED)
find_package(HYPRE REQUIRED HINTS ${AMDIS_DIR})
if(HAVE_HYPRE)
list(APPEND AMDIS_INCLUDE_DIRS ${HYPRE_INCLUDE_DIRECTORIES})
list(APPEND AMDIS_COMPILEFLAGS "-DMTL_HAS_HYPRE")
list(APPEND AMDIS_LIBRARIES ${HYPRE_LIBRARIES} )
if(NOT MPI_FOUND)
find_package(MPI REQUIRED)
if(MPI_FOUND)
......@@ -63,9 +66,6 @@ if(AMDIS_NEED_HYPRE)
list(APPEND AMDIS_INCLUDE_DIRS ${MPI_INCLUDE_PATH})
endif()
endif()
list(APPEND AMDIS_INCLUDE_DIRS ${HYPRE_INCLUDE_DIRECTORIES})
list(APPEND AMDIS_COMPILEFLAGS "-DMTL_HAS_HYPRE")
list(APPEND AMDIS_LIBRARIES ${HYPRE_LIBRARIES} )
endif()
endif(AMDIS_NEED_HYPRE)
......
......@@ -340,10 +340,10 @@ endif(ENABLE_UMFPACK)
if(ENABLE_HYPRE)
include(FindHYPRE.cmake)
message("have hypre: ${HAVE_HYPRE}")
if(HAVE_HYPRE)
if(NOT MPI_FOUND)
if (HAVE_HYPRE)
if (NOT MPI_FOUND)
find_package(MPI REQUIRED)
if(MPI_FOUND)
if (MPI_FOUND)
list(APPEND COMPILEFLAGS "${MPI_COMPILE_FLAGS}")
include_directories(${MPI_INCLUDE_PATH})
endif()
......
......@@ -48,7 +48,7 @@ namespace AMDiS {
return MPI::Wtime() - first_mpi;
#else
time_duration td = microsec_clock::local_time()-first_seq;
return static_cast<double>(td.total_milliseconds())*1.e-3;
return static_cast<double>(td.total_microseconds())*1.e-6;
#endif
}
......
......@@ -29,6 +29,8 @@
#include "parallel/PetscSolver.h"
#include "parallel/StdMpi.h"
#include "tools.h"
namespace AMDiS { namespace Parallel {
PetscSolver::PetscSolver(std::string name)
......@@ -86,10 +88,10 @@ namespace AMDiS { namespace Parallel {
MPI::COMM_WORLD.Barrier();
Timer t;
if (createMatrixData)
fillPetscMatrix(const_cast< Matrix< DOFMatrix* >* >(A.getOriginalMat()));
fillPetscMatrix(const_cast< Matrix< DOFMatrix* >* >(A.getOriginalMat()));
fillPetscRhs(&b);
INFO(info, 8)("creation of parallel data structures needed %.5f seconds\n",
......
......@@ -66,6 +66,8 @@ namespace AMDiS { namespace Parallel {
Parameters::get(name + "->pc_type", precon);
if (!precon.size())
Parameters::get(name + "->left precon", precon);
if (!precon.size())
Parameters::get(name + "->right precon", precon);
if (!matSolverPackage && params.emptyParam.find(precon) == params.emptyParam.end()) {
precon = (params.preconMap.find(precon) != params.preconMap.end() ? params.preconMap[precon] : precon);
PetscOptionsInsertString(("-" + kspPrefix + "pc_type " + precon).c_str());
......
......@@ -24,6 +24,8 @@
#include "parallel/PetscHelper.h"
#include "TransformDOF.h"
#include "tools.h"
namespace AMDiS { namespace Parallel {
using namespace std;
......@@ -114,6 +116,7 @@ namespace AMDiS { namespace Parallel {
{
FUNCNAME("PetscSolverNavierStokes::initSolver()");
printMem("begin(initSolver)");
// Create FGMRES based outer solver
MSG("CREATE POS 1: %p\n", &ksp);
......@@ -128,6 +131,8 @@ namespace AMDiS { namespace Parallel {
// Create null space information.
if (pressureNullSpace)
setConstantNullSpace(ksp, pressureComponent, true);
printMem("end(initSolver)");
}
......@@ -140,6 +145,9 @@ namespace AMDiS { namespace Parallel {
TEST_EXIT(nu)("nu pointer not set!\n");
TEST_EXIT(invTau)("invtau pointer not set!\n");
TEST_EXIT(solution)("solution pointer not set!\n");
printMem("begin(initPreconditioner)");
int dim = componentSpaces[pressureComponent]->getMesh()->getDim();
......@@ -356,6 +364,9 @@ namespace AMDiS { namespace Parallel {
}
setConstantNullSpace(matShellContext.kspLaplace);
printMem("end(initPreconditioner)");
MSG("Setup of Navier-Stokes preconditioner needed %.5f seconds\n",
t.elapsed());
......
......@@ -89,11 +89,17 @@ namespace AMDiS {
TEST_EXIT(preconPair.l != nullptr)("there is no left preconditioner\n");
TEST_EXIT(preconPair.r != nullptr)("there is no right preconditioner\n");
typedef typename mtl::Collection<MatrixType>::value_type value_type;
VectorType r(A * x - b);
VectorType r(num_rows(b)); r = b;
if (two_norm(x) != 0) {
r = A * x ;
r -= b;
}
int error = 0;
if (oem.getInfo() == 0) {
itl::basic_iteration<typename MatrixType::value_type>
itl::basic_iteration<value_type>
iter(r, oem.getMaxIterations(), oem.getRelative(), oem.getTolerance());
error = solver(A, x, b, *(preconPair.l), *(preconPair.r), iter);
......@@ -101,7 +107,7 @@ namespace AMDiS {
oem.setIterations(iter.iterations());
oem.setResidual(iter.resid());
} else {
itl::cyclic_iteration<typename MatrixType::value_type>
itl::cyclic_iteration<value_type>
iter(r, oem.getMaxIterations(), oem.getRelative(), oem.getTolerance(),
oem.getPrint_cycle());
......@@ -122,29 +128,31 @@ namespace AMDiS {
TEST_EXIT(preconPair.l != nullptr)("there is no left preconditioner\n");
TEST_EXIT(preconPair.r != nullptr)("there is no right preconditioner\n");
mtl::matrix::transposed_view<const MatrixType> B(A);
VectorType r(B * x - b);
int error = 0;
if (oem.getInfo() == 0) {
itl::basic_iteration<typename MatrixType::value_type>
iter(r, oem.getMaxIterations(), oem.getRelative(), oem.getTolerance());
error = solver(B, x, b, *(preconPair.l), *(preconPair.r), iter);
oem.setErrorCode(error);
oem.setIterations(iter.iterations());
oem.setResidual(iter.resid());
} else {
itl::cyclic_iteration<typename MatrixType::value_type>
iter(r, oem.getMaxIterations(), oem.getRelative(), oem.getTolerance(),
oem.getPrint_cycle());
error = solver(B, x, b, *(preconPair.l), *(preconPair.r), iter);
oem.setErrorCode(error);
oem.setIterations(iter.iterations());
oem.setResidual(iter.resid());
}
return error;
// typedef typename mtl::Collection<MatrixType>::value_type value_type;
// mtl::matrix::transposed_view<const MatrixType> B(A);
// VectorType r(B * x - b);
// int error = 0;
// if (oem.getInfo() == 0) {
// itl::basic_iteration<value_type>
// iter(r, oem.getMaxIterations(), oem.getRelative(), oem.getTolerance());
//
// error = solver(B, x, b, *(preconPair.l), *(preconPair.r), iter);
// oem.setErrorCode(error);
// oem.setIterations(iter.iterations());
// oem.setResidual(iter.resid());
// } else {
// itl::cyclic_iteration<value_type>
// iter(r, oem.getMaxIterations(), oem.getRelative(), oem.getTolerance(),
// oem.getPrint_cycle());
//
// error = solver(B, x, b, *(preconPair.l), *(preconPair.r), iter);
// oem.setErrorCode(error);
// oem.setIterations(iter.iterations());
// oem.setResidual(iter.resid());
// }
//
// return error;
return 1;
}
......
......@@ -43,8 +43,8 @@ namespace AMDiS {
/// Constructor
MTL4SolverBase(LinearSolver* oem_)
: runner(oem_),
oem(*oem_),
matrix(0,0)
oem(*oem_)/*,
matrix(0,0)*/
{}
protected:
......
......@@ -56,8 +56,9 @@ namespace AMDiS {
mat(mat), mapper(m) {}
};
template< typename MTLMatrix, typename Mapper >
void operator<<(MTLMatrix& matrix, MatMap<const SolverMatrix<Matrix<DOFMatrix* > >, Mapper >& Asolver)
// requires MatrixType to have an inserter
template< typename MatrixType, typename Mapper >
void operator<<(MatrixType& matrix, MatMap<const SolverMatrix<Matrix<DOFMatrix* > >, Mapper >& Asolver)
{
const Matrix< DOFMatrix* >& A = *(Asolver.mat.getOriginalMat());
int ns = A.getSize();
......@@ -72,7 +73,7 @@ namespace AMDiS {
nnz += A[rb][cb]->getBaseMatrix().nnz();
{
typedef mtl::matrix::mapped_inserter< typename Collection< MTLMatrix >::Inserter, Mapper > Inserter;
typedef mtl::matrix::mapped_inserter< typename Collection< MatrixType >::Inserter, Mapper > Inserter;
Inserter ins(matrix, mapper, int(1.2 * nnz / matrix.dim1()));
for (int rb = 0; rb < ns; ++rb) {
mapper.setRow(rb);
......@@ -84,6 +85,37 @@ namespace AMDiS {
}
}
}
// requires MatrixType to have an inserter
template< typename MatrixType, typename Mapper >
void operator<<(MatrixType& matrix, MatMap<Matrix<MatrixType* >, Mapper >& Asolver)
{
Matrix< MatrixType* >& A = *(Asolver.mat);
int ns = A.getSize();
Mapper& mapper(Asolver.mapper);
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]->nnz();
{
typedef mtl::matrix::mapped_inserter< typename Collection< MatrixType >::Inserter, Mapper > Inserter;
Inserter ins(matrix, mapper, int(1.2 * nnz / matrix.dim1()));
for (int rb = 0; rb < ns; ++rb) {
mapper.setRow(rb);
for (int cb = 0; cb < ns; ++cb) {
mapper.setCol(cb);
if (A[rb][cb])
ins << (*A[rb][cb]);
}
}
}
}
template< typename Vector, typename CurMap >
void operator<<(Vector& dest, VecMap<DOFVector<double >, CurMap >& rhs)
......@@ -185,31 +217,31 @@ namespace AMDiS {
/// create and fill \ref PetscMatrix
inline void operator<<(PetscMatrix& mat, const SolverMatrix<Matrix<DOFMatrix*> >& Asolver)
{
const Matrix< DOFMatrix* >& A = *(Asolver.getOriginalMat());
if(mat.isNested) {
std::vector<PetscInt> nnz;
mat.nestMat.resize(A.getNumRows() * A.getNumCols());
for (size_t i = 0; i < static_cast<size_t>(A.getNumRows()); i++) {
for (size_t j = 0; j < static_cast<size_t>(A.getNumCols()); j++) {
size_t idx = i * A.getNumCols() + j;
if (A[i][j] == nullptr
|| num_rows(A[i][j]->getBaseMatrix()) == 0
|| num_cols(A[i][j]->getBaseMatrix()) == 0
|| A[i][j]->getBaseMatrix().nnz() == 0) {
mat.nestMat[idx] = PETSC_NULL;
continue;
}
const Matrix< DOFMatrix* >& A = *(Asolver.getOriginalMat());
if (mat.isNested) {
mat.nestMat[idx] << *(A[i][j]);
}
}
std::vector<PetscInt> nnz;
mat.nestMat.resize(A.getNumRows() * A.getNumCols());
for (size_t i = 0; i < static_cast<size_t>(A.getNumRows()); i++) {
for (size_t j = 0; j < static_cast<size_t>(A.getNumCols()); j++) {
size_t idx = i * A.getNumCols() + j;
if (A[i][j] == nullptr
|| num_rows(A[i][j]->getBaseMatrix()) == 0
|| num_cols(A[i][j]->getBaseMatrix()) == 0
|| A[i][j]->getBaseMatrix().nnz() == 0) {
mat.nestMat[idx] = PETSC_NULL;
continue;
}
// create nested matrix from a vector of block matrices
MatCreateNest(PETSC_COMM_SELF, A.getNumRows(), PETSC_NULL, A.getNumCols(), PETSC_NULL, &(mat.nestMat[0]), &mat.matrix);
}else {
mat.nestMat[idx] << *(A[i][j]);
}
}
// create nested matrix from a vector of block matrices
MatCreateNest(PETSC_COMM_SELF, A.getNumRows(), PETSC_NULL, A.getNumCols(), PETSC_NULL, &(mat.nestMat[0]), &mat.matrix);
} else {
bool initMatrix = false;
unsigned nRows(0);
unsigned nEntries(0);
......
......@@ -63,6 +63,8 @@ namespace AMDiS {
Parameters::get(oem.getName() + "->pc_type", precon);
if (!precon.size())
Parameters::get(oem.getName() + "->left precon", precon);
if (!precon.size())
Parameters::get(oem.getName() + "->right precon", precon);
if (!matSolverPackage && !params.emptyParam[precon]) {
precon = (params.preconMap.find(precon) != params.preconMap.end() ? params.preconMap[precon] : precon);
PetscOptionsInsertString(("-" + kspPrefix + "pc_type " + precon).c_str());
......@@ -138,6 +140,7 @@ namespace AMDiS {
PC pc_;
KSPGetPC(ksp_, &pc_);
PCSetType(pc_, pcType);
PCSetFromOptions(pc_);
}
}
......
......@@ -141,21 +141,23 @@ namespace AMDiS {
return &runner;
}
void setNestedVectors(bool createNestedVector = true)
void setNestedVectors(bool nested = true)
{
vecSol.isNested= createNestedVector;
vecRhs.isNested= createNestedVector;
vecSol.isNested = nested;
vecRhs.isNested = nested;
}
void setNestedMatrix(bool nested=true)
void setNestedMatrix(bool nested = true)
{
petscMat.isNested=nested;
petscMat.isNested=nested;
}
void setNested(bool n)
{
setNestedVectors(n); setNestedMatrix(n);
setNestedVectors(n);
setNestedMatrix(n);
}
protected:
Runner runner;
......@@ -170,10 +172,7 @@ namespace AMDiS {
// transfer matrix and rhs-vector to PETSc data-structures
if (createMatrixData) {
petscMat << A;
if (x.getSize() == 1)
runner.init(A, petscMat.nestMat[0]);
else
runner.init(A, petscMat.matrix);
runner.init(A, petscMat.matrix);
}
vecSol << x;
......@@ -182,10 +181,7 @@ namespace AMDiS {
// solve the linear system using PETSc solvers
t.reset();
if (vecSol.nestVec.size() == 1)
error = runner.solve(petscMat.nestMat[0], vecSol.nestVec[0], vecRhs.nestVec[0]);
else
error = runner.solve(petscMat.matrix, vecSol.vector, vecRhs.vector);
error = runner.solve(petscMat.matrix, vecSol.vector, vecRhs.vector);
// transfer solution from PETSc vector to SystemVector
vecSol.vector >> x;
......
......@@ -31,7 +31,23 @@ namespace itl {
namespace details {
/// Compute the Givens rotation matrix parameters for a and b.
inline void rotmat( double& a, double& b , double& c, double& s)
// template<typename T>
// void rotmat(const T& a, const T& b , T& c, T& s)
// {
// using std::abs; using std::sqrt; using mtl::conj;
//
// const T zero = math::zero(T());
// if (a == zero) {
// c = 0.0;
// s = 1.0;
// } else {
// double temp = abs(a) / sqrt( conj(a)*a + conj(b)*b );
// c = temp;
// s = temp * (b / a);
// }
// }
inline void rotmat(const double& a, const double& b , double& c, double& s)
{
using std::abs; using std::sqrt;
if ( b == 0.0 ) {
......
......@@ -52,28 +52,30 @@ int fgmres_full(const Matrix &A, Vector &x, const Vector &b,
if (size(b) == 0) throw mtl::logic_error("empty rhs vector");
const Scalar zero= math::zero(Scalar()), dbl_tol= 1.e-16, kappa = 10.0;
const Scalar zero= math::zero(Scalar()), breakdown_tol= 1.e-16, kappa = 10.0;
Scalar rho, bnrm2, temp, hr;
Size k, kmax(std::min(size(x), Size(iter.max_iterations() - iter.iterations())));
Vector w(solve(P, b)), r(b - A*x);
Vector w(resource(x)), r(b - A*x);
mtl::matrix::multi_vector<Vector> V(Vector(resource(x), zero), kmax+1);
mtl::matrix::multi_vector<Vector> Z(Vector(resource(x), zero), kmax+1);
mtl::matrix::dense2D<Scalar> H(kmax+1, kmax);
mtl::vector::dense_vector<Scalar> sn(kmax, zero), cs(kmax, zero), s(kmax+1, zero), y(kmax, zero); // replicated in distributed solvers
bnrm2 = two_norm(w);
if (bnrm2 < dbl_tol)
bnrm2 = 1.0;
bnrm2 = two_norm(b);
if (bnrm2 == zero) {
set_to_zero(x);
// b == 0 => solution = 0
return iter.terminate(bnrm2);
}
temp = two_norm(r);
rho= temp * reciprocal(bnrm2);
if (iter.finished(rho))
rho = two_norm(r); // norm of preconditioned residual
if (iter.finished(rho)) // initial guess is good enough solution
return iter;
V.vector(0) = r * reciprocal(temp);
V.vector(0) = r * reciprocal(rho);
H = zero;
s[0] = temp;
s[0] = rho;
// FGMRES iteration
for (k= 0; k < kmax && !iter.finished(rho) ; ++k, ++iter) {
......@@ -104,7 +106,7 @@ int fgmres_full(const Matrix &A, Vector &x, const Vector &b,
}
}
if (H[k+1][k] < dbl_tol)
if (H[k+1][k] < breakdown_tol)
return iter.fail(2, "FGMRES: Singular matrix - nearly hard breakdown");
V.vector(k+1) = w * reciprocal(H[k+1][k]);
......@@ -123,11 +125,11 @@ int fgmres_full(const Matrix &A, Vector &x, const Vector &b,
H[k][k] = cs[k]*H[k][k] + sn[k]*H[k+1][k];
H[k+1][k] = 0.0;
rho = std::abs(s[k+1]) / bnrm2;
rho = std::abs(s[k+1]);
}
// reduce k, to get regular matrix
while (k > 0 && abs(s[k-1]) <= iter.atol()) k--;
// while (k > 0 && abs(s[k-1]) <= iter.atol()) k--;
// iteration is finished -> compute x: solve H*y=g as far as rank of H allows
irange range(k);
......@@ -146,8 +148,8 @@ int fgmres_full(const Matrix &A, Vector &x, const Vector &b,
x += Z.vector(range) * y[range];
r = b - A*x;
return iter.terminate(r);
// r = b - A*x;
return iter.terminate(rho);
}
/// Flexible Generalized Minimal Residual method with restart
......
......@@ -52,7 +52,7 @@ int fgmres_householder_full(const Matrix &A, Vector &x, const Vector &b,
if (size(b) == 0) throw mtl::logic_error("empty rhs vector");
const Scalar zero= math::zero(Scalar()), dbl_tol= 1.e-16;
const Scalar zero= math::zero(Scalar());
Scalar rho, bnrm2, temp, alpha, beta;
Size k, kmax(std::min(size(x), Size(iter.max_iterations() - iter.iterations())));
Vector w(solve(R, b)), r(b - A*x);
......@@ -61,20 +61,19 @@ int fgmres_householder_full(const Matrix &A, Vector &x, const Vector &b,
mtl::vector::dense_vector<Scalar> sn(kmax, zero), cs(kmax, zero), s(kmax+1, zero), y(kmax, zero); // replicated in distributed solvers
mtl::matrix::dense2D<Scalar> H(kmax, kmax);
bnrm2 = two_norm(b);
if (bnrm2 == zero) {
set_to_zero(x);
// b == 0 => solution = 0
return iter.terminate(bnrm2);
}
bnrm2 = two_norm(w);
if (bnrm2 < dbl_tol)
bnrm2 = 1.0;
// TODO: if rhs=0 => solution=0 ??
temp = two_norm(r); // norm of preconditioned residual
rho = temp * reciprocal(bnrm2);
rho = two_norm(r); // norm of preconditioned residual
if (iter.finished(rho)) // initial guess is good enough solution
return iter;
// u = r + sing(r(0))*||r||*e0
beta = signum(r[0])*temp;
beta = signum(r[0])*rho;
w = r;
w[0] += beta;
w *= reciprocal(two_norm(w));
......@@ -140,7 +139,7 @@ int fgmres_householder_full(const Matrix &A, Vector &x, const Vector &b,
irange range(num_rows(H));
H[iall][k] = w[range];
rho = std::abs(s[k+1]) / bnrm2;
rho = std::abs(s[k+1]);
}