1: !
  2: !     Creates a tridiagonal sparse matrix explicitly in Fortran and solves a linear system with it
  3: !
  4: !     The matrix is provided in triples in a way that supports new nonzero values with the same nonzero structure
  5: !
  6: #include <petsc/finclude/petscksp.h>
  7: program main
  8:   use petscksp
  9:   implicit none

 11:   PetscInt i, n
 12:   PetscCount nz
 13:   PetscBool flg
 14:   PetscErrorCode ierr
 15:   PetscScalar, allocatable :: a(:)
 16:   PetscScalar, pointer :: b(:)

 18:   PetscInt, allocatable :: rows(:)
 19:   PetscInt, allocatable :: cols(:)

 21:   Mat J
 22:   Vec rhs, solution
 23:   KSP ksp

 25:   PetscCallA(PetscInitialize(ierr))

 27:   n = 3
 28:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
 29:   nz = 3*n - 4

 31:   PetscCallA(VecCreateSeq(PETSC_COMM_SELF, n, rhs, ierr))
 32:   PetscCallA(VecCreateSeq(PETSC_COMM_SELF, n, solution, ierr))
 33:   allocate (rows(nz), cols(nz), a(nz))

 35:   PetscCallA(VecGetArray(rhs, b, ierr))
 36:   b(1:n) = 1.0
 37:   PetscCallA(VecRestoreArray(rhs, b, ierr))

 39:   rows(1) = 0; cols(1) = 0
 40:   a(1) = 1.0
 41:   do i = 2, n - 1
 42:     rows(2 + 3*(i - 2)) = i - 1; cols(2 + 3*(i - 2)) = i - 2
 43:     a(2 + 3*(i - 2)) = -1.0
 44:     rows(2 + 3*(i - 2) + 1) = i - 1; cols(2 + 3*(i - 2) + 1) = i - 1
 45:     a(2 + 3*(i - 2) + 1) = 2.0
 46:     rows(2 + 3*(i - 2) + 2) = i - 1; cols(2 + 3*(i - 2) + 2) = i
 47:     a(2 + 3*(i - 2) + 2) = -1.0
 48:   end do
 49:   rows(nz) = n - 1; cols(nz) = n - 1
 50:   a(nz) = 1.0

 52:   PetscCallA(MatCreate(PETSC_COMM_SELF, J, ierr))
 53:   PetscCallA(MatSetSizes(J, n, n, n, n, ierr))
 54:   PetscCallA(MatSetType(J, MATSEQAIJ, ierr))
 55:   PetscCallA(MatSetPreallocationCOO(J, nz, rows, cols, ierr))
 56:   PetscCallA(MatSetValuesCOO(J, a, INSERT_VALUES, ierr))

 58:   PetscCallA(KSPCreate(PETSC_COMM_SELF, ksp, ierr))
 59:   PetscCallA(KSPSetErrorIfNotConverged(ksp, PETSC_TRUE, ierr))
 60:   PetscCallA(KSPSetFromOptions(ksp, ierr))
 61:   PetscCallA(KSPSetOperators(ksp, J, J, ierr))

 63:   PetscCallA(KSPSolve(ksp, rhs, solution, ierr))

 65: ! Keep the same size and nonzero structure of the matrix but change its numerical entries
 66:   do i = 2, n - 1
 67:     a(2 + 3*(i - 2) + 1) = 4.0
 68:   end do
 69:   PetscCallA(MatSetValuesCOO(J, a, INSERT_VALUES, ierr))

 71:   PetscCallA(KSPSolve(ksp, rhs, solution, ierr))

 73:   PetscCallA(KSPDestroy(ksp, ierr))
 74:   PetscCallA(VecDestroy(rhs, ierr))
 75:   PetscCallA(VecDestroy(solution, ierr))
 76:   PetscCallA(MatDestroy(J, ierr))

 78:   deallocate (rows, cols, a)

 80:   PetscCallA(PetscFinalize(ierr))
 81: end

 83: !/*TEST
 84: !
 85: !     test:
 86: !       requires: defined(PETSC_USE_SINGLE_LIBRARY)
 87: !       nsize: 3
 88: !       filter: sed 's?ATOL?RTOL?g' | grep -v HERMITIAN | grep -v "shared memory" | grep -v "Mat_0"
 89: !       # use the MPI Linear Solver Server
 90: !       args: -n 20 -mpi_linear_solver_server -mpi_linear_solver_server_view -mpi_linear_solver_server_use_shared_memory false
 91: !       # controls for the use of PCMPI on a particular system
 92: !       args: -mpi_linear_solver_server_minimum_count_per_rank 5 -mpi_linear_solver_server_ksp_view
 93: !       # the usual options for the linear solver (in this case using the server)
 94: !       args: -ksp_monitor -ksp_converged_reason -ksp_view
 95: !
 96: !     test:
 97: !       suffix: 2
 98: !       requires: defined(PETSC_USE_SINGLE_LIBRARY)
 99: !       nsize: 3
100: !       filter: sed 's?ATOL?RTOL?g' | grep -v HERMITIAN | grep -v "shared memory" | grep -v "Mat_0"
101: !       # use the MPI Linear Solver Server
102: !       args: -n 20 -mpi_linear_solver_server -mpi_linear_solver_server_view -mpi_linear_solver_server_use_shared_memory false
103: !       # controls for the use of PCMPI on a particular system
104: !       args: -mpi_linear_solver_server_ksp_view
105: !       # the usual options for the linear solver (in this case using the server)
106: !       args: -ksp_monitor -ksp_converged_reason -ksp_view
107: !
108: !TEST*/