/[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

trunk/paso/src/MKL.c revision 425 by gross, Tue Jan 10 04:10:39 2006 UTC temp/paso/src/MKL.c revision 1387 by trankine, Fri Jan 11 07:45:26 2008 UTC
# Line 1  Line 1 
1    
2  /* $Id$ */  /* $Id$ */
3    
4    /*******************************************************
5     *
6     *           Copyright 2003-2007 by ACceSS MNRF
7     *       Copyright 2007 by University of Queensland
8     *
9     *                http://esscc.uq.edu.au
10     *        Primary Business: Queensland, Australia
11     *  Licensed under the Open Software License version 3.0
12     *     http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16  /**************************************************************/  /**************************************************************/
17    
18  /* Paso: interface to the Intel MKL library */  /* Paso: interface to the Intel MKL library */
# Line 27  void Paso_MKL_free(Paso_SystemMatrix* A) Line 40  void Paso_MKL_free(Paso_SystemMatrix* A)
40        if (A->solver!=NULL) {        if (A->solver!=NULL) {
41             _INTEGER_t mtype = MKL_MTYPE_UNSYM;             _INTEGER_t mtype = MKL_MTYPE_UNSYM;
42             if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;             if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;
43            _INTEGER_t n = A->num_rows;            _INTEGER_t n = A->mainBlock->numRows;
44            _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */            _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
45            _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */            _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
46            _INTEGER_t msglvl=0; /* message level */            _INTEGER_t msglvl=0; /* message level */
# Line 40  void Paso_MKL_free(Paso_SystemMatrix* A) Line 53  void Paso_MKL_free(Paso_SystemMatrix* A)
53    
54            _INTEGER_t phase = MKL_PHASE_RELEASE_MEMORY;            _INTEGER_t phase = MKL_PHASE_RELEASE_MEMORY;
55            PARDISO ((_MKL_DSS_HANDLE_t *)(A->solver), &maxfct, &mnum, &mtype, &phase,            PARDISO ((_MKL_DSS_HANDLE_t *)(A->solver), &maxfct, &mnum, &mtype, &phase,
56                     &n, A->val, A->pattern->ptr, A->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) Paso_setError(TYPE_ERROR,"memory release in paradiso library failed.");
# Line 52  void Paso_MKL_free(Paso_SystemMatrix* A) Line 65  void Paso_MKL_free(Paso_SystemMatrix* A)
65  void Paso_MKL(Paso_SystemMatrix* A,  void Paso_MKL(Paso_SystemMatrix* A,
66                            double* out,                            double* out,
67                            double* in,                            double* in,
68                            Paso_Options* options) {                            Paso_Options* options,
69                              Paso_Performance* pp) {
70  #ifdef MKL  #ifdef MKL
71       double time0;       double time0;
72       index_t i;       index_t i;
# Line 61  void Paso_MKL(Paso_SystemMatrix* A, Line 75  void Paso_MKL(Paso_SystemMatrix* A,
75          Paso_setError(TYPE_ERROR,"Paso_MKL: MKL requires CSR format with index offset 1 and block size 1.");          Paso_setError(TYPE_ERROR,"Paso_MKL: MKL requires CSR format with index offset 1 and block size 1.");
76          return;          return;
77       }       }
78         Performance_startMonitor(pp,PERFORMANCE_ALL);
79       _INTEGER_t mtype = MKL_MTYPE_UNSYM;       _INTEGER_t mtype = MKL_MTYPE_UNSYM;
80       if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;       if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;
81       _INTEGER_t n = A->num_rows;       _INTEGER_t n = A->mainBlock->numRows;
82       _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */       _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
83       _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */       _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
84       _INTEGER_t msglvl=0; /* message level */       _INTEGER_t msglvl=0; /* message level */
# Line 109  void Paso_MKL(Paso_SystemMatrix* A, Line 124  void Paso_MKL(Paso_SystemMatrix* A,
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,
127                   &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,                   &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
128                   iparm, &msglvl, in, out, &error);                   iparm, &msglvl, in, out, &error);
129          if (error != MKL_ERROR_NO) {          if (error != MKL_ERROR_NO) {
130               Paso_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");               Paso_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");
131               Paso_MKL_free(A);               Paso_MKL_free(A);
132               return;          } else {
133          }             /* LDU factorization */
134          /* LDU factorization */             phase = MKL_PHASE_FACTORIZATION;
135          phase = MKL_PHASE_FACTORIZATION;             PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
136          PARDISO(pt, &maxfct, &mnum, &mtype, &phase,                  &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
                 &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,  
137                  iparm, &msglvl, in, out, &error);                  iparm, &msglvl, in, out, &error);
138          if (error != MKL_ERROR_NO) {             if (error != MKL_ERROR_NO) {
139               Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed.");               Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular.");
140               Paso_MKL_free(A);               Paso_MKL_free(A);
141               return;             }
142               if (options->verbose) printf("timing MKL: LDU factorization: %.4e sec.\n",Paso_timer()-time0);
143          }          }
         if (options->verbose) printf("timing MKL: LDU factorization: %.4e sec.\n",Paso_timer()-time0);  
144       }       }
145       /* forward backward substitution\ */       /* forward backward substitution\ */
146       if (Paso_noError())  {       if (Paso_noError())  {
147          time0=Paso_timer();          time0=Paso_timer();
148          phase = MKL_PHASE_SOLVE;          phase = MKL_PHASE_SOLVE;
149          PARDISO (pt, &maxfct, &mnum, &mtype, &phase,          PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
150                   &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,                   &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,
151                   iparm, &msglvl, in, out, &error);                   iparm, &msglvl, in, out, &error);
152          if (options->verbose) printf("timing MKL: solve: %.4e sec\n",Paso_timer()-time0);          if (options->verbose) printf("timing MKL: solve: %.4e sec\n",Paso_timer()-time0);
153          if (error != MKL_ERROR_NO) {          if (error != MKL_ERROR_NO) {
154                Paso_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed.");                Paso_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed. Most likely the matrix is singular.");
               return;  
155          }          }
156       }       }
157         Performance_stopMonitor(pp,PERFORMANCE_ALL);
158  #else  #else
159      Paso_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");      Paso_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");
160  #endif  #endif

Legend:
Removed from v.425  
changed lines
  Added in v.1387

  ViewVC Help
Powered by ViewVC 1.1.26