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 on mat

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

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

135:   PetscObjectTypeCompare((PetscObject)mat, matype, &sametype);
136:   if (sametype) return 0;

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

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

149:   /* should these null spaces be removed? */
150:   MatNullSpaceDestroy(&mat->nullsp);
151:   MatNullSpaceDestroy(&mat->nearnullsp);

153:   PetscMemzero(mat->ops, sizeof(struct _MatOps));
154:   mat->preallocated  = PETSC_FALSE;
155:   mat->assembled     = PETSC_FALSE;
156:   mat->was_assembled = PETSC_FALSE;

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

166:   /* create the new data structure */
167:   (*r)(mat);
168:   return 0;
169: }

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

174:    Not Collective

176:    Input Parameter:
177: .  mat - the matrix

179:    Output Parameter:
180: .  name - name of matrix type

182:    Level: intermediate

184: .seealso: `MatType`, `MatSetType()`
185: @*/
186: PetscErrorCode MatGetType(Mat mat, MatType *type)
187: {
190:   *type = ((PetscObject)mat)->type_name;
191:   return 0;
192: }

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

197:    Not Collective

199:    Input Parameter:
200: .  mat - the matrix

202:    Output Parameter:
203: .  name - name of vector type

205:    Level: intermediate

207: .seealso: `MatType`, `Mat`, `MatSetVecType()`, `VecType`
208: @*/
209: PetscErrorCode MatGetVecType(Mat mat, VecType *vtype)
210: {
213:   *vtype = mat->defaultvectype;
214:   return 0;
215: }

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

220:    Collective on mat

222:    Input Parameters:
223: +  mat   - the matrix object
224: -  vtype - vector type

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

229:   Level: intermediate

231: .seealso: `VecType`, `VecSetType()`, `MatGetVecType()`
232: @*/
233: PetscErrorCode MatSetVecType(Mat mat, VecType vtype)
234: {
236:   PetscFree(mat->defaultvectype);
237:   PetscStrallocpy(vtype, &mat->defaultvectype);
238:   return 0;
239: }

241: /*@C
242:   MatRegister -  - Adds a new matrix type

244:    Not Collective

246:    Input Parameters:
247: +  name - name of a new user-defined matrix type
248: -  routine_create - routine to create method context

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

253:    Sample usage:
254: .vb
255:    MatRegister("my_mat",MyMatCreate);
256: .ve

258:    Then, your solver can be chosen with the procedural interface via
259: $     MatSetType(Mat,"my_mat")
260:    or at runtime via the option
261: $     -mat_type my_mat

263:    Level: advanced

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

274: MatRootName MatRootNameList = NULL;

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

281:   Input Parameters:
282: +     rname - the rootname, for example, `MATAIJ`
283: .     sname - the name of the sequential matrix type, for example, `MATSEQAIJ`
284: -     mname - the name of the parallel matrix type, for example, `MATMPIAIJ`

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

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

295:   Level: developer

297: .seealso: `Mat`, `MatType`, `PetscObjectBaseTypeCompare()`
298: @*/
299: PetscErrorCode MatRegisterRootName(const char rname[], const char sname[], const char mname[])
300: {
301:   MatRootName names;

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