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*/