Actual source code: baijmkl.c

  1: /*
  2:   Defines basic operations for the MATSEQBAIJMKL matrix class.
  3:   Uses sparse BLAS operations from the Intel Math Kernel Library (MKL)
  4:   wherever possible. If used MKL verion is older than 11.3 PETSc default
  5:   code for sparse matrix operations is used.
  6: */

  8: #include <../src/mat/impls/baij/seq/baij.h>
  9: #include <../src/mat/impls/baij/seq/baijmkl/baijmkl.h>
 10: #include <mkl_spblas.h>

 12: static PetscBool PetscSeqBAIJSupportsZeroBased(void)
 13: {
 14:   static PetscBool set = PETSC_FALSE,value;
 15:   int              n=1,ia[1],ja[1];
 16:   float            a[1];
 17:   sparse_status_t  status;
 18:   sparse_matrix_t  A;

 20:   if (!set) {
 21:     status = mkl_sparse_s_create_bsr(&A,SPARSE_INDEX_BASE_ZERO,SPARSE_LAYOUT_COLUMN_MAJOR,n,n,n,ia,ia,ja,a);
 22:     value  = (status != SPARSE_STATUS_NOT_SUPPORTED) ? PETSC_TRUE : PETSC_FALSE;
 23:     (void)   mkl_sparse_destroy(A);
 24:     set    = PETSC_TRUE;
 25:   }
 26:   return value;
 27: }

 29: typedef struct {
 30:   PetscBool           sparse_optimized; /* If PETSC_TRUE, then mkl_sparse_optimize() has been called. */
 31:   sparse_matrix_t     bsrA; /* "Handle" used by SpMV2 inspector-executor routines. */
 32:   struct matrix_descr descr;
 33:   PetscInt            *ai1;
 34:   PetscInt            *aj1;
 35: } Mat_SeqBAIJMKL;

 37: static PetscErrorCode MatAssemblyEnd_SeqBAIJMKL(Mat A, MatAssemblyType mode);
 38: extern PetscErrorCode MatAssemblyEnd_SeqBAIJ(Mat,MatAssemblyType);

 40: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJMKL_SeqBAIJ(Mat A,MatType type,MatReuse reuse,Mat *newmat)
 41: {
 42:   /* This routine is only called to convert a MATBAIJMKL to its base PETSc type, */
 43:   /* so we will ignore 'MatType type'. */
 44:   Mat            B        = *newmat;
 45:   Mat_SeqBAIJMKL *baijmkl = (Mat_SeqBAIJMKL*)A->spptr;

 47:   if (reuse == MAT_INITIAL_MATRIX) MatDuplicate(A,MAT_COPY_VALUES,&B);

 49:   /* Reset the original function pointers. */
 50:   B->ops->duplicate        = MatDuplicate_SeqBAIJ;
 51:   B->ops->assemblyend      = MatAssemblyEnd_SeqBAIJ;
 52:   B->ops->destroy          = MatDestroy_SeqBAIJ;
 53:   B->ops->multtranspose    = MatMultTranspose_SeqBAIJ;
 54:   B->ops->multtransposeadd = MatMultTransposeAdd_SeqBAIJ;
 55:   B->ops->scale            = MatScale_SeqBAIJ;
 56:   B->ops->diagonalscale    = MatDiagonalScale_SeqBAIJ;
 57:   B->ops->axpy             = MatAXPY_SeqBAIJ;

 59:   switch (A->rmap->bs) {
 60:     case 1:
 61:       B->ops->mult    = MatMult_SeqBAIJ_1;
 62:       B->ops->multadd = MatMultAdd_SeqBAIJ_1;
 63:       break;
 64:     case 2:
 65:       B->ops->mult    = MatMult_SeqBAIJ_2;
 66:       B->ops->multadd = MatMultAdd_SeqBAIJ_2;
 67:       break;
 68:     case 3:
 69:       B->ops->mult    = MatMult_SeqBAIJ_3;
 70:       B->ops->multadd = MatMultAdd_SeqBAIJ_3;
 71:       break;
 72:     case 4:
 73:       B->ops->mult    = MatMult_SeqBAIJ_4;
 74:       B->ops->multadd = MatMultAdd_SeqBAIJ_4;
 75:       break;
 76:     case 5:
 77:       B->ops->mult    = MatMult_SeqBAIJ_5;
 78:       B->ops->multadd = MatMultAdd_SeqBAIJ_5;
 79:       break;
 80:     case 6:
 81:       B->ops->mult    = MatMult_SeqBAIJ_6;
 82:       B->ops->multadd = MatMultAdd_SeqBAIJ_6;
 83:       break;
 84:     case 7:
 85:       B->ops->mult    = MatMult_SeqBAIJ_7;
 86:       B->ops->multadd = MatMultAdd_SeqBAIJ_7;
 87:       break;
 88:     case 11:
 89:       B->ops->mult    = MatMult_SeqBAIJ_11;
 90:       B->ops->multadd = MatMultAdd_SeqBAIJ_11;
 91:       break;
 92:     case 15:
 93:       B->ops->mult    = MatMult_SeqBAIJ_15_ver1;
 94:       B->ops->multadd = MatMultAdd_SeqBAIJ_N;
 95:       break;
 96:     default:
 97:       B->ops->mult    = MatMult_SeqBAIJ_N;
 98:       B->ops->multadd = MatMultAdd_SeqBAIJ_N;
 99:       break;
100:   }
101:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqbaijmkl_seqbaij_C",NULL);

103:   /* Free everything in the Mat_SeqBAIJMKL data structure. Currently, this
104:    * simply involves destroying the MKL sparse matrix handle and then freeing
105:    * the spptr pointer. */
106:   if (reuse == MAT_INITIAL_MATRIX) baijmkl = (Mat_SeqBAIJMKL*)B->spptr;

108:   if (baijmkl->sparse_optimized) PetscStackCallStandard(mkl_sparse_destroy,baijmkl->bsrA);
109:   PetscFree2(baijmkl->ai1,baijmkl->aj1);
110:   PetscFree(B->spptr);

112:   /* Change the type of B to MATSEQBAIJ. */
113:   PetscObjectChangeTypeName((PetscObject)B, MATSEQBAIJ);

115:   *newmat = B;
116:   return 0;
117: }

119: static PetscErrorCode MatDestroy_SeqBAIJMKL(Mat A)
120: {
121:   Mat_SeqBAIJMKL *baijmkl = (Mat_SeqBAIJMKL*) A->spptr;

123:   if (baijmkl) {
124:     /* Clean up everything in the Mat_SeqBAIJMKL data structure, then free A->spptr. */
125:     if (baijmkl->sparse_optimized) PetscStackCallStandard(mkl_sparse_destroy,baijmkl->bsrA);
126:     PetscFree2(baijmkl->ai1,baijmkl->aj1);
127:     PetscFree(A->spptr);
128:   }

130:   /* Change the type of A back to SEQBAIJ and use MatDestroy_SeqBAIJ()
131:    * to destroy everything that remains. */
132:   PetscObjectChangeTypeName((PetscObject)A, MATSEQBAIJ);
133:   MatDestroy_SeqBAIJ(A);
134:   return 0;
135: }

137: static PetscErrorCode MatSeqBAIJMKL_create_mkl_handle(Mat A)
138: {
139:   Mat_SeqBAIJ     *a = (Mat_SeqBAIJ*)A->data;
140:   Mat_SeqBAIJMKL  *baijmkl = (Mat_SeqBAIJMKL*)A->spptr;
141:   PetscInt        mbs, nbs, nz, bs;
142:   MatScalar       *aa;
143:   PetscInt        *aj,*ai;
144:   PetscInt        i;

146:   if (baijmkl->sparse_optimized) {
147:     /* Matrix has been previously assembled and optimized. Must destroy old
148:      * matrix handle before running the optimization step again. */
149:     PetscFree2(baijmkl->ai1,baijmkl->aj1);
150:     mkl_sparse_destroy(baijmkl->bsrA);
151:   }
152:   baijmkl->sparse_optimized = PETSC_FALSE;

154:   /* Now perform the SpMV2 setup and matrix optimization. */
155:   baijmkl->descr.type        = SPARSE_MATRIX_TYPE_GENERAL;
156:   baijmkl->descr.mode        = SPARSE_FILL_MODE_LOWER;
157:   baijmkl->descr.diag        = SPARSE_DIAG_NON_UNIT;
158:   mbs = a->mbs;
159:   nbs = a->nbs;
160:   nz  = a->nz;
161:   bs  = A->rmap->bs;
162:   aa  = a->a;

164:   if ((nz!=0) & !(A->structure_only)) {
165:     /* Create a new, optimized sparse matrix handle only if the matrix has nonzero entries.
166:      * The MKL sparse-inspector executor routines don't like being passed an empty matrix. */
167:     if (PetscSeqBAIJSupportsZeroBased()) {
168:       aj   = a->j;
169:       ai   = a->i;
170:       mkl_sparse_x_create_bsr(&(baijmkl->bsrA),SPARSE_INDEX_BASE_ZERO,SPARSE_LAYOUT_COLUMN_MAJOR,mbs,nbs,bs,ai,ai+1,aj,aa);
171:     } else {
172:       PetscMalloc2(mbs+1,&ai,nz,&aj);
173:       for (i=0;i<mbs+1;i++) ai[i] = a->i[i]+1;
174:       for (i=0;i<nz;i++) aj[i] = a->j[i]+1;
175:       aa   = a->a;
176:       mkl_sparse_x_create_bsr(&baijmkl->bsrA,SPARSE_INDEX_BASE_ONE,SPARSE_LAYOUT_COLUMN_MAJOR,mbs,nbs,bs,ai,ai+1,aj,aa);
177:       baijmkl->ai1 = ai;
178:       baijmkl->aj1 = aj;
179:     }
180:     mkl_sparse_set_mv_hint(baijmkl->bsrA,SPARSE_OPERATION_NON_TRANSPOSE,baijmkl->descr,1000);
181:     mkl_sparse_set_memory_hint(baijmkl->bsrA,SPARSE_MEMORY_AGGRESSIVE);
182:     mkl_sparse_optimize(baijmkl->bsrA);
183:     baijmkl->sparse_optimized = PETSC_TRUE;
184:   }
185:   return 0;
186: }

188: static PetscErrorCode MatDuplicate_SeqBAIJMKL(Mat A, MatDuplicateOption op, Mat *M)
189: {
190:   Mat_SeqBAIJMKL *baijmkl;
191:   Mat_SeqBAIJMKL *baijmkl_dest;

193:   MatDuplicate_SeqBAIJ(A,op,M);
194:   baijmkl = (Mat_SeqBAIJMKL*) A->spptr;
195:   PetscNewLog((*M),&baijmkl_dest);
196:   (*M)->spptr = (void*)baijmkl_dest;
197:   PetscMemcpy(baijmkl_dest,baijmkl,sizeof(Mat_SeqBAIJMKL));
198:   baijmkl_dest->sparse_optimized = PETSC_FALSE;
199:   MatSeqBAIJMKL_create_mkl_handle(A);
200:   return 0;
201: }

203: static PetscErrorCode MatMult_SeqBAIJMKL_SpMV2(Mat A,Vec xx,Vec yy)
204: {
205:   Mat_SeqBAIJ        *a = (Mat_SeqBAIJ*)A->data;
206:   Mat_SeqBAIJMKL     *baijmkl=(Mat_SeqBAIJMKL*)A->spptr;
207:   const PetscScalar  *x;
208:   PetscScalar        *y;

210:   /* If there are no nonzero entries, zero yy and return immediately. */
211:   if (!a->nz) {
212:     VecSet(yy,0.0);
213:     return 0;
214:   }

216:   VecGetArrayRead(xx,&x);
217:   VecGetArray(yy,&y);

219:   /* In some cases, we get to this point without mkl_sparse_optimize() having been called, so we check and then call
220:    * it if needed. Eventually, when everything in PETSc is properly updating the matrix state, we should probably
221:    * take a "lazy" approach to creation/updating of the MKL matrix handle and plan to always do it here (when needed). */
222:   if (!baijmkl->sparse_optimized) {
223:     MatSeqBAIJMKL_create_mkl_handle(A);
224:   }

226:   /* Call MKL SpMV2 executor routine to do the MatMult. */
227:   mkl_sparse_x_mv(SPARSE_OPERATION_NON_TRANSPOSE,1.0,baijmkl->bsrA,baijmkl->descr,x,0.0,y);

229:   PetscLogFlops(2.0*a->bs2*a->nz - a->nonzerorowcnt*A->rmap->bs);
230:   VecRestoreArrayRead(xx,&x);
231:   VecRestoreArray(yy,&y);
232:   return 0;
233: }

235: static PetscErrorCode MatMultTranspose_SeqBAIJMKL_SpMV2(Mat A,Vec xx,Vec yy)
236: {
237:   Mat_SeqBAIJ       *a       = (Mat_SeqBAIJ*)A->data;
238:   Mat_SeqBAIJMKL    *baijmkl = (Mat_SeqBAIJMKL*)A->spptr;
239:   const PetscScalar *x;
240:   PetscScalar       *y;

242:   /* If there are no nonzero entries, zero yy and return immediately. */
243:   if (!a->nz) {
244:     VecSet(yy,0.0);
245:     return 0;
246:   }

248:   VecGetArrayRead(xx,&x);
249:   VecGetArray(yy,&y);

251:   /* In some cases, we get to this point without mkl_sparse_optimize() having been called, so we check and then call
252:    * it if needed. Eventually, when everything in PETSc is properly updating the matrix state, we should probably
253:    * take a "lazy" approach to creation/updating of the MKL matrix handle and plan to always do it here (when needed). */
254:   if (!baijmkl->sparse_optimized) {
255:     MatSeqBAIJMKL_create_mkl_handle(A);
256:   }

258:   /* Call MKL SpMV2 executor routine to do the MatMultTranspose. */
259:   mkl_sparse_x_mv(SPARSE_OPERATION_TRANSPOSE,1.0,baijmkl->bsrA,baijmkl->descr,x,0.0,y);

261:   PetscLogFlops(2.0*a->bs2*a->nz - a->nonzerorowcnt*A->rmap->bs);
262:   VecRestoreArrayRead(xx,&x);
263:   VecRestoreArray(yy,&y);
264:   return 0;
265: }

267: static PetscErrorCode MatMultAdd_SeqBAIJMKL_SpMV2(Mat A,Vec xx,Vec yy,Vec zz)
268: {
269:   Mat_SeqBAIJ        *a       = (Mat_SeqBAIJ*)A->data;
270:   Mat_SeqBAIJMKL     *baijmkl = (Mat_SeqBAIJMKL*)A->spptr;
271:   const PetscScalar  *x;
272:   PetscScalar        *y,*z;
273:   PetscInt           m=a->mbs*A->rmap->bs;
274:   PetscInt           i;

276:   /* If there are no nonzero entries, set zz = yy and return immediately. */
277:   if (!a->nz) {
278:     VecCopy(yy,zz);
279:     return 0;
280:   }

282:   VecGetArrayRead(xx,&x);
283:   VecGetArrayPair(yy,zz,&y,&z);

285:   /* In some cases, we get to this point without mkl_sparse_optimize() having been called, so we check and then call
286:    * it if needed. Eventually, when everything in PETSc is properly updating the matrix state, we should probably
287:    * take a "lazy" approach to creation/updating of the MKL matrix handle and plan to always do it here (when needed). */
288:   if (!baijmkl->sparse_optimized) {
289:     MatSeqBAIJMKL_create_mkl_handle(A);
290:   }

292:   /* Call MKL sparse BLAS routine to do the MatMult. */
293:   if (zz == yy) {
294:     /* If zz and yy are the same vector, we can use mkl_sparse_x_mv, which calculates y = alpha*A*x + beta*y,
295:      * with alpha and beta both set to 1.0. */
296:     mkl_sparse_x_mv(SPARSE_OPERATION_NON_TRANSPOSE,1.0,baijmkl->bsrA,baijmkl->descr,x,1.0,z);
297:   } else {
298:     /* zz and yy are different vectors, so we call mkl_sparse_x_mv with alpha=1.0 and beta=0.0, and then
299:      * we add the contents of vector yy to the result; MKL sparse BLAS does not have a MatMultAdd equivalent. */
300:     mkl_sparse_x_mv(SPARSE_OPERATION_NON_TRANSPOSE,1.0,baijmkl->bsrA,baijmkl->descr,x,0.0,z);
301:     for (i=0; i<m; i++) {
302:       z[i] += y[i];
303:     }
304:   }

306:   PetscLogFlops(2.0*a->bs2*a->nz);
307:   VecRestoreArrayRead(xx,&x);
308:   VecRestoreArrayPair(yy,zz,&y,&z);
309:   return 0;
310: }

312: static PetscErrorCode MatMultTransposeAdd_SeqBAIJMKL_SpMV2(Mat A,Vec xx,Vec yy,Vec zz)
313: {
314:   Mat_SeqBAIJ       *a       = (Mat_SeqBAIJ*)A->data;
315:   Mat_SeqBAIJMKL    *baijmkl = (Mat_SeqBAIJMKL*)A->spptr;
316:   const PetscScalar *x;
317:   PetscScalar       *y,*z;
318:   PetscInt          n=a->nbs*A->rmap->bs;
319:   PetscInt          i;
320:   /* Variables not in MatMultTransposeAdd_SeqBAIJ. */

322:   /* If there are no nonzero entries, set zz = yy and return immediately. */
323:   if (!a->nz) {
324:     VecCopy(yy,zz);
325:     return 0;
326:   }

328:   VecGetArrayRead(xx,&x);
329:   VecGetArrayPair(yy,zz,&y,&z);

331:   /* In some cases, we get to this point without mkl_sparse_optimize() having been called, so we check and then call
332:    * it if needed. Eventually, when everything in PETSc is properly updating the matrix state, we should probably
333:    * take a "lazy" approach to creation/updating of the MKL matrix handle and plan to always do it here (when needed). */
334:   if (!baijmkl->sparse_optimized) {
335:     MatSeqBAIJMKL_create_mkl_handle(A);
336:   }

338:   /* Call MKL sparse BLAS routine to do the MatMult. */
339:   if (zz == yy) {
340:     /* If zz and yy are the same vector, we can use mkl_sparse_x_mv, which calculates y = alpha*A*x + beta*y,
341:      * with alpha and beta both set to 1.0. */
342:     mkl_sparse_x_mv(SPARSE_OPERATION_TRANSPOSE,1.0,baijmkl->bsrA,baijmkl->descr,x,1.0,z);
343:   } else {
344:     /* zz and yy are different vectors, so we call mkl_sparse_x_mv with alpha=1.0 and beta=0.0, and then
345:      * we add the contents of vector yy to the result; MKL sparse BLAS does not have a MatMultAdd equivalent. */
346:     mkl_sparse_x_mv(SPARSE_OPERATION_TRANSPOSE,1.0,baijmkl->bsrA,baijmkl->descr,x,0.0,z);
347:     for (i=0; i<n; i++) {
348:       z[i] += y[i];
349:     }
350:   }

352:   PetscLogFlops(2.0*a->bs2*a->nz);
353:   VecRestoreArrayRead(xx,&x);
354:   VecRestoreArrayPair(yy,zz,&y,&z);
355:   return 0;
356: }

358: static PetscErrorCode MatScale_SeqBAIJMKL(Mat inA,PetscScalar alpha)
359: {
360:   MatScale_SeqBAIJ(inA,alpha);
361:   MatSeqBAIJMKL_create_mkl_handle(inA);
362:   return 0;
363: }

365: static PetscErrorCode MatDiagonalScale_SeqBAIJMKL(Mat A,Vec ll,Vec rr)
366: {
367:   MatDiagonalScale_SeqBAIJ(A,ll,rr);
368:   MatSeqBAIJMKL_create_mkl_handle(A);
369:   return 0;
370: }

372: static PetscErrorCode MatAXPY_SeqBAIJMKL(Mat Y,PetscScalar a,Mat X,MatStructure str)
373: {
374:   MatAXPY_SeqBAIJ(Y,a,X,str);
375:   if (str == SAME_NONZERO_PATTERN) {
376:     /* MatAssemblyEnd() is not called if SAME_NONZERO_PATTERN, so we need to force update of the MKL matrix handle. */
377:     MatSeqBAIJMKL_create_mkl_handle(Y);
378:   }
379:   return 0;
380: }
381: /* MatConvert_SeqBAIJ_SeqBAIJMKL converts a SeqBAIJ matrix into a
382:  * SeqBAIJMKL matrix.  This routine is called by the MatCreate_SeqMKLBAIJ()
383:  * routine, but can also be used to convert an assembled SeqBAIJ matrix
384:  * into a SeqBAIJMKL one. */
385: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqBAIJMKL(Mat A,MatType type,MatReuse reuse,Mat *newmat)
386: {
387:   Mat            B = *newmat;
388:   Mat_SeqBAIJMKL *baijmkl;
389:   PetscBool      sametype;

391:   if (reuse == MAT_INITIAL_MATRIX) MatDuplicate(A,MAT_COPY_VALUES,&B);

393:   PetscObjectTypeCompare((PetscObject)A,type,&sametype);
394:   if (sametype) return 0;

396:   PetscNewLog(B,&baijmkl);
397:   B->spptr = (void*)baijmkl;

399:   /* Set function pointers for methods that we inherit from BAIJ but override.
400:    * We also parse some command line options below, since those determine some of the methods we point to. */
401:   B->ops->assemblyend      = MatAssemblyEnd_SeqBAIJMKL;

403:   baijmkl->sparse_optimized = PETSC_FALSE;

405:   PetscObjectComposeFunction((PetscObject)B,"MatScale_SeqBAIJMKL_C",MatScale_SeqBAIJMKL);
406:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqbaijmkl_seqbaij_C",MatConvert_SeqBAIJMKL_SeqBAIJ);

408:   PetscObjectChangeTypeName((PetscObject)B,MATSEQBAIJMKL);
409:   *newmat = B;
410:   return 0;
411: }

413: static PetscErrorCode MatAssemblyEnd_SeqBAIJMKL(Mat A, MatAssemblyType mode)
414: {
415:   if (mode == MAT_FLUSH_ASSEMBLY) return 0;
416:   MatAssemblyEnd_SeqBAIJ(A, mode);
417:   MatSeqBAIJMKL_create_mkl_handle(A);
418:   A->ops->destroy          = MatDestroy_SeqBAIJMKL;
419:   A->ops->mult             = MatMult_SeqBAIJMKL_SpMV2;
420:   A->ops->multtranspose    = MatMultTranspose_SeqBAIJMKL_SpMV2;
421:   A->ops->multadd          = MatMultAdd_SeqBAIJMKL_SpMV2;
422:   A->ops->multtransposeadd = MatMultTransposeAdd_SeqBAIJMKL_SpMV2;
423:   A->ops->scale            = MatScale_SeqBAIJMKL;
424:   A->ops->diagonalscale    = MatDiagonalScale_SeqBAIJMKL;
425:   A->ops->axpy             = MatAXPY_SeqBAIJMKL;
426:   A->ops->duplicate        = MatDuplicate_SeqBAIJMKL;
427:   return 0;
428: }

430: /*@C
431:    MatCreateSeqBAIJMKL - Creates a sparse matrix of type SEQBAIJMKL.
432:    This type inherits from BAIJ and is largely identical, but uses sparse BLAS
433:    routines from Intel MKL whenever possible.
434:    MatMult, MatMultAdd, MatMultTranspose, and MatMultTransposeAdd
435:    operations are currently supported.
436:    If the installed version of MKL supports the "SpMV2" sparse
437:    inspector-executor routines, then those are used by default.
438:    Default PETSc kernels are used otherwise.

440:    Input Parameters:
441: +  comm - MPI communicator, set to PETSC_COMM_SELF
442: .  bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row
443:           blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
444: .  m - number of rows
445: .  n - number of columns
446: .  nz - number of nonzero blocks  per block row (same for all rows)
447: -  nnz - array containing the number of nonzero blocks in the various block rows
448:          (possibly different for each block row) or NULL

450:    Output Parameter:
451: .  A - the matrix

453:    It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
454:    MatXXXXSetPreallocation() paradigm instead of this routine directly.
455:    [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation]

457:    Options Database Keys:
458: +   -mat_no_unroll - uses code that does not unroll the loops in the
459:                      block calculations (much slower)
460: -   -mat_block_size - size of the blocks to use

462:    Level: intermediate

464:    Notes:
465:    The number of rows and columns must be divisible by blocksize.

467:    If the nnz parameter is given then the nz parameter is ignored

469:    A nonzero block is any block that as 1 or more nonzeros in it

471:    The block AIJ format is fully compatible with standard Fortran 77
472:    storage.  That is, the stored row and column indices can begin at
473:    either one (as in Fortran) or zero.  See the users' manual for details.

475:    Specify the preallocated storage with either nz or nnz (not both).
476:    Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
477:    allocation.  See Users-Manual: ch_mat for details.
478:    matrices.

480: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateBAIJ()

482: @*/
483: PetscErrorCode  MatCreateSeqBAIJMKL(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A)
484: {
485:   MatCreate(comm,A);
486:   MatSetSizes(*A,m,n,m,n);
487:   MatSetType(*A,MATSEQBAIJMKL);
488:   MatSeqBAIJSetPreallocation_SeqBAIJ(*A,bs,nz,(PetscInt*)nnz);
489:   return 0;
490: }

492: PETSC_EXTERN PetscErrorCode MatCreate_SeqBAIJMKL(Mat A)
493: {
494:   MatSetType(A,MATSEQBAIJ);
495:   MatConvert_SeqBAIJ_SeqBAIJMKL(A,MATSEQBAIJMKL,MAT_INPLACE_MATRIX,&A);
496:   return 0;
497: }