Actual source code: ex23f.F90

      1: !
      2: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      3: !  SLEPc - Scalable Library for Eigenvalue Problem Computations
      4: !  Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
      5: !
      6: !  This file is part of SLEPc.
      7: !  SLEPc is distributed under a 2-clause BSD license (see LICENSE).
      8: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      9: !
     10: !  Program usage: mpiexec -n <np> ./ex23f [-help] [-t <t>] [-m <m>] [SLEPc opts]
     11: !
     12: !  Description: Computes exp(t*A)*v for a matrix from a Markov model.
     13: !  This is the Fortran90 equivalent to ex23.c
     14: !
     15: !  The command line options are:
     16: !    -t <t>, where <t> = time parameter (multiplies the matrix)
     17: !    -m <m>, where <m> = number of grid subdivisions in each dimension
     18: !
     19: ! ----------------------------------------------------------------------
     20: !
     21: #include <slepc/finclude/slepcmfn.h>
     22: program ex23f
     23:   use slepcmfn
     24:   implicit none

     26: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     27: ! Declarations
     28: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     30:   Mat                :: A      ! problem matrix
     31:   MFN                :: mfn    ! matrix function solver context
     32:   FN                 :: f
     33:   PetscReal          :: tol, norm, cst, pd, pu
     34:   PetscScalar        :: t, z
     35:   Vec                :: v, y
     36:   PetscInt           :: N, m, ncv, maxit, its, ii, jj
     37:   PetscInt           :: i, j, jmax, ix, Istart, Iend
     38:   PetscMPIInt        :: rank
     39:   PetscErrorCode     :: ierr
     40:   PetscBool          :: flg
     41:   MFNConvergedReason :: reason
     42:   PetscScalar, parameter :: one = 1.0

     44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     45: ! Beginning of program
     46: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     48:   PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
     49:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
     50:   m = 15
     51:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-m', m, flg, ierr))
     52:   t = 2.0
     53:   PetscCallA(PetscOptionsGetScalar(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-t', t, flg, ierr))
     54:   N = m*(m + 1)/2
     55:   if (rank == 0) then
     56:     write (*, '(/a,i6,a,i4,a)') 'Markov y=exp(t*A)*e_1, N=', N, ' (m=', m, ')'
     57:   end if

     59: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     60: ! Compute the transition probability matrix, A
     61: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     63:   PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
     64:   PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, N, N, ierr))
     65:   PetscCallA(MatSetFromOptions(A, ierr))
     66:   PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))
     67:   ix = 0
     68:   cst = 0.5/real(m - 1, PETSC_REAL_KIND)
     69:   do i = 1, m
     70:     jmax = m - i + 1
     71:     do j = 1, jmax
     72:       ix = ix + 1
     73:       ii = ix - 1
     74:       if (ix - 1 >= Istart .and. ix <= Iend) then
     75:         if (j /= jmax) then
     76:           pd = cst*(i + j - 1)
     77:           !** north
     78:           if (i == 1) then
     79:             z = 2.0*pd
     80:           else
     81:             z = pd
     82:           end if
     83:           PetscCallA(MatSetValue(A, ii, ix, z, INSERT_VALUES, ierr))
     84:           !** east
     85:           if (j == 1) then
     86:             z = 2.0*pd
     87:           else
     88:             z = pd
     89:           end if
     90:           jj = ix + jmax - 1
     91:           PetscCallA(MatSetValue(A, ii, jj, z, INSERT_VALUES, ierr))
     92:         end if

     94:         !** south
     95:         pu = 0.5 - cst*(i + j - 3)
     96:         z = pu
     97:         if (j > 1) then
     98:           jj = ix - 2
     99:           PetscCallA(MatSetValue(A, ii, jj, z, INSERT_VALUES, ierr))
    100:         end if
    101:         !** west
    102:         if (i > 1) then
    103:           jj = ix - jmax - 2
    104:           PetscCallA(MatSetValue(A, ii, jj, z, INSERT_VALUES, ierr))
    105:         end if
    106:       end if
    107:     end do
    108:   end do
    109:   PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
    110:   PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))

    112: ! ** Set v = e_1
    113:   PetscCallA(MatCreateVecs(A, y, v, ierr))
    114:   ii = 0
    115:   z = 1.0
    116:   PetscCallA(VecSetValue(v, ii, z, INSERT_VALUES, ierr))
    117:   PetscCallA(VecAssemblyBegin(v, ierr))
    118:   PetscCallA(VecAssemblyEnd(v, ierr))

    120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    121: ! Create the solver and set various options
    122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    124: ! ** Create matrix function solver context
    125:   PetscCallA(MFNCreate(PETSC_COMM_WORLD, mfn, ierr))

    127: ! ** Set operator matrix, the function to compute, and other options
    128:   PetscCallA(MFNSetOperator(mfn, A, ierr))
    129:   PetscCallA(MFNGetFN(mfn, f, ierr))
    130:   PetscCallA(FNSetType(f, FNEXP, ierr))
    131:   PetscCallA(FNSetScale(f, t, one, ierr))
    132:   tol = 0.0000001
    133:   PetscCallA(MFNSetTolerances(mfn, tol, PETSC_CURRENT_INTEGER, ierr))

    135: ! ** Set solver parameters at runtime
    136:   PetscCallA(MFNSetFromOptions(mfn, ierr))

    138: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    139: ! Solve the problem, y=exp(t*A)*v
    140: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    142:   PetscCallA(MFNSolve(mfn, v, y, ierr))
    143:   PetscCallA(MFNGetConvergedReason(mfn, reason, ierr))
    144:   PetscCheckA(reason%v >= 0, PETSC_COMM_WORLD, PETSC_ERR_NOT_CONVERGED, 'Solver did not converge')
    145:   PetscCallA(VecNorm(y, NORM_2, norm, ierr))

    147: ! ** Optional: Get some information from the solver and display it
    148:   PetscCallA(MFNGetIterationNumber(mfn, its, ierr))
    149:   if (rank == 0) then
    150:     write (*, '(a,i4)') ' Number of iterations of the method: ', its
    151:   end if
    152:   PetscCallA(MFNGetDimensions(mfn, ncv, ierr))
    153:   if (rank == 0) then
    154:     write (*, '(a,i4)') ' Subspace dimension:', ncv
    155:   end if
    156:   PetscCallA(MFNGetTolerances(mfn, tol, maxit, ierr))
    157:   if (rank == 0) then
    158:     write (*, '(a,f10.7,a,i4)') ' Stopping condition: tol=', tol, ' maxit=', maxit
    159:   end if

    161: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    162: ! Display solution and clean up
    163: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    165:   if (rank == 0) then
    166:     write (*, '(a,f4.1,a,f8.5)') ' Computed vector at time t=', PetscRealPart(t), ' has norm ', norm
    167:   end if

    169:   PetscCallA(MFNDestroy(mfn, ierr))
    170:   PetscCallA(MatDestroy(A, ierr))
    171:   PetscCallA(VecDestroy(v, ierr))
    172:   PetscCallA(VecDestroy(y, ierr))
    173:   PetscCallA(SlepcFinalize(ierr))
    174: end program ex23f

    176: !/*TEST
    177: !
    178: !   test:
    179: !      suffix: 1
    180: !      args: -mfn_ncv 6
    181: !
    182: !TEST*/