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:       program main
     17: #include <slepc/finclude/slepcfn.h>
     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 .eq. 0) then
     47:         write(*,100) n
     48:       endif
     49:  100  format (/'Matrix square root, n =',I3,' (Fortran)')

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

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

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

     71: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     72: !     Scalar evaluation
     73: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

     79:       if (rank .eq. 0) then
     80:         re = PetscRealPart(y)
     81:         im = PetscImaginaryPart(y)
     82:         if (abs(im).lt.1.d-10) then
     83:           write(*,110) 'f', PetscRealPart(x), re
     84:         else
     85:           write(*,120) 'f', PetscRealPart(x), re, im
     86:         endif
     87:         re = PetscRealPart(yp)
     88:         im = PetscImaginaryPart(yp)
     89:         if (abs(im).lt.1.d-10) then
     90:           write(*,110) 'f''', PetscRealPart(x), re
     91:         else
     92:           write(*,120) 'f''', PetscRealPart(x), re, im
     93:         endif
     94:       endif
     95:  110  format (A2,'(',F4.1,') = ',F8.5)
     96:  120  format (A2,'(',F4.1,') = ',F8.5,SP,F8.5,'i')

     98: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     99: !     Compute matrix square root
    100: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

    118: !     *** check error ||S*S-A||_F
    119:       PetscCallA(MatMatMult(S,S,MAT_INITIAL_MATRIX,PETSC_DEFAULT_REAL,R,ierr))
    120:       if (eta .ne. 1.0) then
    121:         alpha = 1.0/(eta*eta)
    122:         PetscCallA(MatScale(R,alpha,ierr))
    123:       endif
    124:       alpha = -tau
    125:       PetscCallA(MatAXPY(R,alpha,A,SAME_NONZERO_PATTERN,ierr))
    126:       PetscCallA(MatNorm(R,NORM_FROBENIUS,nrm,ierr))
    127:       if (nrm<100*PETSC_MACHINE_EPSILON) then
    128:         write (*,*) '||S*S-A||_F < 100*eps'
    129:       else
    130:         write (*,130) nrm
    131:       endif
    132:  130  format ('||S*S-A||_F = ',F8.5)

    134: !     *** Clean up
    135:       PetscCallA(MatDestroy(S,ierr))
    136:       PetscCallA(MatDestroy(R,ierr))
    137:       PetscCallA(MatDestroy(A,ierr))
    138:       PetscCallA(FNDestroy(fn,ierr))
    139:       PetscCallA(SlepcFinalize(ierr))

    141:       contains

    143: ! -----------------------------------------------------------------

    145:       subroutine FillUpMatrix(n,X)
    146:       PetscInt    n,i,j
    147:       PetscScalar X(n,n)

    149:       do i=1,n
    150:         X(i,i) = 2.5
    151:       end do
    152:       do j=1,2
    153:         do i=1,n-j
    154:           X(i,i+j) = 1.d0
    155:           X(i+j,i) = 1.d0
    156:         end do
    157:       end do

    159:       end

    161:       end program main

    163: !/*TEST
    164: !
    165: !   test:
    166: !      suffix: 1
    167: !      nsize: 1
    168: !      args: -fn_scale .13,2 -n 19 -fn_method {{0 1 2 3}shared output}
    169: !      filter: grep -v "computing matrix functions"
    170: !      output_file: output/test7f_1.out
    171: !
    172: !TEST*/