Actual source code: dspai.c


  2: #include <petscmat.h>
  3: #include <petsc/private/petscimpl.h>

  5: /*
  6:      MatDumpSPAI - Dumps a PETSc matrix to a file in an ASCII format
  7:   suitable for the SPAI code of Stephen Barnard to solve. This routine
  8:   is simply here to allow testing of matrices directly with the SPAI
  9:   code, rather then through the PETSc interface.

 11: */
 12: PetscErrorCode MatDumpSPAI(Mat A, FILE *file)
 13: {
 14:   PetscMPIInt size;
 15:   PetscInt    n;
 16:   MPI_Comm    comm;

 18:   PetscFunctionBegin;
 21:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
 22:   PetscCallMPI(MPI_Comm_size(comm, &size));
 23:   PetscCheck(size == 1, comm, PETSC_ERR_SUP, "Only single processor dumps");
 24:   PetscCall(MatGetSize(A, &n, &n));
 25:   /* print the matrix */
 26:   fprintf(file, "%" PetscInt_FMT "\n", n);
 27:   for (PetscInt i = 0; i < n; i++) {
 28:     const PetscInt    *cols;
 29:     const PetscScalar *vals;
 30:     PetscInt           nz;

 32:     PetscCall(MatGetRow(A, i, &nz, &cols, &vals));
 33:     for (PetscInt j = 0; j < nz; j++) fprintf(file, "%" PetscInt_FMT " %d" PetscInt_FMT " %16.14e\n", i + 1, cols[j] + 1, vals[j]);
 34:     PetscCall(MatRestoreRow(A, i, &nz, &cols, &vals));
 35:   }
 36:   PetscFunctionReturn(PETSC_SUCCESS);
 37: }

 39: PetscErrorCode VecDumpSPAI(Vec b, FILE *file)
 40: {
 41:   PetscInt           n;
 42:   const PetscScalar *array;

 44:   PetscFunctionBegin;
 47:   PetscCall(VecGetSize(b, &n));
 48:   PetscCall(VecGetArrayRead(b, &array));
 49:   fprintf(file, "%" PetscInt_FMT "\n", n);
 50:   for (PetscInt i = 0; i < n; i++) fprintf(file, "%" PetscInt_FMT " %16.14e\n", i + 1, array[i]);
 51:   PetscCall(VecRestoreArrayRead(b, &array));
 52:   PetscFunctionReturn(PETSC_SUCCESS);
 53: }