/[escript]/branches/trilinos_from_5897/paso/src/MKL.cpp
ViewVC logotype

Diff of /branches/trilinos_from_5897/paso/src/MKL.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 6008 by caltinay, Fri Feb 5 03:37:49 2016 UTC revision 6009 by caltinay, Wed Mar 2 04:13:26 2016 UTC
# Line 28  Line 28 
28    
29  #include "MKL.h"  #include "MKL.h"
30  #include "Options.h"  #include "Options.h"
31    #include "PasoException.h"
32    
33  namespace paso {  namespace paso {
34    
# Line 57  void MKL_free(SparseMatrix* A) Line 58  void MKL_free(SparseMatrix* A)
58          delete[] A->solver_p;          delete[] A->solver_p;
59          A->solver_p=NULL;          A->solver_p=NULL;
60          if (error != MKL_ERROR_NO)          if (error != MKL_ERROR_NO)
61              Esys_setError(SYSTEM_ERROR,              throw PasoException("Memory release in MKL library failed.");
                           "memory release in MKL library failed.");  
62      }      }
63  #endif  #endif
64  }  }
# Line 68  void MKL_solve(SparseMatrix_ptr A, doubl Line 68  void MKL_solve(SparseMatrix_ptr A, doubl
68  {  {
69  #ifdef MKL  #ifdef MKL
70      if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {      if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
71          Esys_setError(TYPE_ERROR, "Paso: MKL requires CSR format with index offset 1 and block size 1.");          throw PasoException("Paso: MKL requires CSR format with index offset 1 and block size 1.");
         return;  
72      }      }
73    
74      // MKL uses 'long long int' in 64-bit version, escript 'long'. Make sure      // MKL uses 'long long int' in 64-bit version, escript 'long'. Make sure
75      // they are compatible      // they are compatible
76      if (sizeof(ES_MKL_INT) != sizeof(index_t)) {      if (sizeof(ES_MKL_INT) != sizeof(index_t)) {
77          Esys_setError(TYPE_ERROR, "Paso: MKL index type is not compatible with this escript build. Check compile options.");          throw PasoException("Paso: MKL index type is not compatible with this escript build. Check compile options.");
         return;  
78      }      }
79    
80      ES_MKL_INT* ptr = reinterpret_cast<ES_MKL_INT*>(A->pattern->ptr);      ES_MKL_INT* ptr = reinterpret_cast<ES_MKL_INT*>(A->pattern->ptr);
# Line 132  void MKL_solve(SparseMatrix_ptr A, doubl Line 130  void MKL_solve(SparseMatrix_ptr A, doubl
130          A->solver_package = PASO_MKL;          A->solver_package = PASO_MKL;
131          // symbolic factorization          // symbolic factorization
132          phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;          phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
133          time0 = Esys_timer();          time0 = escript::gettime();
134          ES_PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, A->val,          ES_PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, A->val,
135                     ptr, index, &idum, &nrhs, iparm, &msglvl, in, out, &error);                     ptr, index, &idum, &nrhs, iparm, &msglvl, in, out, &error);
136          if (error != MKL_ERROR_NO) {          if (error != MKL_ERROR_NO) {
137               if (verbose)               if (verbose)
138                   printf("MKL: symbolic factorization failed.\n");                   printf("MKL: symbolic factorization failed.\n");
              Esys_setError(VALUE_ERROR,"symbolic factorization in MKL library failed.");  
139               MKL_free(A.get());               MKL_free(A.get());
140                 throw PasoException("symbolic factorization in MKL library failed.");
141          } else {          } else {
142              // LDU factorization              // LDU factorization
143              phase = MKL_PHASE_FACTORIZATION;              phase = MKL_PHASE_FACTORIZATION;
# Line 148  void MKL_solve(SparseMatrix_ptr A, doubl Line 146  void MKL_solve(SparseMatrix_ptr A, doubl
146              if (error != MKL_ERROR_NO) {              if (error != MKL_ERROR_NO) {
147                  if (verbose)                  if (verbose)
148                      printf("MKL: LDU factorization failed.\n");                      printf("MKL: LDU factorization failed.\n");
                 Esys_setError(ZERO_DIVISION_ERROR, "factorization in MKL library failed. Most likely the matrix is singular.");  
149                  MKL_free(A.get());                  MKL_free(A.get());
150                    throw PasoException("factorization in MKL library failed. Most likely the matrix is singular.");
151             }             }
152             if (verbose)             if (verbose)
153                 printf("MKL: LDU factorization completed (time = %e).\n", Esys_timer()-time0);                 printf("MKL: LDU factorization completed (time = %e).\n", escript::gettime()-time0);
154          }          }
155      }      }
156      // forward backward substitution      // forward backward substitution
157      if (Esys_noError())  {      time0 = escript::gettime();
158          time0 = Esys_timer();      phase = MKL_PHASE_SOLVE;
159          phase = MKL_PHASE_SOLVE;      ES_PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, A->val,
160          ES_PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, A->val,                 ptr, index, &idum, &nrhs, iparm, &msglvl, in, out, &error);
161                     ptr, index, &idum, &nrhs, iparm, &msglvl, in, out, &error);      if (verbose) printf("MKL: solve completed.\n");
162          if (verbose) printf("MKL: solve completed.\n");      if (error != MKL_ERROR_NO) {
163          if (error != MKL_ERROR_NO) {          if (verbose)
164              if (verbose)              printf("MKL: forward/backward substitution failed.\n");
165                  printf("MKL: forward/backward substitution failed.\n");          throw PasoException("forward/backward substitution in MKL library failed. Most likely the matrix is singular.");
166              Esys_setError(ZERO_DIVISION_ERROR, "forward/backward substitution in MKL library failed. Most likely the matrix is singular.");      } else {
167          } else {          if (verbose)
168              if (verbose)              printf("MKL: forward/backward substitution completed (time = %e).\n", escript::gettime()-time0);
                 printf("MKL: forward/backward substitution completed (time = %e).\n", Esys_timer()-time0);  
         }  
169      }      }
170  #else  #else
171      Esys_setError(SYSTEM_ERROR, "Paso: MKL is not available.");      throw PasoException("Paso: MKL is not available.");
172  #endif  #endif
173  }  }
174    

Legend:
Removed from v.6008  
changed lines
  Added in v.6009

  ViewVC Help
Powered by ViewVC 1.1.26