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