* Pardiso Solver

AMDiS/src/PardisoSolver.cc

@@ -4,6 +4,7 @@
 
 #ifdef HAVE_MKL
 
+#include <mkl.h>
 #include <mkl_pardiso.h>
 
 namespace AMDiS {
@@ -120,7 +121,7 @@ namespace AMDiS {
 
     iparm[0] = 1; // No solver default
     iparm[1] = 2; // Fill-in reordering from METIS
-    iparm[2] = 1; // Number of threads
+    iparm[2] = mkl_get_max_threads(); // 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
@@ -151,7 +152,52 @@ namespace AMDiS {
     PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nRhs, 
 	    iparm, &msglvl, &ddum, &ddum, &error);
     
+    TEST_EXIT(error == 0)("Intel MKL Pardiso error during symbolic factorization: %d\n", error);
+
+    // Numerical factorization
+    phase = 22;
+
+    PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nRhs, 
+	    iparm, &msglvl, &ddum, &ddum, &error);
+    
+    TEST_EXIT(error == 0)("Intel MKL Pardiso error during numerical factorization: %d\n", error);
+
+    // Back substitution and iterative refinement
+    phase = 33;
+    iparm[7] = 2; // Max numbers of iterative refinement steps
+
+    PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nRhs,
+	    iparm, &msglvl, bvec, xvec, &error);
+
+    TEST_EXIT(error == 0)("Intel MKL Pardiso error during solution: %d\n", error);
+
+    // Termination and release of memory
+    phase = -1;
+
+    PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nRhs, 
+	    iparm, &msglvl, &ddum, &ddum, &error);    
+
+
+    // Copy and resort solution.
+    for (int i = 0; i < x->getSize(); i++) {
+      DOFVector<double>::Iterator it(x->getDOFVector(i), USED_DOFS);
+      
+      int counter = 0;
+      for (it.reset(); !it.end(); it++, counter++) {
+	*it = xvec[counter * nComponents + i];
+      }
+    }
+
+    // Calculate and print the residual.
+    *p = *x;
+    *p *= -1.0;
+    matVec->matVec(NoTranspose, *p, *r);
+    *r += *b;
     
+    this->residual = norm(r);
+
+    MSG("Residual: %e\n", this->residual);
+
 
     free(a);
     free(ja);