1 #include <petscmat.h> 2 #include <petsc/private/petscimpl.h> 3 4 /* 5 MatDumpSPAI - Dumps a PETSc matrix to a file in an ASCII format 6 suitable for the SPAI code of Stephen Barnard to solve. This routine 7 is simply here to allow testing of matrices directly with the SPAI 8 code, rather than through the PETSc interface. 9 10 */ 11 PetscErrorCode MatDumpSPAI(Mat A, FILE *file) 12 { 13 PetscMPIInt size; 14 PetscInt n; 15 MPI_Comm comm; 16 17 PetscFunctionBegin; 18 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 19 PetscAssertPointer(file, 2); 20 PetscCall(PetscObjectGetComm((PetscObject)A, &comm)); 21 PetscCallMPI(MPI_Comm_size(comm, &size)); 22 PetscCheck(size == 1, comm, PETSC_ERR_SUP, "Only single processor dumps"); 23 PetscCall(MatGetSize(A, &n, &n)); 24 /* print the matrix */ 25 fprintf(file, "%" PetscInt_FMT "\n", n); 26 for (PetscInt i = 0; i < n; i++) { 27 const PetscInt *cols; 28 const PetscScalar *vals; 29 PetscInt nz; 30 31 PetscCall(MatGetRow(A, i, &nz, &cols, &vals)); 32 for (PetscInt j = 0; j < nz; j++) fprintf(file, "%" PetscInt_FMT " %d" PetscInt_FMT " %16.14e\n", i + 1, cols[j] + 1, vals[j]); 33 PetscCall(MatRestoreRow(A, i, &nz, &cols, &vals)); 34 } 35 PetscFunctionReturn(PETSC_SUCCESS); 36 } 37 38 PetscErrorCode VecDumpSPAI(Vec b, FILE *file) 39 { 40 PetscInt n; 41 const PetscScalar *array; 42 43 PetscFunctionBegin; 44 PetscValidHeaderSpecific(b, VEC_CLASSID, 1); 45 PetscAssertPointer(file, 2); 46 PetscCall(VecGetSize(b, &n)); 47 PetscCall(VecGetArrayRead(b, &array)); 48 fprintf(file, "%" PetscInt_FMT "\n", n); 49 for (PetscInt i = 0; i < n; i++) fprintf(file, "%" PetscInt_FMT " %16.14e\n", i + 1, array[i]); 50 PetscCall(VecRestoreArrayRead(b, &array)); 51 PetscFunctionReturn(PETSC_SUCCESS); 52 } 53