1: #include <../src/mat/impls/aij/mpi/mpiaij.h>
  2: /*@C
  3:   MatCreateMPIAIJMKL - Creates a sparse parallel matrix whose local
  4:   portions are stored as `MATSEQAIJMKL` matrices (a matrix class that inherits
  5:   from `MATSEQAIJ` but uses some operations provided by Intel MKL).

  7:   Collective

  9:   Input Parameters:
 10: + comm  - MPI communicator
 11: . m     - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
 12:            This value should be the same as the local size used in creating the
 13:            y vector for the matrix-vector product y = Ax.
 14: . n     - This value should be the same as the local size used in creating the
 15:        x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have
 16:        calculated if N is given) For square matrices n is almost always `m`.
 17: . M     - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given)
 18: . N     - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given)
 19: . d_nz  - number of nonzeros per row in DIAGONAL portion of local submatrix
 20:            (same value is used for all local rows)
 21: . d_nnz - array containing the number of nonzeros in the various rows of the
 22:            DIAGONAL portion of the local submatrix (possibly different for each row)
 23:            or `NULL`, if `d_nz` is used to specify the nonzero structure.
 24:            The size of this array is equal to the number of local rows, i.e `m`.
 25:            For matrices you plan to factor you must leave room for the diagonal entry and
 26:            put in the entry even if it is zero.
 27: . o_nz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
 28:            submatrix (same value is used for all local rows).
 29: - o_nnz - array containing the number of nonzeros in the various rows of the
 30:            OFF-DIAGONAL portion of the local submatrix (possibly different for
 31:            each row) or `NULL`, if `o_nz` is used to specify the nonzero
 32:            structure. The size of this array is equal to the number
 33:            of local rows, i.e `m`.

 35:   Output Parameter:
 36: . A - the matrix

 38:   Options Database Key:
 39: . -mat_aijmkl_no_spmv2 - disables use of the SpMV2 inspector-executor routines

 41:   Level: intermediate

 43:   Notes:
 44:   If the *_nnz parameter is given then the *_nz parameter is ignored

 46:   `m`,`n`,`M`,`N` parameters specify the size of the matrix, and its partitioning across
 47:   processors, while `d_nz`,`d_nnz`,`o_nz`,`o_nnz` parameters specify the approximate
 48:   storage requirements for this matrix.

 50:   If `PETSC_DECIDE` or `PETSC_DETERMINE` is used for a particular argument on one
 51:   processor than it must be used on all processors that share the object for
 52:   that argument.

 54:   The user MUST specify either the local or global matrix dimensions
 55:   (possibly both).

 57:   The parallel matrix is partitioned such that the first m0 rows belong to
 58:   process 0, the next m1 rows belong to process 1, the next m2 rows belong
 59:   to process 2 etc.. where m0,m1,m2... are the input parameter `m`.

 61:   The DIAGONAL portion of the local submatrix of a processor can be defined
 62:   as the submatrix which is obtained by extraction the part corresponding
 63:   to the rows r1-r2 and columns r1-r2 of the global matrix, where r1 is the
 64:   first row that belongs to the processor, and r2 is the last row belonging
 65:   to the this processor. This is a square mxm matrix. The remaining portion
 66:   of the local submatrix (mxN) constitute the OFF-DIAGONAL portion.

 68:   If `o_nnz`, `d_nnz` are specified, then `o_nz`, and `d_nz` are ignored.

 70:   When calling this routine with a single process communicator, a matrix of
 71:   type `MATSEQAIJMKL` is returned.  If a matrix of type `MATMPIAIJMKL` is desired
 72:   for this type of communicator, use the construction mechanism
 73: .vb
 74:   MatCreate(...,&A);
 75:   MatSetType(A,MPIAIJMKL);
 76:   MatMPIAIJSetPreallocation(A,...);
 77: .ve

 79: .seealso: [](ch_matrices), `Mat`, [Sparse Matrix Creation](sec_matsparse), `MATMPIAIJMKL`, `MatCreate()`, `MatCreateSeqAIJMKL()`, `MatSetValues()`
 80: @*/
 81: PetscErrorCode MatCreateMPIAIJMKL(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[], Mat *A)
 82: {
 83:   PetscMPIInt size;

 85:   PetscFunctionBegin;
 86:   PetscCall(MatCreate(comm, A));
 87:   PetscCall(MatSetSizes(*A, m, n, M, N));
 88:   PetscCallMPI(MPI_Comm_size(comm, &size));
 89:   if (size > 1) {
 90:     PetscCall(MatSetType(*A, MATMPIAIJMKL));
 91:     PetscCall(MatMPIAIJSetPreallocation(*A, d_nz, d_nnz, o_nz, o_nnz));
 92:   } else {
 93:     PetscCall(MatSetType(*A, MATSEQAIJMKL));
 94:     PetscCall(MatSeqAIJSetPreallocation(*A, d_nz, d_nnz));
 95:   }
 96:   PetscFunctionReturn(PETSC_SUCCESS);
 97: }

 99: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJMKL(Mat, MatType, MatReuse, Mat *);

101: static PetscErrorCode MatMPIAIJSetPreallocation_MPIAIJMKL(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
102: {
103:   Mat_MPIAIJ *b = (Mat_MPIAIJ *)B->data;

105:   PetscFunctionBegin;
106:   PetscCall(MatMPIAIJSetPreallocation_MPIAIJ(B, d_nz, d_nnz, o_nz, o_nnz));
107:   PetscCall(MatConvert_SeqAIJ_SeqAIJMKL(b->A, MATSEQAIJMKL, MAT_INPLACE_MATRIX, &b->A));
108:   PetscCall(MatConvert_SeqAIJ_SeqAIJMKL(b->B, MATSEQAIJMKL, MAT_INPLACE_MATRIX, &b->B));
109:   PetscFunctionReturn(PETSC_SUCCESS);
110: }

112: PETSC_INTERN PetscErrorCode MatConvert_MPIAIJ_MPIAIJMKL(Mat A, MatType type, MatReuse reuse, Mat *newmat)
113: {
114:   Mat B = *newmat;

116:   PetscFunctionBegin;
117:   if (reuse == MAT_INITIAL_MATRIX) PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
118:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATMPIAIJMKL));
119:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMPIAIJSetPreallocation_C", MatMPIAIJSetPreallocation_MPIAIJMKL));
120:   *newmat = B;
121:   PetscFunctionReturn(PETSC_SUCCESS);
122: }

124: PETSC_EXTERN PetscErrorCode MatCreate_MPIAIJMKL(Mat A)
125: {
126:   PetscFunctionBegin;
127:   PetscCall(MatSetType(A, MATMPIAIJ));
128:   PetscCall(MatConvert_MPIAIJ_MPIAIJMKL(A, MATMPIAIJMKL, MAT_INPLACE_MATRIX, &A));
129:   PetscFunctionReturn(PETSC_SUCCESS);
130: }

132: /*MC
133:    MATAIJMKL - MATAIJMKL = "AIJMKL" - A matrix type to be used for sparse matrices.

135:    This matrix type is identical to `MATSEQAIJMKL` when constructed with a single process communicator,
136:    and `MATMPIAIJMKL` otherwise.  As a result, for single process communicators,
137:   MatSeqAIJSetPreallocation() is supported, and similarly `MatMPIAIJSetPreallocation()` is supported
138:   for communicators controlling multiple processes.  It is recommended that you call both of
139:   the above preallocation routines for simplicity.

141:    Options Database Key:
142: . -mat_type aijmkl - sets the matrix type to `MATAIJMKL` during a call to `MatSetFromOptions()`

144:   Level: beginner

146: .seealso: [](ch_matrices), `Mat`, `MATMPIAIJMKL`, `MATSEQAIJMKL`, `MatCreateMPIAIJMKL()`, `MATSEQAIJMKL`, `MATMPIAIJMKL`, `MATSEQAIJSELL`, `MATMPIAIJSELL`, `MATSEQAIJPERM`, `MATMPIAIJPERM`
147: M*/