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: }