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 8 program main 9#include "petsc/finclude/petsc.h" 10 use petsc 11 implicit none 12 13 PetscInt:: IR(1),IC(1),I,J,DMS=4 ! Set DMS=3 for a 3x3 squared system 14 PetscErrorCode ierr; 15 PetscReal :: MV(12),X(3),B(4),BI(1) 16 Mat:: MTX 17 Vec:: PTCB,PTCX 18 KSP:: KK 19 PetscInt one,three 20 21 call PetscInitialize(PETSC_NULL_CHARACTER,ierr) 22 if (ierr .ne. 0) then 23 print*,'Unable to initialize PETSc' 24 stop 25 endif 26 27 one=1 28 three=3 29 call MatCreate(PETSC_COMM_WORLD,mtx,ierr) 30 call MatSetSizes(mtx,PETSC_DECIDE,PETSC_DECIDE,DMS,three,ierr) 31 call MatSetFromOptions(mtx,ierr) 32 call MatSetUp(mtx,ierr) 33 call MatSetOption(mtx,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,ierr) 34 MV=(/1.,1.,1.,1.,-1.,1.,1.,1.,-1.,1.,1.,2./) 35 36 do J=1,3 37 do I=1,DMS 38 IR(1)=I-1; IC(1)=J-1; X(1)=MV(J+(I-1)*3) 39 call MatSetValues(MTX,one,IR,one,IC,X,INSERT_VALUES,ierr) 40 end do 41 end do 42 43 call MatAssemblyBegin(MTX,MAT_FINAL_ASSEMBLY,ierr) 44 call MatAssemblyEnd(MTX,MAT_FINAL_ASSEMBLY,ierr) 45 46 X=0.; B=(/6.,2.,0.,9./) 47 call VecCreate(PETSC_COMM_WORLD,PTCB,ierr) ! RHS vector 48 call VecSetSizes(PTCB,PETSC_DECIDE,DMS,ierr) 49 call VecSetFromOptions(PTCB,ierr) 50 51 do I=1,DMS 52 IR(1)=I-1 53 BI(1)=B(i) 54 call VecSetValues(PTCB,one,IR,BI,INSERT_VALUES,ierr) 55 end do 56 57 call vecAssemblyBegin(PTCB,ierr); 58 call vecAssemblyEnd(PTCB,ierr) 59 60 call VecCreate(PETSC_COMM_WORLD,PTCX,ierr) ! Solution vector 61 call VecSetSizes(PTCX,PETSC_DECIDE,three,ierr) 62 call VecSetFromOptions(PTCX,ierr) 63 call vecAssemblyBegin(PTCX,ierr); 64 call vecAssemblyEnd(PTCX,ierr) 65 66 call KSPCreate(PETSC_COMM_WORLD,KK,ierr) 67 call KSPSetOperators(KK,MTX,MTX,ierr) 68 call KSPSetFromOptions(KK,ierr) 69 call KSPSetUp(KK,ierr);CHKERRA(ierr) 70 call KSPSolve(KK,PTCB,PTCX,ierr) 71 call VecView(PTCX,PETSC_VIEWER_STDOUT_WORLD,ierr) 72 73 call MatDestroy(MTX,ierr) 74 call KSPDestroy(KK,ierr) 75 call VecDestroy(PTCB,ierr) 76 call VecDestroy(PTCX,ierr) 77 call PetscFinalize(ierr) 78 end program main 79 80!/*TEST 81! build: 82! requires: !complex 83! test: 84! args: -ksp_type cgls -pc_type none 85! 86!TEST*/ 87