/[escript]/trunk/paso/src/MKL.c
ViewVC logotype

Diff of /trunk/paso/src/MKL.c

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

revision 3258 by artak, Wed Apr 28 02:21:23 2010 UTC revision 3259 by jfenwick, Mon Oct 11 01:48:14 2010 UTC
# Line 56  void Paso_MKL_free(Paso_SystemMatrix* A) Line 56  void Paso_MKL_free(Paso_SystemMatrix* A)
56                     &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,                     &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
57                     iparm, &msglvl,&ddum, &ddum, &error);                     iparm, &msglvl,&ddum, &ddum, &error);
58            MEMFREE(A->solver);            MEMFREE(A->solver);
59            if (error != MKL_ERROR_NO) Paso_setError(TYPE_ERROR,"memory release in paradiso library failed.");            if (error != MKL_ERROR_NO) Esys_setError(TYPE_ERROR,"memory release in paradiso library failed.");
60       }       }
61  #endif  #endif
62  }  }
# Line 72  void Paso_MKL(Paso_SystemMatrix* A, Line 72  void Paso_MKL(Paso_SystemMatrix* A,
72       index_t i;       index_t i;
73    
74       if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {       if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
75          Paso_setError(TYPE_ERROR,"Paso_MKL: MKL requires CSR format with index offset 1 and block size 1.");          Esys_setError(TYPE_ERROR,"Paso_MKL: MKL requires CSR format with index offset 1 and block size 1.");
76          return;          return;
77       }       }
78       options->converged=FALSE;       options->converged=FALSE;
# Line 117  void Paso_MKL(Paso_SystemMatrix* A, Line 117  void Paso_MKL(Paso_SystemMatrix* A,
117       if (pt==NULL) {       if (pt==NULL) {
118          /* allocate address pointer */          /* allocate address pointer */
119          pt=MEMALLOC(64,_MKL_DSS_HANDLE_t);          pt=MEMALLOC(64,_MKL_DSS_HANDLE_t);
120          if (Paso_checkPtr(pt)) return;          if (Esys_checkPtr(pt)) return;
121          A->solver=(void*) pt;          A->solver=(void*) pt;
122          for (i=0;i<64;++i) pt[i]=NULL;          for (i=0;i<64;++i) pt[i]=NULL;
123          time0=Paso_timer();          time0=Esys_timer();
124          /* symbolic factorization */          /* symbolic factorization */
125          phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;          phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
126          PARDISO (pt, &maxfct, &mnum, &mtype, &phase,          PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
# Line 128  void Paso_MKL(Paso_SystemMatrix* A, Line 128  void Paso_MKL(Paso_SystemMatrix* A,
128                   iparm, &msglvl, in, out, &error);                   iparm, &msglvl, in, out, &error);
129          if (error != MKL_ERROR_NO) {          if (error != MKL_ERROR_NO) {
130               if (options->verbose) printf("MKL: symbolic factorization factorization failed.\n");               if (options->verbose) printf("MKL: symbolic factorization factorization failed.\n");
131               Paso_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");               Esys_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");
132               Paso_MKL_free(A);               Paso_MKL_free(A);
133          } else {          } else {
134             /* LDU factorization */             /* LDU factorization */
# Line 138  void Paso_MKL(Paso_SystemMatrix* A, Line 138  void Paso_MKL(Paso_SystemMatrix* A,
138                  iparm, &msglvl, in, out, &error);                  iparm, &msglvl, in, out, &error);
139             if (error != MKL_ERROR_NO) {             if (error != MKL_ERROR_NO) {
140               if (options->verbose) printf("MKL: LDU factorization failed.\n");               if (options->verbose) printf("MKL: LDU factorization failed.\n");
141               Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular.");               Esys_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular.");
142               Paso_MKL_free(A);               Paso_MKL_free(A);
143             }             }
144             if (options->verbose) printf("MKL: LDU factorization completed.\n");             if (options->verbose) printf("MKL: LDU factorization completed.\n");
145          }          }
146          options->set_up_time=Paso_timer()-time0;          options->set_up_time=Esys_timer()-time0;
147       } else {       } else {
148          options->set_up_time=0;          options->set_up_time=0;
149       }       }
150       /* forward backward substitution\ */       /* forward backward substitution\ */
151       if (Paso_noError())  {       if (Esys_noError())  {
152          time0=Paso_timer();          time0=Esys_timer();
153          phase = MKL_PHASE_SOLVE;          phase = MKL_PHASE_SOLVE;
154          PARDISO (pt, &maxfct, &mnum, &mtype, &phase,          PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
155                   &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,                   &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
# Line 157  void Paso_MKL(Paso_SystemMatrix* A, Line 157  void Paso_MKL(Paso_SystemMatrix* A,
157          if (options->verbose) printf("MKL: solve completed.\n");          if (options->verbose) printf("MKL: solve completed.\n");
158          if (error != MKL_ERROR_NO) {          if (error != MKL_ERROR_NO) {
159                if (options->verbose) printf("MKL: forward/backward substitution failed.\n");                if (options->verbose) printf("MKL: forward/backward substitution failed.\n");
160                Paso_setError(VALUE_ERROR,"forward/backward substitution in paradiso library failed. Most likely the matrix is singular.");                Esys_setError(VALUE_ERROR,"forward/backward substitution in paradiso library failed. Most likely the matrix is singular.");
161          } else {          } else {
162              if (options->verbose) printf("MKL: forward/backward substitution completed.\n");              if (options->verbose) printf("MKL: forward/backward substitution completed.\n");
163              options->residual_norm=0.;              options->residual_norm=0.;
# Line 166  void Paso_MKL(Paso_SystemMatrix* A, Line 166  void Paso_MKL(Paso_SystemMatrix* A,
166              options->num_inner_iter=0;              options->num_inner_iter=0;
167              options->converged=TRUE;              options->converged=TRUE;
168          }          }
169          options->time=Paso_timer()-time0 + options->set_up_time;          options->time=Esys_timer()-time0 + options->set_up_time;
170       }       }
171       Performance_stopMonitor(pp,PERFORMANCE_ALL);       Performance_stopMonitor(pp,PERFORMANCE_ALL);
172  #else  #else
173      Paso_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");      Esys_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");
174  #endif  #endif
175  }  }
176    
# Line 196  void Paso_MKL_free1(Paso_SparseMatrix* A Line 196  void Paso_MKL_free1(Paso_SparseMatrix* A
196                     &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,                     &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
197                     iparm, &msglvl,&ddum, &ddum, &error);                     iparm, &msglvl,&ddum, &ddum, &error);
198            MEMFREE(A->solver);            MEMFREE(A->solver);
199            if (error != MKL_ERROR_NO) Paso_setError(TYPE_ERROR,"memory release in paradiso library failed.");            if (error != MKL_ERROR_NO) Esys_setError(TYPE_ERROR,"memory release in paradiso library failed.");
200       }       }
201  #endif  #endif
202  }  }
# Line 210  void Paso_MKL1(Paso_SparseMatrix* A, Line 210  void Paso_MKL1(Paso_SparseMatrix* A,
210       index_t i;       index_t i;
211    
212       if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {       if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
213          Paso_setError(TYPE_ERROR,"Paso_MKL: MKL requires CSR format with index offset 1 and block size 1.");          Esys_setError(TYPE_ERROR,"Paso_MKL: MKL requires CSR format with index offset 1 and block size 1.");
214          return;          return;
215       }       }
216       _INTEGER_t mtype = MKL_MTYPE_UNSYM;       _INTEGER_t mtype = MKL_MTYPE_UNSYM;
# Line 248  void Paso_MKL1(Paso_SparseMatrix* A, Line 248  void Paso_MKL1(Paso_SparseMatrix* A,
248       if (pt==NULL) {       if (pt==NULL) {
249          /* allocate address pointer */          /* allocate address pointer */
250          pt=MEMALLOC(64,_MKL_DSS_HANDLE_t);          pt=MEMALLOC(64,_MKL_DSS_HANDLE_t);
251          if (Paso_checkPtr(pt)) return;          if (Esys_checkPtr(pt)) return;
252          A->solver=(void*) pt;          A->solver=(void*) pt;
253          for (i=0;i<64;++i) pt[i]=NULL;          for (i=0;i<64;++i) pt[i]=NULL;
254          /* symbolic factorization */          /* symbolic factorization */
# Line 257  void Paso_MKL1(Paso_SparseMatrix* A, Line 257  void Paso_MKL1(Paso_SparseMatrix* A,
257                   &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,                   &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
258                   iparm, &msglvl, in, out, &error);                   iparm, &msglvl, in, out, &error);
259          if (error != MKL_ERROR_NO) {          if (error != MKL_ERROR_NO) {
260               Paso_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");               Esys_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");
261               Paso_MKL_free1(A);               Paso_MKL_free1(A);
262          } else {          } else {
263             /* LDU factorization */             /* LDU factorization */
# Line 266  void Paso_MKL1(Paso_SparseMatrix* A, Line 266  void Paso_MKL1(Paso_SparseMatrix* A,
266                  &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,                  &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
267                  iparm, &msglvl, in, out, &error);                  iparm, &msglvl, in, out, &error);
268             if (error != MKL_ERROR_NO) {             if (error != MKL_ERROR_NO) {
269               Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular.");               Esys_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular.");
270               Paso_MKL_free1(A);               Paso_MKL_free1(A);
271             }             }
272             if (verbose) printf("MKL: LDU factorization completed.\n");             if (verbose) printf("MKL: LDU factorization completed.\n");
273          }          }
274       }       }
275       /* forward backward substitution\ */       /* forward backward substitution\ */
276       if (Paso_noError())  {       if (Esys_noError())  {
277          phase = MKL_PHASE_SOLVE;          phase = MKL_PHASE_SOLVE;
278          PARDISO (pt, &maxfct, &mnum, &mtype, &phase,          PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
279                   &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,                   &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
280                   iparm, &msglvl, in, out, &error);                   iparm, &msglvl, in, out, &error);
281          if (verbose) printf("MKL: solve completed.\n");          if (verbose) printf("MKL: solve completed.\n");
282          if (error != MKL_ERROR_NO) {          if (error != MKL_ERROR_NO) {
283                Paso_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed. Most likely the matrix is singular.");                Esys_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed. Most likely the matrix is singular.");
284          }          }
285       }       }
286  #else  #else
287      Paso_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");      Esys_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");
288  #endif  #endif
289  }  }
290  /*  /*

Legend:
Removed from v.3258  
changed lines
  Added in v.3259

  ViewVC Help
Powered by ViewVC 1.1.26