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> ./test14f [-help] [-n <n>] [all SLEPc options]
 11: !
 12: !  Description: Simple example that tests solving a DSNHEP problem.
 13: !
 14: !  The command line options are:
 15: !    -n <n>, where <n> = matrix size
 16: !
 17: ! ----------------------------------------------------------------------
 18: !
 19: #include <slepc/finclude/slepcds.h>
 20: program test14f
 21:   use slepcds
 22:   implicit none

 24: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 25: ! Declarations
 26: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 28:   Mat                  :: A   ! problem matrix
 29:   DS                   :: ds  ! dense solver context
 30:   PetscInt             :: n, i, ld
 31:   PetscMPIInt          :: rank
 32:   PetscErrorCode       :: ierr
 33:   PetscBool            :: flg
 34:   PetscScalar          :: wr(100), wi(100)
 35:   PetscReal            :: re, im
 36:   PetscScalar, pointer :: aa(:, :)

 38: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 39: ! Beginning of program
 40: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 42:   PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
 43:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
 44:   n = 10
 45:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
 46:   PetscCheckA(n <= 100, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, 'Program currently limited to n=100')

 48:   if (rank == 0) then
 49:     write (*, '(/a,i3,a)') 'Solve a Dense System of type NHEP, n =', n, ' (Fortran)'
 50:   end if

 52: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 53: ! Create DS object
 54: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 56:   PetscCallA(DSCreate(PETSC_COMM_WORLD, ds, ierr))
 57:   PetscCallA(DSSetType(ds, DSNHEP, ierr))
 58:   PetscCallA(DSSetFromOptions(ds, ierr))
 59:   ld = n
 60:   PetscCallA(DSAllocate(ds, ld, ierr))
 61:   PetscCallA(DSSetDimensions(ds, n, 0_PETSC_INT_KIND, 0_PETSC_INT_KIND, ierr))

 63: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 64: ! Fill with Grcar matrix
 65: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 67:   PetscCallA(DSGetMat(ds, DS_MAT_A, A, ierr))
 68:   PetscCallA(MatDenseGetArray(A, aa, ierr))
 69:   call FillUpMatrix(n, aa)
 70:   PetscCallA(MatDenseRestoreArray(A, aa, ierr))
 71:   PetscCallA(DSRestoreMat(ds, DS_MAT_A, A, ierr))
 72:   PetscCallA(DSSetState(ds, DS_STATE_INTERMEDIATE, ierr))

 74: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 75: ! Solve the problem and show eigenvalues
 76: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 78:   PetscCallA(DSSolve(ds, wr, wi, ierr))
 79: ! PetscCallA(DSSort(ds,wr,wi,PETSC_NULL_SCALAR,PETSC_NULL_SCALAR,PETSC_NULL_INTEGER,ierr))

 81:   if (rank == 0) then
 82:     write (*, *) 'Computed eigenvalues ='
 83:     do i = 1, n
 84: #if defined(PETSC_USE_COMPLEX)
 85:       re = PetscRealPart(wr(i))
 86:       im = PetscImaginaryPart(wr(i))
 87: #else
 88:       re = wr(i)
 89:       im = wi(i)
 90: #endif
 91:       if (abs(im) < 1.d-10) then
 92:         write (*, '(a,f8.5)') '  ', re
 93:       else
 94:         write (*, '(a,f8.5,sp,f8.5,a)') '  ', re, im, 'i'
 95:       end if
 96:     end do
 97:   end if

 99: ! *** Clean up
100:   PetscCallA(DSDestroy(ds, ierr))
101:   PetscCallA(SlepcFinalize(ierr))

103: contains

105:   subroutine FillUpMatrix(n, X)
106:     PetscInt    :: n, i, j
107:     PetscScalar :: X(n, n)

109:     do i = 2, n
110:       X(i, i - 1) = -1.d0
111:     end do
112:     do j = 0, 3
113:       do i = 1, n - j
114:         X(i, i + j) = 1.d0
115:       end do
116:     end do

118:   end

120: end program test14f

122: !/*TEST
123: !
124: !   test:
125: !      suffix: 1
126: !      requires: !complex
127: !
128: !TEST*/