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 program main 7#include <petsc/finclude/petscksp.h> 8 use petscksp 9 implicit none 10 11 PetscInt i,n 12 PetscCount nz 13 PetscBool flg 14 PetscErrorCode ierr 15 PetscScalar,ALLOCATABLE :: a(:) 16 PetscScalar,pointer :: b(:) 17 18 PetscInt,ALLOCATABLE :: rows(:) 19 PetscInt,ALLOCATABLE :: cols(:) 20 21 Mat J 22 Vec rhs,solution 23 KSP ksp 24 25 PetscCallA(PetscInitialize(ierr)) 26 27 n = 3 28 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr)) 29 nz = 3*n - 4 30 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)) 34 35 PetscCallA(VecGetArray(rhs,b,ierr)) 36 do i=1,n 37 b(i) = 1.0 38 enddo 39 PetscCallA(VecRestoreArray(rhs,b,ierr)) 40 41 rows(1) = 0; cols(1) = 0 42 a(1) = 1.0 43 do i=2,n-1 44 rows(2+3*(i-2)) = i-1; cols(2+3*(i-2)) = i-2 45 a(2+3*(i-2)) = -1.0 46 rows(2+3*(i-2)+1) = i-1; cols(2+3*(i-2)+1) = i-1 47 a(2+3*(i-2)+1) = 2.0 48 rows(2+3*(i-2)+2) = i-1; cols(2+3*(i-2)+2) = i 49 a(2+3*(i-2)+2) = -1.0 50 enddo 51 rows(nz) = n-1; cols(nz) = n-1 52 a(nz) = 1.0 53 54 PetscCallA(MatCreate(PETSC_COMM_SELF,J,ierr)) 55 PetscCallA(MatSetSizes(J,n,n,n,n,ierr)) 56 PetscCallA(MatSetType(J,MATSEQAIJ,ierr)) 57 PetscCallA(MatSetPreallocationCOO(J,nz,rows,cols,ierr)) 58 PetscCallA(MatSetValuesCOO(J,a,INSERT_VALUES,ierr)) 59 60 PetscCallA(KSPCreate(PETSC_COMM_SELF,ksp,ierr)) 61 PetscCallA(KSPSetErrorIfNotConverged(ksp,PETSC_TRUE,ierr)) 62 PetscCallA(KSPSetFromOptions(ksp,ierr)) 63 PetscCallA(KSPSetOperators(ksp,J,J,ierr)) 64 65 PetscCallA(KSPSolve(ksp,rhs,solution,ierr)) 66 67! Keep the same size and nonzero structure of the matrix but change its numerical entries 68 do i=2,n-1 69 a(2+3*(i-2)+1) = 4.0 70 enddo 71 PetscCallA(MatSetValuesCOO(J,a,INSERT_VALUES,ierr)) 72 73 PetscCallA(KSPSolve(ksp,rhs,solution,ierr)) 74 75 PetscCallA(KSPDestroy(ksp,ierr)) 76 PetscCallA(VecDestroy(rhs,ierr)) 77 PetscCallA(VecDestroy(solution,ierr)) 78 PetscCallA(MatDestroy(J,ierr)) 79 80 DEALLOCATE (rows, cols, a) 81 82 PetscCallA(PetscFinalize(ierr)) 83 end 84 85!/*TEST 86! 87! test: 88! requires: defined(PETSC_USE_SINGLE_LIBRARY) 89! nsize: 3 90! filter: sed 's?ATOL?RTOL?g' | grep -v HERMITIAN | grep -v "shared memory" | grep -v "Mat_0" 91! # use the MPI Linear Solver Server 92! args: -n 20 -mpi_linear_solver_server -mpi_linear_solver_server_view -mpi_linear_solver_server_use_shared_memory false 93! # controls for the use of PCMPI on a particular system 94! args: -mpi_linear_solver_server_minimum_count_per_rank 5 -mpi_linear_solver_server_ksp_view 95! # the usual options for the linear solver (in this case using the server) 96! args: -ksp_monitor -ksp_converged_reason -ksp_view 97! 98! test: 99! suffix: 2 100! requires: defined(PETSC_USE_SINGLE_LIBRARY) 101! nsize: 3 102! filter: sed 's?ATOL?RTOL?g' | grep -v HERMITIAN | grep -v "shared memory" | grep -v "Mat_0" 103! # use the MPI Linear Solver Server 104! args: -n 20 -mpi_linear_solver_server -mpi_linear_solver_server_view -mpi_linear_solver_server_use_shared_memory false 105! # controls for the use of PCMPI on a particular system 106! args: -mpi_linear_solver_server_ksp_view 107! # the usual options for the linear solver (in this case using the server) 108! args: -ksp_monitor -ksp_converged_reason -ksp_view 109! 110!TEST*/ 111