/[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 3282 by jfenwick, Mon Oct 11 01:48:14 2010 UTC revision 3283 by gross, Mon Oct 18 22:39:28 2010 UTC
# Line 34  Line 34 
34    
35  /*  free any extra stuff possibly used by the MKL library */  /*  free any extra stuff possibly used by the MKL library */
36    
37  void Paso_MKL_free(Paso_SystemMatrix* A) {  void Paso_MKL_free(Paso_SparseMatrix* A) {
38  #ifdef MKL  #ifdef MKL
39        index_t i;        index_t i;
40        if (A->solver!=NULL) {        if ((A->solver_p!=NULL) && (A->solver_package=PASO_MKL) ) {
41             _INTEGER_t mtype = MKL_MTYPE_UNSYM;            _INTEGER_t mtype = MKL_MTYPE_UNSYM;
42             if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;            _INTEGER_t n = A->numRows;
           _INTEGER_t n = A->mainBlock->numRows;  
43            _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */            _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
44            _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */            _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
45            _INTEGER_t msglvl=0; /* message level */            _INTEGER_t msglvl=0; /* message level */
# Line 52  void Paso_MKL_free(Paso_SystemMatrix* A) Line 51  void Paso_MKL_free(Paso_SystemMatrix* A)
51            for (i=0;i<64;++i) iparm[i]=0;            for (i=0;i<64;++i) iparm[i]=0;
52    
53            _INTEGER_t phase = MKL_PHASE_RELEASE_MEMORY;            _INTEGER_t phase = MKL_PHASE_RELEASE_MEMORY;
54            PARDISO ((_MKL_DSS_HANDLE_t *)(A->solver), &maxfct, &mnum, &mtype, &phase,            PARDISO ((_MKL_DSS_HANDLE_t *)(A->solver_p), &maxfct, &mnum, &mtype, &phase,
55                     &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,
56                     iparm, &msglvl,&ddum, &ddum, &error);                     iparm, &msglvl,&ddum, &ddum, &error);
57            MEMFREE(A->solver);            MEMFREE(A->solver_p);
58            if (error != MKL_ERROR_NO) Esys_setError(TYPE_ERROR,"memory release in paradiso library failed.");            if (error != MKL_ERROR_NO) Paso_setError(TYPE_ERROR,"memory release in PARDISO library failed.");
59       }       }
60  #endif  #endif
61  }  }
62    
63  /*  call the solver: */  /*  call the solver: */
64    
65  void Paso_MKL(Paso_SystemMatrix* A,  void Paso_MKL(Paso_SparseMatrix* A,
66                            double* out,                double* out,
67                            double* in,                double* in,
68                            Paso_Options* options,            index_t reordering,
69                            Paso_Performance* pp) {            dim_t numRefinements,
70              bool_t verbose)
71    {        
72    
73  #ifdef MKL  #ifdef MKL
74       double time0, time1;       double time0;
75       index_t i;       index_t i;
76    
77       if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {       if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
78          Esys_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.");
79          return;          return;
80       }       }
81       options->converged=FALSE;      
      Performance_startMonitor(pp,PERFORMANCE_ALL);  
82       _INTEGER_t mtype = MKL_MTYPE_UNSYM;       _INTEGER_t mtype = MKL_MTYPE_UNSYM;
83       if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;       _INTEGER_t n = A->numRows;
      _INTEGER_t n = A->mainBlock->numRows;  
84       _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */       _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */
85       _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */       _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */
86       _INTEGER_t msglvl=0; /* message level */       _INTEGER_t msglvl=0; /* message level */
# Line 90  void Paso_MKL(Paso_SystemMatrix* A, Line 91  void Paso_MKL(Paso_SystemMatrix* A,
91    
92       _INTEGER_t error=MKL_ERROR_NO;  /* error code */       _INTEGER_t error=MKL_ERROR_NO;  /* error code */
93       _INTEGER_t iparm[64]; /* parameters */       _INTEGER_t iparm[64]; /* parameters */
94       _MKL_DSS_HANDLE_t* pt = (_MKL_DSS_HANDLE_t *)(A->solver);       _MKL_DSS_HANDLE_t* pt = (_MKL_DSS_HANDLE_t *)(A->solver_p);
95       /* set iparm */       /* set iparm */
96       for (i=0;i<64;++i) iparm[i]=0;       for (i=0;i<64;++i) iparm[i]=0;
97       iparm[0] = 1; /* no default settings*/       iparm[0] = 1; /* no default settings*/
98       switch (options->reordering) {       switch (reordering) {
99              case PASO_MINIMUM_FILL_IN:              case PASO_MINIMUM_FILL_IN:
100                 iparm[1]=MKL_REORDERING_MINIMUM_DEGREE;                 iparm[1]=MKL_REORDERING_MINIMUM_DEGREE;
101                 break;                 break;
# Line 108  void Paso_MKL(Paso_SystemMatrix* A, Line 109  void Paso_MKL(Paso_SystemMatrix* A,
109       iparm[2] = 1;       iparm[2] = 1;
110       #endif       #endif
111       iparm[5] = 0; /* store solution into output array */       iparm[5] = 0; /* store solution into output array */
112       iparm[7] = 2; /* maximum number of refinements */       iparm[7] = numRefinements; /* maximum number of refinements */
113       iparm[9] = 13; /* 10**(-iparm[9]) preturbation of pivot elements */       iparm[9] = 13; /* 10**(-iparm[9]) preturbation of pivot elements */
114       iparm[10] = 1; /* rescaling the matrix before factorization started */       iparm[10] = 1; /* rescaling the matrix before factorization started */
115       iparm[17] =0; /* =-1 report number of non-zeroes */       iparm[17] =0; /* =-1 report number of non-zeroes */
# Line 118  void Paso_MKL(Paso_SystemMatrix* A, Line 119  void Paso_MKL(Paso_SystemMatrix* A,
119          /* allocate address pointer */          /* allocate address pointer */
120          pt=MEMALLOC(64,_MKL_DSS_HANDLE_t);          pt=MEMALLOC(64,_MKL_DSS_HANDLE_t);
121          if (Esys_checkPtr(pt)) return;          if (Esys_checkPtr(pt)) return;
122          A->solver=(void*) pt;          A->solver_p=(void*) pt;
123        A->solver_package=PASO_MKL;
124          for (i=0;i<64;++i) pt[i]=NULL;          for (i=0;i<64;++i) pt[i]=NULL;
125          time0=Esys_timer();          time0=Esys_timer();
126          /* symbolic factorization */          /* symbolic factorization */
127          phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;          phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;
128          PARDISO (pt, &maxfct, &mnum, &mtype, &phase,          PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
129                   &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,
130                   iparm, &msglvl, in, out, &error);                   iparm, &msglvl, in, out, &error);
131          if (error != MKL_ERROR_NO) {          if (error != MKL_ERROR_NO) {
132               if (options->verbose) printf("MKL: symbolic factorization factorization failed.\n");               if (verbose) printf("MKL: symbolic factorization factorization failed.\n");
133               Esys_setError(VALUE_ERROR,"symbolic factorization in paradiso library failed.");               Paso_setError(VALUE_ERROR,"symbolic factorization in PARDISO library failed.");
134               Paso_MKL_free(A);               Paso_MKL_free(A);
135          } else {          } else {
136             /* LDU factorization */             /* LDU factorization */
137             phase = MKL_PHASE_FACTORIZATION;             phase = MKL_PHASE_FACTORIZATION;
138             PARDISO(pt, &maxfct, &mnum, &mtype, &phase,             PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
139                  &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,
140                  iparm, &msglvl, in, out, &error);                  iparm, &msglvl, in, out, &error);
141             if (error != MKL_ERROR_NO) {             if (error != MKL_ERROR_NO) {
142               if (options->verbose) printf("MKL: LDU factorization failed.\n");               if (verbose) printf("MKL: LDU factorization failed.\n");
143               Esys_setError(ZERO_DIVISION_ERROR,"factorization in paradiso library failed. Most likely the matrix is singular.");               Paso_setError(ZERO_DIVISION_ERROR,"factorization in PARDISO library failed. Most likely the matrix is singular.");
144               Paso_MKL_free(A);               Paso_MKL_free(A);
145             }             }
146             if (options->verbose) printf("MKL: LDU factorization completed.\n");             if (verbose) printf("MKL: LDU factorization completed (time = %e).\n",Paso_timer()-time0);
         }  
         options->set_up_time=Esys_timer()-time0;  
      } else {  
         options->set_up_time=0;  
      }  
      /* forward backward substitution\ */  
      if (Esys_noError())  {  
         time0=Esys_timer();  
         phase = MKL_PHASE_SOLVE;  
         PARDISO (pt, &maxfct, &mnum, &mtype, &phase,  
                  &n, A->mainBlock->val, A->mainBlock->pattern->ptr, A->mainBlock->pattern->index, &idum, &nrhs,  
                  iparm, &msglvl, in, out, &error);  
         if (options->verbose) printf("MKL: solve completed.\n");  
         if (error != MKL_ERROR_NO) {  
               if (options->verbose) printf("MKL: forward/backward substitution failed.\n");  
               Esys_setError(VALUE_ERROR,"forward/backward substitution in paradiso library failed. Most likely the matrix is singular.");  
         } else {  
             if (options->verbose) printf("MKL: forward/backward substitution completed.\n");  
             options->residual_norm=0.;  
             options->num_iter=0;  
             options->num_level=0;  
             options->num_inner_iter=0;  
             options->converged=TRUE;  
147          }          }
         options->time=Esys_timer()-time0 + options->set_up_time;  
      }  
      Performance_stopMonitor(pp,PERFORMANCE_ALL);  
 #else  
     Esys_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");  
 #endif  
 }  
   
 void Paso_MKL_free1(Paso_SparseMatrix* A) {  
 #ifdef MKL  
       index_t i;  
       if (A->solver!=NULL) {  
            _INTEGER_t mtype = MKL_MTYPE_UNSYM;  
            if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;  
           _INTEGER_t n = A->numRows;  
           _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */  
           _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */  
           _INTEGER_t msglvl=0; /* message level */  
           _INTEGER_t nrhs=1; /* number of right hand sides */  
           _INTEGER_t idum; /* dummy integer */  
           _DOUBLE_PRECISION_t ddum; /* dummy float */  
           _INTEGER_t error=MKL_ERROR_NO;  /* error code */  
           _INTEGER_t iparm[64]; /* parameters */  
           for (i=0;i<64;++i) iparm[i]=0;  
   
           _INTEGER_t phase = MKL_PHASE_RELEASE_MEMORY;  
           PARDISO ((_MKL_DSS_HANDLE_t *)(A->solver), &maxfct, &mnum, &mtype, &phase,  
                    &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,  
                    iparm, &msglvl,&ddum, &ddum, &error);  
           MEMFREE(A->solver);  
           if (error != MKL_ERROR_NO) Esys_setError(TYPE_ERROR,"memory release in paradiso library failed.");  
      }  
 #endif  
 }  
 /*  call the solver: */  
   
 void Paso_MKL1(Paso_SparseMatrix* A,  
                           double* out,  
                           double* in,  
                           bool_t verbose) {  
 #ifdef MKL  
      index_t i;  
   
      if (! (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {  
         Esys_setError(TYPE_ERROR,"Paso_MKL: MKL requires CSR format with index offset 1 and block size 1.");  
         return;  
      }  
      _INTEGER_t mtype = MKL_MTYPE_UNSYM;  
      if (A->type & MATRIX_FORMAT_SYM) mtype=MKL_MTYPE_SYM;  
      _INTEGER_t n = A->numRows;  
      _INTEGER_t maxfct=1; /* number of factorizations on the same pattern */  
      _INTEGER_t mnum =1; /* factoriztion to be handeled in this call */  
      _INTEGER_t msglvl=0; /* message level */  
      _INTEGER_t nrhs=1; /* number of right hand sides */  
      _INTEGER_t idum; /* dummy integer */  
      _DOUBLE_PRECISION_t ddum; /* dummy float */  
      _INTEGER_t phase = MKL_PHASE_SYMBOLIC_FACTORIZATION;  
   
      _INTEGER_t error=MKL_ERROR_NO;  /* error code */  
      _INTEGER_t iparm[64]; /* parameters */  
      _MKL_DSS_HANDLE_t* pt = (_MKL_DSS_HANDLE_t *)(A->solver);  
      /* set iparm */  
      for (i=0;i<64;++i) iparm[i]=0;  
      iparm[0] = 1; /* no default settings*/  
      iparm[1]=MKL_REORDERING_MINIMUM_DEGREE;  
      #ifdef _OPENMP  
      iparm[2] = omp_get_max_threads();  
      #else  
      iparm[2] = 1;  
      #endif  
       
      iparm[5] = 0; /* store solution into output array */  
      iparm[7] = 2; /* maximum number of refinements */  
      iparm[9] = 13; /* 10**(-iparm[9]) preturbation of pivot elements */  
      iparm[10] = 1; /* rescaling the matrix before factorization started */  
      iparm[17] =0; /* =-1 report number of non-zeroes */  
      iparm[18] =0; /* =-1 report flops */  
148    
   
      if (pt==NULL) {  
         /* allocate address pointer */  
         pt=MEMALLOC(64,_MKL_DSS_HANDLE_t);  
         if (Esys_checkPtr(pt)) return;  
         A->solver=(void*) pt;  
         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) {  
              Esys_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) {  
              Esys_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");  
         }  
      }  
149       /* forward backward substitution\ */       /* forward backward substitution\ */
150       if (Esys_noError())  {       if (Esys_noError())  {
151            time0=Esys_timer();
152          phase = MKL_PHASE_SOLVE;          phase = MKL_PHASE_SOLVE;
153          PARDISO (pt, &maxfct, &mnum, &mtype, &phase,          PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
154                   &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,                   &n, A->val, A->pattern->ptr, A->pattern->index, &idum, &nrhs,
155                   iparm, &msglvl, in, out, &error);                   iparm, &msglvl, in, out, &error);
156          if (verbose) printf("MKL: solve completed.\n");          if (verbose) printf("MKL: solve completed.\n");
157          if (error != MKL_ERROR_NO) {          if (error != MKL_ERROR_NO) {
158                Esys_setError(VALUE_ERROR,"forward/backward substition in paradiso library failed. Most likely the matrix is singular.");                if (verbose) printf("MKL: forward/backward substitution failed.\n");
159                  Paso_setError(VALUE_ERROR,"forward/backward substitution in PARDISO library failed. Most likely the matrix is singular.");
160            } else {
161           if (verbose) printf("MKL: forward/backward substitution completed (time = %e).\n",Paso_timer()-time0);
162          }          }
163       }       }
164  #else  #else
165      Esys_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");      Esys_setError(SYSTEM_ERROR,"Paso_MKL:MKL is not avialble.");
166  #endif  #endif
167  }  }
 /*  
  * $Log$  
  *  
  */  
Legend:
Removed from v.3282  
changed lines
  Added in v.3283
  ViewVC Help
Powered by ViewVC 1.1.26