Actual source code: ex16f.F

  1: !
  2:       program main
  3:        implicit none

 5:  #include finclude/petscsys.h
 6:  #include finclude/petscvec.h
 7:  #include finclude/petscmat.h
 8:  #include finclude/petscpc.h
 9:  #include finclude/petscksp.h
 10:  #include finclude/petscviewer.h
 11:  #include finclude/petscis.h
 12: !
 13: !  This example is a modified Fortran version of ex6.c.  It tests the use of
 14: !  options prefixes in PETSc. Two linear problems are solved in this program.
 15: !  The first problem is read from a file. The second problem is constructed
 16: !  from the first, by eliminating some of the entries of the linear matrix 'A'.

 18: !  Each solve is distinguished by a unique prefix - 'a' for the first, 'b'
 19: !  for the second.  With the prefix the user can distinguish between the various
 20: !  options (command line, from .petscrc file, etc.) for each of the solvers.
 21: !  Input arguments are:
 22: !     -f <input_file> : file to load.  For a 5X5 example of the 5-pt. stencil
 23: !                       use the file petsc/src/mat/examples/mat.ex.binary

 25:       PetscErrorCode  ierr
 26:       PetscInt its,ione,ifive,izero
 27:       PetscErrorCode flg
 28:       PetscScalar      norm,none,five
 29:       Vec              x,b,u
 30:       Mat              A
 31:       KSP             ksp1,ksp2
 32:       character*(128)  f
 33:       PetscViewer      fd
 34:       IS               isrow
 35:       none  = -1.0
 36:       five  = 5.0
 37:       ifive = 5
 38:       ione  = 1
 39:       izero = 0

 41:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 42: 
 43: ! Read in matrix and RHS
 44:       call PetscOptionsGetString(PETSC_NULL_CHARACTER,'-f',f,flg,ierr)
 45:       call PetscViewerBinaryOpen(PETSC_COMM_WORLD,f,FILE_MODE_READ,            &
 46:      &     fd,ierr)
 47:       if (ierr .ne. 0) then
 48:         print*, 'Unable to open file ',f
 49:         SETERRQ(1,' ',ierr)
 50:       endif

 52:       call MatLoad(fd,MATSEQAIJ,A,ierr)
 53:       if (ierr .ne. 0) then
 54:         print*, 'Unable to load matrix '
 55:         SETERRQ(1,' ',ierr)
 56:       endif

 58:       call VecLoad(fd,PETSC_NULL_CHARACTER,b,ierr)
 59:       call PetscViewerDestroy(fd,ierr)

 61: ! Set up solution
 62:       call VecDuplicate(b,x,ierr)
 63:       call VecDuplicate(b,u,ierr)

 65: ! Solve system-1
 66:       call KSPCreate(PETSC_COMM_WORLD,ksp1,ierr)
 67:       call KSPSetOptionsPrefix(ksp1,'a',ierr)
 68:       call KSPAppendOptionsPrefix(ksp1,'_',ierr)
 69:       call KSPSetOperators(ksp1,A,A,DIFFERENT_NONZERO_PATTERN,ierr)
 70:       call KSPSetFromOptions(ksp1,ierr)
 71:       call KSPSolve(ksp1,b,x,ierr)

 73: ! Show result
 74:       call MatMult(A,x,u,ierr)
 75:       call VecAXPY(u,none,b,ierr)
 76:       call VecNorm(u,NORM_2,norm,ierr)
 77:       call KSPGetIterationNumber(ksp1,its,ierr)


 80:       write(6,100) norm,its
 81:   100 format('Residual norm ',e10.4,' iterations ',i5)

 83: ! Create system 2 by striping off some rows of the matrix
 84:       call ISCreateStride(PETSC_COMM_SELF,ifive,izero,ione,isrow,ierr)
 85:       call MatZeroRowsIS(A,isrow,five,ierr)

 87: ! Solve system-2
 88:       call KSPCreate(PETSC_COMM_WORLD,ksp2,ierr)
 89:       call KSPSetOptionsPrefix(ksp2,'b',ierr)
 90:       call KSPAppendOptionsPrefix(ksp2,'_',ierr)
 91:       call KSPSetOperators(ksp2,A,A,DIFFERENT_NONZERO_PATTERN,ierr)
 92:       call KSPSetFromOptions(ksp2,ierr)
 93:       call KSPSolve(ksp2,b,x,ierr)

 95: ! Show result
 96:       call MatMult(A,x,u,ierr)
 97:       call VecAXPY(u,none,b,ierr)
 98:       call VecNorm(u,NORM_2,norm,ierr)
 99:       call KSPGetIterationNumber(ksp2,its,ierr)
100:       write(6,100) norm,its

102: ! Cleanup
103:       call KSPDestroy(ksp1,ierr)
104:       call KSPDestroy(ksp2,ierr)
105:       call VecDestroy(b,ierr)
106:       call VecDestroy(x,ierr)
107:       call VecDestroy(u,ierr)
108:       call MatDestroy(A,ierr)
109:       call ISDestroy(isrow,ierr)

111:       call PetscFinalize(ierr)
112:       end