Actual source code: matreg.c

  1: /*
  2:      Mechanism for register PETSc matrix types
  3: */
  4: #include <petsc/private/matimpl.h>

  6: PetscBool MatRegisterAllCalled = PETSC_FALSE;

  8: /*
  9:    Contains the list of registered Mat routines
 10: */
 11: PetscFunctionList MatList = NULL;

 13: /* MatGetRootType_Private - Gets the root type of the input matrix's type (e.g., MATAIJ for MATSEQAIJ)

 15:    Not Collective

 17:    Input Parameter:
 18: .  mat      - the input matrix, could be sequential or MPI

 20:    Output Parameter:
 21: .  rootType  - the root matrix type

 23:    Level: developer

 25: .seealso: `MatGetType()`, `MatSetType()`, `MatType`, `Mat`
 26: */
 27: PetscErrorCode MatGetRootType_Private(Mat mat, MatType *rootType)
 28: {
 29:   PetscBool   found = PETSC_FALSE;
 30:   MatRootName names = MatRootNameList;
 31:   MatType     inType;

 33:   PetscFunctionBegin;
 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: }

 50: /* MatGetMPIMatType_Private - Gets the MPI type corresponding to the input matrix's type (e.g., MATMPIAIJ for MATSEQAIJ)

 52:    Not Collective

 54:    Input Parameter:
 55: .  mat      - the input matrix, could be sequential or MPI

 57:    Output Parameter:
 58: .  MPIType  - the parallel (MPI) matrix type

 60:    Level: developer

 62: .seealso: `MatGetType()`, `MatSetType()`, `MatType`, `Mat`
 63: */
 64: PetscErrorCode MatGetMPIMatType_Private(Mat mat, MatType *MPIType)
 65: {
 66:   PetscBool   found = PETSC_FALSE;
 67:   MatRootName names = MatRootNameList;
 68:   MatType     inType;

 70:   PetscFunctionBegin;
 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: }

 88: /*@
 89:   MatSetType - Builds matrix object for a particular matrix type

 91:   Collective

 93:   Input Parameters:
 94: + mat    - the matrix object
 95: - matype - matrix type

 97:   Options Database Key:
 98: . -mat_type  <method> - Sets the type; see `MatType`

100:   Level: intermediate

102:   Note:
103:   See `MatType` for possible values

105: .seealso: [](ch_matrices), `Mat`, `PCSetType()`, `VecSetType()`, `MatCreate()`, `MatType`
106: @*/
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);

113:   PetscFunctionBegin;

116:   /* make this special case fast */
117:   PetscCall(PetscObjectTypeCompare((PetscObject)mat, matype, &sametype));
118:   if (sametype) PetscFunctionReturn(PETSC_SUCCESS);

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);

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

140:   PetscCall(PetscObjectTypeCompare((PetscObject)mat, matype, &sametype));
141:   if (sametype) PetscFunctionReturn(PETSC_SUCCESS);

143:   PetscCall(PetscFunctionListFind(MatList, matype, &r));
144:   PetscCheck(r, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown Mat type given: %s", matype);

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;

162:   /* should these null spaces be removed? */
163:   PetscCall(MatNullSpaceDestroy(&mat->nullsp));
164:   PetscCall(MatNullSpaceDestroy(&mat->nearnullsp));

166:   PetscCall(PetscMemzero(mat->ops, sizeof(struct _MatOps)));
167:   mat->preallocated  = PETSC_FALSE;
168:   mat->assembled     = PETSC_FALSE;
169:   mat->was_assembled = PETSC_FALSE;

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));

179:   /* create the new data structure */
180:   PetscCall((*r)(mat));
181:   PetscFunctionReturn(PETSC_SUCCESS);
182: }

184: /*@
185:   MatGetType - Gets the matrix type as a string from the matrix object.

187:   Not Collective

189:   Input Parameter:
190: . mat - the matrix

192:   Output Parameter:
193: . type - name of matrix type

195:   Level: intermediate

197: .seealso: [](ch_matrices), `Mat`, `MatType`, `MatSetType()`
198: @*/
199: PetscErrorCode MatGetType(Mat mat, MatType *type)
200: {
201:   PetscFunctionBegin;
203:   PetscAssertPointer(type, 2);
204:   *type = ((PetscObject)mat)->type_name;
205:   PetscFunctionReturn(PETSC_SUCCESS);
206: }

208: /*@
209:   MatGetVecType - Gets the vector type the matrix will return with `MatCreateVecs()`

211:   Not Collective

213:   Input Parameter:
214: . mat - the matrix

216:   Output Parameter:
217: . vtype - name of vector type

219:   Level: intermediate

221: .seealso: [](ch_matrices), `Mat`, `MatType`, `MatSetVecType()`, `VecType`
222: @*/
223: PetscErrorCode MatGetVecType(Mat mat, VecType *vtype)
224: {
225:   PetscFunctionBegin;
227:   PetscAssertPointer(vtype, 2);
228:   *vtype = mat->defaultvectype;
229:   PetscFunctionReturn(PETSC_SUCCESS);
230: }

232: /*@
233:   MatSetVecType - Set the vector type the matrix will return with `MatCreateVecs()`

235:   Collective

237:   Input Parameters:
238: + mat   - the matrix object
239: - vtype - vector type

241:   Level: advanced

243:   Note:
244:   This is rarely needed in practice since each matrix object internally sets the proper vector type.

246: .seealso: [](ch_matrices), `Mat`, `VecType`, `VecSetType()`, `MatGetVecType()`
247: @*/
248: PetscErrorCode MatSetVecType(Mat mat, VecType vtype)
249: {
250:   PetscFunctionBegin;
252:   PetscCall(PetscFree(mat->defaultvectype));
253:   PetscCall(PetscStrallocpy(vtype, &mat->defaultvectype));
254:   PetscFunctionReturn(PETSC_SUCCESS);
255: }

257: /*@C
258:   MatRegister -  - Adds a new matrix type implementation

260:   Not Collective, No Fortran Support

262:   Input Parameters:
263: + sname    - name of a new user-defined matrix type
264: - function - routine to create method context

266:   Level: advanced

268:   Note:
269:   `MatRegister()` may be called multiple times to add several user-defined solvers.

271:   Example Usage:
272: .vb
273:    MatRegister("my_mat", MyMatCreate);
274: .ve

276:   Then, your solver can be chosen with the procedural interface via `MatSetType(Mat, "my_mat")` or at runtime via the option
277:   `-mat_type my_mat`

279: .seealso: [](ch_matrices), `Mat`, `MatType`, `MatSetType()`, `MatRegisterAll()`
280: @*/
281: PetscErrorCode MatRegister(const char sname[], PetscErrorCode (*function)(Mat))
282: {
283:   PetscFunctionBegin;
284:   PetscCall(MatInitializePackage());
285:   PetscCall(PetscFunctionListAdd(&MatList, sname, function));
286:   PetscFunctionReturn(PETSC_SUCCESS);
287: }

289: MatRootName MatRootNameList = NULL;

291: /*@
292:   MatRegisterRootName - Registers a name that can be used for either a sequential or its corresponding parallel matrix type.

294:   Input Parameters:
295: + rname - the rootname, for example, `MATAIJ`
296: . sname - the name of the sequential matrix type, for example, `MATSEQAIJ`
297: - mname - the name of the parallel matrix type, for example, `MATMPIAIJ`

299:   Level: developer

301:   Notes:
302:   `MatSetType()` and `-mat_type name` will automatically use the sequential or parallel version
303:   based on the size of the MPI communicator associated with the matrix.

305:   The matrix rootname should not be confused with the base type of the function
306:   `PetscObjectBaseTypeCompare()`

308:   Developer Notes:
309:   PETSc vectors have a similar rootname that indicates PETSc should automatically select the appropriate `VecType` based on the
310:   size of the communicator but it is implemented by simply having additional `VecCreate_RootName()` registerer routines that dispatch to the
311:   appropriate creation routine. Why have two different ways of implementing the same functionality for different types of objects? It is
312:   confusing.

314: .seealso: [](ch_matrices), `Mat`, `MatType`, `PetscObjectBaseTypeCompare()`
315: @*/
316: PetscErrorCode MatRegisterRootName(const char rname[], const char sname[], const char mname[])
317: {
318:   MatRootName names;

320:   PetscFunctionBegin;
321:   PetscCall(PetscNew(&names));
322:   PetscCall(PetscStrallocpy(rname, &names->rname));
323:   PetscCall(PetscStrallocpy(sname, &names->sname));
324:   PetscCall(PetscStrallocpy(mname, &names->mname));
325:   if (!MatRootNameList) {
326:     MatRootNameList = names;
327:   } else {
328:     MatRootName next = MatRootNameList;
329:     while (next->next) next = next->next;
330:     next->next = names;
331:   }
332:   PetscFunctionReturn(PETSC_SUCCESS);
333: }