/[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 3009 by jfenwick, Thu Jan 28 02:03:15 2010 UTC revision 3010 by artak, Tue Apr 27 05:10:46 2010 UTC
# Line 209  void Paso_MKL1(Paso_SparseMatrix* A, Line 209  void Paso_MKL1(Paso_SparseMatrix* A,
209                            bool_t verbose) {                            bool_t verbose) {
210  #ifdef MKL  #ifdef MKL
211       index_t i;       index_t i;
212        
213         /***** FOR AMG **/
214         double **xx;
215         double **bb;
216         Paso_SparseMatrix *block[MAX_BLOCK_SIZE];
217        
218         xx=MEMALLOC(A->row_block_size,double*);
219         bb=MEMALLOC(A->row_block_size,double*);
220         if (Paso_checkPtr(xx) || Paso_checkPtr(bb)) return;
221         /****/
222    
223    
224       if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {       if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
225          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.");
# Line 246  void Paso_MKL1(Paso_SparseMatrix* A, Line 257  void Paso_MKL1(Paso_SparseMatrix* A,
257       iparm[18] =0; /* =-1 report flops */       iparm[18] =0; /* =-1 report flops */
258    
259    
260       if (pt==NULL) {       /******* FOR AMG ****/
261          /* allocate address pointer */       for (i=0;i<A->row_block_size;i++) {
262          pt=MEMALLOC(64,_MKL_DSS_HANDLE_t);            xx[i]=MEMALLOC(A->numRows,double);
263          if (Paso_checkPtr(pt)) return;            bb[i]=MEMALLOC(A->numRows,double);
264          A->solver=(void*) pt;            if (Paso_checkPtr(xx[i]) && Paso_checkPtr(bb[i])) return;
         for (i=0;i<64;++i) pt[i]=NULL;  
         /* symbolic factorization */  
         phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;  
         PARDISO (pt, &maxfct, &mnum, &mtype, &phase,  
                  &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,  
                  iparm, &msglvl, in, out, &error);  
         if (error != MKL_ERROR_NO) {  
              Paso_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");  
              Paso_MKL_free1(A);  
         } else {  
            /* LDU factorization */  
            phase = MKL_PHASE_FACTORIZATION;  
            PARDISO(pt, &maxfct, &mnum, &mtype, &phase,  
                 &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,  
                 iparm, &msglvl, in, out, &error);  
            if (error != MKL_ERROR_NO) {  
              Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular.");  
              Paso_MKL_free1(A);  
            }  
            if (verbose) printf("MKL: LDU factorization completed.\n");  
         }  
265       }       }
266       /* forward backward substitution\ */      
267       if (Paso_noError())  {       /*#pragma omp parallel for private(i,j) schedule(static)*/
268          phase = MKL_PHASE_SOLVE;       for (i=0;i<A->numRows;i++) {
269          PARDISO (pt, &maxfct, &mnum, &mtype, &phase,           for (j=0;j<A->row_block_size;j++) {
270                   &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,            bb[j][i]=in[A->row_block_size*i+j];
271                   iparm, &msglvl, in, out, &error);            xx[j][i]=0;  
272          if (verbose) printf("MKL: solve completed.\n");           }
273          if (error != MKL_ERROR_NO) {        }
274                Paso_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed. Most likely the matrix is singular.");      
275          }       for (i=0;i<MAX_BLOCK_SIZE;++i) {
276              block[i]=NULL;
277       }       }
278         /*****************/
279    
280         for (i=0;i<A->row_block_size;++i) {
281            block[i]=Paso_SparseMatrix_getBlock(A,i+1);
282            
283              if (pt==NULL) {
284                 /* allocate address pointer */
285                 pt=MEMALLOC(64,_MKL_DSS_HANDLE_t);
286                 if (Paso_checkPtr(pt)) return;
287                 block[i]->solver=(void*) pt;
288                 for (i=0;i<64;++i) pt[i]=NULL;
289                 /* symbolic factorization */
290                 phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
291                 PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
292                          &n, block[i]->val, block[i]->pattern->ptr, block[i]->pattern->index, &idum, &nrhs,
293                          iparm, &msglvl, bb[i], xx[i], &error);
294                 if (error != MKL_ERROR_NO) {
295                      Paso_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");
296                      Paso_MKL_free1(block[i]);
297                 } else {
298                    /* LDU factorization */
299                    phase = MKL_PHASE_FACTORIZATION;
300                    PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
301                         &n, block[i]->val, block[i]->pattern->ptr, block[i]->pattern->index, &idum, &nrhs,
302                         iparm, &msglvl, bb[i], xx[i], &error);
303                    if (error != MKL_ERROR_NO) {
304                      Paso_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular.");
305                      Paso_MKL_free1(block[i]);
306                    }
307                    if (verbose) printf("MKL: LDU factorization completed.\n");
308                 }
309              }
310              /* forward backward substitution\ */
311              if (Paso_noError())  {
312                 phase = MKL_PHASE_SOLVE;
313                 PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
314                          &n, block[i]->val, block[i]->pattern->ptr, block[i]->pattern->index, &idum, &nrhs,
315                          iparm, &msglvl, bb[i], xx[i], &error);
316                 if (verbose) printf("MKL: solve completed.\n");
317                 if (error != MKL_ERROR_NO) {
318                       Paso_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed. Most likely the matrix is singular.");
319                 }
320              }
321         }
322        
323         /***** FOR AMG ********/
324         /*#pragma omp parallel for private(i,j) schedule(static)*/
325          for (i=0;i<A->numRows;i++) {
326              for (j=0;j<A->row_block_size;j++) {
327              out[A->row_block_size*i+j]=xx[j][i];
328              }
329           }
330          
331         for (i=0;i<A->row_block_size;i++) {
332                    MEMFREE(xx[i]);
333                    MEMFREE(bb[i]);
334                    Paso_SparseMatrix_free(block[i]);
335          }
336          /*****************/
337  #else  #else
338      Paso_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");      Paso_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");
339  #endif  #endif

Legend:
Removed from v.3009  
changed lines
  Added in v.3010

  ViewVC Help
Powered by ViewVC 1.1.26