1: !
  2: !  Tests MatGetNullSpace() and PetscObjectIsNull()
  3: !
  4: #include <petsc/finclude/petscmat.h>
  5: program main

  7:   use petscmat
  8:   implicit none

 10:   PetscErrorCode ierr
 11:   Mat A, B
 12:   Mat C, SC
 13:   MatNullSpace sp, sp1
 14:   PetscInt rend
 15:   PetscScalar, parameter :: sone = 1.0
 16:   Vec x, y

 18:   PetscCallA(PetscInitialize(ierr))

 20:   PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
 21:   PetscCallA(MatCreate(PETSC_COMM_WORLD, B, ierr))

 23:   PetscCallA(MatGetNullSpace(A, sp, ierr))
 24:   PetscCheckA(PetscObjectIsNull(sp), PETSC_COMM_SELF, PETSC_ERR_PLIB, 'Matrix null space should not exist')

 26:   PetscCallA(MatSetNullSpace(A, PETSC_NULL_MATNULLSPACE, ierr))
 27:   PetscCallA(MatGetNullSpace(A, sp, ierr))
 28:   PetscCheckA(PetscObjectIsNull(sp), PETSC_COMM_SELF, PETSC_ERR_PLIB, 'Matrix null space should not exist')

 30:   PetscCallA(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0_PETSC_INT_KIND, PETSC_NULL_VEC_ARRAY, sp, ierr))
 31:   PetscCallA(MatSetNullSpace(A, sp, ierr))
 32:   PetscCallA(MatGetNullSpace(A, sp1, ierr))
 33:   PetscCheckA(.not. PetscObjectIsNull(sp), PETSC_COMM_SELF, PETSC_ERR_PLIB, 'Matrix null space should not exist')
 34:   PetscCallA(MatNullSpaceDestroy(sp, ierr))

 36:   PetscCallA(MatCreateSeqDense(PETSC_COMM_WORLD, 1_PETSC_INT_KIND, 1_PETSC_INT_KIND, PETSC_NULL_SCALAR_ARRAY, C, ierr))
 37:   PetscCallA(MatSetValues(C, 1_PETSC_INT_KIND, [0_PETSC_INT_KIND], 1_PETSC_INT_KIND, [0_PETSC_INT_KIND], [sone], INSERT_VALUES, ierr))
 38:   PetscCallA(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY, ierr))
 39:   PetscCallA(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY, ierr))
 40:   PetscCallA(MatCreateSchurComplement(C, C, C, C, PETSC_NULL_MAT, SC, ierr))
 41:   PetscCallA(MatGetOwnershipRange(SC, PETSC_NULL_INTEGER, rend, ierr))
 42:   PetscCallA(VecCreateSeq(PETSC_COMM_SELF, 1_PETSC_INT_KIND, x, ierr))
 43:   PetscCallA(VecDuplicate(x, y, ierr))
 44:   PetscCallA(VecSetValues(x, 1_PETSC_INT_KIND, [0_PETSC_INT_KIND], [sone], INSERT_VALUES, ierr))
 45:   PetscCallA(VecAssemblyBegin(x, ierr))
 46:   PetscCallA(VecAssemblyEnd(x, ierr))
 47:   PetscCallA(MatMult(SC, x, y, ierr))
 48:   PetscCallA(VecView(y, PETSC_VIEWER_STDOUT_SELF, ierr))
 49:   PetscCallA(VecSetRandom(x, PETSC_NULL_RANDOM, ierr))
 50:   PetscCallA(VecView(x, PETSC_VIEWER_STDOUT_SELF, ierr))

 52:   PetscCallA(MatDestroy(SC, ierr))
 53:   PetscCallA(MatDestroy(C, ierr))
 54:   PetscCallA(VecDestroy(x, ierr))
 55:   PetscCallA(VecDestroy(y, ierr))
 56:   PetscCallA(MatDestroy(A, ierr))
 57:   PetscCallA(MatDestroy(B, ierr))
 58:   PetscCallA(PetscFinalize(ierr))
 59: end

 61: !/*TEST
 62: !
 63: !   test:
 64: !      requires: !complex
 65: !
 66: !TEST*/