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