Actual source code: ex54f.F90

petsc-3.12.2 2019-11-22
Report Typos and Errors
  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
 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:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 22:       if (ierr .ne. 0) then
 23:         print*,'Unable to initialize PETSc'
 24:         stop
 25:       endif

 27:       one=1
 28:       three=3
 29:       call MatCreate(PETSC_COMM_WORLD,mtx,ierr)
 30:       call MatSetSizes(mtx,PETSC_DECIDE,PETSC_DECIDE,DMS,three,ierr)
 31:       call MatSetFromOptions(mtx,ierr)
 32:       call MatSetUp(mtx,ierr)
 33:       call MatSetOption(mtx,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,ierr)
 34:       MV=(/1.,1.,1.,1.,-1.,1.,1.,1.,-1.,1.,1.,2./)

 36:       do J=1,3
 37:          do I=1,DMS
 38:             IR(1)=I-1; IC(1)=J-1; X(1)=MV(J+(I-1)*3)
 39:             call MatSetValues(MTX,one,IR,one,IC,X,INSERT_VALUES,ierr)
 40:          end do
 41:       end do

 43:       call MatAssemblyBegin(MTX,MAT_FINAL_ASSEMBLY,ierr)
 44:       call MatAssemblyEnd(MTX,MAT_FINAL_ASSEMBLY,ierr)

 46:       X=0.; B=(/6.,2.,0.,9./)
 47:       call VecCreate(PETSC_COMM_WORLD,PTCB,ierr)   ! RHS vector
 48:       call VecSetSizes(PTCB,PETSC_DECIDE,DMS,ierr)
 49:       call VecSetFromOptions(PTCB,ierr)

 51:       do I=1,DMS
 52:          IR(1)=I-1
 53:          BI(1)=B(i)
 54:          call VecSetValues(PTCB,one,IR,BI,INSERT_VALUES,ierr)
 55:       end do

 57:       call vecAssemblyBegin(PTCB,ierr);
 58:       call vecAssemblyEnd(PTCB,ierr)

 60:       call VecCreate(PETSC_COMM_WORLD,PTCX,ierr)   ! Solution vector
 61:       call VecSetSizes(PTCX,PETSC_DECIDE,three,ierr)
 62:       call VecSetFromOptions(PTCX,ierr)
 63:       call vecAssemblyBegin(PTCX,ierr);
 64:       call vecAssemblyEnd(PTCX,ierr)

 66:       call KSPCreate(PETSC_COMM_WORLD,KK,ierr)
 67:       call KSPSetOperators(KK,MTX,MTX,ierr)
 68:       call KSPSetFromOptions(KK,ierr)
 69:       call KSPSetUp(KK,ierr);CHKERRA(ierr)
 70:       call KSPSolve(KK,PTCB,PTCX,ierr)
 71:       call VecView(PTCX,PETSC_VIEWER_STDOUT_WORLD,ierr)

 73:       call MatDestroy(MTX,ierr)
 74:       call KSPDestroy(KK,ierr)
 75:       call VecDestroy(PTCB,ierr)
 76:       call VecDestroy(PTCX,ierr)
 77:       call PetscFinalize(ierr)
 78:       end program main

 80: !/*TEST
 81: !     build:
 82: !       requires: !complex
 83: !     test:
 84: !       args: -ksp_type cgls -pc_type none
 85: !
 86: !TEST*/