Actual source code: ex6f.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> ./ex6f [-help] [-m <m>] [all SLEPc options]
     11: !
     12: !  Description: Eigensystem from the Ising model for ferromagnetic materials.
     13: !  Information about the model can be found at the following
     14: !  site https://math.nist.gov/MatrixMarket/data/NEP
     15: !
     16: !  The command line options are:
     17: !    -m <m>, where <m> is the number of 2x2 blocks, i.e. matrix size N=2*m
     18: !
     19: ! ----------------------------------------------------------------------
     20: !
     21: #include <slepc/finclude/slepceps.h>

     23: module ex6fmodule
     24:   use slepceps
     25:   implicit none

     27: contains
     28:   ! -------------------------------------------------------------------
     29:   ! The actual routine for the matrix-vector product
     30:   ! See https://math.nist.gov/MatrixMarket/data/NEP/mvmisg/mvmisg.html
     31:   !
     32:   ! Computes Y(:,1:M) = op(A)*X(:,1:M)
     33:   ! where op(A) is A or A' and A is the Ising matrix.
     34:   !
     35:   !   trans   (input) integer
     36:   !           If trans = 0, compute y(:,1:M) = A*x(:,1:M)
     37:   !           If trans = 1, compute y(:,1:M) = A'*x(:,1:M)
     38:   !
     39:   !   n       (input) integer
     40:   !           The order of the matrix, it has to be an even number.
     41:   !
     42:   !   m       (input) integeR
     43:   !           The number of columns of x to multiply.
     44:   !
     45:   !   x       (input) double precision array, dimension (ldx, m)
     46:   !           x contains the matrix (vectors) x.
     47:   !
     48:   !   ldx     (input) integer
     49:   !           The leading dimension of array x, ldx >= max(1, n)
     50:   !
     51:   !   y       (output) double precision array, dimension (ldx, m)
     52:   !           contains the product of the matrix op(A) with x.
     53:   !
     54:   !   ldy     (input) integer
     55:   !           The leading dimension of array y, ldy >= max(1, n)
     56:   !
     57:   subroutine mvmisg(trans, n, m, x, ldx, y, ldy)
     58:     use petscsys
     59:     implicit none

     61:     PetscInt    :: ldy, ldx, m, n, trans
     62:     PetscScalar :: y(ldy, *), x(ldx, *)
     63:     PetscInt    :: i, k
     64:     PetscReal   :: alpha, beta, cosa, cosb, sina, sinb
     65:     PetscScalar :: temp, temp1

     67:     alpha = PETSC_PI/4
     68:     beta = PETSC_PI/4
     69:     cosa = cos(alpha)
     70:     sina = sin(alpha)
     71:     cosb = cos(beta)
     72:     sinb = sin(beta)

     74:     if (trans == 0) then

     76: !     ** Compute y(:,1:m) = A*x(:,1:m)

     78:       do k = 1, m
     79:         y(1, k) = cosb*x(1, k) - sinb*x(n, k)
     80:         do i = 2, n - 1, 2
     81:           y(i, k) = cosb*x(i, k) + sinb*x(i + 1, k)
     82:           y(i + 1, k) = -sinb*x(i, k) + cosb*x(i + 1, k)
     83:         end do
     84:         y(n, k) = sinb*x(1, k) + cosb*x(n, k)
     85:         do i = 1, n, 2
     86:           temp = cosa*y(i, k) + sina*y(i + 1, k)
     87:           y(i + 1, k) = -sina*y(i, k) + cosa*y(i + 1, k)
     88:           y(i, k) = temp
     89:         end do
     90:       end do

     92:     else if (trans == 1) then

     94: !     ** Compute y(:1:m) = A'*x(:,1:m)

     96:       do k = 1, m
     97:         do i = 1, n, 2
     98:           y(i, k) = cosa*x(i, k) - sina*x(i + 1, k)
     99:           y(i + 1, k) = sina*x(i, k) + cosa*x(i + 1, k)
    100:         end do
    101:         temp = cosb*y(1, k) + sinb*y(n, k)
    102:         do i = 2, n - 1, 2
    103:           temp1 = cosb*y(i, k) - sinb*y(i + 1, k)
    104:           y(i + 1, k) = sinb*y(i, k) + cosb*y(i + 1, k)
    105:           y(i, k) = temp1
    106:         end do
    107:         y(n, k) = -sinb*y(1, k) + cosb*y(n, k)
    108:         y(1, k) = temp
    109:       end do

    111:     end if
    112:   end subroutine

    114:   ! -------------------------------------------------------------------
    115:   ! MatMult_Ising - user provided matrix-vector multiply
    116:   !
    117:   ! Input Parameters:
    118:   !   A - matrix
    119:   !   x - input vector
    120:   !
    121:   ! Output Parameter:
    122:   !   y - output vector
    123:   !
    124:   subroutine MatMult_Ising(A, x, y, ierr)
    125:     use petscmat
    126:     implicit none

    128:     Mat                  :: A
    129:     Vec                  :: x, y
    130:     PetscInt             :: trans, N
    131:     PetscScalar, pointer :: xx(:), yy(:)
    132:     PetscErrorCode, intent(out) :: ierr

    134:     PetscCall(MatGetSize(A, N, PETSC_NULL_INTEGER, ierr))
    135:     PetscCall(VecGetArrayRead(x, xx, ierr))
    136:     PetscCall(VecGetArray(y, yy, ierr))

    138:     trans = 0
    139:     call mvmisg(trans, N, 1_PETSC_INT_KIND, xx, N, yy, N)

    141:     PetscCall(VecRestoreArrayRead(x, xx, ierr))
    142:     PetscCall(VecRestoreArray(y, yy, ierr))
    143:   end subroutine

    145:   ! -------------------------------------------------------------------
    146:   ! MatMultTranspose_Ising - user provided transpose matrix-vector multiply
    147:   !
    148:   ! Input Parameters:
    149:   !   A - matrix
    150:   !   x - input vector
    151:   !
    152:   ! Output Parameter:
    153:   !   y - output vector
    154:   !
    155:   subroutine MatMultTranspose_Ising(A, x, y, ierr)
    156:     use petscmat
    157:     implicit none

    159:     Mat                  :: A
    160:     Vec                  :: x, y
    161:     PetscInt             :: trans, N
    162:     PetscScalar, pointer :: xx(:), yy(:)
    163:     PetscErrorCode, intent(out) :: ierr

    165:     PetscCall(MatGetSize(A, N, PETSC_NULL_INTEGER, ierr))
    166:     PetscCall(VecGetArrayRead(x, xx, ierr))
    167:     PetscCall(VecGetArray(y, yy, ierr))

    169:     trans = 1
    170:     call mvmisg(trans, N, 1_PETSC_INT_KIND, xx, N, yy, N)

    172:     PetscCall(VecRestoreArrayRead(x, xx, ierr))
    173:     PetscCall(VecRestoreArray(y, yy, ierr))
    174:   end subroutine

    176: end module ex6fmodule

    178: program ex6f
    179:   use slepceps
    180:   use ex6fmodule
    181:   implicit none

    183: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    184: ! Declarations
    185: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    187:   Mat            :: A    ! operator matrix
    188:   EPS            :: eps  ! aeigenproblem solver context
    189:   EPSType        :: tname
    190:   PetscReal      :: tol
    191:   PetscInt       :: N, m
    192:   PetscInt       :: nev, maxit, its
    193:   PetscMPIInt    :: sz, rank
    194:   PetscErrorCode :: ierr
    195:   PetscBool      :: flg, terse

    197: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    198: ! Beginning of program
    199: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    201:   PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
    202: #if defined(PETSC_USE_COMPLEX)
    203:   SETERRA(PETSC_COMM_SELF, PETSC_ERR_SUP, 'This example requires real numbers')
    204: #endif
    205:   PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, sz, ierr))
    206:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
    207:   PetscCheckA(sz == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, 'This is a uniprocessor example only!')
    208:   m = 30
    209:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-m', m, flg, ierr))
    210:   N = 2*m

    212:   if (rank == 0) then
    213:     write (*, '(/a,i6,a/)') 'Ising Model Eigenproblem, m=', m, ', (N=2*m)'
    214:   end if

    216: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    217: !  Register the matrix-vector subroutine for the operator that defines
    218: !  the eigensystem, Ax=kx
    219: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    221:   PetscCallA(MatCreateShell(PETSC_COMM_WORLD, N, N, N, N, PETSC_NULL_INTEGER, A, ierr))
    222:   PetscCallA(MatShellSetOperation(A, MATOP_MULT, MatMult_Ising, ierr))
    223:   PetscCallA(MatShellSetOperation(A, MATOP_MULT_TRANSPOSE, MatMultTranspose_Ising, ierr))

    225: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    226: ! Create the eigensolver and display info
    227: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    229: ! ** Create eigensolver context
    230:   PetscCallA(EPSCreate(PETSC_COMM_WORLD, eps, ierr))

    232: ! ** Set operators. In this case, it is a standard eigenvalue problem
    233:   PetscCallA(EPSSetOperators(eps, A, PETSC_NULL_MAT, ierr))
    234:   PetscCallA(EPSSetProblemType(eps, EPS_NHEP, ierr))

    236: ! ** Set solver parameters at runtime
    237:   PetscCallA(EPSSetFromOptions(eps, ierr))

    239: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    240: ! Solve the eigensystem
    241: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    243:   PetscCallA(EPSSolve(eps, ierr))
    244:   PetscCallA(EPSGetIterationNumber(eps, its, ierr))
    245:   if (rank == 0) then
    246:     write (*, '(a,i4)') ' Number of iterations of the method: ', its
    247:   end if

    249: ! ** Optional: Get some information from the solver and display it
    250:   PetscCallA(EPSGetType(eps, tname, ierr))
    251:   if (rank == 0) then
    252:     write (*, '(a,a)') ' Solution method: ', tname
    253:   end if
    254:   PetscCallA(EPSGetDimensions(eps, nev, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
    255:   if (rank == 0) then
    256:     write (*, '(a,i2)') ' Number of requested eigenvalues:', nev
    257:   end if
    258:   PetscCallA(EPSGetTolerances(eps, tol, maxit, ierr))
    259:   if (rank == 0) then
    260:     write (*, '(a,1pe11.4,a,i6)') ' Stopping condition: tol=', tol, ', maxit=', maxit
    261:   end if

    263: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    264: ! Display solution and clean up
    265: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    267: ! ** show detailed info unless -terse option is given by user
    268:   PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-terse', terse, ierr))
    269:   if (terse) then
    270:     PetscCallA(EPSErrorView(eps, EPS_ERROR_RELATIVE, PETSC_NULL_VIEWER, ierr))
    271:   else
    272:     PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL, ierr))
    273:     PetscCallA(EPSConvergedReasonView(eps, PETSC_VIEWER_STDOUT_WORLD, ierr))
    274:     PetscCallA(EPSErrorView(eps, EPS_ERROR_RELATIVE, PETSC_VIEWER_STDOUT_WORLD, ierr))
    275:     PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD, ierr))
    276:   end if
    277:   PetscCallA(EPSDestroy(eps, ierr))
    278:   PetscCallA(MatDestroy(A, ierr))
    279:   PetscCallA(SlepcFinalize(ierr))
    280: end program ex6f

    282: !/*TEST
    283: !
    284: !   testset:
    285: !      args: -eps_max_it 1000 -eps_ncv 12 -eps_tol 1e-5 -eps_nev 4 -eps_largest_imaginary -terse
    286: !      requires: !complex
    287: !      output_file: output/ex6f_1.out
    288: !      filter: grep -v iterations | sed -e 's/-0.00000/0.00000/g'
    289: !      test:
    290: !         suffix: 1
    291: !      test:
    292: !         suffix: 1_ts
    293: !         args: -eps_two_sided
    294: !         requires: !single
    295: !
    296: !TEST*/