Actual source code: ex16f.F90

petsc-3.12.2 2019-11-22
Report Typos and Errors
  1: !
  2:       program main
  3:  #include <petsc/finclude/petscksp.h>
  4:       use petscksp
  5:       implicit none

  7: !
  8: !  This example is a modified Fortran version of ex6.c.  It tests the use of
  9: !  options prefixes in PETSc. Two linear problems are solved in this program.
 10: !  The first problem is read from a file. The second problem is constructed
 11: !  from the first, by eliminating some of the entries of the linear matrix 'A'.

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

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

 37:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 38:       if (ierr .ne. 0) then
 39:         print*,'Unable to initialize PETSc'
 40:         stop
 41:       endif

 43: ! Read in matrix and RHS
 44:       call PetscOptionsGetString(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-f',f,flg,ierr);CHKERRA(ierr)
 45:       call PetscViewerBinaryOpen(PETSC_COMM_WORLD,f,FILE_MODE_READ,fd,ierr);CHKERRA(ierr)

 47:       call MatCreate(PETSC_COMM_WORLD,A,ierr)
 48:       call MatSetType(A, MATSEQAIJ,ierr)
 49:       call MatLoad(A,fd,ierr)

 51:       call VecCreate(PETSC_COMM_WORLD,b,ierr)
 52:       call VecLoad(b,fd,ierr)
 53:       call PetscViewerDestroy(fd,ierr)

 55: ! Set up solution
 56:       call VecDuplicate(b,x,ierr)
 57:       call VecDuplicate(b,u,ierr)

 59: ! Solve system-1
 60:       call KSPCreate(PETSC_COMM_WORLD,ksp1,ierr)
 61:       call KSPSetOptionsPrefix(ksp1,'a',ierr)
 62:       call KSPAppendOptionsPrefix(ksp1,'_',ierr)
 63:       call KSPSetOperators(ksp1,A,A,ierr)
 64:       call KSPSetFromOptions(ksp1,ierr)
 65:       call KSPSolve(ksp1,b,x,ierr)

 67: ! Show result
 68:       call MatMult(A,x,u,ierr)
 69:       call VecAXPY(u,none,b,ierr)
 70:       call VecNorm(u,NORM_2,norm,ierr)
 71:       call KSPGetIterationNumber(ksp1,its,ierr)


 74:       write(6,100) norm,its
 75:   100 format('Residual norm ',e11.4,' iterations ',i5)

 77: ! Create system 2 by striping off some rows of the matrix
 78:       call ISCreateStride(PETSC_COMM_SELF,ifive,izero,ione,isrow,ierr)
 79:       call MatZeroRowsIS(A,isrow,five,PETSC_NULL_VEC,                   &
 80:      &                   PETSC_NULL_VEC,ierr)

 82: ! Solve system-2
 83:       call KSPCreate(PETSC_COMM_WORLD,ksp2,ierr)
 84:       call KSPSetOptionsPrefix(ksp2,'b',ierr)
 85:       call KSPAppendOptionsPrefix(ksp2,'_',ierr)
 86:       call KSPSetOperators(ksp2,A,A,ierr)
 87:       call KSPSetFromOptions(ksp2,ierr)
 88:       call KSPSolve(ksp2,b,x,ierr)

 90: ! Show result
 91:       call MatMult(A,x,u,ierr)
 92:       call VecAXPY(u,none,b,ierr)
 93:       call VecNorm(u,NORM_2,norm,ierr)
 94:       call KSPGetIterationNumber(ksp2,its,ierr)
 95:       write(6,100) norm,its

 97: ! Cleanup
 98:       call KSPDestroy(ksp1,ierr)
 99:       call KSPDestroy(ksp2,ierr)
100:       call VecDestroy(b,ierr)
101:       call VecDestroy(x,ierr)
102:       call VecDestroy(u,ierr)
103:       call MatDestroy(A,ierr)
104:       call ISDestroy(isrow,ierr)

106:       call PetscFinalize(ierr)
107:       end

109: !/*TEST
110: !
111: !    test:
112: !      args: -f ${DATAFILESPATH}/matrices/arco1 -options_left no
113: !      requires: datafilespath double  !complex !define(PETSC_USE_64BIT_INDICES)
114: !
115: !TEST*/