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

  PetscInt:: IR(1), IC(1), I, J, DMS = 4 ! Set DMS=3 for a 3x3 squared system
  PetscErrorCode ierr
  PetscReal :: MV(12), X(3), B(4), BI(1)
  Mat:: MTX
  Vec:: PTCB, PTCX
  KSP:: KK
  PetscInt one, three

  PetscCallA(PetscInitialize(ierr))
  one = 1
  three = 3
  PetscCallA(MatCreate(PETSC_COMM_WORLD, mtx, ierr))
  PetscCallA(MatSetSizes(mtx, PETSC_DECIDE, PETSC_DECIDE, DMS, three, ierr))
  PetscCallA(MatSetFromOptions(mtx, ierr))
  PetscCallA(MatSetUp(mtx, ierr))
  PetscCallA(MatSetOption(mtx, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE, ierr))
  MV = (/1., 1., 1., 1., -1., 1., 1., 1., -1., 1., 1., 2./)

  do J = 1, 3
    do I = 1, DMS
      IR(1) = I - 1; IC(1) = J - 1; X(1) = MV(J + (I - 1)*3)
      PetscCallA(MatSetValues(MTX, one, IR, one, IC, X, INSERT_VALUES, ierr))
    end do
  end do

  PetscCallA(MatAssemblyBegin(MTX, MAT_FINAL_ASSEMBLY, ierr))
  PetscCallA(MatAssemblyEnd(MTX, MAT_FINAL_ASSEMBLY, ierr))

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

  do I = 1, DMS
    IR(1) = I - 1
    BI(1) = B(i)
    PetscCallA(VecSetValues(PTCB, one, IR, BI, INSERT_VALUES, ierr))
  end do

  PetscCallA(vecAssemblyBegin(PTCB, ierr))
  PetscCallA(vecAssemblyEnd(PTCB, ierr))

  PetscCallA(VecCreate(PETSC_COMM_WORLD, PTCX, ierr))   ! Solution vector
  PetscCallA(VecSetSizes(PTCX, PETSC_DECIDE, three, ierr))
  PetscCallA(VecSetFromOptions(PTCX, ierr))
  PetscCallA(vecAssemblyBegin(PTCX, ierr))
  PetscCallA(vecAssemblyEnd(PTCX, ierr))

  PetscCallA(KSPCreate(PETSC_COMM_WORLD, KK, ierr))
  PetscCallA(KSPSetOperators(KK, MTX, MTX, ierr))
  PetscCallA(KSPSetFromOptions(KK, ierr))
  PetscCallA(KSPSetUp(KK, ierr))
  PetscCallA(KSPSolve(KK, PTCB, PTCX, ierr))
  PetscCallA(VecView(PTCX, PETSC_VIEWER_STDOUT_WORLD, ierr))

  PetscCallA(MatDestroy(MTX, ierr))
  PetscCallA(KSPDestroy(KK, ierr))
  PetscCallA(VecDestroy(PTCB, ierr))
  PetscCallA(VecDestroy(PTCX, ierr))
  PetscCallA(PetscFinalize(ierr))
end program main

!/*TEST
!     build:
!       requires: !complex
!     test:
!       args: -ksp_type cgls -pc_type none
!
!TEST*/