Actual source code: ex6f.F
petsc-3.4.2 2013-07-02

      1: !
      2: !  Description: This example demonstrates repeated linear solves as
      3: !  well as the use of different preconditioner and linear system
      4: !  matrices.  This example also illustrates how to save PETSc objects
      5: !  in common blocks.
      6: !
      7: !/*T
      8: !  Concepts: KSP^repeatedly solving linear systems;
      9: !  Concepts: KSP^different matrices for linear system and preconditioner;
     10: !  Processors: n
     11: !T*/
     12: !
     13: !  The following include statements are required for KSP Fortran programs:
     14: !     petscsys.h       - base PETSc routines
     15: !     petscvec.h    - vectors
     16: !     petscmat.h    - matrices
     17: !     petscpc.h     - preconditioners
     18: !     petscksp.h    - Krylov subspace methods
     19: !  Other include statements may be needed if using additional PETSc
     20: !  routines in a Fortran program, e.g.,
     21: !     petscviewer.h - viewers
     22: !     petscis.h     - index sets
     23: !
     24:       program main
     25: #include <finclude/petscsys.h>
     26: #include <finclude/petscvec.h>
     27: #include <finclude/petscmat.h>
     28: #include <finclude/petscpc.h>
     29: #include <finclude/petscksp.h>

     31: !  Variables:
     32: !
     33: !  A       - matrix that defines linear system
     34: !  ksp    - KSP context
     35: !  ksp     - KSP context
     36: !  x, b, u - approx solution, RHS, exact solution vectors
     37: !
     38:       Vec     x,u,b
     39:       Mat     A
     40:       KSP    ksp
     41:       PetscInt i,j,II,JJ,m,n
     42:       PetscInt Istart,Iend
     43:       PetscInt nsteps,one
     44:       PetscErrorCode ierr
     45:       PetscBool  flg
     46:       PetscScalar  v

     49:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
     50:       m      = 3
     51:       n      = 3
     52:       nsteps = 2
     53:       one    = 1
     54:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-m',m,flg,ierr)
     55:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg,ierr)
     56:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-nsteps',nsteps,    &
     57:      &     flg,ierr)

     59: !  Create parallel matrix, specifying only its global dimensions.
     60: !  When using MatCreate(), the matrix format can be specified at
     61: !  runtime. Also, the parallel partitioning of the matrix is
     62: !  determined by PETSc at runtime.

     64:       call MatCreate(PETSC_COMM_WORLD,A,ierr)
     65:       call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,ierr)
     66:       call MatSetFromOptions(A,ierr)
     67:       call MatSetUp(A,ierr)

     69: !  The matrix is partitioned by contiguous chunks of rows across the
     70: !  processors.  Determine which rows of the matrix are locally owned.

     72:       call MatGetOwnershipRange(A,Istart,Iend,ierr)

     74: !  Set matrix elements.
     75: !   - Each processor needs to insert only elements that it owns
     76: !     locally (but any non-local elements will be sent to the
     77: !     appropriate processor during matrix assembly).
     78: !   - Always specify global rows and columns of matrix entries.

     80:       do 10, II=Istart,Iend-1
     81:         v = -1.0
     82:         i = II/n
     83:         j = II - i*n
     84:         if (i.gt.0) then
     85:           JJ = II - n
     86:           call MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr)
     87:         endif
     88:         if (i.lt.m-1) then
     89:           JJ = II + n
     90:           call MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr)
     91:         endif
     92:         if (j.gt.0) then
     93:           JJ = II - 1
     94:           call MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr)
     95:         endif
     96:         if (j.lt.n-1) then
     97:           JJ = II + 1
     98:           call MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr)
     99:         endif
    100:         v = 4.0
    101:         call  MatSetValues(A,one,II,one,II,v,ADD_VALUES,ierr)
    102:  10   continue

    104: !  Assemble matrix, using the 2-step process:
    105: !       MatAssemblyBegin(), MatAssemblyEnd()
    106: !  Computations can be done while messages are in transition
    107: !  by placing code between these two statements.

    109:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
    110:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)

    112: !  Create parallel vectors.
    113: !   - When using VecCreate(), the parallel partitioning of the vector
    114: !     is determined by PETSc at runtime.
    115: !   - Note: We form 1 vector from scratch and then duplicate as needed.

    117:       call VecCreate(PETSC_COMM_WORLD,u,ierr)
    118:       call VecSetSizes(u,PETSC_DECIDE,m*n,ierr)
    119:       call VecSetFromOptions(u,ierr)
    120:       call VecDuplicate(u,b,ierr)
    121:       call VecDuplicate(b,x,ierr)

    123: !  Create linear solver context

    125:       call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)

    127: !  Set runtime options (e.g., -ksp_type <type> -pc_type <type>)

    129:       call KSPSetFromOptions(ksp,ierr)

    131: !  Solve several linear systems in succession

    133:       do 100 i=1,nsteps
    134:          call solve1(ksp,A,x,b,u,i,nsteps,ierr)
    135:  100  continue

    137: !  Free work space.  All PETSc objects should be destroyed when they
    138: !  are no longer needed.

    140:       call VecDestroy(u,ierr)
    141:       call VecDestroy(x,ierr)
    142:       call VecDestroy(b,ierr)
    143:       call MatDestroy(A,ierr)
    144:       call KSPDestroy(ksp,ierr)

    146:       call PetscFinalize(ierr)
    147:       end

    149: ! -----------------------------------------------------------------------
    150: !
    151:       subroutine solve1(ksp,A,x,b,u,count,nsteps,ierr)

    153: #include <finclude/petscsys.h>
    154: #include <finclude/petscvec.h>
    155: #include <finclude/petscmat.h>
    156: #include <finclude/petscpc.h>
    157: #include <finclude/petscksp.h>

    159: !
    160: !   solve1 - This routine is used for repeated linear system solves.
    161: !   We update the linear system matrix each time, but retain the same
    162: !   preconditioning matrix for all linear solves.
    163: !
    164: !      A - linear system matrix
    165: !      A2 - preconditioning matrix
    166: !
    167:       PetscScalar  v,val
    168:       PetscInt II,Istart,Iend
    169:       PetscInt count,nsteps,one
    170:       PetscErrorCode ierr
    171:       Mat     A
    172:       KSP     ksp
    173:       Vec     x,b,u

    175: ! Use common block to retain matrix between successive subroutine calls
    176:       Mat              A2
    177:       PetscMPIInt      rank
    178:       PetscBool        pflag
    179:       common /my_data/ A2,pflag,rank

    181:       one = 1
    182: ! First time thorough: Create new matrix to define the linear system
    183:       if (count .eq. 1) then
    184:         call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
    185:         pflag = .false.
    186:         call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-mat_view',       &
    187:      &       pflag,ierr)
    188:         if (pflag) then
    189:           if (rank .eq. 0) write(6,100)
    190:           call flush(6)
    191:         endif
    192:         call MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,A2,ierr)
    193: ! All other times: Set previous solution as initial guess for next solve.
    194:       else
    195:         call KSPSetInitialGuessNonzero(ksp,PETSC_TRUE,ierr)
    196:       endif

    198: ! Alter the matrix A a bit
    199:       call MatGetOwnershipRange(A,Istart,Iend,ierr)
    200:       do 20, II=Istart,Iend-1
    201:         v = 2.0
    202:         call MatSetValues(A,one,II,one,II,v,ADD_VALUES,ierr)
    203:  20   continue
    204:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
    205:       if (pflag) then
    206:         if (rank .eq. 0) write(6,110)
    207:         call flush(6)
    208:       endif
    209:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)

    211: ! Set the exact solution; compute the right-hand-side vector
    212:       val = 1.0*count
    213:       call VecSet(u,val,ierr)
    214:       call MatMult(A,u,b,ierr)

    216: ! Set operators, keeping the identical preconditioner matrix for
    217: ! all linear solves.  This approach is often effective when the
    218: ! linear systems do not change very much between successive steps.
    219:       call KSPSetOperators(ksp,A,A2,SAME_PRECONDITIONER,ierr)

    221: ! Solve linear system
    222:       call KSPSolve(ksp,b,x,ierr)

    224: ! Destroy the preconditioner matrix on the last time through
    225:       if (count .eq. nsteps) call MatDestroy(A2,ierr)

    227:  100  format('previous matrix: preconditioning')
    228:  110  format('next matrix: defines linear system')

    230:       end