Actual source code: ex54f.F90
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
8: program main
9: #include "petsc/finclude/petsc.h"
10: use petsc
11: implicit none
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
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./)
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
38: PetscCallA(MatAssemblyBegin(MTX,MAT_FINAL_ASSEMBLY,ierr))
39: PetscCallA(MatAssemblyEnd(MTX,MAT_FINAL_ASSEMBLY,ierr))
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))
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
52: PetscCallA(vecAssemblyBegin(PTCB,ierr))
53: PetscCallA(vecAssemblyEnd(PTCB,ierr))
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))
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))
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
75: !/*TEST
76: ! build:
77: ! requires: !complex
78: ! test:
79: ! args: -ksp_type cgls -pc_type none
80: !
81: !TEST*/