Actual source code: matreg.c


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

  7: PetscBool MatRegisterAllCalled = PETSC_FALSE;

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

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

 16:    Not Collective

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

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

 24:    Level: developer

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

 34:   PetscFunctionBegin;
 36:   PetscCall(MatGetType(mat, &inType));
 37:   while (names) {
 38:     PetscCall(PetscStrcmp(inType, names->mname, &found));
 39:     if (!found) PetscCall(PetscStrcmp(inType, names->sname, &found));
 40:     if (found) {
 41:       found     = PETSC_TRUE;
 42:       *rootType = names->rname;
 43:       break;
 44:     }
 45:     names = names->next;
 46:   }
 47:   if (!found) *rootType = inType;
 48:   PetscFunctionReturn(PETSC_SUCCESS);
 49: }

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

 53:    Not Collective

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

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

 61:    Level: developer

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

 71:   PetscFunctionBegin;
 73:   PetscCall(MatGetType(mat, &inType));
 74:   while (names) {
 75:     PetscCall(PetscStrcmp(inType, names->sname, &found));
 76:     if (!found) PetscCall(PetscStrcmp(inType, names->mname, &found));
 77:     if (!found) PetscCall(PetscStrcmp(inType, names->rname, &found));
 78:     if (found) {
 79:       found    = PETSC_TRUE;
 80:       *MPIType = names->mname;
 81:       break;
 82:     }
 83:     names = names->next;
 84:   }
 85:   PetscCheck(found, PETSC_COMM_SELF, PETSC_ERR_SUP, "No corresponding parallel (MPI) type for this matrix");
 86:   PetscFunctionReturn(PETSC_SUCCESS);
 87: }

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

 92:    Collective

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

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

101:   Level: intermediate

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

106: .seealso: [](ch_matrices), `Mat`, `PCSetType()`, `VecSetType()`, `MatCreate()`, `MatType`, `Mat`
107: @*/
108: PetscErrorCode MatSetType(Mat mat, MatType matype)
109: {
110:   PetscBool   sametype, found, subclass = PETSC_FALSE, matMPI = PETSC_FALSE, requestSeq = PETSC_FALSE;
111:   MatRootName names = MatRootNameList;
112:   PetscErrorCode (*r)(Mat);

114:   PetscFunctionBegin;

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

121:   /* suppose with one MPI process, one created an MPIAIJ (mpiaij) matrix with MatCreateMPIAIJWithArrays(), and later tried
122:      to change its type via '-mat_type aijcusparse'. Even there is only one MPI rank, we need to adapt matype to
123:      'mpiaijcusparse' so that it will be treated as a subclass of MPIAIJ and proper MatCovert() will be called.
124:   */
125:   if (((PetscObject)mat)->type_name) PetscCall(PetscStrbeginswith(((PetscObject)mat)->type_name, "mpi", &matMPI)); /* mat claims itself is an 'mpi' matrix */
126:   if (matype) PetscCall(PetscStrbeginswith(matype, "seq", &requestSeq));                                           /* user is requesting a 'seq' matrix */
127:   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);

129:   while (names) {
130:     PetscCall(PetscStrcmp(matype, names->rname, &found));
131:     if (found) {
132:       PetscMPIInt size;
133:       PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)mat), &size));
134:       if (size == 1 && !matMPI) matype = names->sname; /* try to align the requested type (matype) with the existing type per seq/mpi */
135:       else matype = names->mname;
136:       break;
137:     }
138:     names = names->next;
139:   }

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

144:   PetscCall(PetscFunctionListFind(MatList, matype, &r));
145:   PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown Mat type given: %s", matype);

147:   if (mat->assembled && ((PetscObject)mat)->type_name) PetscCall(PetscStrbeginswith(matype, ((PetscObject)mat)->type_name, &subclass));
148:   if (subclass) { /* mat is a subclass of the requested 'matype'? */
149:     PetscCall(MatConvert(mat, matype, MAT_INPLACE_MATRIX, &mat));
150:     PetscFunctionReturn(PETSC_SUCCESS);
151:   }
152:   PetscTryTypeMethod(mat, destroy);
153:   mat->ops->destroy = NULL;

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

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

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

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

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

180:    Not Collective

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

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

188:    Level: intermediate

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

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

204:    Not Collective

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

209:    Output Parameter:
210: .  name - name of vector type

212:    Level: intermediate

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

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

228:    Collective

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

234:   Level: advanced

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

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

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

253:    Not Collective

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

259:    Level: advanced

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

264:    Sample usage:
265: .vb
266:    MatRegister("my_mat", MyMatCreate);
267: .ve

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

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

284: MatRootName MatRootNameList = NULL;

286: /*@C
287:       MatRegisterRootName - Registers a name that can be used for either a sequential or its corresponding parallel matrix type. `MatSetType()`
288:         and `-mat_type name` will automatically use the sequential or parallel version based on the size of the MPI communicator associated with the
289:         matrix.

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

296:   Level: developer

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

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

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

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