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 version 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: #if defined(PETSC_HAVE_MKL_INTEL_ILP64)
 11:   #define MKL_ILP64
 12: #endif
 13: #include <mkl_spblas.h>

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

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

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

 40: static PetscErrorCode MatAssemblyEnd_SeqBAIJMKL(Mat A, MatAssemblyType mode);
 41: extern PetscErrorCode MatAssemblyEnd_SeqBAIJ(Mat, MatAssemblyType);

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

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

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

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

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

111:   if (baijmkl->sparse_optimized) PetscCallExternal(mkl_sparse_destroy, baijmkl->bsrA);
112:   PetscFree2(baijmkl->ai1, baijmkl->aj1);
113:   PetscFree(B->spptr);

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

118:   *newmat = B;
119:   return 0;
120: }

122: static PetscErrorCode MatDestroy_SeqBAIJMKL(Mat A)
123: {
124:   Mat_SeqBAIJMKL *baijmkl = (Mat_SeqBAIJMKL *)A->spptr;

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

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

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

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

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

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

191: static PetscErrorCode MatDuplicate_SeqBAIJMKL(Mat A, MatDuplicateOption op, Mat *M)
192: {
193:   Mat_SeqBAIJMKL *baijmkl;
194:   Mat_SeqBAIJMKL *baijmkl_dest;

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

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

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

219:   VecGetArrayRead(xx, &x);
220:   VecGetArray(yy, &y);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

301:   PetscLogFlops(2.0 * a->bs2 * a->nz);
302:   VecRestoreArrayRead(xx, &x);
303:   VecRestoreArrayPair(yy, zz, &y, &z);
304:   return 0;
305: }

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

317:   /* If there are no nonzero entries, set zz = yy and return immediately. */
318:   if (!a->nz) {
319:     VecCopy(yy, zz);
320:     return 0;
321:   }

323:   VecGetArrayRead(xx, &x);
324:   VecGetArrayPair(yy, zz, &y, &z);

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

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

343:   PetscLogFlops(2.0 * a->bs2 * a->nz);
344:   VecRestoreArrayRead(xx, &x);
345:   VecRestoreArrayPair(yy, zz, &y, &z);
346:   return 0;
347: }

349: static PetscErrorCode MatScale_SeqBAIJMKL(Mat inA, PetscScalar alpha)
350: {
351:   MatScale_SeqBAIJ(inA, alpha);
352:   MatSeqBAIJMKL_create_mkl_handle(inA);
353:   return 0;
354: }

356: static PetscErrorCode MatDiagonalScale_SeqBAIJMKL(Mat A, Vec ll, Vec rr)
357: {
358:   MatDiagonalScale_SeqBAIJ(A, ll, rr);
359:   MatSeqBAIJMKL_create_mkl_handle(A);
360:   return 0;
361: }

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

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

384:   PetscObjectTypeCompare((PetscObject)A, type, &sametype);
385:   if (sametype) return 0;

387:   PetscNew(&baijmkl);
388:   B->spptr = (void *)baijmkl;

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

394:   baijmkl->sparse_optimized = PETSC_FALSE;

396:   PetscObjectComposeFunction((PetscObject)B, "MatScale_SeqBAIJMKL_C", MatScale_SeqBAIJMKL);
397:   PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqbaijmkl_seqbaij_C", MatConvert_SeqBAIJMKL_SeqBAIJ);

399:   PetscObjectChangeTypeName((PetscObject)B, MATSEQBAIJMKL);
400:   *newmat = B;
401:   return 0;
402: }

404: static PetscErrorCode MatAssemblyEnd_SeqBAIJMKL(Mat A, MatAssemblyType mode)
405: {
406:   if (mode == MAT_FLUSH_ASSEMBLY) return 0;
407:   MatAssemblyEnd_SeqBAIJ(A, mode);
408:   MatSeqBAIJMKL_create_mkl_handle(A);
409:   A->ops->destroy          = MatDestroy_SeqBAIJMKL;
410:   A->ops->mult             = MatMult_SeqBAIJMKL_SpMV2;
411:   A->ops->multtranspose    = MatMultTranspose_SeqBAIJMKL_SpMV2;
412:   A->ops->multadd          = MatMultAdd_SeqBAIJMKL_SpMV2;
413:   A->ops->multtransposeadd = MatMultTransposeAdd_SeqBAIJMKL_SpMV2;
414:   A->ops->scale            = MatScale_SeqBAIJMKL;
415:   A->ops->diagonalscale    = MatDiagonalScale_SeqBAIJMKL;
416:   A->ops->axpy             = MatAXPY_SeqBAIJMKL;
417:   A->ops->duplicate        = MatDuplicate_SeqBAIJMKL;
418:   return 0;
419: }

421: /*@C
422:    MatCreateSeqBAIJMKL - Creates a sparse matrix of type `MATSEQBAIJMKL`.
423:    This type inherits from `MATSEQBAIJ` and is largely identical, but uses sparse BLAS
424:    routines from Intel MKL whenever possible.
425:    `MatMult()`, `MatMultAdd()`, `MatMultTranspose()`, and `MatMultTransposeAdd()`
426:    operations are currently supported.
427:    If the installed version of MKL supports the "SpMV2" sparse
428:    inspector-executor routines, then those are used by default.
429:    Default PETSc kernels are used otherwise.

431:    Input Parameters:
432: +  comm - MPI communicator, set to `PETSC_COMM_SELF`
433: .  bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row
434:           blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
435: .  m - number of rows
436: .  n - number of columns
437: .  nz - number of nonzero blocks  per block row (same for all rows)
438: -  nnz - array containing the number of nonzero blocks in the various block rows
439:          (possibly different for each block row) or NULL

441:    Output Parameter:
442: .  A - the matrix

444:    It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`,
445:    MatXXXXSetPreallocation() paradigm instead of this routine directly.
446:    [MatXXXXSetPreallocation() is, for example, `MatSeqBAIJSetPreallocation()`]

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

452:    Level: intermediate

454:    Notes:
455:    The number of rows and columns must be divisible by blocksize.

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

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

461:    The `MATSEQBAIJ` format is fully compatible with standard Fortran 77
462:    storage.  That is, the stored row and column indices can begin at
463:    either one (as in Fortran) or zero.  See the users' manual for details.

465:    Specify the preallocated storage with either nz or nnz (not both).
466:    Set nz = `PETSC_DEFAULT` and nnz = NULL for PETSc to control dynamic memory
467:    allocation.  See [Sparse Matrices](sec_matsparse) for details.
468:    matrices.

470: .seealso: [Sparse Matrices](sec_matsparse), `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateBAIJ()`
471: @*/
472: PetscErrorCode MatCreateSeqBAIJMKL(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt nz, const PetscInt nnz[], Mat *A)
473: {
474:   MatCreate(comm, A);
475:   MatSetSizes(*A, m, n, m, n);
476:   MatSetType(*A, MATSEQBAIJMKL);
477:   MatSeqBAIJSetPreallocation_SeqBAIJ(*A, bs, nz, (PetscInt *)nnz);
478:   return 0;
479: }

481: PETSC_EXTERN PetscErrorCode MatCreate_SeqBAIJMKL(Mat A)
482: {
483:   MatSetType(A, MATSEQBAIJ);
484:   MatConvert_SeqBAIJ_SeqBAIJMKL(A, MATSEQBAIJMKL, MAT_INPLACE_MATRIX, &A);
485:   return 0;
486: }