From a9c641bc3d879c49a62e76c84f18a9982a9e8b51 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Tue, 10 Jun 2008 10:36:47 +0000
Subject: [PATCH] * Further work on PardisoSolver
---
AMDiS/bin/Makefile.am | 2 +-
AMDiS/bin/Makefile.in | 2 +-
AMDiS/src/PardisoSolver.cc | 149 +++++++++++++++++++++++++++++++++++++
3 files changed, 151 insertions(+), 2 deletions(-)
diff --git a/AMDiS/bin/Makefile.am b/AMDiS/bin/Makefile.am
index 5fed9ecd..ffa328d0 100644
--- a/AMDiS/bin/Makefile.am
+++ b/AMDiS/bin/Makefile.am
@@ -34,7 +34,7 @@ if ENABLE_UMFPACK
endif
if ENABLE_MKL
- libamdis_la_CXXFLAGS += -DHAVE_MKL=1
+ libamdis_la_CXXFLAGS += -DHAVE_MKL=1 -I${MKL_INC}
endif
INCLUDES = $(AMDIS_INCLUDES) $(PARALLEL_INCLUDES)
diff --git a/AMDiS/bin/Makefile.in b/AMDiS/bin/Makefile.in
index 0b37740b..6ea1ff01 100644
--- a/AMDiS/bin/Makefile.in
+++ b/AMDiS/bin/Makefile.in
@@ -41,7 +41,7 @@ host_triplet = @host@
@ENABLE_UMFPACK_TRUE@ -I$(LIB_DIR)/AMD/Include \
@ENABLE_UMFPACK_TRUE@ -I$(LIB_DIR)/UMFPACK/Include
-@ENABLE_MKL_TRUE@am__append_3 = -DHAVE_MKL=1
+@ENABLE_MKL_TRUE@am__append_3 = -DHAVE_MKL=1 -I${MKL_INC}
@AMDIS_DEBUG_TRUE@am__append_4 = -g -O0 -Wall -DDEBUG=1 $(OPENMP_FLAG) -ftemplate-depth-30 $(INCLUDES) #-pedantic
@AMDIS_DEBUG_FALSE@am__append_5 = -O2 -Wall -DDEBUG=0 $(OPENMP_FLAG) -ftemplate-depth-30 $(INCLUDES) #-pedantic
subdir = bin
diff --git a/AMDiS/src/PardisoSolver.cc b/AMDiS/src/PardisoSolver.cc
index c6597e99..98f749ca 100644
--- a/AMDiS/src/PardisoSolver.cc
+++ b/AMDiS/src/PardisoSolver.cc
@@ -4,12 +4,161 @@
#ifdef HAVE_MKL
+#include <mkl_pardiso.h>
+
namespace AMDiS {
template<>
int PardisoSolver<SystemVector>::solveSystem(MatVecMultiplier<SystemVector> *matVec,
SystemVector *x, SystemVector *b)
{
+ FUNCNAME("PardisoSolver::solveSystem()");
+
+ TEST_EXIT(x->getSize() == b->getSize())("Vectors x and b must have the same size!");
+
+ // Extract the matrix of DOF-matrices.
+ StandardMatVec<Matrix<DOFMatrix*>, SystemVector> *stdMatVec =
+ dynamic_cast<StandardMatVec<Matrix<DOFMatrix*>, SystemVector> *>(matVec);
+ Matrix<DOFMatrix*> *m = stdMatVec->getMatrix();
+
+ // Number of systems.
+ int nComponents = m->getSize();
+ // Size of the new composed matrix.
+ int newMatrixSize = ((*m)[0][0])->getFESpace()->getAdmin()->getUsedSize() * nComponents;
+
+ // The new matrix has to be stored in compressed row format, therefore
+ // the rows are collected.
+ ::std::vector< ::std::vector< MatEntry > > rows(newMatrixSize, ::std::vector<MatEntry>(0));
+
+ // Counter for the number of non-zero elements in the new matrix.
+ int nElements = 0;
+
+ for (int stencilRow = 0; stencilRow < nComponents; stencilRow++) {
+ for (int stencilCol = 0; stencilCol < nComponents; stencilCol++) {
+
+ if (!(*m)[stencilRow][stencilCol]) {
+ continue;
+ }
+
+ DOFMatrix::Iterator matrixRow((*m)[stencilRow][stencilCol], USED_DOFS);
+ int rowIndex = 0;
+ for (matrixRow.reset(); !matrixRow.end(); matrixRow++, rowIndex++) {
+ for (int i = 0; i < static_cast<int>((*matrixRow).size()); i++) {
+ if ((*matrixRow)[i].col >= 0) {
+ MatEntry me;
+ me.entry = (*matrixRow)[i].entry;
+ // The col field is used to store the row number of the new element.
+ me.col = ((*matrixRow)[i].col * nComponents) + stencilCol;
+
+ rows[(rowIndex * nComponents) + stencilRow].push_back(me);
+
+ nElements++;
+ }
+ }
+ }
+
+ }
+ }
+
+ double *a = (double*)malloc(sizeof(double) * nElements);
+ int *ja = (int*)malloc(sizeof(int) * nElements);
+ int *ia = (int*)malloc(sizeof(int) * (newMatrixSize + 1));
+ double *bvec = (double*)malloc(sizeof(double) * newMatrixSize);
+ double *xvec = (double*)malloc(sizeof(double) * newMatrixSize);
+
+ int elCounter = 0;
+ int rowCounter = 0;
+ ia[0] = 1;
+
+ for (::std::vector< ::std::vector< MatEntry > >::iterator rowsIt = rows.begin();
+ rowsIt != rows.end();
+ ++rowsIt) {
+
+ sort((*rowsIt).begin(), (*rowsIt).end(), CmpMatEntryCol());
+
+ ia[rowCounter + 1] = ia[rowCounter] + (*rowsIt).size();
+
+ for (::std::vector<MatEntry>::iterator rowIt = (*rowsIt).begin();
+ rowIt != (*rowsIt).end();
+ rowIt++) {
+ a[elCounter] = (*rowIt).entry;
+ ja[elCounter] = (*rowIt).col + 1;
+
+ elCounter++;
+ }
+
+ rowCounter++;
+ }
+
+ // Resort the right hand side of the linear system.
+ for (int i = 0; i < b->getSize(); i++) {
+ DOFVector<double>::Iterator it(b->getDOFVector(i), USED_DOFS);
+
+ int counter = 0;
+ for (it.reset(); !it.end(); ++it, counter++) {
+ bvec[counter * nComponents + i] = *it;
+ }
+ }
+
+ // real unsymmetric matrix
+ int mtype = 11;
+
+ // number of right hand sides
+ int nRhs = 1;
+
+ // Pardiso internal memory
+ void *pt[64];
+ for (int i = 0; i < 64; i++) {
+ pt[i] = 0;
+ }
+
+ // Pardiso control parameters
+ int iparm[64];
+ for (int i = 0; i < 64; i++) {
+ iparm[i] = 0;
+ }
+
+ iparm[0] = 1; // No solver default
+ iparm[1] = 2; // Fill-in reordering from METIS
+ iparm[2] = 1; // Number of threads
+ iparm[7] = 2; // Max numbers of iterative refinement steps
+ iparm[9] = 13; // Perturb the pivot elements with 1e-13
+ iparm[10] = 1; // Use nonsymmetric permutation and scaling MPS
+ iparm[17] = -1; // Output: Number of nonzeros in the factor LU
+ iparm[18] = -1; // Output: Mflops for LU factorization
+
+ // Maximum number of numerical factorizations
+ int maxfct = 1;
+
+ // Which factorization to use
+ int mnum = 1;
+
+ // Print statistical information in file
+ int msglvl = 1;
+
+ // Error flag
+ int error = 0;
+
+ int n = newMatrixSize;
+
+ // Reordering and symbolic factorization
+ int phase = 11;
+
+ double ddum;
+
+ int idum;
+
+ PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nRhs,
+ iparm, &msglvl, &ddum, &ddum, &error);
+
+
+
+ free(a);
+ free(ja);
+ free(ia);
+ free(b);
+ free(x);
+
return(1);
}
}
--
GitLab