Actual source code: test7f.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> ./test7f [-help] [-n <n>] [-verbose] [-inplace]
     11: !
     12: !  Description: Simple example that tests the matrix square root.
     13: !
     14: ! ----------------------------------------------------------------------
     15: !
     16: #include <slepc/finclude/slepcfn.h>
     17: program test7f
     18:   use slepcfn
     19:   implicit none

     21: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     22: ! Declarations
     23: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     25:   Mat                  :: A, S, R
     26:   FN                   :: fn
     27:   PetscInt             :: n
     28:   PetscMPIInt          :: rank
     29:   PetscErrorCode       :: ierr
     30:   PetscBool            :: flg, verbose, inplace
     31:   PetscReal            :: re, im, nrm
     32:   PetscScalar          :: tau, eta, alpha, x, y, yp
     33:   PetscScalar, pointer :: aa(:, :)

     35: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     36: ! Beginning of program
     37: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     39:   PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
     40:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
     41:   n = 10
     42:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
     43:   PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-verbose', verbose, ierr))
     44:   PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-inplace', inplace, ierr))

     46:   if (rank == 0) then
     47:     write (*, '(/a,i3,a)') 'Matrix square root, n =', n, ' (Fortran)'
     48:   end if

     50: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     51: ! Create FN object and matrix
     52: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     54:   PetscCallA(FNCreate(PETSC_COMM_WORLD, fn, ierr))
     55:   PetscCallA(FNSetType(fn, FNSQRT, ierr))
     56:   tau = 0.15
     57:   eta = 1.0
     58:   PetscCallA(FNSetScale(fn, tau, eta, ierr))
     59:   PetscCallA(FNSetFromOptions(fn, ierr))
     60:   PetscCallA(FNGetScale(fn, tau, eta, ierr))
     61:   PetscCallA(FNView(fn, PETSC_NULL_VIEWER, ierr))

     63:   PetscCallA(MatCreateSeqDense(PETSC_COMM_SELF, n, n, PETSC_NULL_SCALAR_ARRAY, A, ierr))
     64:   PetscCallA(PetscObjectSetName(A, 'A', ierr))
     65:   PetscCallA(MatDenseGetArray(A, aa, ierr))
     66:   call FillUpMatrix(n, aa)
     67:   PetscCallA(MatDenseRestoreArray(A, aa, ierr))
     68:   PetscCallA(MatSetOption(A, MAT_HERMITIAN, PETSC_TRUE, ierr))

     70: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     71: ! Scalar evaluation
     72: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     74:   x = 2.2
     75:   PetscCallA(FNEvaluateFunction(fn, x, y, ierr))
     76:   PetscCallA(FNEvaluateDerivative(fn, x, yp, ierr))

     78:   if (rank == 0) then
     79:     re = PetscRealPart(y)
     80:     im = PetscImaginaryPart(y)
     81:     if (abs(im) < 1.d-10) then
     82:       write (*, '(a3,f3.1,a,f8.5)') 'f(', PetscRealPart(x), ') = ', re
     83:     else
     84:       write (*, '(a3,f3.1,a,f10.5,sp,f9.5,a)') 'f(', PetscRealPart(x), ') = ', re, im, 'i'
     85:     end if
     86:     re = PetscRealPart(yp)
     87:     im = PetscImaginaryPart(yp)
     88:     if (abs(im) < 1.d-10) then
     89:       write (*, '(a3,f3.1,a,f8.5)') 'f''(', PetscRealPart(x), ') = ', re
     90:     else
     91:       write (*, '(a3,f3.1,a,f8.5,sp,f8.5,a)') 'f''(', PetscRealPart(x), ') = ', re, im, 'i'
     92:     end if
     93:   end if

     95: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     96: ! Compute matrix square root
     97: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     99:   PetscCallA(MatCreateSeqDense(PETSC_COMM_SELF, n, n, PETSC_NULL_SCALAR_ARRAY, S, ierr))
    100:   PetscCallA(PetscObjectSetName(S, 'S', ierr))
    101:   if (inplace) then
    102:     PetscCallA(MatCopy(A, S, SAME_NONZERO_PATTERN, ierr))
    103:     PetscCallA(MatSetOption(S, MAT_HERMITIAN, PETSC_TRUE, ierr))
    104:     PetscCallA(FNEvaluateFunctionMat(fn, S, PETSC_NULL_MAT, ierr))
    105:   else
    106:     PetscCallA(FNEvaluateFunctionMat(fn, A, S, ierr))
    107:   end if
    108:   if (verbose) then
    109:     if (rank == 0) write (*, *) 'Matrix A - - - - - - - -'
    110:     PetscCallA(MatView(A, PETSC_NULL_VIEWER, ierr))
    111:     if (rank == 0) write (*, *) 'Computed sqrtm(A) - - - - - - - -'
    112:     PetscCallA(MatView(S, PETSC_NULL_VIEWER, ierr))
    113:   end if

    115: ! *** check error ||S*S-A||_F
    116:   PetscCallA(MatMatMult(S, S, MAT_INITIAL_MATRIX, PETSC_DEFAULT_REAL, R, ierr))
    117:   if (eta /= 1.0) then
    118:     alpha = 1.0/(eta*eta)
    119:     PetscCallA(MatScale(R, alpha, ierr))
    120:   end if
    121:   alpha = -tau
    122:   PetscCallA(MatAXPY(R, alpha, A, SAME_NONZERO_PATTERN, ierr))
    123:   PetscCallA(MatNorm(R, NORM_FROBENIUS, nrm, ierr))
    124:   if (nrm < 100*PETSC_MACHINE_EPSILON) then
    125:     write (*, *) '||S*S-A||_F < 100*eps'
    126:   else
    127:     write (*, '(a,f8.5)') '||S*S-A||_F = ', nrm
    128:   end if

    130: ! *** Clean up
    131:   PetscCallA(MatDestroy(S, ierr))
    132:   PetscCallA(MatDestroy(R, ierr))
    133:   PetscCallA(MatDestroy(A, ierr))
    134:   PetscCallA(FNDestroy(fn, ierr))
    135:   PetscCallA(SlepcFinalize(ierr))

    137: contains

    139:   subroutine FillUpMatrix(n, X)
    140:     PetscInt    :: n, i, j
    141:     PetscScalar :: X(n, n)

    143:     do i = 1, n
    144:       X(i, i) = 2.5
    145:     end do
    146:     do j = 1, 2
    147:       do i = 1, n - j
    148:         X(i, i + j) = 1.d0
    149:         X(i + j, i) = 1.d0
    150:       end do
    151:     end do

    153:   end

    155: end program test7f

    157: !/*TEST
    158: !
    159: !   test:
    160: !      suffix: 1
    161: !      nsize: 1
    162: !      args: -fn_scale .13,2 -n 19 -fn_method {{0 1 2 3}} -inplace {{0 1}}
    163: !      filter: grep -v "computing matrix functions"
    164: !      output_file: output/test7f_1.out
    165: !
    166: !TEST*/