Actual source code: mkl_cpardiso.c

petsc-master 2019-09-20
Report Typos and Errors

  2: #include <petscsys.h>
  3: #include <../src/mat/impls/aij/mpi/mpiaij.h> /*I  "petscmat.h"  I*/
  4: #include <../src/mat/impls/sbaij/mpi/mpisbaij.h>

  6: #if defined(PETSC_HAVE_MKL_INTEL_ILP64)
  7: #define MKL_ILP64
  8: #endif
  9: #include <mkl.h>
 10: #include <mkl_cluster_sparse_solver.h>

 12: /*
 13:  *  Possible mkl_cpardiso phases that controls the execution of the solver.
 14:  *  For more information check mkl_cpardiso manual.
 15:  */
 16: #define JOB_ANALYSIS 11
 17: #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION 12
 18: #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 13
 19: #define JOB_NUMERICAL_FACTORIZATION 22
 20: #define JOB_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 23
 21: #define JOB_SOLVE_ITERATIVE_REFINEMENT 33
 22: #define JOB_SOLVE_FORWARD_SUBSTITUTION 331
 23: #define JOB_SOLVE_DIAGONAL_SUBSTITUTION 332
 24: #define JOB_SOLVE_BACKWARD_SUBSTITUTION 333
 25: #define JOB_RELEASE_OF_LU_MEMORY 0
 26: #define JOB_RELEASE_OF_ALL_MEMORY -1

 28: #define IPARM_SIZE 64
 29: #define INT_TYPE MKL_INT

 31: static const char *Err_MSG_CPardiso(int errNo) {
 32:   switch (errNo) {
 33:     case -1:
 34:       return "input inconsistent"; break;
 35:     case -2:
 36:       return "not enough memory"; break;
 37:     case -3:
 38:       return "reordering problem"; break;
 39:     case -4:
 40:       return "zero pivot, numerical factorization or iterative refinement problem"; break;
 41:     case -5:
 42:       return "unclassified (internal) error"; break;
 43:     case -6:
 44:       return "preordering failed (matrix types 11, 13 only)"; break;
 45:     case -7:
 46:       return "diagonal matrix problem"; break;
 47:     case -8:
 48:       return "32-bit integer overflow problem"; break;
 49:     case -9:
 50:       return "not enough memory for OOC"; break;
 51:     case -10:
 52:       return "problems with opening OOC temporary files"; break;
 53:     case -11:
 54:       return "read/write problems with the OOC data file"; break;
 55:     default :
 56:       return "unknown error";
 57:   }
 58: }

 60: /*
 61:  *  Internal data structure.
 62:  *  For more information check mkl_cpardiso manual.
 63:  */

 65: typedef struct {

 67:   /* Configuration vector */
 68:   INT_TYPE     iparm[IPARM_SIZE];

 70:   /*
 71:    * Internal mkl_cpardiso memory location.
 72:    * After the first call to mkl_cpardiso do not modify pt, as that could cause a serious memory leak.
 73:    */
 74:   void         *pt[IPARM_SIZE];

 76:   MPI_Comm     comm_mkl_cpardiso;

 78:   /* Basic mkl_cpardiso info*/
 79:   INT_TYPE     phase, maxfct, mnum, mtype, n, nrhs, msglvl, err;

 81:   /* Matrix structure */
 82:   PetscScalar  *a;

 84:   INT_TYPE     *ia, *ja;

 86:   /* Number of non-zero elements */
 87:   INT_TYPE     nz;

 89:   /* Row permutaton vector*/
 90:   INT_TYPE     *perm;

 92:   /* Define is matrix preserve sparce structure. */
 93:   MatStructure matstruc;

 95:   PetscErrorCode (*ConvertToTriples)(Mat, MatReuse, PetscInt*, PetscInt**, PetscInt**, PetscScalar**);

 97:   /* True if mkl_cpardiso function have been used. */
 98:   PetscBool CleanUp;
 99: } Mat_MKL_CPARDISO;

101: /*
102:  * Copy the elements of matrix A.
103:  * Input:
104:  *   - Mat A: MATSEQAIJ matrix
105:  *   - int shift: matrix index.
106:  *     - 0 for c representation
107:  *     - 1 for fortran representation
108:  *   - MatReuse reuse:
109:  *     - MAT_INITIAL_MATRIX: Create a new aij representation
110:  *     - MAT_REUSE_MATRIX: Reuse all aij representation and just change values
111:  * Output:
112:  *   - int *nnz: Number of nonzero-elements.
113:  *   - int **r pointer to i index
114:  *   - int **c pointer to j elements
115:  *   - MATRIXTYPE **v: Non-zero elements
116:  */
117: PetscErrorCode MatCopy_seqaij_seqaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v)
118: {
119:   Mat_SeqAIJ *aa=(Mat_SeqAIJ*)A->data;

122:   *v=aa->a;
123:   if (reuse == MAT_INITIAL_MATRIX) {
124:     *r   = (INT_TYPE*)aa->i;
125:     *c   = (INT_TYPE*)aa->j;
126:     *nnz = aa->nz;
127:   }
128:   return(0);
129: }

131: PetscErrorCode MatConvertToTriples_mpiaij_mpiaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v)
132: {
133:   const PetscInt    *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj;
134:   PetscErrorCode    ierr;
135:   PetscInt          rstart,nz,i,j,countA,countB;
136:   PetscInt          *row,*col;
137:   const PetscScalar *av, *bv;
138:   PetscScalar       *val;
139:   Mat_MPIAIJ        *mat = (Mat_MPIAIJ*)A->data;
140:   Mat_SeqAIJ        *aa  = (Mat_SeqAIJ*)(mat->A)->data;
141:   Mat_SeqAIJ        *bb  = (Mat_SeqAIJ*)(mat->B)->data;
142:   PetscInt          colA_start,jB,jcol;

145:   ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart=A->rmap->rstart;
146:   av=aa->a; bv=bb->a;

148:   garray = mat->garray;

150:   if (reuse == MAT_INITIAL_MATRIX) {
151:     nz   = aa->nz + bb->nz;
152:     *nnz = nz;
153:     PetscMalloc3(m+1,&row,nz,&col,nz,&val);
154:     *r = row; *c = col; *v = val;
155:   } else {
156:     row = *r; col = *c; val = *v;
157:   }

159:   nz = 0;
160:   for (i=0; i<m; i++) {
161:     row[i] = nz;
162:     countA     = ai[i+1] - ai[i];
163:     countB     = bi[i+1] - bi[i];
164:     ajj        = aj + ai[i]; /* ptr to the beginning of this row */
165:     bjj        = bj + bi[i];

167:     /* B part, smaller col index */
168:     colA_start = rstart + ajj[0]; /* the smallest global col index of A */
169:     jB         = 0;
170:     for (j=0; j<countB; j++) {
171:       jcol = garray[bjj[j]];
172:       if (jcol > colA_start) break;
173:       col[nz]   = jcol;
174:       val[nz++] = *bv++;
175:     }
176:     jB = j;

178:     /* A part */
179:     for (j=0; j<countA; j++) {
180:       col[nz]   = rstart + ajj[j];
181:       val[nz++] = *av++;
182:     }

184:     /* B part, larger col index */
185:     for (j=jB; j<countB; j++) {
186:       col[nz]   = garray[bjj[j]];
187:       val[nz++] = *bv++;
188:     }
189:   }
190:   row[m] = nz;

192:   return(0);
193: }

195: PetscErrorCode MatConvertToTriples_mpibaij_mpibaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v)
196: {
197:   const PetscInt    *ai, *aj, *bi, *bj,*garray,bs=A->rmap->bs,bs2=bs*bs,m=A->rmap->n/bs,*ajj,*bjj;
198:   PetscErrorCode    ierr;
199:   PetscInt          rstart,nz,i,j,countA,countB;
200:   PetscInt          *row,*col;
201:   const PetscScalar *av, *bv;
202:   PetscScalar       *val;
203:   Mat_MPIBAIJ       *mat = (Mat_MPIBAIJ*)A->data;
204:   Mat_SeqBAIJ       *aa  = (Mat_SeqBAIJ*)(mat->A)->data;
205:   Mat_SeqBAIJ       *bb  = (Mat_SeqBAIJ*)(mat->B)->data;
206:   PetscInt          colA_start,jB,jcol;

209:   ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart=A->rmap->rstart/bs;
210:   av=aa->a; bv=bb->a;

212:   garray = mat->garray;

214:   if (reuse == MAT_INITIAL_MATRIX) {
215:     nz   = aa->nz + bb->nz;
216:     *nnz = nz;
217:     PetscMalloc3(m+1,&row,nz,&col,nz*bs2,&val);
218:     *r = row; *c = col; *v = val;
219:   } else {
220:     row = *r; col = *c; val = *v;
221:   }

223:   nz = 0;
224:   for (i=0; i<m; i++) {
225:     row[i]     = nz+1;
226:     countA     = ai[i+1] - ai[i];
227:     countB     = bi[i+1] - bi[i];
228:     ajj        = aj + ai[i]; /* ptr to the beginning of this row */
229:     bjj        = bj + bi[i];

231:     /* B part, smaller col index */
232:     colA_start = rstart + (countA > 0 ? ajj[0] : 0); /* the smallest global col index of A */
233:     jB         = 0;
234:     for (j=0; j<countB; j++) {
235:       jcol = garray[bjj[j]];
236:       if (jcol > colA_start) break;
237:       col[nz++] = jcol + 1;
238:     }
239:     jB = j;
240:     PetscArraycpy(val,bv,jB*bs2);
241:     val += jB*bs2;
242:     bv  += jB*bs2;

244:     /* A part */
245:     for (j=0; j<countA; j++) col[nz++] = rstart + ajj[j] + 1;
246:     PetscArraycpy(val,av,countA*bs2);
247:     val += countA*bs2;
248:     av  += countA*bs2;

250:     /* B part, larger col index */
251:     for (j=jB; j<countB; j++) col[nz++] = garray[bjj[j]] + 1;
252:     PetscArraycpy(val,bv,(countB-jB)*bs2);
253:     val += (countB-jB)*bs2;
254:     bv  += (countB-jB)*bs2;
255:   }
256:   row[m] = nz+1;

258:   return(0);
259: }

261: PetscErrorCode MatConvertToTriples_mpisbaij_mpisbaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v)
262: {
263:   const PetscInt    *ai, *aj, *bi, *bj,*garray,bs=A->rmap->bs,bs2=bs*bs,m=A->rmap->n/bs,*ajj,*bjj;
264:   PetscErrorCode    ierr;
265:   PetscInt          rstart,nz,i,j,countA,countB;
266:   PetscInt          *row,*col;
267:   const PetscScalar *av, *bv;
268:   PetscScalar       *val;
269:   Mat_MPISBAIJ      *mat = (Mat_MPISBAIJ*)A->data;
270:   Mat_SeqSBAIJ      *aa  = (Mat_SeqSBAIJ*)(mat->A)->data;
271:   Mat_SeqBAIJ       *bb  = (Mat_SeqBAIJ*)(mat->B)->data;

274:   ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart=A->rmap->rstart/bs;
275:   av=aa->a; bv=bb->a;

277:   garray = mat->garray;

279:   if (reuse == MAT_INITIAL_MATRIX) {
280:     nz   = aa->nz + bb->nz;
281:     *nnz = nz;
282:     PetscMalloc3(m+1,&row,nz,&col,nz*bs2,&val);
283:     *r = row; *c = col; *v = val;
284:   } else {
285:     row = *r; col = *c; val = *v;
286:   }

288:   nz = 0;
289:   for (i=0; i<m; i++) {
290:     row[i]     = nz+1;
291:     countA     = ai[i+1] - ai[i];
292:     countB     = bi[i+1] - bi[i];
293:     ajj        = aj + ai[i]; /* ptr to the beginning of this row */
294:     bjj        = bj + bi[i];

296:     /* A part */
297:     for (j=0; j<countA; j++) col[nz++] = rstart + ajj[j] + 1;
298:     PetscArraycpy(val,av,countA*bs2);
299:     val += countA*bs2;
300:     av  += countA*bs2;

302:     /* B part, larger col index */
303:     for (j=0; j<countB; j++) col[nz++] = garray[bjj[j]] + 1;
304:     PetscArraycpy(val,bv,countB*bs2);
305:     val += countB*bs2;
306:     bv  += countB*bs2;
307:   }
308:   row[m] = nz+1;

310:   return(0);
311: }

313: /*
314:  * Free memory for Mat_MKL_CPARDISO structure and pointers to objects.
315:  */
316: PetscErrorCode MatDestroy_MKL_CPARDISO(Mat A)
317: {
318:   Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data;
319:   PetscErrorCode   ierr;

322:   /* Terminate instance, deallocate memories */
323:   if (mat_mkl_cpardiso->CleanUp) {
324:     mat_mkl_cpardiso->phase = JOB_RELEASE_OF_ALL_MEMORY;

326:     cluster_sparse_solver (
327:       mat_mkl_cpardiso->pt,
328:       &mat_mkl_cpardiso->maxfct,
329:       &mat_mkl_cpardiso->mnum,
330:       &mat_mkl_cpardiso->mtype,
331:       &mat_mkl_cpardiso->phase,
332:       &mat_mkl_cpardiso->n,
333:       NULL,
334:       NULL,
335:       NULL,
336:       mat_mkl_cpardiso->perm,
337:       &mat_mkl_cpardiso->nrhs,
338:       mat_mkl_cpardiso->iparm,
339:       &mat_mkl_cpardiso->msglvl,
340:       NULL,
341:       NULL,
342:       &mat_mkl_cpardiso->comm_mkl_cpardiso,
343:       (PetscInt*)&mat_mkl_cpardiso->err);
344:   }

346:   if (mat_mkl_cpardiso->ConvertToTriples != MatCopy_seqaij_seqaij_MKL_CPARDISO) {
347:     PetscFree3(mat_mkl_cpardiso->ia,mat_mkl_cpardiso->ja,mat_mkl_cpardiso->a);
348:   }
349:   MPI_Comm_free(&(mat_mkl_cpardiso->comm_mkl_cpardiso));
350:   PetscFree(A->data);

352:   /* clear composed functions */
353:   PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL);
354:   PetscObjectComposeFunction((PetscObject)A,"MatMkl_CPardisoSetCntl_C",NULL);
355:   return(0);
356: }

358: /*
359:  * Computes Ax = b
360:  */
361: PetscErrorCode MatSolve_MKL_CPARDISO(Mat A,Vec b,Vec x)
362: {
363:   Mat_MKL_CPARDISO   *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)(A)->data;
364:   PetscErrorCode    ierr;
365:   PetscScalar       *xarray;
366:   const PetscScalar *barray;

369:   mat_mkl_cpardiso->nrhs = 1;
370:   VecGetArray(x,&xarray);
371:   VecGetArrayRead(b,&barray);

373:   /* solve phase */
374:   /*-------------*/
375:   mat_mkl_cpardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
376:   cluster_sparse_solver (
377:     mat_mkl_cpardiso->pt,
378:     &mat_mkl_cpardiso->maxfct,
379:     &mat_mkl_cpardiso->mnum,
380:     &mat_mkl_cpardiso->mtype,
381:     &mat_mkl_cpardiso->phase,
382:     &mat_mkl_cpardiso->n,
383:     mat_mkl_cpardiso->a,
384:     mat_mkl_cpardiso->ia,
385:     mat_mkl_cpardiso->ja,
386:     mat_mkl_cpardiso->perm,
387:     &mat_mkl_cpardiso->nrhs,
388:     mat_mkl_cpardiso->iparm,
389:     &mat_mkl_cpardiso->msglvl,
390:     (void*)barray,
391:     (void*)xarray,
392:     &mat_mkl_cpardiso->comm_mkl_cpardiso,
393:     (PetscInt*)&mat_mkl_cpardiso->err);

395:   if (mat_mkl_cpardiso->err < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_CPARDISO: err=%d, msg = \"%s\". Please check manual\n",mat_mkl_cpardiso->err,Err_MSG_CPardiso(mat_mkl_cpardiso->err));

397:   VecRestoreArray(x,&xarray);
398:   VecRestoreArrayRead(b,&barray);
399:   mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
400:   return(0);
401: }

403: PetscErrorCode MatSolveTranspose_MKL_CPARDISO(Mat A,Vec b,Vec x)
404: {
405:   Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data;
406:   PetscErrorCode   ierr;

409: #if defined(PETSC_USE_COMPLEX)
410:   mat_mkl_cpardiso->iparm[12 - 1] = 1;
411: #else
412:   mat_mkl_cpardiso->iparm[12 - 1] = 2;
413: #endif
414:   MatSolve_MKL_CPARDISO(A,b,x);
415:   mat_mkl_cpardiso->iparm[12 - 1] = 0;
416:   return(0);
417: }

419: PetscErrorCode MatMatSolve_MKL_CPARDISO(Mat A,Mat B,Mat X)
420: {
421:   Mat_MKL_CPARDISO  *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)(A)->data;
422:   PetscErrorCode    ierr;
423:   PetscScalar       *xarray;
424:   const PetscScalar *barray;

427:   MatGetSize(B,NULL,(PetscInt*)&mat_mkl_cpardiso->nrhs);

429:   if (mat_mkl_cpardiso->nrhs > 0) {
430:     MatDenseGetArrayRead(B,&barray);
431:     MatDenseGetArray(X,&xarray);

433:     /* solve phase */
434:     /*-------------*/
435:     mat_mkl_cpardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
436:     cluster_sparse_solver (
437:       mat_mkl_cpardiso->pt,
438:       &mat_mkl_cpardiso->maxfct,
439:       &mat_mkl_cpardiso->mnum,
440:       &mat_mkl_cpardiso->mtype,
441:       &mat_mkl_cpardiso->phase,
442:       &mat_mkl_cpardiso->n,
443:       mat_mkl_cpardiso->a,
444:       mat_mkl_cpardiso->ia,
445:       mat_mkl_cpardiso->ja,
446:       mat_mkl_cpardiso->perm,
447:       &mat_mkl_cpardiso->nrhs,
448:       mat_mkl_cpardiso->iparm,
449:       &mat_mkl_cpardiso->msglvl,
450:       (void*)barray,
451:       (void*)xarray,
452:       &mat_mkl_cpardiso->comm_mkl_cpardiso,
453:       (PetscInt*)&mat_mkl_cpardiso->err);
454:     if (mat_mkl_cpardiso->err < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_CPARDISO: err=%d, msg = \"%s\". Please check manual\n",mat_mkl_cpardiso->err,Err_MSG_CPardiso(mat_mkl_cpardiso->err));
455:     MatDenseRestoreArrayRead(B,&barray);
456:     MatDenseRestoreArray(X,&xarray);

458:   }
459:   mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
460:   return(0);

462: }

464: /*
465:  * LU Decomposition
466:  */
467: PetscErrorCode MatFactorNumeric_MKL_CPARDISO(Mat F,Mat A,const MatFactorInfo *info)
468: {
469:   Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)(F)->data;
470:   PetscErrorCode   ierr;

473:   mat_mkl_cpardiso->matstruc = SAME_NONZERO_PATTERN;
474:   (*mat_mkl_cpardiso->ConvertToTriples)(A, MAT_REUSE_MATRIX,&mat_mkl_cpardiso->nz,&mat_mkl_cpardiso->ia,&mat_mkl_cpardiso->ja,&mat_mkl_cpardiso->a);

476:   mat_mkl_cpardiso->phase = JOB_NUMERICAL_FACTORIZATION;
477:   cluster_sparse_solver (
478:     mat_mkl_cpardiso->pt,
479:     &mat_mkl_cpardiso->maxfct,
480:     &mat_mkl_cpardiso->mnum,
481:     &mat_mkl_cpardiso->mtype,
482:     &mat_mkl_cpardiso->phase,
483:     &mat_mkl_cpardiso->n,
484:     mat_mkl_cpardiso->a,
485:     mat_mkl_cpardiso->ia,
486:     mat_mkl_cpardiso->ja,
487:     mat_mkl_cpardiso->perm,
488:     &mat_mkl_cpardiso->nrhs,
489:     mat_mkl_cpardiso->iparm,
490:     &mat_mkl_cpardiso->msglvl,
491:     NULL,
492:     NULL,
493:     &mat_mkl_cpardiso->comm_mkl_cpardiso,
494:     &mat_mkl_cpardiso->err);
495:   if (mat_mkl_cpardiso->err < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_CPARDISO: err=%d, msg = \"%s\". Please check manual\n",mat_mkl_cpardiso->err,Err_MSG_CPardiso(mat_mkl_cpardiso->err));

497:   mat_mkl_cpardiso->matstruc = SAME_NONZERO_PATTERN;
498:   mat_mkl_cpardiso->CleanUp  = PETSC_TRUE;
499:   return(0);
500: }

502: /* Sets mkl_cpardiso options from the options database */
503: PetscErrorCode PetscSetMKL_CPARDISOFromOptions(Mat F, Mat A)
504: {
505:   Mat_MKL_CPARDISO    *mat_mkl_cpardiso = (Mat_MKL_CPARDISO*)F->data;
506:   PetscErrorCode      ierr;
507:   PetscInt            icntl,threads;
508:   PetscBool           flg;

511:   PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MKL_CPARDISO Options","Mat");
512:   PetscOptionsInt("-mat_mkl_cpardiso_65","Number of threads to use","None",threads,&threads,&flg);
513:   if (flg) mkl_set_num_threads((int)threads);

515:   PetscOptionsInt("-mat_mkl_cpardiso_66","Maximum number of factors with identical sparsity structure that must be kept in memory at the same time","None",mat_mkl_cpardiso->maxfct,&icntl,&flg);
516:   if (flg) mat_mkl_cpardiso->maxfct = icntl;

518:   PetscOptionsInt("-mat_mkl_cpardiso_67","Indicates the actual matrix for the solution phase","None",mat_mkl_cpardiso->mnum,&icntl,&flg);
519:   if (flg) mat_mkl_cpardiso->mnum = icntl;

521:   PetscOptionsInt("-mat_mkl_cpardiso_68","Message level information","None",mat_mkl_cpardiso->msglvl,&icntl,&flg);
522:   if (flg) mat_mkl_cpardiso->msglvl = icntl;

524:   PetscOptionsInt("-mat_mkl_cpardiso_69","Defines the matrix type","None",mat_mkl_cpardiso->mtype,&icntl,&flg);
525:   if (flg) {
526:     mat_mkl_cpardiso->mtype = icntl;
527: #if defined(PETSC_USE_REAL_SINGLE)
528:     mat_mkl_cpardiso->iparm[27] = 1;
529: #else
530:     mat_mkl_cpardiso->iparm[27] = 0;
531: #endif
532:     mat_mkl_cpardiso->iparm[34] = 1;
533:   }
534:   PetscOptionsInt("-mat_mkl_cpardiso_1","Use default values","None",mat_mkl_cpardiso->iparm[0],&icntl,&flg);

536:   if (flg && icntl != 0) {
537:     PetscOptionsInt("-mat_mkl_cpardiso_2","Fill-in reducing ordering for the input matrix","None",mat_mkl_cpardiso->iparm[1],&icntl,&flg);
538:     if (flg) mat_mkl_cpardiso->iparm[1] = icntl;

540:     PetscOptionsInt("-mat_mkl_cpardiso_4","Preconditioned CGS/CG","None",mat_mkl_cpardiso->iparm[3],&icntl,&flg);
541:     if (flg) mat_mkl_cpardiso->iparm[3] = icntl;

543:     PetscOptionsInt("-mat_mkl_cpardiso_5","User permutation","None",mat_mkl_cpardiso->iparm[4],&icntl,&flg);
544:     if (flg) mat_mkl_cpardiso->iparm[4] = icntl;

546:     PetscOptionsInt("-mat_mkl_cpardiso_6","Write solution on x","None",mat_mkl_cpardiso->iparm[5],&icntl,&flg);
547:     if (flg) mat_mkl_cpardiso->iparm[5] = icntl;

549:     PetscOptionsInt("-mat_mkl_cpardiso_8","Iterative refinement step","None",mat_mkl_cpardiso->iparm[7],&icntl,&flg);
550:     if (flg) mat_mkl_cpardiso->iparm[7] = icntl;

552:     PetscOptionsInt("-mat_mkl_cpardiso_10","Pivoting perturbation","None",mat_mkl_cpardiso->iparm[9],&icntl,&flg);
553:     if (flg) mat_mkl_cpardiso->iparm[9] = icntl;

555:     PetscOptionsInt("-mat_mkl_cpardiso_11","Scaling vectors","None",mat_mkl_cpardiso->iparm[10],&icntl,&flg);
556:     if (flg) mat_mkl_cpardiso->iparm[10] = icntl;

558:     PetscOptionsInt("-mat_mkl_cpardiso_12","Solve with transposed or conjugate transposed matrix A","None",mat_mkl_cpardiso->iparm[11],&icntl,&flg);
559:     if (flg) mat_mkl_cpardiso->iparm[11] = icntl;

561:     PetscOptionsInt("-mat_mkl_cpardiso_13","Improved accuracy using (non-) symmetric weighted matching","None",mat_mkl_cpardiso->iparm[12],&icntl,
562:       &flg);
563:     if (flg) mat_mkl_cpardiso->iparm[12] = icntl;

565:     PetscOptionsInt("-mat_mkl_cpardiso_18","Numbers of non-zero elements","None",mat_mkl_cpardiso->iparm[17],&icntl,
566:       &flg);
567:     if (flg) mat_mkl_cpardiso->iparm[17] = icntl;

569:     PetscOptionsInt("-mat_mkl_cpardiso_19","Report number of floating point operations","None",mat_mkl_cpardiso->iparm[18],&icntl,&flg);
570:     if (flg) mat_mkl_cpardiso->iparm[18] = icntl;

572:     PetscOptionsInt("-mat_mkl_cpardiso_21","Pivoting for symmetric indefinite matrices","None",mat_mkl_cpardiso->iparm[20],&icntl,&flg);
573:     if (flg) mat_mkl_cpardiso->iparm[20] = icntl;

575:     PetscOptionsInt("-mat_mkl_cpardiso_24","Parallel factorization control","None",mat_mkl_cpardiso->iparm[23],&icntl,&flg);
576:     if (flg) mat_mkl_cpardiso->iparm[23] = icntl;

578:     PetscOptionsInt("-mat_mkl_cpardiso_25","Parallel forward/backward solve control","None",mat_mkl_cpardiso->iparm[24],&icntl,&flg);
579:     if (flg) mat_mkl_cpardiso->iparm[24] = icntl;

581:     PetscOptionsInt("-mat_mkl_cpardiso_27","Matrix checker","None",mat_mkl_cpardiso->iparm[26],&icntl,&flg);
582:     if (flg) mat_mkl_cpardiso->iparm[26] = icntl;

584:     PetscOptionsInt("-mat_mkl_cpardiso_31","Partial solve and computing selected components of the solution vectors","None",mat_mkl_cpardiso->iparm[30],&icntl,&flg);
585:     if (flg) mat_mkl_cpardiso->iparm[30] = icntl;

587:     PetscOptionsInt("-mat_mkl_cpardiso_34","Optimal number of threads for conditional numerical reproducibility (CNR) mode","None",mat_mkl_cpardiso->iparm[33],&icntl,&flg);
588:     if (flg) mat_mkl_cpardiso->iparm[33] = icntl;

590:     PetscOptionsInt("-mat_mkl_cpardiso_60","Intel MKL_CPARDISO mode","None",mat_mkl_cpardiso->iparm[59],&icntl,&flg);
591:     if (flg) mat_mkl_cpardiso->iparm[59] = icntl;
592:   }

594:   PetscOptionsEnd();
595:   return(0);
596: }

598: PetscErrorCode PetscInitialize_MKL_CPARDISO(Mat A, Mat_MKL_CPARDISO *mat_mkl_cpardiso)
599: {
600:   PetscErrorCode  ierr;
601:   PetscInt        bs;
602:   PetscBool       match;
603:   PetscMPIInt     size;


607:   MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(mat_mkl_cpardiso->comm_mkl_cpardiso));
608:   MPI_Comm_size(mat_mkl_cpardiso->comm_mkl_cpardiso, &size);

610:   mat_mkl_cpardiso->CleanUp = PETSC_FALSE;
611:   mat_mkl_cpardiso->maxfct = 1;
612:   mat_mkl_cpardiso->mnum = 1;
613:   mat_mkl_cpardiso->n = A->rmap->N;
614:   if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36];
615:   mat_mkl_cpardiso->msglvl = 0;
616:   mat_mkl_cpardiso->nrhs = 1;
617:   mat_mkl_cpardiso->err = 0;
618:   mat_mkl_cpardiso->phase = -1;
619: #if defined(PETSC_USE_COMPLEX)
620:   mat_mkl_cpardiso->mtype = 13;
621: #else
622:   mat_mkl_cpardiso->mtype = 11;
623: #endif

625: #if defined(PETSC_USE_REAL_SINGLE)
626:   mat_mkl_cpardiso->iparm[27] = 1;
627: #else
628:   mat_mkl_cpardiso->iparm[27] = 0;
629: #endif

631:   mat_mkl_cpardiso->iparm[ 0] =  1; /* Solver default parameters overriden with provided by iparm */
632:   mat_mkl_cpardiso->iparm[ 1] =  2; /* Use METIS for fill-in reordering */
633:   mat_mkl_cpardiso->iparm[ 5] =  0; /* Write solution into x */
634:   mat_mkl_cpardiso->iparm[ 7] =  2; /* Max number of iterative refinement steps */
635:   mat_mkl_cpardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */
636:   mat_mkl_cpardiso->iparm[10] =  1; /* Use nonsymmetric permutation and scaling MPS */
637:   mat_mkl_cpardiso->iparm[12] =  1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
638:   mat_mkl_cpardiso->iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
639:   mat_mkl_cpardiso->iparm[18] = -1; /* Output: Mflops for LU factorization */
640:   mat_mkl_cpardiso->iparm[26] =  1; /* Check input data for correctness */

642:   mat_mkl_cpardiso->iparm[39] = 0;
643:   if (size > 1) {
644:     mat_mkl_cpardiso->iparm[39] = 2;
645:     mat_mkl_cpardiso->iparm[40] = A->rmap->rstart;
646:     mat_mkl_cpardiso->iparm[41] = A->rmap->rend-1;
647:   }
648:   PetscObjectTypeCompareAny((PetscObject)A,&match,MATMPIBAIJ,MATMPISBAIJ,"");
649:   if (match) {
650:     MatGetBlockSize(A,&bs);
651:     mat_mkl_cpardiso->iparm[36] = bs;
652:     mat_mkl_cpardiso->iparm[40] /= bs;
653:     mat_mkl_cpardiso->iparm[41] /= bs;
654:     mat_mkl_cpardiso->iparm[40]++;
655:     mat_mkl_cpardiso->iparm[41]++;
656:     mat_mkl_cpardiso->iparm[34] = 0;  /* Fortran style */
657:   } else {
658:     mat_mkl_cpardiso->iparm[34] = 1;  /* C style */
659:   }

661:   mat_mkl_cpardiso->perm = 0;
662:   return(0);
663: }

665: /*
666:  * Symbolic decomposition. Mkl_Pardiso analysis phase.
667:  */
668: PetscErrorCode MatLUFactorSymbolic_AIJMKL_CPARDISO(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
669: {
670:   Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO*)F->data;
671:   PetscErrorCode  ierr;

674:   mat_mkl_cpardiso->matstruc = DIFFERENT_NONZERO_PATTERN;

676:   /* Set MKL_CPARDISO options from the options database */
677:   PetscSetMKL_CPARDISOFromOptions(F,A);
678:   (*mat_mkl_cpardiso->ConvertToTriples)(A,MAT_INITIAL_MATRIX,&mat_mkl_cpardiso->nz,&mat_mkl_cpardiso->ia,&mat_mkl_cpardiso->ja,&mat_mkl_cpardiso->a);

680:   mat_mkl_cpardiso->n = A->rmap->N;
681:   if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36];

683:   /* analysis phase */
684:   /*----------------*/
685:   mat_mkl_cpardiso->phase = JOB_ANALYSIS;

687:   cluster_sparse_solver (
688:     mat_mkl_cpardiso->pt,
689:     &mat_mkl_cpardiso->maxfct,
690:     &mat_mkl_cpardiso->mnum,
691:     &mat_mkl_cpardiso->mtype,
692:     &mat_mkl_cpardiso->phase,
693:     &mat_mkl_cpardiso->n,
694:     mat_mkl_cpardiso->a,
695:     mat_mkl_cpardiso->ia,
696:     mat_mkl_cpardiso->ja,
697:     mat_mkl_cpardiso->perm,
698:     &mat_mkl_cpardiso->nrhs,
699:     mat_mkl_cpardiso->iparm,
700:     &mat_mkl_cpardiso->msglvl,
701:     NULL,
702:     NULL,
703:     &mat_mkl_cpardiso->comm_mkl_cpardiso,
704:     (PetscInt*)&mat_mkl_cpardiso->err);

706:   if (mat_mkl_cpardiso->err < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_CPARDISO: err=%d, msg = \"%s\".Check manual\n",mat_mkl_cpardiso->err,Err_MSG_CPardiso(mat_mkl_cpardiso->err));

708:   mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
709:   F->ops->lufactornumeric = MatFactorNumeric_MKL_CPARDISO;
710:   F->ops->solve           = MatSolve_MKL_CPARDISO;
711:   F->ops->solvetranspose  = MatSolveTranspose_MKL_CPARDISO;
712:   F->ops->matsolve        = MatMatSolve_MKL_CPARDISO;
713:   return(0);
714: }

716: PetscErrorCode MatCholeskyFactorSymbolic_AIJMKL_CPARDISO(Mat F,Mat A,IS perm,const MatFactorInfo *info)
717: {
718:   Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO*)F->data;
719:   PetscErrorCode  ierr;

722:   mat_mkl_cpardiso->matstruc = DIFFERENT_NONZERO_PATTERN;


725:   /* Set MKL_CPARDISO options from the options database */
726:   PetscSetMKL_CPARDISOFromOptions(F,A);
727:   (*mat_mkl_cpardiso->ConvertToTriples)(A,MAT_INITIAL_MATRIX,&mat_mkl_cpardiso->nz,&mat_mkl_cpardiso->ia,&mat_mkl_cpardiso->ja,&mat_mkl_cpardiso->a);

729:   mat_mkl_cpardiso->n = A->rmap->N;
730:   if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36];
731: #if defined(PETSC_USE_COMPLEX)
732:   SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO CHOLESKY with complex scalars! Use MAT_FACTOR_LU instead",((PetscObject)A)->type_name);
733: #endif
734:   if (A->spd_set && A->spd) mat_mkl_cpardiso->mtype = 2;
735:   else                      mat_mkl_cpardiso->mtype = -2;

737:   /* analysis phase */
738:   /*----------------*/
739:   mat_mkl_cpardiso->phase = JOB_ANALYSIS;

741:   cluster_sparse_solver (
742:     mat_mkl_cpardiso->pt,
743:     &mat_mkl_cpardiso->maxfct,
744:     &mat_mkl_cpardiso->mnum,
745:     &mat_mkl_cpardiso->mtype,
746:     &mat_mkl_cpardiso->phase,
747:     &mat_mkl_cpardiso->n,
748:     mat_mkl_cpardiso->a,
749:     mat_mkl_cpardiso->ia,
750:     mat_mkl_cpardiso->ja,
751:     mat_mkl_cpardiso->perm,
752:     &mat_mkl_cpardiso->nrhs,
753:     mat_mkl_cpardiso->iparm,
754:     &mat_mkl_cpardiso->msglvl,
755:     NULL,
756:     NULL,
757:     &mat_mkl_cpardiso->comm_mkl_cpardiso,
758:     (PetscInt*)&mat_mkl_cpardiso->err);

760:   if (mat_mkl_cpardiso->err < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_CPARDISO: err=%d, msg = \"%s\".Check manual\n",mat_mkl_cpardiso->err,Err_MSG_CPardiso(mat_mkl_cpardiso->err));

762:   mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
763:   F->ops->choleskyfactornumeric = MatFactorNumeric_MKL_CPARDISO;
764:   F->ops->solve                 = MatSolve_MKL_CPARDISO;
765:   F->ops->solvetranspose        = MatSolveTranspose_MKL_CPARDISO;
766:   F->ops->matsolve              = MatMatSolve_MKL_CPARDISO;
767:   return(0);
768: }

770: PetscErrorCode MatView_MKL_CPARDISO(Mat A, PetscViewer viewer)
771: {
772:   PetscErrorCode    ierr;
773:   PetscBool         iascii;
774:   PetscViewerFormat format;
775:   Mat_MKL_CPARDISO  *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data;
776:   PetscInt          i;

779:   /* check if matrix is mkl_cpardiso type */
780:   if (A->ops->solve != MatSolve_MKL_CPARDISO) return(0);

782:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
783:   if (iascii) {
784:     PetscViewerGetFormat(viewer,&format);
785:     if (format == PETSC_VIEWER_ASCII_INFO) {
786:       PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO run parameters:\n");
787:       PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO phase:             %d \n",mat_mkl_cpardiso->phase);
788:       for (i = 1; i <= 64; i++) {
789:         PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO iparm[%d]:     %d \n",i, mat_mkl_cpardiso->iparm[i - 1]);
790:       }
791:       PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO maxfct:     %d \n", mat_mkl_cpardiso->maxfct);
792:       PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO mnum:     %d \n", mat_mkl_cpardiso->mnum);
793:       PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO mtype:     %d \n", mat_mkl_cpardiso->mtype);
794:       PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO n:     %d \n", mat_mkl_cpardiso->n);
795:       PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO nrhs:     %d \n", mat_mkl_cpardiso->nrhs);
796:       PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO msglvl:     %d \n", mat_mkl_cpardiso->msglvl);
797:     }
798:   }
799:   return(0);
800: }

802: PetscErrorCode MatGetInfo_MKL_CPARDISO(Mat A, MatInfoType flag, MatInfo *info)
803: {
804:   Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data;

807:   info->block_size        = 1.0;
808:   info->nz_allocated      = mat_mkl_cpardiso->nz + 0.0;
809:   info->nz_unneeded       = 0.0;
810:   info->assemblies        = 0.0;
811:   info->mallocs           = 0.0;
812:   info->memory            = 0.0;
813:   info->fill_ratio_given  = 0;
814:   info->fill_ratio_needed = 0;
815:   info->factor_mallocs    = 0;
816:   return(0);
817: }

819: PetscErrorCode MatMkl_CPardisoSetCntl_MKL_CPARDISO(Mat F,PetscInt icntl,PetscInt ival)
820: {
821:   Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)F->data;

824:   if (icntl <= 64) {
825:     mat_mkl_cpardiso->iparm[icntl - 1] = ival;
826:   } else {
827:     if (icntl == 65) mkl_set_num_threads((int)ival);
828:     else if (icntl == 66) mat_mkl_cpardiso->maxfct = ival;
829:     else if (icntl == 67) mat_mkl_cpardiso->mnum = ival;
830:     else if (icntl == 68) mat_mkl_cpardiso->msglvl = ival;
831:     else if (icntl == 69) {
832:       mat_mkl_cpardiso->mtype = ival;
833: #if defined(PETSC_USE_REAL_SINGLE)
834:       mat_mkl_cpardiso->iparm[27] = 1;
835: #else
836:       mat_mkl_cpardiso->iparm[27] = 0;
837: #endif
838:       mat_mkl_cpardiso->iparm[34] = 1;
839:     }
840:   }
841:   return(0);
842: }

844: /*@
845:   MatMkl_CPardisoSetCntl - Set Mkl_Pardiso parameters

847:    Logically Collective on Mat

849:    Input Parameters:
850: +  F - the factored matrix obtained by calling MatGetFactor()
851: .  icntl - index of Mkl_Pardiso parameter
852: -  ival - value of Mkl_Pardiso parameter

854:   Options Database:
855: .   -mat_mkl_cpardiso_<icntl> <ival>

857:    Level: Intermediate

859:    Notes:
860:     This routine cannot be used if you are solving the linear system with TS, SNES, or KSP, only if you directly call MatGetFactor() so use the options 
861:           database approach when working with TS, SNES, or KSP.

863:    References:
864: .      Mkl_Pardiso Users' Guide

866: .seealso: MatGetFactor()
867: @*/
868: PetscErrorCode MatMkl_CPardisoSetCntl(Mat F,PetscInt icntl,PetscInt ival)
869: {

873:   PetscTryMethod(F,"MatMkl_CPardisoSetCntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));
874:   return(0);
875: }

877: static PetscErrorCode MatFactorGetSolverType_mkl_cpardiso(Mat A, MatSolverType *type)
878: {
880:   *type = MATSOLVERMKL_CPARDISO;
881:   return(0);
882: }

884: /* MatGetFactor for MPI AIJ matrices */
885: static PetscErrorCode MatGetFactor_mpiaij_mkl_cpardiso(Mat A,MatFactorType ftype,Mat *F)
886: {
887:   Mat              B;
888:   PetscErrorCode   ierr;
889:   Mat_MKL_CPARDISO *mat_mkl_cpardiso;
890:   PetscBool        isSeqAIJ,isMPIBAIJ,isMPISBAIJ;

893:   /* Create the factorization matrix */

895:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);
896:   PetscObjectTypeCompare((PetscObject)A,MATMPIBAIJ,&isMPIBAIJ);
897:   PetscObjectTypeCompare((PetscObject)A,MATMPISBAIJ,&isMPISBAIJ);
898:   MatCreate(PetscObjectComm((PetscObject)A),&B);
899:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
900:   PetscStrallocpy("mkl_cpardiso",&((PetscObject)B)->type_name);
901:   MatSetUp(B);

903:   PetscNewLog(B,&mat_mkl_cpardiso);

905:   if (isSeqAIJ) mat_mkl_cpardiso->ConvertToTriples = MatCopy_seqaij_seqaij_MKL_CPARDISO;
906:   else if (isMPIBAIJ) mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpibaij_mpibaij_MKL_CPARDISO;
907:   else if (isMPISBAIJ) mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij_MKL_CPARDISO;
908:   else          mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij_MKL_CPARDISO;

910:   if (ftype == MAT_FACTOR_LU) B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_CPARDISO;
911:   else B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJMKL_CPARDISO;
912:   B->ops->destroy = MatDestroy_MKL_CPARDISO;

914:   B->ops->view    = MatView_MKL_CPARDISO;
915:   B->ops->getinfo = MatGetInfo_MKL_CPARDISO;

917:   B->factortype   = ftype;
918:   B->assembled    = PETSC_TRUE;           /* required by -ksp_view */

920:   B->data = mat_mkl_cpardiso;

922:   /* set solvertype */
923:   PetscFree(B->solvertype);
924:   PetscStrallocpy(MATSOLVERMKL_CPARDISO,&B->solvertype);

926:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_mkl_cpardiso);
927:   PetscObjectComposeFunction((PetscObject)B,"MatMkl_CPardisoSetCntl_C",MatMkl_CPardisoSetCntl_MKL_CPARDISO);
928:   PetscInitialize_MKL_CPARDISO(A, mat_mkl_cpardiso);

930:   *F = B;
931:   return(0);
932: }

934: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_MKL_CPardiso(void)
935: {
937: 
939:   MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATMPIAIJ,MAT_FACTOR_LU,MatGetFactor_mpiaij_mkl_cpardiso);
940:   MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_mpiaij_mkl_cpardiso);
941:   MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATMPIBAIJ,MAT_FACTOR_LU,MatGetFactor_mpiaij_mkl_cpardiso);
942:   MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATMPISBAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_mpiaij_mkl_cpardiso);
943:   return(0);
944: }