Actual source code: ex16f.F90

  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*(PETSC_MAX_PATH_LEN) 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:       PetscCallA(PetscInitialize(ierr))

 39: !     Read in matrix and RHS
 40:       PetscCallA(PetscOptionsGetString(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-f',f,flg,ierr))
 41:       PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_WORLD,f,FILE_MODE_READ,fd,ierr))

 43:       PetscCallA(MatCreate(PETSC_COMM_WORLD,A,ierr))
 44:       PetscCallA(MatSetType(A, MATSEQAIJ,ierr))
 45:       PetscCallA(MatLoad(A,fd,ierr))

 47:       PetscCallA(VecCreate(PETSC_COMM_WORLD,b,ierr))
 48:       PetscCallA(VecLoad(b,fd,ierr))
 49:       PetscCallA(PetscViewerDestroy(fd,ierr))

 51: ! Set up solution
 52:       PetscCallA(VecDuplicate(b,x,ierr))
 53:       PetscCallA(VecDuplicate(b,u,ierr))

 55: ! Solve system-1
 56:       PetscCallA(KSPCreate(PETSC_COMM_WORLD,ksp1,ierr))
 57:       PetscCallA(KSPSetOptionsPrefix(ksp1,'a',ierr))
 58:       PetscCallA(KSPAppendOptionsPrefix(ksp1,'_',ierr))
 59:       PetscCallA(KSPSetOperators(ksp1,A,A,ierr))
 60:       PetscCallA(KSPSetFromOptions(ksp1,ierr))
 61:       PetscCallA(KSPSolve(ksp1,b,x,ierr))

 63: ! Show result
 64:       PetscCallA(MatMult(A,x,u,ierr))
 65:       PetscCallA(VecAXPY(u,none,b,ierr))
 66:       PetscCallA(VecNorm(u,NORM_2,norm,ierr))
 67:       PetscCallA(KSPGetIterationNumber(ksp1,its,ierr))

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

 72: ! Create system 2 by striping off some rows of the matrix
 73:       PetscCallA(ISCreateStride(PETSC_COMM_SELF,ifive,izero,ione,isrow,ierr))
 74:       PetscCallA(MatZeroRowsIS(A,isrow,five,PETSC_NULL_VEC,PETSC_NULL_VEC,ierr))

 76: ! Solve system-2
 77:       PetscCallA(KSPCreate(PETSC_COMM_WORLD,ksp2,ierr))
 78:       PetscCallA(KSPSetOptionsPrefix(ksp2,'b',ierr))
 79:       PetscCallA(KSPAppendOptionsPrefix(ksp2,'_',ierr))
 80:       PetscCallA(KSPSetOperators(ksp2,A,A,ierr))
 81:       PetscCallA(KSPSetFromOptions(ksp2,ierr))
 82:       PetscCallA(KSPSolve(ksp2,b,x,ierr))

 84: ! Show result
 85:       PetscCallA(MatMult(A,x,u,ierr))
 86:       PetscCallA(VecAXPY(u,none,b,ierr))
 87:       PetscCallA(VecNorm(u,NORM_2,norm,ierr))
 88:       PetscCallA(KSPGetIterationNumber(ksp2,its,ierr))
 89:       write(6,100) norm,its

 91: ! Cleanup
 92:       PetscCallA(KSPDestroy(ksp1,ierr))
 93:       PetscCallA(KSPDestroy(ksp2,ierr))
 94:       PetscCallA(VecDestroy(b,ierr))
 95:       PetscCallA(VecDestroy(x,ierr))
 96:       PetscCallA(VecDestroy(u,ierr))
 97:       PetscCallA(MatDestroy(A,ierr))
 98:       PetscCallA(ISDestroy(isrow,ierr))

100:       PetscCallA(PetscFinalize(ierr))
101:       end

103: !/*TEST
104: !
105: !    test:
106: !      args: -f ${DATAFILESPATH}/matrices/arco1 -options_left no
107: !      requires: datafilespath double  !complex !defined(PETSC_USE_64BIT_INDICES)
108: !
109: !TEST*/