1: static char help[] = "Test MatNullSpaceView()/MatNullSpaceLoad() with a binary viewer.\n\n";

  3: #include <petscmat.h>

  5: int main(int argc, char **argv)
  6: {
  7:   MatNullSpace nsp, loaded;
  8:   PetscViewer  viewer;
  9:   Vec          v;
 10:   const Vec   *vecs;
 11:   PetscInt     rstart, rend, n;
 12:   PetscBool    has_cnst, equal;
 13:   PetscScalar  value;

 15:   PetscFunctionBeginUser;
 16:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 17:   PetscCall(VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, 2, &v));
 18:   PetscCall(VecGetOwnershipRange(v, &rstart, &rend));
 19:   for (PetscInt i = rstart; i < rend; i++) {
 20:     value = i ? -1.0 : 1.0;
 21:     PetscCall(VecSetValue(v, i, value, INSERT_VALUES));
 22:   }
 23:   PetscCall(VecAssemblyBegin(v));
 24:   PetscCall(VecAssemblyEnd(v));
 25:   PetscCall(VecNormalize(v, NULL));
 26:   PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 1, &v, &nsp));

 28:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "nullspace.dat", FILE_MODE_WRITE, &viewer));
 29:   PetscCall(MatNullSpaceView(nsp, viewer));
 30:   PetscCall(PetscViewerDestroy(&viewer));
 31:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "nullspace.dat", FILE_MODE_READ, &viewer));
 32:   PetscCall(MatNullSpaceLoad(viewer, &loaded));
 33:   PetscCall(PetscViewerDestroy(&viewer));

 35:   PetscCall(MatNullSpaceGetVecs(loaded, &has_cnst, &n, &vecs));
 36:   PetscCheck(has_cnst, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Loaded null space does not contain the constant vector");
 37:   PetscCheck(n == 1, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Loaded null space contains %" PetscInt_FMT " vectors instead of 1", n);
 38:   PetscCall(VecEqual(v, vecs[0], &equal));
 39:   PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Loaded null space vector differs from the original");

 41:   PetscCall(MatNullSpaceDestroy(&loaded));
 42:   PetscCall(MatNullSpaceDestroy(&nsp));
 43:   PetscCall(VecDestroy(&v));
 44:   PetscCall(PetscFinalize());
 45:   return 0;
 46: }

 48: /*TEST

 50:   testset:
 51:     output_file: output/empty.out
 52:     temporaries: nullspace.dat nullspace.dat.info

 54:     test:
 55:       suffix: 1
 56:       nsize: {{1 2}}
 57:       args: -viewer_binary_skip_header {{0 1}}

 59: TEST*/