Actual source code: convert.c
1: #include <petsc/private/matimpl.h>
3: /*
4: MatConvert_Basic - Converts from any input format to another format.
5: Does not do preallocation so in general will be slow
6: */
7: PETSC_INTERN PetscErrorCode MatConvert_Basic(Mat mat, MatType newtype, MatReuse reuse, Mat *newmat)
8: {
9: Mat M;
10: const PetscScalar *vwork;
11: PetscInt i, rstart, rend, nz;
12: const PetscInt *cwork;
13: PetscBool isSBAIJ;
15: PetscFunctionBegin;
16: if (!mat->ops->getrow) { /* missing get row, use matvecs */
17: PetscCall(MatConvert_Shell(mat, newtype, reuse, newmat));
18: PetscFunctionReturn(PETSC_SUCCESS);
19: }
20: PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATSEQSBAIJ, &isSBAIJ));
21: if (!isSBAIJ) PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATMPISBAIJ, &isSBAIJ));
22: PetscCheck(!isSBAIJ, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot convert from SBAIJ matrix since cannot obtain entire rows of matrix");
24: if (reuse == MAT_REUSE_MATRIX) {
25: M = *newmat;
26: } else {
27: PetscInt m, n, lm, ln;
28: PetscCall(MatGetSize(mat, &m, &n));
29: PetscCall(MatGetLocalSize(mat, &lm, &ln));
30: PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &M));
31: PetscCall(MatSetSizes(M, lm, ln, m, n));
32: PetscCall(MatSetBlockSizesFromMats(M, mat, mat));
33: PetscCall(MatSetType(M, newtype));
34: PetscCall(MatSetUp(M));
36: PetscCall(MatSetOption(M, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE));
37: PetscCall(MatSetOption(M, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
38: PetscCall(PetscObjectTypeCompare((PetscObject)M, MATSEQSBAIJ, &isSBAIJ));
39: if (!isSBAIJ) PetscCall(PetscObjectTypeCompare((PetscObject)M, MATMPISBAIJ, &isSBAIJ));
40: if (isSBAIJ) PetscCall(MatSetOption(M, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE));
41: }
43: PetscCall(MatGetOwnershipRange(mat, &rstart, &rend));
44: for (i = rstart; i < rend; i++) {
45: PetscCall(MatGetRow(mat, i, &nz, &cwork, &vwork));
46: PetscCall(MatSetValues(M, 1, &i, nz, cwork, vwork, INSERT_VALUES));
47: PetscCall(MatRestoreRow(mat, i, &nz, &cwork, &vwork));
48: }
49: PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY));
50: PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY));
52: if (reuse == MAT_INPLACE_MATRIX) {
53: PetscCall(MatHeaderReplace(mat, &M));
54: } else {
55: *newmat = M;
56: }
57: PetscFunctionReturn(PETSC_SUCCESS);
58: }