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 Parameters:
 19: .  mat      - the input matrix, could be sequential or MPI

 21:    Output Parameters:
 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;

 35:   MatGetType(mat, &inType);
 36:   while (names) {
 37:     PetscStrcmp(inType, names->mname, &found);
 38:     if (!found) 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:   return 0;
 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 Parameters:
 55: .  mat      - the input matrix, could be sequential or MPI

 57:    Output Parameters:
 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;

 71:   MatGetType(mat, &inType);
 72:   while (names) {
 73:     PetscStrcmp(inType, names->sname, &found);
 74:     if (!found) PetscStrcmp(inType, names->mname, &found);
 75:     if (!found) PetscStrcmp(inType, names->rname, &found);
 76:     if (found) {
 77:       found    = PETSC_TRUE;
 78:       *MPIType = names->mname;
 79:       break;
 80:     }
 81:     names = names->next;
 82:   }
 84:   return 0;
 85: }

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

 90:    Collective

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

 96:    Options Database Key:
 97: .  -mat_type  <method> - Sets the type; use -help for a list
 98:     of available methods (for instance, seqaij)

100:    Note:
101:    See "${PETSC_DIR}/include/petscmat.h" for available methods

103:   Level: intermediate

105: .seealso: `PCSetType()`, `VecSetType()`, `MatCreate()`, `MatType`, `Mat`
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);


115:   /* make this special case fast */
116:   PetscObjectTypeCompare((PetscObject)mat, matype, &sametype);
117:   if (sametype) return 0;

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

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

139:   PetscObjectTypeCompare((PetscObject)mat, matype, &sametype);
140:   if (sametype) return 0;

142:   PetscFunctionListFind(MatList, matype, &r);

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

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

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

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

170:   /* create the new data structure */
171:   (*r)(mat);
172:   return 0;
173: }

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

178:    Not Collective

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

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

186:    Level: intermediate

188: .seealso: `MatType`, `MatSetType()`
189: @*/
190: PetscErrorCode MatGetType(Mat mat, MatType *type)
191: {
194:   *type = ((PetscObject)mat)->type_name;
195:   return 0;
196: }

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

201:    Not Collective

203:    Input Parameter:
204: .  mat - the matrix

206:    Output Parameter:
207: .  name - name of vector type

209:    Level: intermediate

211: .seealso: `MatType`, `Mat`, `MatSetVecType()`, `VecType`
212: @*/
213: PetscErrorCode MatGetVecType(Mat mat, VecType *vtype)
214: {
217:   *vtype = mat->defaultvectype;
218:   return 0;
219: }

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

224:    Collective

226:    Input Parameters:
227: +  mat   - the matrix object
228: -  vtype - vector type

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

233:   Level: intermediate

235: .seealso: `VecType`, `VecSetType()`, `MatGetVecType()`
236: @*/
237: PetscErrorCode MatSetVecType(Mat mat, VecType vtype)
238: {
240:   PetscFree(mat->defaultvectype);
241:   PetscStrallocpy(vtype, &mat->defaultvectype);
242:   return 0;
243: }

245: /*@C
246:   MatRegister -  - Adds a new matrix type

248:    Not Collective

250:    Input Parameters:
251: +  name - name of a new user-defined matrix type
252: -  routine_create - routine to create method context

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

257:    Sample usage:
258: .vb
259:    MatRegister("my_mat",MyMatCreate);
260: .ve

262:    Then, your solver can be chosen with the procedural interface via
263: $     MatSetType(Mat,"my_mat")
264:    or at runtime via the option
265: $     -mat_type my_mat

267:    Level: advanced

269: .seealso: `Mat`, `MatType`, `MatSetType()`, `MatRegisterAll()`
270: @*/
271: PetscErrorCode MatRegister(const char sname[], PetscErrorCode (*function)(Mat))
272: {
273:   MatInitializePackage();
274:   PetscFunctionListAdd(&MatList, sname, function);
275:   return 0;
276: }

278: MatRootName MatRootNameList = NULL;

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

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

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

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

299:   Level: developer

301: .seealso: `Mat`, `MatType`, `PetscObjectBaseTypeCompare()`
302: @*/
303: PetscErrorCode MatRegisterRootName(const char rname[], const char sname[], const char mname[])
304: {
305:   MatRootName names;

307:   PetscNew(&names);
308:   PetscStrallocpy(rname, &names->rname);
309:   PetscStrallocpy(sname, &names->sname);
310:   PetscStrallocpy(mname, &names->mname);
311:   if (!MatRootNameList) {
312:     MatRootNameList = names;
313:   } else {
314:     MatRootName next = MatRootNameList;
315:     while (next->next) next = next->next;
316:     next->next = names;
317:   }
318:   return 0;
319: }