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: /*@C
 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:   PetscTryTypeMethod(mat, destroy);
152:   mat->ops->destroy = NULL;

154:   /* should these null spaces be removed? */
155:   PetscCall(MatNullSpaceDestroy(&mat->nullsp));
156:   PetscCall(MatNullSpaceDestroy(&mat->nearnullsp));

158:   PetscCall(PetscMemzero(mat->ops, sizeof(struct _MatOps)));
159:   mat->preallocated  = PETSC_FALSE;
160:   mat->assembled     = PETSC_FALSE;
161:   mat->was_assembled = PETSC_FALSE;

163:   /*
164:    Increment, rather than reset these: the object is logically the same, so its logging and
165:    state is inherited.  Furthermore, resetting makes it possible for the same state to be
166:    obtained with a different structure, confusing the PC.
167:   */
168:   mat->nonzerostate++;
169:   PetscCall(PetscObjectStateIncrease((PetscObject)mat));

171:   /* create the new data structure */
172:   PetscCall((*r)(mat));
173:   PetscFunctionReturn(PETSC_SUCCESS);
174: }

176: /*@C
177:   MatGetType - Gets the matrix type as a string from the matrix object.

179:   Not Collective

181:   Input Parameter:
182: . mat - the matrix

184:   Output Parameter:
185: . type - name of matrix type

187:   Level: intermediate

189: .seealso: [](ch_matrices), `Mat`, `MatType`, `MatSetType()`
190: @*/
191: PetscErrorCode MatGetType(Mat mat, MatType *type)
192: {
193:   PetscFunctionBegin;
195:   PetscAssertPointer(type, 2);
196:   *type = ((PetscObject)mat)->type_name;
197:   PetscFunctionReturn(PETSC_SUCCESS);
198: }

200: /*@C
201:   MatGetVecType - Gets the vector type the matrix will return with `MatCreateVecs()`

203:   Not Collective

205:   Input Parameter:
206: . mat - the matrix

208:   Output Parameter:
209: . vtype - name of vector type

211:   Level: intermediate

213: .seealso: [](ch_matrices), `Mat`, `MatType`, `MatSetVecType()`, `VecType`
214: @*/
215: PetscErrorCode MatGetVecType(Mat mat, VecType *vtype)
216: {
217:   PetscFunctionBegin;
219:   PetscAssertPointer(vtype, 2);
220:   *vtype = mat->defaultvectype;
221:   PetscFunctionReturn(PETSC_SUCCESS);
222: }

224: /*@C
225:   MatSetVecType - Set the vector type the matrix will return with `MatCreateVecs()`

227:   Collective

229:   Input Parameters:
230: + mat   - the matrix object
231: - vtype - vector type

233:   Level: advanced

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

238: .seealso: [](ch_matrices), `Mat`, `VecType`, `VecSetType()`, `MatGetVecType()`
239: @*/
240: PetscErrorCode MatSetVecType(Mat mat, VecType vtype)
241: {
242:   PetscFunctionBegin;
244:   PetscCall(PetscFree(mat->defaultvectype));
245:   PetscCall(PetscStrallocpy(vtype, &mat->defaultvectype));
246:   PetscFunctionReturn(PETSC_SUCCESS);
247: }

249: /*@C
250:   MatRegister -  - Adds a new matrix type implementation

252:   Not Collective

254:   Input Parameters:
255: + sname    - name of a new user-defined matrix type
256: - function - routine to create method context

258:   Level: advanced

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

263:   Example Usage:
264: .vb
265:    MatRegister("my_mat", MyMatCreate);
266: .ve

268:   Then, your solver can be chosen with the procedural interface via
269: $     MatSetType(Mat, "my_mat")
270:   or at runtime via the option
271: $     -mat_type my_mat

273: .seealso: [](ch_matrices), `Mat`, `MatType`, `MatSetType()`, `MatRegisterAll()`
274: @*/
275: PetscErrorCode MatRegister(const char sname[], PetscErrorCode (*function)(Mat))
276: {
277:   PetscFunctionBegin;
278:   PetscCall(MatInitializePackage());
279:   PetscCall(PetscFunctionListAdd(&MatList, sname, function));
280:   PetscFunctionReturn(PETSC_SUCCESS);
281: }

283: MatRootName MatRootNameList = NULL;

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

288:   Input Parameters:
289: + rname - the rootname, for example, `MATAIJ`
290: . sname - the name of the sequential matrix type, for example, `MATSEQAIJ`
291: - mname - the name of the parallel matrix type, for example, `MATMPIAIJ`

293:   Level: developer

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

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

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

308: .seealso: [](ch_matrices), `Mat`, `MatType`, `PetscObjectBaseTypeCompare()`
309: @*/
310: PetscErrorCode MatRegisterRootName(const char rname[], const char sname[], const char mname[])
311: {
312:   MatRootName names;

314:   PetscFunctionBegin;
315:   PetscCall(PetscNew(&names));
316:   PetscCall(PetscStrallocpy(rname, &names->rname));
317:   PetscCall(PetscStrallocpy(sname, &names->sname));
318:   PetscCall(PetscStrallocpy(mname, &names->mname));
319:   if (!MatRootNameList) {
320:     MatRootNameList = names;
321:   } else {
322:     MatRootName next = MatRootNameList;
323:     while (next->next) next = next->next;
324:     next->next = names;
325:   }
326:   PetscFunctionReturn(PETSC_SUCCESS);
327: }