1: ! Solve the system for (x,y,z):
  2: !   x + y + z = 6
  3: !   x - y + z = 2
  4: !   x + y - z = 0
  5: !   x + y + 2*z = 9    This equation is used if DMS=4 (else set DMS=3)
  6: ! => x=1 , y=2 , z=3
  7: #include "petsc/finclude/petsc.h"
  8: program main
  9:   use petsc
 10:   implicit none

 12:   PetscInt:: IR(1), IC(1), I, J, DMS = 4 ! Set DMS=3 for a 3x3 squared system
 13:   PetscErrorCode ierr
 14:   PetscReal, parameter :: MV(12) = [1., 1., 1., 1., -1., 1., 1., 1., -1., 1., 1., 2.]
 15:   PetscReal, parameter :: B(4) = [6., 2., 0., 9.]
 16:   PetscReal :: X(3), BI(1)
 17:   Mat:: MTX
 18:   Vec:: PTCB, PTCX
 19:   KSP:: KK

 21:   PetscCallA(PetscInitialize(ierr))

 23:   PetscCallA(MatCreate(PETSC_COMM_WORLD, mtx, ierr))
 24:   PetscCallA(MatSetSizes(mtx, PETSC_DECIDE, PETSC_DECIDE, DMS, 3_PETSC_INT_KIND, ierr))
 25:   PetscCallA(MatSetFromOptions(mtx, ierr))
 26:   PetscCallA(MatSetUp(mtx, ierr))
 27:   PetscCallA(MatSetOption(mtx, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE, ierr))

 29:   do J = 1, 3
 30:     do I = 1, DMS
 31:       IR(1) = I - 1
 32:       IC(1) = J - 1
 33:       X(1) = MV(J + (I - 1)*3)
 34:       PetscCallA(MatSetValues(MTX, 1_PETSC_INT_KIND, IR, 1_PETSC_INT_KIND, IC, X, INSERT_VALUES, ierr))
 35:     end do
 36:   end do

 38:   PetscCallA(MatAssemblyBegin(MTX, MAT_FINAL_ASSEMBLY, ierr))
 39:   PetscCallA(MatAssemblyEnd(MTX, MAT_FINAL_ASSEMBLY, ierr))

 41:   X = 0.
 42:   PetscCallA(VecCreate(PETSC_COMM_WORLD, PTCB, ierr))   ! RHS vector
 43:   PetscCallA(VecSetSizes(PTCB, PETSC_DECIDE, DMS, ierr))
 44:   PetscCallA(VecSetFromOptions(PTCB, ierr))

 46:   do I = 1, DMS
 47:     IR(1) = I - 1
 48:     BI(1) = B(i)
 49:     PetscCallA(VecSetValues(PTCB, 1_PETSC_INT_KIND, IR, BI, INSERT_VALUES, ierr))
 50:   end do

 52:   PetscCallA(vecAssemblyBegin(PTCB, ierr))
 53:   PetscCallA(vecAssemblyEnd(PTCB, ierr))

 55:   PetscCallA(VecCreate(PETSC_COMM_WORLD, PTCX, ierr))   ! Solution vector
 56:   PetscCallA(VecSetSizes(PTCX, PETSC_DECIDE, 3_PETSC_INT_KIND, ierr))
 57:   PetscCallA(VecSetFromOptions(PTCX, ierr))
 58:   PetscCallA(vecAssemblyBegin(PTCX, ierr))
 59:   PetscCallA(vecAssemblyEnd(PTCX, ierr))

 61:   PetscCallA(KSPCreate(PETSC_COMM_WORLD, KK, ierr))
 62:   PetscCallA(KSPSetOperators(KK, MTX, MTX, ierr))
 63:   PetscCallA(KSPSetFromOptions(KK, ierr))
 64:   PetscCallA(KSPSetUp(KK, ierr))
 65:   PetscCallA(KSPSolve(KK, PTCB, PTCX, ierr))
 66:   PetscCallA(VecView(PTCX, PETSC_VIEWER_STDOUT_WORLD, ierr))

 68:   PetscCallA(MatDestroy(MTX, ierr))
 69:   PetscCallA(KSPDestroy(KK, ierr))
 70:   PetscCallA(VecDestroy(PTCB, ierr))
 71:   PetscCallA(VecDestroy(PTCX, ierr))
 72:   PetscCallA(PetscFinalize(ierr))
 73: end program main

 75: !/*TEST
 76: !
 77: !     build:
 78: !       requires: !complex
 79: !     test:
 80: !       args: -ksp_type cgls -pc_type none
 81: !
 82: !TEST*/