Actual source code: test7f.F

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: !
  2: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3: !  SLEPc - Scalable Library for Eigenvalue Problem Computations
  4: !  Copyright (c) 2002-2019, 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    aa(1),tau,eta,alpha,x,y,yp
 33:       PetscOffset    ia

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

 39:       call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
 40:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
 41:       n = 10
 42:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,  &
 43:      &                        '-n',n,flg,ierr)
 44:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
 45:      &                        '-verbose',verbose,ierr)
 46:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
 47:      &                        '-inplace',inplace,ierr)

 49:       if (rank .eq. 0) then
 50:         write(*,100) n
 51:       endif
 52:  100  format (/'Matrix square root, n =',I3,' (Fortran)')

 54: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 55: !     Create FN object and matrix
 56: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 58:       call FNCreate(PETSC_COMM_WORLD,fn,ierr)
 59:       call FNSetType(fn,FNSQRT,ierr)
 60:       tau = 0.15
 61:       eta = 1.0
 62:       call FNSetScale(fn,tau,eta,ierr)
 63:       call FNSetFromOptions(fn,ierr)
 64:       call FNGetScale(fn,tau,eta,ierr)
 65:       call FNView(fn,PETSC_NULL_VIEWER,ierr)

 67:       call MatCreateSeqDense(PETSC_COMM_SELF,n,n,PETSC_NULL_SCALAR,A,   &
 68:      &     ierr)
 69:       call PetscObjectSetName(A,'A',ierr)
 70:       call MatDenseGetArray(A,aa,ia,ierr)
 71:       call FillUpMatrix(n,aa(ia+1))
 72:       call MatDenseRestoreArray(A,aa,ia,ierr)
 73:       call MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE,ierr)

 75: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 76: !     Scalar evaluation
 77: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 79:       x = 2.2
 80:       call FNEvaluateFunction(fn,x,y,ierr)
 81:       call FNEvaluateDerivative(fn,x,yp,ierr)

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

102: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: !     Compute matrix square root
104: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

106:       call MatCreateSeqDense(PETSC_COMM_SELF,n,n,PETSC_NULL_SCALAR,S,   &
107:      &     ierr)
108:       call PetscObjectSetName(S,'S',ierr)
109:       call MatCreateSeqDense(PETSC_COMM_SELF,n,n,PETSC_NULL_SCALAR,R,   &
110:      &     ierr)
111:       call PetscObjectSetName(R,'R',ierr)
112:       if (inplace) then
113:         call MatCopy(A,S,SAME_NONZERO_PATTERN,ierr)
114:         call MatSetOption(S,MAT_HERMITIAN,PETSC_TRUE,ierr)
115:         call FNEvaluateFunctionMat(fn,S,PETSC_NULL_MAT,ierr)
116:       else
117:         call FNEvaluateFunctionMat(fn,A,S,ierr)
118:       endif
119:       if (verbose) then
120:         if (rank .eq. 0) write (*,*) 'Matrix A - - - - - - - -'
121:         call MatView(A,PETSC_NULL_VIEWER,ierr)
122:         if (rank .eq. 0) write (*,*) 'Computed sqrtm(A) - - - - - - - -'
123:         call MatView(S,PETSC_NULL_VIEWER,ierr)
124:       endif

126: !     *** check error ||S*S-A||_F
127:       call MatMatMult(S,S,MAT_REUSE_MATRIX,PETSC_DEFAULT_REAL,R,ierr)
128:       if (eta .ne. 1.0) then
129:         alpha = 1.0/(eta*eta)
130:         call MatScale(R,alpha,ierr)
131:       endif
132:       alpha = -tau
133:       call MatAXPY(R,alpha,A,SAME_NONZERO_PATTERN,ierr)
134:       call MatNorm(R,NORM_FROBENIUS,nrm,ierr)
135:       if (nrm<100*PETSC_MACHINE_EPSILON) then
136:         write (*,*) '||S*S-A||_F < 100*eps'
137:       else
138:         write (*,130) nrm
139:       endif
140:  130  format ('||S*S-A||_F = ',F8.5)

142: !     *** Clean up
143:       call MatDestroy(S,ierr)
144:       call MatDestroy(R,ierr)
145:       call MatDestroy(A,ierr)
146:       call FNDestroy(fn,ierr)
147:       call SlepcFinalize(ierr)
148:       end

150: ! -----------------------------------------------------------------

152:       subroutine FillUpMatrix(n,X)
153:       PetscInt    n,i,j
154:       PetscScalar X(n,n)

156:       do i=1,n
157:         X(i,i) = 2.5
158:       end do
159:       do j=1,2
160:         do i=1,n-j
161:           X(i,i+j) = 1.d0
162:           X(i+j,i) = 1.d0
163:         end do
164:       end do
165:       return
166:       end

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