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 PetscCallA(PetscInitialize(ierr)) 22 one=1 23 three=3 24 PetscCallA(MatCreate(PETSC_COMM_WORLD,mtx,ierr)) 25 PetscCallA(MatSetSizes(mtx,PETSC_DECIDE,PETSC_DECIDE,DMS,three,ierr)) 26 PetscCallA(MatSetFromOptions(mtx,ierr)) 27 PetscCallA(MatSetUp(mtx,ierr)) 28 PetscCallA(MatSetOption(mtx,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,ierr)) 29 MV=(/1.,1.,1.,1.,-1.,1.,1.,1.,-1.,1.,1.,2./) 30 31 do J=1,3 32 do I=1,DMS 33 IR(1)=I-1; IC(1)=J-1; X(1)=MV(J+(I-1)*3) 34 PetscCallA(MatSetValues(MTX,one,IR,one,IC,X,INSERT_VALUES,ierr)) 35 end do 36 end do 37 38 PetscCallA(MatAssemblyBegin(MTX,MAT_FINAL_ASSEMBLY,ierr)) 39 PetscCallA(MatAssemblyEnd(MTX,MAT_FINAL_ASSEMBLY,ierr)) 40 41 X=0.; B=(/6.,2.,0.,9./) 42 PetscCallA(VecCreate(PETSC_COMM_WORLD,PTCB,ierr)) ! RHS vector 43 PetscCallA(VecSetSizes(PTCB,PETSC_DECIDE,DMS,ierr)) 44 PetscCallA(VecSetFromOptions(PTCB,ierr)) 45 46 do I=1,DMS 47 IR(1)=I-1 48 BI(1)=B(i) 49 PetscCallA(VecSetValues(PTCB,one,IR,BI,INSERT_VALUES,ierr)) 50 end do 51 52 PetscCallA(vecAssemblyBegin(PTCB,ierr)) 53 PetscCallA(vecAssemblyEnd(PTCB,ierr)) 54 55 PetscCallA(VecCreate(PETSC_COMM_WORLD,PTCX,ierr)) ! Solution vector 56 PetscCallA(VecSetSizes(PTCX,PETSC_DECIDE,three,ierr)) 57 PetscCallA(VecSetFromOptions(PTCX,ierr)) 58 PetscCallA(vecAssemblyBegin(PTCX,ierr)) 59 PetscCallA(vecAssemblyEnd(PTCX,ierr)) 60 61 PetscCallA(KSPCreate(PETSC_COMM_WORLD,KK,ierr)) 62 PetscCallA(KSPSetOperators(KK,MTX,MTX,ierr)) 63 PetscCallA(KSPSetFromOptions(KK,ierr)) 64 PetscCallA(KSPSetUp(KK,ierr)) 65 PetscCallA(KSPSolve(KK,PTCB,PTCX,ierr)) 66 PetscCallA(VecView(PTCX,PETSC_VIEWER_STDOUT_WORLD,ierr)) 67 68 PetscCallA(MatDestroy(MTX,ierr)) 69 PetscCallA(KSPDestroy(KK,ierr)) 70 PetscCallA(VecDestroy(PTCB,ierr)) 71 PetscCallA(VecDestroy(PTCX,ierr)) 72 PetscCallA(PetscFinalize(ierr)) 73 end program main 74 75!/*TEST 76! build: 77! requires: !complex 78! test: 79! args: -ksp_type cgls -pc_type none 80! 81!TEST*/ 82