1 #include <petsc/private/matimpl.h> 2 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; 14 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"); 23 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)); 35 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 } 42 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)); 51 52 if (reuse == MAT_INPLACE_MATRIX) { 53 PetscCall(MatHeaderReplace(mat, &M)); 54 } else { 55 *newmat = M; 56 } 57 PetscFunctionReturn(PETSC_SUCCESS); 58 } 59