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 :: MV(12), X(3), B(4), BI(1)
 15:   Mat:: MTX
 16:   Vec:: PTCB, PTCX
 17:   KSP:: KK
 18:   PetscInt one, three

 20:   PetscCallA(PetscInitialize(ierr))
 21:   one = 1
 22:   three = 3
 23:   PetscCallA(MatCreate(PETSC_COMM_WORLD, mtx, ierr))
 24:   PetscCallA(MatSetSizes(mtx, PETSC_DECIDE, PETSC_DECIDE, DMS, three, ierr))
 25:   PetscCallA(MatSetFromOptions(mtx, ierr))
 26:   PetscCallA(MatSetUp(mtx, ierr))
 27:   PetscCallA(MatSetOption(mtx, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE, ierr))
 28:   MV = (/1., 1., 1., 1., -1., 1., 1., 1., -1., 1., 1., 2./)

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

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

 40:   X = 0.; B = (/6., 2., 0., 9./)
 41:   PetscCallA(VecCreate(PETSC_COMM_WORLD, PTCB, ierr))   ! RHS vector
 42:   PetscCallA(VecSetSizes(PTCB, PETSC_DECIDE, DMS, ierr))
 43:   PetscCallA(VecSetFromOptions(PTCB, ierr))

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

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

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

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

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

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