1 /*
2 Mechanism for register PETSc matrix types
3 */
4 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
5
6 PetscBool MatRegisterAllCalled = PETSC_FALSE;
7
8 /*
9 Contains the list of registered Mat routines
10 */
11 PetscFunctionList MatList = NULL;
12
13 /* MatGetRootType_Private - Gets the root type of the input matrix's type (e.g., MATAIJ for MATSEQAIJ)
14
15 Not Collective
16
17 Input Parameter:
18 . mat - the input matrix, could be sequential or MPI
19
20 Output Parameter:
21 . rootType - the root matrix type
22
23 Level: developer
24
25 .seealso: `MatGetType()`, `MatSetType()`, `MatType`, `Mat`
26 */
MatGetRootType_Private(Mat mat,MatType * rootType)27 PetscErrorCode MatGetRootType_Private(Mat mat, MatType *rootType)
28 {
29 PetscBool found = PETSC_FALSE;
30 MatRootName names = MatRootNameList;
31 MatType inType;
32
33 PetscFunctionBegin;
34 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1);
35 PetscCall(MatGetType(mat, &inType));
36 while (names) {
37 PetscCall(PetscStrcmp(inType, names->mname, &found));
38 if (!found) PetscCall(PetscStrcmp(inType, names->sname, &found));
39 if (found) {
40 found = PETSC_TRUE;
41 *rootType = names->rname;
42 break;
43 }
44 names = names->next;
45 }
46 if (!found) *rootType = inType;
47 PetscFunctionReturn(PETSC_SUCCESS);
48 }
49
50 /* MatGetMPIMatType_Private - Gets the MPI type corresponding to the input matrix's type (e.g., MATMPIAIJ for MATSEQAIJ)
51
52 Not Collective
53
54 Input Parameter:
55 . mat - the input matrix, could be sequential or MPI
56
57 Output Parameter:
58 . MPIType - the parallel (MPI) matrix type
59
60 Level: developer
61
62 .seealso: `MatGetType()`, `MatSetType()`, `MatType`, `Mat`
63 */
MatGetMPIMatType_Private(Mat mat,MatType * MPIType)64 PetscErrorCode MatGetMPIMatType_Private(Mat mat, MatType *MPIType)
65 {
66 PetscBool found = PETSC_FALSE;
67 MatRootName names = MatRootNameList;
68 MatType inType;
69
70 PetscFunctionBegin;
71 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1);
72 PetscCall(MatGetType(mat, &inType));
73 while (names) {
74 PetscCall(PetscStrcmp(inType, names->sname, &found));
75 if (!found) PetscCall(PetscStrcmp(inType, names->mname, &found));
76 if (!found) PetscCall(PetscStrcmp(inType, names->rname, &found));
77 if (found) {
78 found = PETSC_TRUE;
79 *MPIType = names->mname;
80 break;
81 }
82 names = names->next;
83 }
84 PetscCheck(found, PETSC_COMM_SELF, PETSC_ERR_SUP, "No corresponding parallel (MPI) type for this matrix");
85 PetscFunctionReturn(PETSC_SUCCESS);
86 }
87
88 /*@
89 MatSetType - Builds matrix object for a particular matrix type
90
91 Collective
92
93 Input Parameters:
94 + mat - the matrix object
95 - matype - matrix type
96
97 Options Database Key:
98 . -mat_type <method> - Sets the type; see `MatType`
99
100 Level: intermediate
101
102 Note:
103 See `MatType` for possible values
104
105 .seealso: [](ch_matrices), `Mat`, `PCSetType()`, `VecSetType()`, `MatCreate()`, `MatType`
106 @*/
MatSetType(Mat mat,MatType matype)107 PetscErrorCode MatSetType(Mat mat, MatType matype)
108 {
109 PetscBool sametype, found, subclass = PETSC_FALSE, matMPI = PETSC_FALSE, requestSeq = PETSC_FALSE;
110 MatRootName names = MatRootNameList;
111 PetscErrorCode (*r)(Mat);
112
113 PetscFunctionBegin;
114 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1);
115
116 /* make this special case fast */
117 PetscCall(PetscObjectTypeCompare((PetscObject)mat, matype, &sametype));
118 if (sametype) PetscFunctionReturn(PETSC_SUCCESS);
119
120 /* suppose with one MPI process, one created an MPIAIJ (mpiaij) matrix with MatCreateMPIAIJWithArrays(), and later tried
121 to change its type via '-mat_type aijcusparse'. Even there is only one MPI rank, we need to adapt matype to
122 'mpiaijcusparse' so that it will be treated as a subclass of MPIAIJ and proper MatCovert() will be called.
123 */
124 if (((PetscObject)mat)->type_name) PetscCall(PetscStrbeginswith(((PetscObject)mat)->type_name, "mpi", &matMPI)); /* mat claims itself is an 'mpi' matrix */
125 if (matype) PetscCall(PetscStrbeginswith(matype, "seq", &requestSeq)); /* user is requesting a 'seq' matrix */
126 PetscCheck(!(matMPI && requestSeq), PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Changing an MPI matrix (%s) to a sequential one (%s) is not allowed. Please remove the 'seq' prefix from your matrix type.", ((PetscObject)mat)->type_name, matype);
127
128 while (names) {
129 PetscCall(PetscStrcmp(matype, names->rname, &found));
130 if (found) {
131 PetscMPIInt size;
132 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)mat), &size));
133 if (size == 1 && !matMPI) matype = names->sname; /* try to align the requested type (matype) with the existing type per seq/mpi */
134 else matype = names->mname;
135 break;
136 }
137 names = names->next;
138 }
139
140 PetscCall(PetscObjectTypeCompare((PetscObject)mat, matype, &sametype));
141 if (sametype) PetscFunctionReturn(PETSC_SUCCESS);
142
143 PetscCall(PetscFunctionListFind(MatList, matype, &r));
144 PetscCheck(r, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown Mat type given: %s", matype);
145
146 if (mat->assembled && ((PetscObject)mat)->type_name) PetscCall(PetscStrbeginswith(matype, ((PetscObject)mat)->type_name, &subclass));
147 if (subclass) { /* mat is a subclass of the requested 'matype'? */
148 PetscCall(MatConvert(mat, matype, MAT_INPLACE_MATRIX, &mat));
149 PetscFunctionReturn(PETSC_SUCCESS);
150 }
151 if (names && mat->assembled) {
152 PetscCall(PetscStrbeginswith(names->rname, "sell", &sametype));
153 if (sametype) { /* mattype is MATSELL or its subclass */
154 PetscCall(MatConvert(mat, MATSELL, MAT_INPLACE_MATRIX, &mat)); /* convert to matsell first */
155 PetscCall(MatConvert(mat, matype, MAT_INPLACE_MATRIX, &mat));
156 PetscFunctionReturn(PETSC_SUCCESS);
157 }
158 }
159 PetscTryTypeMethod(mat, destroy);
160 mat->ops->destroy = NULL;
161
162 /* should these null spaces be removed? */
163 PetscCall(MatNullSpaceDestroy(&mat->nullsp));
164 PetscCall(MatNullSpaceDestroy(&mat->nearnullsp));
165
166 PetscCall(PetscMemzero(mat->ops, sizeof(struct _MatOps)));
167 mat->preallocated = PETSC_FALSE;
168 mat->assembled = PETSC_FALSE;
169 mat->was_assembled = PETSC_FALSE;
170
171 /*
172 Increment, rather than reset these: the object is logically the same, so its logging and
173 state is inherited. Furthermore, resetting makes it possible for the same state to be
174 obtained with a different structure, confusing the PC.
175 */
176 mat->nonzerostate++;
177 PetscCall(PetscObjectStateIncrease((PetscObject)mat));
178
179 /* create the new data structure */
180 PetscCall((*r)(mat));
181 PetscFunctionReturn(PETSC_SUCCESS);
182 }
183
184 /*@
185 MatGetType - Gets the matrix type as a string from the matrix object.
186
187 Not Collective
188
189 Input Parameter:
190 . mat - the matrix
191
192 Output Parameter:
193 . type - name of matrix type
194
195 Level: intermediate
196
197 .seealso: [](ch_matrices), `Mat`, `MatType`, `MatSetType()`
198 @*/
MatGetType(Mat mat,MatType * type)199 PetscErrorCode MatGetType(Mat mat, MatType *type)
200 {
201 PetscFunctionBegin;
202 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1);
203 PetscAssertPointer(type, 2);
204 *type = ((PetscObject)mat)->type_name;
205 PetscFunctionReturn(PETSC_SUCCESS);
206 }
207
208 /*@
209 MatGetVecType - Gets the vector type the matrix will return with `MatCreateVecs()`
210
211 Not Collective
212
213 Input Parameter:
214 . mat - the matrix
215
216 Output Parameter:
217 . vtype - name of vector type
218
219 Level: intermediate
220
221 .seealso: [](ch_matrices), `Mat`, `MatType`, `MatSetVecType()`, `VecType`
222 @*/
MatGetVecType(Mat mat,VecType * vtype)223 PetscErrorCode MatGetVecType(Mat mat, VecType *vtype)
224 {
225 PetscFunctionBegin;
226 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1);
227 PetscAssertPointer(vtype, 2);
228 *vtype = mat->defaultvectype;
229 PetscFunctionReturn(PETSC_SUCCESS);
230 }
231
232 /*@
233 MatSetVecType - Set the vector type the matrix will return with `MatCreateVecs()`
234
235 Collective
236
237 Input Parameters:
238 + mat - the matrix object
239 - vtype - vector type
240
241 Level: advanced
242
243 Note:
244 This is rarely needed in practice since each matrix object internally sets the proper vector type.
245
246 .seealso: [](ch_matrices), `Mat`, `VecType`, `VecSetType()`, `MatGetVecType()`
247 @*/
MatSetVecType(Mat mat,VecType vtype)248 PetscErrorCode MatSetVecType(Mat mat, VecType vtype)
249 {
250 PetscFunctionBegin;
251 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1);
252 PetscCall(PetscFree(mat->defaultvectype));
253 PetscCall(PetscStrallocpy(vtype, &mat->defaultvectype));
254 PetscFunctionReturn(PETSC_SUCCESS);
255 }
256
257 /*@C
258 MatRegister - Adds a new matrix type implementation that is usable as a `Mat` in PETSc
259
260 Not Collective, No Fortran Support
261
262 Input Parameters:
263 + sname - name of a new user-defined matrix type
264 - function - routine to create method context
265
266 Level: advanced
267
268 Note:
269 `MatRegister()` may be called multiple times to add several user-defined solvers.
270
271 A simpler alternative to using `MatRegister()` for an application specific matrix format is to use `MatCreateShell()`, which
272 generates a `Mat` of `MatType` `MATSHELL`. One can then use `MatShellSetContext()` and `MatShellSetOperation()` to provide
273 the data structures and routines customized for their matrix format.
274
275 Example Usage:
276 .vb
277 MatRegister("my_mat", MyMatCreate);
278 .ve
279
280 Then, your solver can be chosen with the procedural interface via `MatSetType(Mat, "my_mat")` or at runtime via the option
281 `-mat_type my_mat`
282
283 .seealso: [](ch_matrices), `Mat`, `MatType`, `MatSetType()`, `MatRegisterAll()`
284 @*/
MatRegister(const char sname[],PetscErrorCode (* function)(Mat))285 PetscErrorCode MatRegister(const char sname[], PetscErrorCode (*function)(Mat))
286 {
287 PetscFunctionBegin;
288 PetscCall(MatInitializePackage());
289 PetscCall(PetscFunctionListAdd(&MatList, sname, function));
290 PetscFunctionReturn(PETSC_SUCCESS);
291 }
292
293 MatRootName MatRootNameList = NULL;
294
295 /*@
296 MatRegisterRootName - Registers a name that can be used for either a sequential or its corresponding parallel matrix type.
297
298 Input Parameters:
299 + rname - the rootname, for example, `MATAIJ`
300 . sname - the name of the sequential matrix type, for example, `MATSEQAIJ`
301 - mname - the name of the parallel matrix type, for example, `MATMPIAIJ`
302
303 Level: developer
304
305 Notes:
306 `MatSetType()` and `-mat_type name` will automatically use the sequential or parallel version
307 based on the size of the MPI communicator associated with the matrix.
308
309 The matrix rootname should not be confused with the base type of the function
310 `PetscObjectBaseTypeCompare()`
311
312 Developer Notes:
313 PETSc vectors have a similar rootname that indicates PETSc should automatically select the appropriate `VecType` based on the
314 size of the communicator but it is implemented by simply having additional `VecCreate_RootName()` registerer routines that dispatch to the
315 appropriate creation routine. Why have two different ways of implementing the same functionality for different types of objects? It is
316 confusing.
317
318 .seealso: [](ch_matrices), `Mat`, `MatType`, `PetscObjectBaseTypeCompare()`
319 @*/
MatRegisterRootName(const char rname[],const char sname[],const char mname[])320 PetscErrorCode MatRegisterRootName(const char rname[], const char sname[], const char mname[])
321 {
322 MatRootName names;
323
324 PetscFunctionBegin;
325 PetscCall(PetscNew(&names));
326 PetscCall(PetscStrallocpy(rname, &names->rname));
327 PetscCall(PetscStrallocpy(sname, &names->sname));
328 PetscCall(PetscStrallocpy(mname, &names->mname));
329 if (!MatRootNameList) {
330 MatRootNameList = names;
331 } else {
332 MatRootName next = MatRootNameList;
333 while (next->next) next = next->next;
334 next->next = names;
335 }
336 PetscFunctionReturn(PETSC_SUCCESS);
337 }
338