1c6db04a5SJed Brown #include <petscmat.h>
2af0996ceSBarry Smith #include <petsc/private/petscimpl.h>
301a81e61SBarry Smith
401a81e61SBarry Smith /*
501a81e61SBarry Smith MatDumpSPAI - Dumps a PETSc matrix to a file in an ASCII format
601a81e61SBarry Smith suitable for the SPAI code of Stephen Barnard to solve. This routine
701a81e61SBarry Smith is simply here to allow testing of matrices directly with the SPAI
8*467446fbSPierre Jolivet code, rather than through the PETSc interface.
901a81e61SBarry Smith
1001a81e61SBarry Smith */
MatDumpSPAI(Mat A,FILE * file)11d71ae5a4SJacob Faibussowitsch PetscErrorCode MatDumpSPAI(Mat A, FILE *file)
12d71ae5a4SJacob Faibussowitsch {
133ba16761SJacob Faibussowitsch PetscMPIInt size;
143ba16761SJacob Faibussowitsch PetscInt n;
1501a81e61SBarry Smith MPI_Comm comm;
1601a81e61SBarry Smith
175f80ce2aSJacob Faibussowitsch PetscFunctionBegin;
183ba16761SJacob Faibussowitsch PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
194f572ea9SToby Isaac PetscAssertPointer(file, 2);
203ba16761SJacob Faibussowitsch PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
213ba16761SJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(comm, &size));
223ba16761SJacob Faibussowitsch PetscCheck(size == 1, comm, PETSC_ERR_SUP, "Only single processor dumps");
239566063dSJacob Faibussowitsch PetscCall(MatGetSize(A, &n, &n));
2401a81e61SBarry Smith /* print the matrix */
253ba16761SJacob Faibussowitsch fprintf(file, "%" PetscInt_FMT "\n", n);
263ba16761SJacob Faibussowitsch for (PetscInt i = 0; i < n; i++) {
273ba16761SJacob Faibussowitsch const PetscInt *cols;
283ba16761SJacob Faibussowitsch const PetscScalar *vals;
293ba16761SJacob Faibussowitsch PetscInt nz;
303ba16761SJacob Faibussowitsch
319566063dSJacob Faibussowitsch PetscCall(MatGetRow(A, i, &nz, &cols, &vals));
323ba16761SJacob Faibussowitsch for (PetscInt j = 0; j < nz; j++) fprintf(file, "%" PetscInt_FMT " %d" PetscInt_FMT " %16.14e\n", i + 1, cols[j] + 1, vals[j]);
339566063dSJacob Faibussowitsch PetscCall(MatRestoreRow(A, i, &nz, &cols, &vals));
3401a81e61SBarry Smith }
353ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
3601a81e61SBarry Smith }
3701a81e61SBarry Smith
VecDumpSPAI(Vec b,FILE * file)38d71ae5a4SJacob Faibussowitsch PetscErrorCode VecDumpSPAI(Vec b, FILE *file)
39d71ae5a4SJacob Faibussowitsch {
403ba16761SJacob Faibussowitsch PetscInt n;
413ba16761SJacob Faibussowitsch const PetscScalar *array;
4201a81e61SBarry Smith
435f80ce2aSJacob Faibussowitsch PetscFunctionBegin;
443ba16761SJacob Faibussowitsch PetscValidHeaderSpecific(b, VEC_CLASSID, 1);
454f572ea9SToby Isaac PetscAssertPointer(file, 2);
469566063dSJacob Faibussowitsch PetscCall(VecGetSize(b, &n));
473ba16761SJacob Faibussowitsch PetscCall(VecGetArrayRead(b, &array));
483ba16761SJacob Faibussowitsch fprintf(file, "%" PetscInt_FMT "\n", n);
493ba16761SJacob Faibussowitsch for (PetscInt i = 0; i < n; i++) fprintf(file, "%" PetscInt_FMT " %16.14e\n", i + 1, array[i]);
503ba16761SJacob Faibussowitsch PetscCall(VecRestoreArrayRead(b, &array));
513ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
5201a81e61SBarry Smith }
53