Actual source code: test7f.F
slepc-3.11.2 2019-07-30
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*/