Actual source code: aijmkl.c
1: /*
2: Defines basic operations for the MATSEQAIJMKL matrix class.
3: This class is derived from the MATSEQAIJ class and retains the
4: compressed row storage (aka Yale sparse matrix format) but uses
5: sparse BLAS operations from the Intel Math Kernel Library (MKL)
6: wherever possible.
7: */
9: #include <../src/mat/impls/aij/seq/aij.h>
10: #include <../src/mat/impls/aij/seq/aijmkl/aijmkl.h>
11: #if defined(PETSC_HAVE_MKL_INTEL_ILP64)
12: #define MKL_ILP64
13: #endif
14: #include <mkl_spblas.h>
16: typedef struct {
17: PetscBool no_SpMV2; /* If PETSC_TRUE, then don't use the MKL SpMV2 inspector-executor routines. */
18: PetscBool eager_inspection; /* If PETSC_TRUE, then call mkl_sparse_optimize() in MatDuplicate()/MatAssemblyEnd(). */
19: PetscBool sparse_optimized; /* If PETSC_TRUE, then mkl_sparse_optimize() has been called. */
20: PetscObjectState state;
21: #if defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE)
22: sparse_matrix_t csrA; /* "Handle" used by SpMV2 inspector-executor routines. */
23: struct matrix_descr descr;
24: #endif
25: } Mat_SeqAIJMKL;
27: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJMKL_SeqAIJ(Mat A, MatType type, MatReuse reuse, Mat *newmat)
28: {
29: /* This routine is only called to convert a MATAIJMKL to its base PETSc type, */
30: /* so we will ignore 'MatType type'. */
31: Mat B = *newmat;
32: #if defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE)
33: Mat_SeqAIJMKL *aijmkl = (Mat_SeqAIJMKL *)A->spptr;
34: #endif
36: PetscFunctionBegin;
37: if (reuse == MAT_INITIAL_MATRIX) PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
39: /* Reset the original function pointers. */
40: B->ops->duplicate = MatDuplicate_SeqAIJ;
41: B->ops->assemblyend = MatAssemblyEnd_SeqAIJ;
42: B->ops->destroy = MatDestroy_SeqAIJ;
43: B->ops->mult = MatMult_SeqAIJ;
44: B->ops->multtranspose = MatMultTranspose_SeqAIJ;
45: B->ops->multadd = MatMultAdd_SeqAIJ;
46: B->ops->multtransposeadd = MatMultTransposeAdd_SeqAIJ;
47: B->ops->productsetfromoptions = MatProductSetFromOptions_SeqAIJ;
48: B->ops->matmultsymbolic = MatMatMultSymbolic_SeqAIJ_SeqAIJ;
49: B->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ;
50: B->ops->mattransposemultnumeric = MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ;
51: B->ops->transposematmultnumeric = MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ;
52: B->ops->ptapnumeric = MatPtAPNumeric_SeqAIJ_SeqAIJ;
54: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaijmkl_seqaij_C", NULL));
56: #if defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE)
57: /* Free everything in the Mat_SeqAIJMKL data structure. Currently, this
58: * simply involves destroying the MKL sparse matrix handle and then freeing
59: * the spptr pointer. */
60: if (reuse == MAT_INITIAL_MATRIX) aijmkl = (Mat_SeqAIJMKL *)B->spptr;
62: if (aijmkl->sparse_optimized) PetscCallExternal(mkl_sparse_destroy, aijmkl->csrA);
63: #endif /* PETSC_HAVE_MKL_SPARSE_OPTIMIZE */
64: PetscCall(PetscFree(B->spptr));
66: /* Change the type of B to MATSEQAIJ. */
67: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQAIJ));
69: *newmat = B;
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: static PetscErrorCode MatDestroy_SeqAIJMKL(Mat A)
74: {
75: Mat_SeqAIJMKL *aijmkl = (Mat_SeqAIJMKL *)A->spptr;
77: PetscFunctionBegin;
78: /* If MatHeaderMerge() was used, then this SeqAIJMKL matrix will not have an spptr pointer. */
79: if (aijmkl) {
80: /* Clean up everything in the Mat_SeqAIJMKL data structure, then free A->spptr. */
81: #if defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE)
82: if (aijmkl->sparse_optimized) PetscCallExternal(mkl_sparse_destroy, aijmkl->csrA);
83: #endif /* PETSC_HAVE_MKL_SPARSE_OPTIMIZE */
84: PetscCall(PetscFree(A->spptr));
85: }
87: /* Change the type of A back to SEQAIJ and use MatDestroy_SeqAIJ()
88: * to destroy everything that remains. */
89: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATSEQAIJ));
90: /* I don't call MatSetType(). I believe this is because that
91: * is only to be called when *building* a matrix. I could be wrong, but
92: * that is how things work for the SuperLU matrix class. */
93: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaijmkl_seqaij_C", NULL));
94: PetscCall(MatDestroy_SeqAIJ(A));
95: PetscFunctionReturn(PETSC_SUCCESS);
96: }
98: /* MatSeqAIJMKL_create_mkl_handle(), if called with an AIJMKL matrix that has not had mkl_sparse_optimize() called for it,
99: * creates an MKL sparse matrix handle from the AIJ arrays and calls mkl_sparse_optimize().
100: * If called with an AIJMKL matrix for which aijmkl->sparse_optimized == PETSC_TRUE, then it destroys the old matrix
101: * handle, creates a new one, and then calls mkl_sparse_optimize().
102: * Although in normal MKL usage it is possible to have a valid matrix handle on which mkl_sparse_optimize() has not been
103: * called, for AIJMKL the handle creation and optimization step always occur together, so we don't handle the case of
104: * an unoptimized matrix handle here. */
105: PETSC_INTERN PetscErrorCode MatSeqAIJMKL_create_mkl_handle(Mat A)
106: {
107: #if !defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE)
108: /* If the MKL library does not have mkl_sparse_optimize(), then this routine
109: * does nothing. We make it callable anyway in this case because it cuts
110: * down on littering the code with #ifdefs. */
111: PetscFunctionBegin;
112: PetscFunctionReturn(PETSC_SUCCESS);
113: #else
114: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
115: Mat_SeqAIJMKL *aijmkl = (Mat_SeqAIJMKL *)A->spptr;
116: PetscInt m, n;
117: MatScalar *aa;
118: PetscInt *aj, *ai;
120: PetscFunctionBegin;
121: #if !defined(PETSC_MKL_SPBLAS_DEPRECATED)
122: /* For MKL versions that still support the old, non-inspector-executor interfaces versions, we simply exit here if the no_SpMV2
123: * option has been specified. For versions that have deprecated the old interfaces (version 18, update 2 and later), we must
124: * use the new inspector-executor interfaces, but we can still use the old, non-inspector-executor code by not calling
125: * mkl_sparse_optimize() later. */
126: if (aijmkl->no_SpMV2) PetscFunctionReturn(PETSC_SUCCESS);
127: #endif
129: if (aijmkl->sparse_optimized) {
130: /* Matrix has been previously assembled and optimized. Must destroy old
131: * matrix handle before running the optimization step again. */
132: PetscCallExternal(mkl_sparse_destroy, aijmkl->csrA);
133: }
134: aijmkl->sparse_optimized = PETSC_FALSE;
136: /* Now perform the SpMV2 setup and matrix optimization. */
137: aijmkl->descr.type = SPARSE_MATRIX_TYPE_GENERAL;
138: aijmkl->descr.mode = SPARSE_FILL_MODE_LOWER;
139: aijmkl->descr.diag = SPARSE_DIAG_NON_UNIT;
140: m = A->rmap->n;
141: n = A->cmap->n;
142: aj = a->j; /* aj[k] gives column index for element aa[k]. */
143: aa = a->a; /* Nonzero elements stored row-by-row. */
144: ai = a->i; /* ai[k] is the position in aa and aj where row k starts. */
145: if (a->nz && aa && !A->structure_only) {
146: /* Create a new, optimized sparse matrix handle only if the matrix has nonzero entries.
147: * The MKL sparse-inspector executor routines don't like being passed an empty matrix. */
148: PetscCallExternal(mkl_sparse_x_create_csr, &aijmkl->csrA, SPARSE_INDEX_BASE_ZERO, (MKL_INT)m, (MKL_INT)n, (MKL_INT *)ai, (MKL_INT *)(ai + 1), (MKL_INT *)aj, aa);
149: PetscCallExternal(mkl_sparse_set_mv_hint, aijmkl->csrA, SPARSE_OPERATION_NON_TRANSPOSE, aijmkl->descr, 1000);
150: PetscCallExternal(mkl_sparse_set_memory_hint, aijmkl->csrA, SPARSE_MEMORY_AGGRESSIVE);
151: if (!aijmkl->no_SpMV2) PetscCallExternal(mkl_sparse_optimize, aijmkl->csrA);
152: aijmkl->sparse_optimized = PETSC_TRUE;
153: PetscCall(PetscObjectStateGet((PetscObject)A, &aijmkl->state));
154: } else {
155: aijmkl->csrA = NULL;
156: }
157: PetscFunctionReturn(PETSC_SUCCESS);
158: #endif
159: }
161: #if defined(PETSC_HAVE_MKL_SPARSE_SP2M_FEATURE)
162: /* Take an already created but empty matrix and set up the nonzero structure from an MKL sparse matrix handle. */
163: static PetscErrorCode MatSeqAIJMKL_setup_structure_from_mkl_handle(MPI_Comm comm, sparse_matrix_t csrA, PetscInt nrows, PetscInt ncols, Mat A)
164: {
165: sparse_index_base_t indexing;
166: PetscInt m, n;
167: PetscInt *aj, *ai, *unused;
168: MatScalar *aa;
169: Mat_SeqAIJMKL *aijmkl;
171: PetscFunctionBegin;
172: if (csrA) {
173: /* Note: Must pass in &unused below since MKL can't accept NULL for this output array we don't actually want. */
174: PetscCallExternal(mkl_sparse_x_export_csr, csrA, &indexing, (MKL_INT *)&m, (MKL_INT *)&n, (MKL_INT **)&ai, (MKL_INT **)&unused, (MKL_INT **)&aj, &aa);
175: PetscCheck((m == nrows) && (n == ncols), PETSC_COMM_SELF, PETSC_ERR_LIB, "Number of rows/columns does not match those from mkl_sparse_x_export_csr()");
176: } else {
177: aj = ai = NULL;
178: aa = NULL;
179: }
181: PetscCall(MatSetType(A, MATSEQAIJ));
182: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, nrows, ncols));
183: /* We use MatSeqAIJSetPreallocationCSR() instead of MatCreateSeqAIJWithArrays() because we must copy the arrays exported
184: * from MKL; MKL developers tell us that modifying the arrays may cause unexpected results when using the MKL handle, and
185: * they will be destroyed when the MKL handle is destroyed.
186: * (In the interest of reducing memory consumption in future, can we figure out good ways to deal with this?) */
187: if (csrA) {
188: PetscCall(MatSeqAIJSetPreallocationCSR(A, ai, aj, NULL));
189: } else {
190: /* Since MatSeqAIJSetPreallocationCSR does initial set up and assembly begin/end, we must do that ourselves here. */
191: PetscCall(MatSetUp(A));
192: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
193: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
194: }
196: /* We now have an assembled sequential AIJ matrix created from copies of the exported arrays from the MKL matrix handle.
197: * Now turn it into a MATSEQAIJMKL. */
198: PetscCall(MatConvert_SeqAIJ_SeqAIJMKL(A, MATSEQAIJMKL, MAT_INPLACE_MATRIX, &A));
200: aijmkl = (Mat_SeqAIJMKL *)A->spptr;
201: aijmkl->csrA = csrA;
203: /* The below code duplicates much of what is in MatSeqAIJKL_create_mkl_handle(). I dislike this code duplication, but
204: * MatSeqAIJMKL_create_mkl_handle() cannot be used because we don't need to create a handle -- we've already got one,
205: * and just need to be able to run the MKL optimization step. */
206: aijmkl->descr.type = SPARSE_MATRIX_TYPE_GENERAL;
207: aijmkl->descr.mode = SPARSE_FILL_MODE_LOWER;
208: aijmkl->descr.diag = SPARSE_DIAG_NON_UNIT;
209: if (csrA) {
210: PetscCallExternal(mkl_sparse_set_mv_hint, aijmkl->csrA, SPARSE_OPERATION_NON_TRANSPOSE, aijmkl->descr, 1000);
211: PetscCallExternal(mkl_sparse_set_memory_hint, aijmkl->csrA, SPARSE_MEMORY_AGGRESSIVE);
212: }
213: PetscCall(PetscObjectStateGet((PetscObject)A, &aijmkl->state));
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
216: #endif /* PETSC_HAVE_MKL_SPARSE_SP2M_FEATURE */
218: /* MatSeqAIJMKL_update_from_mkl_handle() updates the matrix values array from the contents of the associated MKL sparse matrix handle.
219: * This is needed after mkl_sparse_sp2m() with SPARSE_STAGE_FINALIZE_MULT has been used to compute new values of the matrix in
220: * MatMatMultNumeric(). */
221: #if defined(PETSC_HAVE_MKL_SPARSE_SP2M_FEATURE)
222: static PetscErrorCode MatSeqAIJMKL_update_from_mkl_handle(Mat A)
223: {
224: PetscInt i;
225: PetscInt nrows, ncols;
226: PetscInt nz;
227: PetscInt *ai, *aj, *unused;
228: PetscScalar *aa;
229: Mat_SeqAIJMKL *aijmkl = (Mat_SeqAIJMKL *)A->spptr;
230: sparse_index_base_t indexing;
232: PetscFunctionBegin;
233: /* Exit immediately in case of the MKL matrix handle being NULL; this will be the case for empty matrices (zero rows or columns). */
234: if (!aijmkl->csrA) PetscFunctionReturn(PETSC_SUCCESS);
236: /* Note: Must pass in &unused below since MKL can't accept NULL for this output array we don't actually want. */
237: PetscCallExternal(mkl_sparse_x_export_csr, aijmkl->csrA, &indexing, (MKL_INT *)&nrows, (MKL_INT *)&ncols, (MKL_INT **)&ai, (MKL_INT **)&unused, (MKL_INT **)&aj, &aa);
239: /* We can't just do a copy from the arrays exported by MKL to those used for the PETSc AIJ storage, because the MKL and PETSc
240: * representations differ in small ways (e.g., more explicit nonzeros per row due to preallocation). */
241: for (i = 0; i < nrows; i++) {
242: nz = ai[i + 1] - ai[i];
243: PetscCall(MatSetValues_SeqAIJ(A, 1, &i, nz, aj + ai[i], aa + ai[i], INSERT_VALUES));
244: }
246: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
247: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
249: PetscCall(PetscObjectStateGet((PetscObject)A, &aijmkl->state));
250: /* At this point our matrix has a valid MKL handle, the contents of which match the PETSc AIJ representation.
251: * The MKL handle has *not* had mkl_sparse_optimize() called on it, though -- the MKL developers have confirmed
252: * that the matrix inspection/optimization step is not performed when matrix-matrix multiplication is finalized. */
253: aijmkl->sparse_optimized = PETSC_FALSE;
254: PetscFunctionReturn(PETSC_SUCCESS);
255: }
256: #endif /* PETSC_HAVE_MKL_SPARSE_SP2M_FEATURE */
258: #if defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE)
259: PETSC_INTERN PetscErrorCode MatSeqAIJMKL_view_mkl_handle(Mat A, PetscViewer viewer)
260: {
261: PetscInt i, j, k;
262: PetscInt nrows, ncols;
263: PetscInt nz;
264: PetscInt *ai, *aj, *unused;
265: PetscScalar *aa;
266: Mat_SeqAIJMKL *aijmkl = (Mat_SeqAIJMKL *)A->spptr;
267: sparse_index_base_t indexing;
269: PetscFunctionBegin;
270: PetscCall(PetscViewerASCIIPrintf(viewer, "Contents of MKL sparse matrix handle for MATSEQAIJMKL object:\n"));
272: /* Exit immediately in case of the MKL matrix handle being NULL; this will be the case for empty matrices (zero rows or columns). */
273: if (!aijmkl->csrA) {
274: PetscCall(PetscViewerASCIIPrintf(viewer, "MKL matrix handle is NULL\n"));
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
278: /* Note: Must pass in &unused below since MKL can't accept NULL for this output array we don't actually want. */
279: PetscCallExternal(mkl_sparse_x_export_csr, aijmkl->csrA, &indexing, (MKL_INT *)&nrows, (MKL_INT *)&ncols, (MKL_INT **)&ai, (MKL_INT **)&unused, (MKL_INT **)&aj, &aa);
281: k = 0;
282: for (i = 0; i < nrows; i++) {
283: PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ": ", i));
284: nz = ai[i + 1] - ai[i];
285: for (j = 0; j < nz; j++) {
286: if (aa) {
287: PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ", %g) ", aj[k], PetscRealPart(aa[k])));
288: } else {
289: PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ", NULL)", aj[k]));
290: }
291: k++;
292: }
293: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
294: }
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
297: #endif /* PETSC_HAVE_MKL_SPARSE_OPTIMIZE */
299: static PetscErrorCode MatDuplicate_SeqAIJMKL(Mat A, MatDuplicateOption op, Mat *M)
300: {
301: Mat_SeqAIJMKL *aijmkl = (Mat_SeqAIJMKL *)A->spptr;
302: Mat_SeqAIJMKL *aijmkl_dest;
304: PetscFunctionBegin;
305: PetscCall(MatDuplicate_SeqAIJ(A, op, M));
306: aijmkl_dest = (Mat_SeqAIJMKL *)(*M)->spptr;
307: PetscCall(PetscArraycpy(aijmkl_dest, aijmkl, 1));
308: aijmkl_dest->sparse_optimized = PETSC_FALSE;
309: if (aijmkl->eager_inspection) PetscCall(MatSeqAIJMKL_create_mkl_handle(A));
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: static PetscErrorCode MatAssemblyEnd_SeqAIJMKL(Mat A, MatAssemblyType mode)
314: {
315: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
316: Mat_SeqAIJMKL *aijmkl;
318: PetscFunctionBegin;
319: if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(PETSC_SUCCESS);
321: /* Since a MATSEQAIJMKL matrix is really just a MATSEQAIJ with some
322: * extra information and some different methods, call the AssemblyEnd
323: * routine for a MATSEQAIJ.
324: * I'm not sure if this is the best way to do this, but it avoids
325: * a lot of code duplication. */
326: a->inode.use = PETSC_FALSE; /* Must disable: otherwise the MKL routines won't get used. */
327: PetscCall(MatAssemblyEnd_SeqAIJ(A, mode));
329: /* If the user has requested "eager" inspection, create the optimized MKL sparse handle (if needed; the function checks).
330: * (The default is to do "lazy" inspection, deferring this until something like MatMult() is called.) */
331: aijmkl = (Mat_SeqAIJMKL *)A->spptr;
332: if (aijmkl->eager_inspection) PetscCall(MatSeqAIJMKL_create_mkl_handle(A));
333: PetscFunctionReturn(PETSC_SUCCESS);
334: }
336: #if !defined(PETSC_MKL_SPBLAS_DEPRECATED)
337: static PetscErrorCode MatMult_SeqAIJMKL(Mat A, Vec xx, Vec yy)
338: {
339: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
340: const PetscScalar *x;
341: PetscScalar *y;
342: const MatScalar *aa;
343: PetscInt m = A->rmap->n;
344: PetscInt n = A->cmap->n;
345: PetscScalar alpha = 1.0;
346: PetscScalar beta = 0.0;
347: const PetscInt *aj, *ai;
348: char matdescra[6];
350: /* Variables not in MatMult_SeqAIJ. */
351: char transa = 'n'; /* Used to indicate to MKL that we are not computing the transpose product. */
353: PetscFunctionBegin;
354: matdescra[0] = 'g'; /* Indicates to MKL that we using a general CSR matrix. */
355: matdescra[3] = 'c'; /* Indicates to MKL that we use C-style (0-based) indexing. */
356: PetscCall(VecGetArrayRead(xx, &x));
357: PetscCall(VecGetArray(yy, &y));
358: aj = a->j; /* aj[k] gives column index for element aa[k]. */
359: aa = a->a; /* Nonzero elements stored row-by-row. */
360: ai = a->i; /* ai[k] is the position in aa and aj where row k starts. */
362: /* Call MKL sparse BLAS routine to do the MatMult. */
363: mkl_xcsrmv(&transa, &m, &n, &alpha, matdescra, aa, aj, ai, ai + 1, x, &beta, y);
365: PetscCall(PetscLogFlops(2.0 * a->nz - a->nonzerorowcnt));
366: PetscCall(VecRestoreArrayRead(xx, &x));
367: PetscCall(VecRestoreArray(yy, &y));
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
370: #endif
372: #if defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE)
373: PetscErrorCode MatMult_SeqAIJMKL_SpMV2(Mat A, Vec xx, Vec yy)
374: {
375: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
376: Mat_SeqAIJMKL *aijmkl = (Mat_SeqAIJMKL *)A->spptr;
377: const PetscScalar *x;
378: PetscScalar *y;
379: PetscObjectState state;
381: PetscFunctionBegin;
382: /* If there are no nonzero entries, zero yy and return immediately. */
383: if (!a->nz) {
384: PetscCall(VecGetArray(yy, &y));
385: PetscCall(PetscArrayzero(y, A->rmap->n));
386: PetscCall(VecRestoreArray(yy, &y));
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }
390: PetscCall(VecGetArrayRead(xx, &x));
391: PetscCall(VecGetArray(yy, &y));
393: /* In some cases, we get to this point without mkl_sparse_optimize() having been called, so we check and then call
394: * it if needed. Eventually, when everything in PETSc is properly updating the matrix state, we should probably
395: * take a "lazy" approach to creation/updating of the MKL matrix handle and plan to always do it here (when needed). */
396: PetscCall(PetscObjectStateGet((PetscObject)A, &state));
397: if (!aijmkl->sparse_optimized || aijmkl->state != state) PetscCall(MatSeqAIJMKL_create_mkl_handle(A));
399: /* Call MKL SpMV2 executor routine to do the MatMult. */
400: PetscCallExternal(mkl_sparse_x_mv, SPARSE_OPERATION_NON_TRANSPOSE, 1.0, aijmkl->csrA, aijmkl->descr, x, 0.0, y);
402: PetscCall(PetscLogFlops(2.0 * a->nz - a->nonzerorowcnt));
403: PetscCall(VecRestoreArrayRead(xx, &x));
404: PetscCall(VecRestoreArray(yy, &y));
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
407: #endif /* PETSC_HAVE_MKL_SPARSE_OPTIMIZE */
409: #if !defined(PETSC_MKL_SPBLAS_DEPRECATED)
410: static PetscErrorCode MatMultTranspose_SeqAIJMKL(Mat A, Vec xx, Vec yy)
411: {
412: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
413: const PetscScalar *x;
414: PetscScalar *y;
415: const MatScalar *aa;
416: PetscInt m = A->rmap->n;
417: PetscInt n = A->cmap->n;
418: PetscScalar alpha = 1.0;
419: PetscScalar beta = 0.0;
420: const PetscInt *aj, *ai;
421: char matdescra[6];
423: /* Variables not in MatMultTranspose_SeqAIJ. */
424: char transa = 't'; /* Used to indicate to MKL that we are computing the transpose product. */
426: PetscFunctionBegin;
427: matdescra[0] = 'g'; /* Indicates to MKL that we using a general CSR matrix. */
428: matdescra[3] = 'c'; /* Indicates to MKL that we use C-style (0-based) indexing. */
429: PetscCall(VecGetArrayRead(xx, &x));
430: PetscCall(VecGetArray(yy, &y));
431: aj = a->j; /* aj[k] gives column index for element aa[k]. */
432: aa = a->a; /* Nonzero elements stored row-by-row. */
433: ai = a->i; /* ai[k] is the position in aa and aj where row k starts. */
435: /* Call MKL sparse BLAS routine to do the MatMult. */
436: mkl_xcsrmv(&transa, &m, &n, &alpha, matdescra, aa, aj, ai, ai + 1, x, &beta, y);
438: PetscCall(PetscLogFlops(2.0 * a->nz - a->nonzerorowcnt));
439: PetscCall(VecRestoreArrayRead(xx, &x));
440: PetscCall(VecRestoreArray(yy, &y));
441: PetscFunctionReturn(PETSC_SUCCESS);
442: }
443: #endif
445: #if defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE)
446: PetscErrorCode MatMultTranspose_SeqAIJMKL_SpMV2(Mat A, Vec xx, Vec yy)
447: {
448: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
449: Mat_SeqAIJMKL *aijmkl = (Mat_SeqAIJMKL *)A->spptr;
450: const PetscScalar *x;
451: PetscScalar *y;
452: PetscObjectState state;
454: PetscFunctionBegin;
455: /* If there are no nonzero entries, zero yy and return immediately. */
456: if (!a->nz) {
457: PetscCall(VecGetArray(yy, &y));
458: PetscCall(PetscArrayzero(y, A->cmap->n));
459: PetscCall(VecRestoreArray(yy, &y));
460: PetscFunctionReturn(PETSC_SUCCESS);
461: }
463: PetscCall(VecGetArrayRead(xx, &x));
464: PetscCall(VecGetArray(yy, &y));
466: /* In some cases, we get to this point without mkl_sparse_optimize() having been called, so we check and then call
467: * it if needed. Eventually, when everything in PETSc is properly updating the matrix state, we should probably
468: * take a "lazy" approach to creation/updating of the MKL matrix handle and plan to always do it here (when needed). */
469: PetscCall(PetscObjectStateGet((PetscObject)A, &state));
470: if (!aijmkl->sparse_optimized || aijmkl->state != state) PetscCall(MatSeqAIJMKL_create_mkl_handle(A));
472: /* Call MKL SpMV2 executor routine to do the MatMultTranspose. */
473: PetscCallExternal(mkl_sparse_x_mv, SPARSE_OPERATION_TRANSPOSE, 1.0, aijmkl->csrA, aijmkl->descr, x, 0.0, y);
475: PetscCall(PetscLogFlops(2.0 * a->nz - a->nonzerorowcnt));
476: PetscCall(VecRestoreArrayRead(xx, &x));
477: PetscCall(VecRestoreArray(yy, &y));
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
480: #endif /* PETSC_HAVE_MKL_SPARSE_OPTIMIZE */
482: #if !defined(PETSC_MKL_SPBLAS_DEPRECATED)
483: static PetscErrorCode MatMultAdd_SeqAIJMKL(Mat A, Vec xx, Vec yy, Vec zz)
484: {
485: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
486: const PetscScalar *x;
487: PetscScalar *y, *z;
488: const MatScalar *aa;
489: PetscInt m = A->rmap->n;
490: PetscInt n = A->cmap->n;
491: const PetscInt *aj, *ai;
492: PetscInt i;
494: /* Variables not in MatMultAdd_SeqAIJ. */
495: char transa = 'n'; /* Used to indicate to MKL that we are not computing the transpose product. */
496: PetscScalar alpha = 1.0;
497: PetscScalar beta;
498: char matdescra[6];
500: PetscFunctionBegin;
501: matdescra[0] = 'g'; /* Indicates to MKL that we using a general CSR matrix. */
502: matdescra[3] = 'c'; /* Indicates to MKL that we use C-style (0-based) indexing. */
504: PetscCall(VecGetArrayRead(xx, &x));
505: PetscCall(VecGetArrayPair(yy, zz, &y, &z));
506: aj = a->j; /* aj[k] gives column index for element aa[k]. */
507: aa = a->a; /* Nonzero elements stored row-by-row. */
508: ai = a->i; /* ai[k] is the position in aa and aj where row k starts. */
510: /* Call MKL sparse BLAS routine to do the MatMult. */
511: if (zz == yy) {
512: /* If zz and yy are the same vector, we can use MKL's mkl_xcsrmv(), which calculates y = alpha*A*x + beta*y. */
513: beta = 1.0;
514: mkl_xcsrmv(&transa, &m, &n, &alpha, matdescra, aa, aj, ai, ai + 1, x, &beta, z);
515: } else {
516: /* zz and yy are different vectors, so call MKL's mkl_xcsrmv() with beta=0, then add the result to z.
517: * MKL sparse BLAS does not have a MatMultAdd equivalent. */
518: beta = 0.0;
519: mkl_xcsrmv(&transa, &m, &n, &alpha, matdescra, aa, aj, ai, ai + 1, x, &beta, z);
520: for (i = 0; i < m; i++) z[i] += y[i];
521: }
523: PetscCall(PetscLogFlops(2.0 * a->nz));
524: PetscCall(VecRestoreArrayRead(xx, &x));
525: PetscCall(VecRestoreArrayPair(yy, zz, &y, &z));
526: PetscFunctionReturn(PETSC_SUCCESS);
527: }
528: #endif
530: #if defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE)
531: PetscErrorCode MatMultAdd_SeqAIJMKL_SpMV2(Mat A, Vec xx, Vec yy, Vec zz)
532: {
533: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
534: Mat_SeqAIJMKL *aijmkl = (Mat_SeqAIJMKL *)A->spptr;
535: const PetscScalar *x;
536: PetscScalar *y, *z;
537: PetscInt m = A->rmap->n;
538: PetscInt i;
540: /* Variables not in MatMultAdd_SeqAIJ. */
541: PetscObjectState state;
543: PetscFunctionBegin;
544: /* If there are no nonzero entries, set zz = yy and return immediately. */
545: if (!a->nz) {
546: PetscCall(VecGetArrayPair(yy, zz, &y, &z));
547: PetscCall(PetscArraycpy(z, y, m));
548: PetscCall(VecRestoreArrayPair(yy, zz, &y, &z));
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: PetscCall(VecGetArrayRead(xx, &x));
553: PetscCall(VecGetArrayPair(yy, zz, &y, &z));
555: /* In some cases, we get to this point without mkl_sparse_optimize() having been called, so we check and then call
556: * it if needed. Eventually, when everything in PETSc is properly updating the matrix state, we should probably
557: * take a "lazy" approach to creation/updating of the MKL matrix handle and plan to always do it here (when needed). */
558: PetscCall(PetscObjectStateGet((PetscObject)A, &state));
559: if (!aijmkl->sparse_optimized || aijmkl->state != state) PetscCall(MatSeqAIJMKL_create_mkl_handle(A));
561: /* Call MKL sparse BLAS routine to do the MatMult. */
562: if (zz == yy) {
563: /* If zz and yy are the same vector, we can use mkl_sparse_x_mv, which calculates y = alpha*A*x + beta*y,
564: * with alpha and beta both set to 1.0. */
565: PetscCallExternal(mkl_sparse_x_mv, SPARSE_OPERATION_NON_TRANSPOSE, 1.0, aijmkl->csrA, aijmkl->descr, x, 1.0, z);
566: } else {
567: /* zz and yy are different vectors, so we call mkl_sparse_x_mv with alpha=1.0 and beta=0.0, and then
568: * we add the contents of vector yy to the result; MKL sparse BLAS does not have a MatMultAdd equivalent. */
569: PetscCallExternal(mkl_sparse_x_mv, SPARSE_OPERATION_NON_TRANSPOSE, 1.0, aijmkl->csrA, aijmkl->descr, x, 0.0, z);
570: for (i = 0; i < m; i++) z[i] += y[i];
571: }
573: PetscCall(PetscLogFlops(2.0 * a->nz));
574: PetscCall(VecRestoreArrayRead(xx, &x));
575: PetscCall(VecRestoreArrayPair(yy, zz, &y, &z));
576: PetscFunctionReturn(PETSC_SUCCESS);
577: }
578: #endif /* PETSC_HAVE_MKL_SPARSE_OPTIMIZE */
580: #if !defined(PETSC_MKL_SPBLAS_DEPRECATED)
581: static PetscErrorCode MatMultTransposeAdd_SeqAIJMKL(Mat A, Vec xx, Vec yy, Vec zz)
582: {
583: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
584: const PetscScalar *x;
585: PetscScalar *y, *z;
586: const MatScalar *aa;
587: PetscInt m = A->rmap->n;
588: PetscInt n = A->cmap->n;
589: const PetscInt *aj, *ai;
590: PetscInt i;
592: /* Variables not in MatMultTransposeAdd_SeqAIJ. */
593: char transa = 't'; /* Used to indicate to MKL that we are computing the transpose product. */
594: PetscScalar alpha = 1.0;
595: PetscScalar beta;
596: char matdescra[6];
598: PetscFunctionBegin;
599: matdescra[0] = 'g'; /* Indicates to MKL that we using a general CSR matrix. */
600: matdescra[3] = 'c'; /* Indicates to MKL that we use C-style (0-based) indexing. */
602: PetscCall(VecGetArrayRead(xx, &x));
603: PetscCall(VecGetArrayPair(yy, zz, &y, &z));
604: aj = a->j; /* aj[k] gives column index for element aa[k]. */
605: aa = a->a; /* Nonzero elements stored row-by-row. */
606: ai = a->i; /* ai[k] is the position in aa and aj where row k starts. */
608: /* Call MKL sparse BLAS routine to do the MatMult. */
609: if (zz == yy) {
610: /* If zz and yy are the same vector, we can use MKL's mkl_xcsrmv(), which calculates y = alpha*A*x + beta*y. */
611: beta = 1.0;
612: mkl_xcsrmv(&transa, &m, &n, &alpha, matdescra, aa, aj, ai, ai + 1, x, &beta, z);
613: } else {
614: /* zz and yy are different vectors, so call MKL's mkl_xcsrmv() with beta=0, then add the result to z.
615: * MKL sparse BLAS does not have a MatMultAdd equivalent. */
616: beta = 0.0;
617: mkl_xcsrmv(&transa, &m, &n, &alpha, matdescra, aa, aj, ai, ai + 1, x, &beta, z);
618: for (i = 0; i < n; i++) z[i] += y[i];
619: }
621: PetscCall(PetscLogFlops(2.0 * a->nz));
622: PetscCall(VecRestoreArrayRead(xx, &x));
623: PetscCall(VecRestoreArrayPair(yy, zz, &y, &z));
624: PetscFunctionReturn(PETSC_SUCCESS);
625: }
626: #endif
628: #if defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE)
629: PetscErrorCode MatMultTransposeAdd_SeqAIJMKL_SpMV2(Mat A, Vec xx, Vec yy, Vec zz)
630: {
631: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
632: Mat_SeqAIJMKL *aijmkl = (Mat_SeqAIJMKL *)A->spptr;
633: const PetscScalar *x;
634: PetscScalar *y, *z;
635: PetscInt n = A->cmap->n;
636: PetscObjectState state;
638: /* Variables not in MatMultTransposeAdd_SeqAIJ. */
640: PetscFunctionBegin;
641: /* If there are no nonzero entries, set zz = yy and return immediately. */
642: if (!a->nz) {
643: PetscCall(VecGetArrayPair(yy, zz, &y, &z));
644: PetscCall(PetscArraycpy(z, y, n));
645: PetscCall(VecRestoreArrayPair(yy, zz, &y, &z));
646: PetscFunctionReturn(PETSC_SUCCESS);
647: }
649: PetscCall(VecGetArrayRead(xx, &x));
650: PetscCall(VecGetArrayPair(yy, zz, &y, &z));
652: /* In some cases, we get to this point without mkl_sparse_optimize() having been called, so we check and then call
653: * it if needed. Eventually, when everything in PETSc is properly updating the matrix state, we should probably
654: * take a "lazy" approach to creation/updating of the MKL matrix handle and plan to always do it here (when needed). */
655: PetscCall(PetscObjectStateGet((PetscObject)A, &state));
656: if (!aijmkl->sparse_optimized || aijmkl->state != state) PetscCall(MatSeqAIJMKL_create_mkl_handle(A));
658: /* Call MKL sparse BLAS routine to do the MatMult. */
659: if (zz == yy) {
660: /* If zz and yy are the same vector, we can use mkl_sparse_x_mv, which calculates y = alpha*A*x + beta*y,
661: * with alpha and beta both set to 1.0. */
662: PetscCallExternal(mkl_sparse_x_mv, SPARSE_OPERATION_TRANSPOSE, 1.0, aijmkl->csrA, aijmkl->descr, x, 1.0, z);
663: } else {
664: /* zz and yy are different vectors, so we call mkl_sparse_x_mv with alpha=1.0 and beta=0.0, and then
665: * we add the contents of vector yy to the result; MKL sparse BLAS does not have a MatMultAdd equivalent. */
666: PetscCallExternal(mkl_sparse_x_mv, SPARSE_OPERATION_TRANSPOSE, 1.0, aijmkl->csrA, aijmkl->descr, x, 0.0, z);
667: for (PetscInt i = 0; i < n; i++) z[i] += y[i];
668: }
670: PetscCall(PetscLogFlops(2.0 * a->nz));
671: PetscCall(VecRestoreArrayRead(xx, &x));
672: PetscCall(VecRestoreArrayPair(yy, zz, &y, &z));
673: PetscFunctionReturn(PETSC_SUCCESS);
674: }
675: #endif /* PETSC_HAVE_MKL_SPARSE_OPTIMIZE */
677: #if defined(PETSC_HAVE_MKL_SPARSE_SP2M_FEATURE)
678: static PetscErrorCode MatMatMultSymbolic_SeqAIJMKL_SeqAIJMKL_Private(Mat A, const sparse_operation_t transA, Mat B, const sparse_operation_t transB, Mat C)
679: {
680: Mat_SeqAIJMKL *a = (Mat_SeqAIJMKL *)A->spptr, *b = (Mat_SeqAIJMKL *)B->spptr;
681: sparse_matrix_t csrA, csrB, csrC;
682: PetscInt nrows, ncols;
683: struct matrix_descr descr_type_gen;
684: PetscObjectState state;
686: PetscFunctionBegin;
687: /* Determine the number of rows and columns that the result matrix C will have. We have to do this ourselves because MKL does
688: * not handle sparse matrices with zero rows or columns. */
689: if (transA == SPARSE_OPERATION_NON_TRANSPOSE) nrows = A->rmap->N;
690: else nrows = A->cmap->N;
691: if (transB == SPARSE_OPERATION_NON_TRANSPOSE) ncols = B->cmap->N;
692: else ncols = B->rmap->N;
694: PetscCall(PetscObjectStateGet((PetscObject)A, &state));
695: if (!a->sparse_optimized || a->state != state) PetscCall(MatSeqAIJMKL_create_mkl_handle(A));
696: PetscCall(PetscObjectStateGet((PetscObject)B, &state));
697: if (!b->sparse_optimized || b->state != state) PetscCall(MatSeqAIJMKL_create_mkl_handle(B));
698: csrA = a->csrA;
699: csrB = b->csrA;
700: descr_type_gen.type = SPARSE_MATRIX_TYPE_GENERAL;
702: if (csrA && csrB) {
703: PetscCallExternal(mkl_sparse_sp2m, transA, descr_type_gen, csrA, transB, descr_type_gen, csrB, SPARSE_STAGE_FULL_MULT_NO_VAL, &csrC);
704: } else {
705: csrC = NULL;
706: }
708: PetscCall(MatSeqAIJMKL_setup_structure_from_mkl_handle(PETSC_COMM_SELF, csrC, nrows, ncols, C));
709: PetscFunctionReturn(PETSC_SUCCESS);
710: }
712: PetscErrorCode MatMatMultNumeric_SeqAIJMKL_SeqAIJMKL_Private(Mat A, const sparse_operation_t transA, Mat B, const sparse_operation_t transB, Mat C)
713: {
714: Mat_SeqAIJMKL *a = (Mat_SeqAIJMKL *)A->spptr, *b = (Mat_SeqAIJMKL *)B->spptr, *c = (Mat_SeqAIJMKL *)C->spptr;
715: sparse_matrix_t csrA, csrB, csrC;
716: struct matrix_descr descr_type_gen;
717: PetscObjectState state;
719: PetscFunctionBegin;
720: PetscCall(PetscObjectStateGet((PetscObject)A, &state));
721: if (!a->sparse_optimized || a->state != state) PetscCall(MatSeqAIJMKL_create_mkl_handle(A));
722: PetscCall(PetscObjectStateGet((PetscObject)B, &state));
723: if (!b->sparse_optimized || b->state != state) PetscCall(MatSeqAIJMKL_create_mkl_handle(B));
724: csrA = a->csrA;
725: csrB = b->csrA;
726: csrC = c->csrA;
727: descr_type_gen.type = SPARSE_MATRIX_TYPE_GENERAL;
729: if (csrA && csrB) {
730: PetscCallExternal(mkl_sparse_sp2m, transA, descr_type_gen, csrA, transB, descr_type_gen, csrB, SPARSE_STAGE_FINALIZE_MULT, &csrC);
731: } else {
732: csrC = NULL;
733: }
735: /* Have to update the PETSc AIJ representation for matrix C from contents of MKL handle. */
736: PetscCall(MatSeqAIJMKL_update_from_mkl_handle(C));
737: PetscFunctionReturn(PETSC_SUCCESS);
738: }
740: PetscErrorCode MatMatMultSymbolic_SeqAIJMKL_SeqAIJMKL(Mat A, Mat B, PetscReal fill, Mat C)
741: {
742: PetscFunctionBegin;
743: PetscCall(MatMatMultSymbolic_SeqAIJMKL_SeqAIJMKL_Private(A, SPARSE_OPERATION_NON_TRANSPOSE, B, SPARSE_OPERATION_NON_TRANSPOSE, C));
744: PetscFunctionReturn(PETSC_SUCCESS);
745: }
747: PetscErrorCode MatMatMultNumeric_SeqAIJMKL_SeqAIJMKL(Mat A, Mat B, Mat C)
748: {
749: PetscFunctionBegin;
750: PetscCall(MatMatMultNumeric_SeqAIJMKL_SeqAIJMKL_Private(A, SPARSE_OPERATION_NON_TRANSPOSE, B, SPARSE_OPERATION_NON_TRANSPOSE, C));
751: PetscFunctionReturn(PETSC_SUCCESS);
752: }
754: PetscErrorCode MatTransposeMatMultNumeric_SeqAIJMKL_SeqAIJMKL(Mat A, Mat B, Mat C)
755: {
756: PetscFunctionBegin;
757: PetscCall(MatMatMultNumeric_SeqAIJMKL_SeqAIJMKL_Private(A, SPARSE_OPERATION_TRANSPOSE, B, SPARSE_OPERATION_NON_TRANSPOSE, C));
758: PetscFunctionReturn(PETSC_SUCCESS);
759: }
761: PetscErrorCode MatTransposeMatMultSymbolic_SeqAIJMKL_SeqAIJMKL(Mat A, Mat B, PetscReal fill, Mat C)
762: {
763: PetscFunctionBegin;
764: PetscCall(MatMatMultSymbolic_SeqAIJMKL_SeqAIJMKL_Private(A, SPARSE_OPERATION_TRANSPOSE, B, SPARSE_OPERATION_NON_TRANSPOSE, C));
765: PetscFunctionReturn(PETSC_SUCCESS);
766: }
768: PetscErrorCode MatMatTransposeMultSymbolic_SeqAIJMKL_SeqAIJMKL(Mat A, Mat B, PetscReal fill, Mat C)
769: {
770: PetscFunctionBegin;
771: PetscCall(MatMatMultSymbolic_SeqAIJMKL_SeqAIJMKL_Private(A, SPARSE_OPERATION_NON_TRANSPOSE, B, SPARSE_OPERATION_TRANSPOSE, C));
772: PetscFunctionReturn(PETSC_SUCCESS);
773: }
775: PetscErrorCode MatMatTransposeMultNumeric_SeqAIJMKL_SeqAIJMKL(Mat A, Mat B, Mat C)
776: {
777: PetscFunctionBegin;
778: PetscCall(MatMatMultNumeric_SeqAIJMKL_SeqAIJMKL_Private(A, SPARSE_OPERATION_NON_TRANSPOSE, B, SPARSE_OPERATION_TRANSPOSE, C));
779: PetscFunctionReturn(PETSC_SUCCESS);
780: }
782: static PetscErrorCode MatProductNumeric_AtB_SeqAIJMKL_SeqAIJMKL(Mat C)
783: {
784: Mat_Product *product = C->product;
785: Mat A = product->A, B = product->B;
787: PetscFunctionBegin;
788: PetscCall(MatTransposeMatMultNumeric_SeqAIJMKL_SeqAIJMKL(A, B, C));
789: PetscFunctionReturn(PETSC_SUCCESS);
790: }
792: static PetscErrorCode MatProductSymbolic_AtB_SeqAIJMKL_SeqAIJMKL(Mat C)
793: {
794: Mat_Product *product = C->product;
795: Mat A = product->A, B = product->B;
796: PetscReal fill = product->fill;
798: PetscFunctionBegin;
799: PetscCall(MatTransposeMatMultSymbolic_SeqAIJMKL_SeqAIJMKL(A, B, fill, C));
800: C->ops->productnumeric = MatProductNumeric_AtB_SeqAIJMKL_SeqAIJMKL;
801: PetscFunctionReturn(PETSC_SUCCESS);
802: }
804: PetscErrorCode MatPtAPNumeric_SeqAIJMKL_SeqAIJMKL_SymmetricReal(Mat A, Mat P, Mat C)
805: {
806: Mat Ct;
807: Vec zeros;
808: Mat_SeqAIJMKL *a = (Mat_SeqAIJMKL *)A->spptr, *p = (Mat_SeqAIJMKL *)P->spptr, *c = (Mat_SeqAIJMKL *)C->spptr;
809: sparse_matrix_t csrA, csrP, csrC;
810: PetscBool set, flag;
811: struct matrix_descr descr_type_sym;
812: PetscObjectState state;
814: PetscFunctionBegin;
815: PetscCall(MatIsSymmetricKnown(A, &set, &flag));
816: PetscCheck(set && flag, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatPtAPNumeric_SeqAIJMKL_SeqAIJMKL_SymmetricReal() called on matrix A not marked as symmetric");
818: PetscCall(PetscObjectStateGet((PetscObject)A, &state));
819: if (!a->sparse_optimized || a->state != state) PetscCall(MatSeqAIJMKL_create_mkl_handle(A));
820: PetscCall(PetscObjectStateGet((PetscObject)P, &state));
821: if (!p->sparse_optimized || p->state != state) PetscCall(MatSeqAIJMKL_create_mkl_handle(P));
822: csrA = a->csrA;
823: csrP = p->csrA;
824: csrC = c->csrA;
825: descr_type_sym.type = SPARSE_MATRIX_TYPE_SYMMETRIC;
826: descr_type_sym.mode = SPARSE_FILL_MODE_UPPER;
827: descr_type_sym.diag = SPARSE_DIAG_NON_UNIT;
829: /* the call below won't work for complex matrices. (We protect this when pointers are assigned in MatConvert.) */
830: PetscCallExternal(mkl_sparse_sypr, SPARSE_OPERATION_TRANSPOSE, csrP, csrA, descr_type_sym, &csrC, SPARSE_STAGE_FINALIZE_MULT);
832: /* Update the PETSc AIJ representation for matrix C from contents of MKL handle.
833: * This is more complicated than it should be: it turns out that, though mkl_sparse_sypr() will accept a full AIJ/CSR matrix,
834: * the output matrix only contains the upper or lower triangle (we arbitrarily have chosen upper) of the symmetric matrix.
835: * We have to fill in the missing portion, which we currently do below by forming the transpose and performing at MatAXPY
836: * operation. This may kill any performance benefit of using the optimized mkl_sparse_sypr() routine. Performance might
837: * improve if we come up with a more efficient way to do this, or we can convince the MKL team to provide an option to output
838: * the full matrix. */
839: PetscCall(MatSeqAIJMKL_update_from_mkl_handle(C));
840: PetscCall(MatTranspose(C, MAT_INITIAL_MATRIX, &Ct));
841: PetscCall(MatCreateVecs(C, &zeros, NULL));
842: PetscCall(VecSetFromOptions(zeros));
843: PetscCall(MatDiagonalSet(Ct, zeros, INSERT_VALUES));
844: PetscCall(MatAXPY(C, 1.0, Ct, DIFFERENT_NONZERO_PATTERN));
845: /* Note: The MatAXPY() call destroys the MatProduct, so we must recreate it. */
846: PetscCall(MatProductCreateWithMat(A, P, NULL, C));
847: PetscCall(MatProductSetType(C, MATPRODUCT_PtAP));
848: PetscCall(MatSeqAIJMKL_create_mkl_handle(C));
849: PetscCall(VecDestroy(&zeros));
850: PetscCall(MatDestroy(&Ct));
851: PetscFunctionReturn(PETSC_SUCCESS);
852: }
854: PetscErrorCode MatProductSymbolic_PtAP_SeqAIJMKL_SeqAIJMKL_SymmetricReal(Mat C)
855: {
856: Mat_Product *product = C->product;
857: Mat A = product->A, P = product->B;
858: Mat_SeqAIJMKL *a = (Mat_SeqAIJMKL *)A->spptr, *p = (Mat_SeqAIJMKL *)P->spptr;
859: sparse_matrix_t csrA, csrP, csrC;
860: struct matrix_descr descr_type_sym;
861: PetscObjectState state;
863: PetscFunctionBegin;
864: PetscCall(PetscObjectStateGet((PetscObject)A, &state));
865: if (!a->sparse_optimized || a->state != state) PetscCall(MatSeqAIJMKL_create_mkl_handle(A));
866: PetscCall(PetscObjectStateGet((PetscObject)P, &state));
867: if (!p->sparse_optimized || p->state != state) PetscCall(MatSeqAIJMKL_create_mkl_handle(P));
868: csrA = a->csrA;
869: csrP = p->csrA;
870: descr_type_sym.type = SPARSE_MATRIX_TYPE_SYMMETRIC;
871: descr_type_sym.mode = SPARSE_FILL_MODE_UPPER;
872: descr_type_sym.diag = SPARSE_DIAG_NON_UNIT;
874: /* the call below won't work for complex matrices. (We protect this when pointers are assigned in MatConvert.) */
875: if (csrP && csrA) {
876: PetscCallExternal(mkl_sparse_sypr, SPARSE_OPERATION_TRANSPOSE, csrP, csrA, descr_type_sym, &csrC, SPARSE_STAGE_FULL_MULT_NO_VAL);
877: } else {
878: csrC = NULL;
879: }
881: /* Update the I and J arrays of the PETSc AIJ representation for matrix C from contents of MKL handle.
882: * Note that, because mkl_sparse_sypr() only computes one triangle of the symmetric matrix, this representation will only contain
883: * the upper triangle of the symmetric matrix. We fix this in MatPtAPNumeric_SeqAIJMKL_SeqAIJMKL_SymmetricReal(). I believe that
884: * leaving things in this incomplete state is OK because the numeric product should follow soon after, but am not certain if this
885: * is guaranteed. */
886: PetscCall(MatSeqAIJMKL_setup_structure_from_mkl_handle(PETSC_COMM_SELF, csrC, P->cmap->N, P->cmap->N, C));
888: C->ops->productnumeric = MatProductNumeric_PtAP;
889: PetscFunctionReturn(PETSC_SUCCESS);
890: }
892: static PetscErrorCode MatProductSetFromOptions_SeqAIJMKL_AB(Mat C)
893: {
894: PetscFunctionBegin;
895: C->ops->productsymbolic = MatProductSymbolic_AB;
896: C->ops->matmultsymbolic = MatMatMultSymbolic_SeqAIJMKL_SeqAIJMKL;
897: PetscFunctionReturn(PETSC_SUCCESS);
898: }
900: static PetscErrorCode MatProductSetFromOptions_SeqAIJMKL_AtB(Mat C)
901: {
902: PetscFunctionBegin;
903: C->ops->productsymbolic = MatProductSymbolic_AtB_SeqAIJMKL_SeqAIJMKL;
904: PetscFunctionReturn(PETSC_SUCCESS);
905: }
907: static PetscErrorCode MatProductSetFromOptions_SeqAIJMKL_ABt(Mat C)
908: {
909: PetscFunctionBegin;
910: C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ;
911: C->ops->productsymbolic = MatProductSymbolic_ABt;
912: PetscFunctionReturn(PETSC_SUCCESS);
913: }
915: static PetscErrorCode MatProductSetFromOptions_SeqAIJMKL_PtAP(Mat C)
916: {
917: Mat_Product *product = C->product;
918: Mat A = product->A;
919: PetscBool set, flag;
921: PetscFunctionBegin;
922: /* AIJMKL only has an optimized routine for PtAP when A is symmetric and real. */
923: PetscCall(MatIsSymmetricKnown(A, &set, &flag));
924: PetscCheck(!PetscDefined(USE_COMPLEX) && set && flag, PETSC_COMM_SELF, PETSC_ERR_SUP, "MATPRODUCT_PtAP not supported for type SeqAIJMKL");
925: /* we don't set C->ops->productnumeric here, as this must happen in MatProductSymbolic_PtAP_SeqAIJMKL_SeqAIJMKL_SymmetricReal() */
926: C->ops->productsymbolic = MatProductSymbolic_PtAP_SeqAIJMKL_SeqAIJMKL_SymmetricReal;
927: PetscFunctionReturn(PETSC_SUCCESS);
928: }
930: static PetscErrorCode MatProductSetFromOptions_SeqAIJMKL_RARt(Mat C)
931: {
932: PetscFunctionBegin;
933: C->ops->productsymbolic = NULL; /* MatProductSymbolic_Unsafe() will be used. */
934: PetscFunctionReturn(PETSC_SUCCESS);
935: }
937: static PetscErrorCode MatProductSetFromOptions_SeqAIJMKL_ABC(Mat C)
938: {
939: PetscFunctionBegin;
940: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MATPRODUCT_ABC not supported for type SeqAIJMKL");
941: }
943: PetscErrorCode MatProductSetFromOptions_SeqAIJMKL(Mat C)
944: {
945: Mat_Product *product = C->product;
947: PetscFunctionBegin;
948: switch (product->type) {
949: case MATPRODUCT_AB:
950: PetscCall(MatProductSetFromOptions_SeqAIJMKL_AB(C));
951: break;
952: case MATPRODUCT_AtB:
953: PetscCall(MatProductSetFromOptions_SeqAIJMKL_AtB(C));
954: break;
955: case MATPRODUCT_ABt:
956: PetscCall(MatProductSetFromOptions_SeqAIJMKL_ABt(C));
957: break;
958: case MATPRODUCT_PtAP:
959: PetscCall(MatProductSetFromOptions_SeqAIJMKL_PtAP(C));
960: break;
961: case MATPRODUCT_RARt:
962: PetscCall(MatProductSetFromOptions_SeqAIJMKL_RARt(C));
963: break;
964: case MATPRODUCT_ABC:
965: PetscCall(MatProductSetFromOptions_SeqAIJMKL_ABC(C));
966: break;
967: default:
968: break;
969: }
970: PetscFunctionReturn(PETSC_SUCCESS);
971: }
972: #endif /* PETSC_HAVE_MKL_SPARSE_SP2M_FEATURE */
974: /* MatConvert_SeqAIJ_SeqAIJMKL converts a SeqAIJ matrix into a
975: * SeqAIJMKL matrix. This routine is called by the MatCreate_SeqAIJMKL()
976: * routine, but can also be used to convert an assembled SeqAIJ matrix
977: * into a SeqAIJMKL one. */
978: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJMKL(Mat A, MatType type, MatReuse reuse, Mat *newmat)
979: {
980: Mat B = *newmat;
981: Mat_SeqAIJMKL *aijmkl;
982: PetscBool set;
983: PetscBool sametype;
985: PetscFunctionBegin;
986: if (reuse == MAT_INITIAL_MATRIX) PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
988: PetscCall(PetscObjectTypeCompare((PetscObject)A, type, &sametype));
989: if (sametype) PetscFunctionReturn(PETSC_SUCCESS);
991: PetscCall(PetscNew(&aijmkl));
992: B->spptr = (void *)aijmkl;
994: /* Set function pointers for methods that we inherit from AIJ but override.
995: * We also parse some command line options below, since those determine some of the methods we point to. */
996: B->ops->duplicate = MatDuplicate_SeqAIJMKL;
997: B->ops->assemblyend = MatAssemblyEnd_SeqAIJMKL;
998: B->ops->destroy = MatDestroy_SeqAIJMKL;
1000: aijmkl->sparse_optimized = PETSC_FALSE;
1001: aijmkl->no_SpMV2 = PetscDefined(HAVE_MKL_SPARSE_OPTIMIZE) ? PETSC_FALSE : PETSC_TRUE; /* Default to using the SpMV2 routines if our MKL supports them. */
1002: aijmkl->eager_inspection = PETSC_FALSE;
1004: /* Parse command line options. */
1005: PetscOptionsBegin(PetscObjectComm((PetscObject)A), ((PetscObject)A)->prefix, "AIJMKL Options", "Mat");
1006: PetscCall(PetscOptionsBool("-mat_aijmkl_no_spmv2", "Disable use of inspector-executor (SpMV 2) routines", "None", (PetscBool)aijmkl->no_SpMV2, (PetscBool *)&aijmkl->no_SpMV2, &set));
1007: PetscCall(PetscOptionsBool("-mat_aijmkl_eager_inspection", "Run inspection at matrix assembly time, instead of waiting until needed by an operation", "None", (PetscBool)aijmkl->eager_inspection, (PetscBool *)&aijmkl->eager_inspection, &set));
1008: PetscOptionsEnd();
1009: #if !defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE)
1010: if (!aijmkl->no_SpMV2) {
1011: PetscCall(PetscInfo(B, "User requested use of MKL SpMV2 routines, but MKL version does not support mkl_sparse_optimize(); defaulting to non-SpMV2 routines.\n"));
1012: aijmkl->no_SpMV2 = PETSC_TRUE;
1013: }
1014: #endif
1016: #if defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE)
1017: B->ops->mult = MatMult_SeqAIJMKL_SpMV2;
1018: B->ops->multtranspose = MatMultTranspose_SeqAIJMKL_SpMV2;
1019: B->ops->multadd = MatMultAdd_SeqAIJMKL_SpMV2;
1020: B->ops->multtransposeadd = MatMultTransposeAdd_SeqAIJMKL_SpMV2;
1021: #if defined(PETSC_HAVE_MKL_SPARSE_SP2M_FEATURE)
1022: B->ops->productsetfromoptions = MatProductSetFromOptions_SeqAIJMKL;
1023: B->ops->matmultsymbolic = MatMatMultSymbolic_SeqAIJMKL_SeqAIJMKL;
1024: B->ops->matmultnumeric = MatMatMultNumeric_SeqAIJMKL_SeqAIJMKL;
1025: B->ops->mattransposemultnumeric = MatMatTransposeMultNumeric_SeqAIJMKL_SeqAIJMKL;
1026: B->ops->transposematmultnumeric = MatTransposeMatMultNumeric_SeqAIJMKL_SeqAIJMKL;
1027: #if !defined(PETSC_USE_COMPLEX)
1028: B->ops->ptapnumeric = MatPtAPNumeric_SeqAIJMKL_SeqAIJMKL_SymmetricReal;
1029: #else
1030: B->ops->ptapnumeric = NULL;
1031: #endif
1032: #endif
1033: #endif /* PETSC_HAVE_MKL_SPARSE_OPTIMIZE */
1035: #if !defined(PETSC_MKL_SPBLAS_DEPRECATED)
1036: /* In MKL version 18, update 2, the old sparse BLAS interfaces were marked as deprecated. If "no_SpMV2" has been specified by the
1037: * user and the old SpBLAS interfaces are deprecated in our MKL version, we use the new _SpMV2 routines (set above), but do not
1038: * call mkl_sparse_optimize(), which results in the old numerical kernels (without the inspector-executor model) being used. For
1039: * versions in which the older interface has not been deprecated, we use the old interface. */
1040: if (aijmkl->no_SpMV2) {
1041: B->ops->mult = MatMult_SeqAIJMKL;
1042: B->ops->multtranspose = MatMultTranspose_SeqAIJMKL;
1043: B->ops->multadd = MatMultAdd_SeqAIJMKL;
1044: B->ops->multtransposeadd = MatMultTransposeAdd_SeqAIJMKL;
1045: }
1046: #endif
1048: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaijmkl_seqaij_C", MatConvert_SeqAIJMKL_SeqAIJ));
1050: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQAIJMKL));
1051: *newmat = B;
1052: PetscFunctionReturn(PETSC_SUCCESS);
1053: }
1055: /*@C
1056: MatCreateSeqAIJMKL - Creates a sparse matrix of type `MATSEQAIJMKL`.
1058: Collective
1060: Input Parameters:
1061: + comm - MPI communicator, set to `PETSC_COMM_SELF`
1062: . m - number of rows
1063: . n - number of columns
1064: . nz - number of nonzeros per row (same for all rows)
1065: - nnz - array containing the number of nonzeros in the various rows
1066: (possibly different for each row) or `NULL`
1068: Output Parameter:
1069: . A - the matrix
1071: Options Database Keys:
1072: + -mat_aijmkl_no_spmv2 - disable use of the SpMV2 inspector-executor routines
1073: - -mat_aijmkl_eager_inspection - perform MKL "inspection" phase upon matrix assembly; default is to do "lazy" inspection,
1074: performing this step the first time the matrix is applied
1076: Level: intermediate
1078: Notes:
1079: If `nnz` is given then `nz` is ignored
1081: This type inherits from `MATSEQAIJ` and is largely identical, but uses sparse BLAS
1082: routines from Intel MKL whenever possible.
1084: If the installed version of MKL supports the "SpMV2" sparse
1085: inspector-executor routines, then those are used by default.
1087: `MatMult()`, `MatMultAdd()`, `MatMultTranspose()`, `MatMultTransposeAdd()`, `MatMatMult()`, `MatTransposeMatMult()`, and `MatPtAP()`
1088: (for symmetric A) operations are currently supported.
1089: MKL version 18, update 2 or later is required for `MatPtAP()`, `MatPtAPNumeric()` and `MatMatMultNumeric()`.
1091: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateMPIAIJMKL()`, `MatSetValues()`
1092: @*/
1093: PetscErrorCode MatCreateSeqAIJMKL(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt nz, const PetscInt nnz[], Mat *A)
1094: {
1095: PetscFunctionBegin;
1096: PetscCall(MatCreate(comm, A));
1097: PetscCall(MatSetSizes(*A, m, n, m, n));
1098: PetscCall(MatSetType(*A, MATSEQAIJMKL));
1099: PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*A, nz, nnz));
1100: PetscFunctionReturn(PETSC_SUCCESS);
1101: }
1103: PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJMKL(Mat A)
1104: {
1105: PetscFunctionBegin;
1106: PetscCall(MatSetType(A, MATSEQAIJ));
1107: PetscCall(MatConvert_SeqAIJ_SeqAIJMKL(A, MATSEQAIJMKL, MAT_INPLACE_MATRIX, &A));
1108: PetscFunctionReturn(PETSC_SUCCESS);
1109: }