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 8program 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)) 73end program main 74 75!/*TEST 76! build: 77! requires: !complex 78! test: 79! args: -ksp_type cgls -pc_type none 80! 81!TEST*/ 82