Actual source code: test7f.F
slepc-3.20.2 2024-03-15
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: 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 MatDenseGetArrayF90(A,aa,ierr)
71: call FillUpMatrix(n,aa)
72: call MatDenseRestoreArrayF90(A,aa,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: if (inplace) then
110: call MatCopy(A,S,SAME_NONZERO_PATTERN,ierr)
111: call MatSetOption(S,MAT_HERMITIAN,PETSC_TRUE,ierr)
112: call FNEvaluateFunctionMat(fn,S,PETSC_NULL_MAT,ierr)
113: else
114: call FNEvaluateFunctionMat(fn,A,S,ierr)
115: endif
116: if (verbose) then
117: if (rank .eq. 0) write (*,*) 'Matrix A - - - - - - - -'
118: call MatView(A,PETSC_NULL_VIEWER,ierr)
119: if (rank .eq. 0) write (*,*) 'Computed sqrtm(A) - - - - - - - -'
120: call MatView(S,PETSC_NULL_VIEWER,ierr)
121: endif
123: ! *** check error ||S*S-A||_F
124: call MatMatMult(S,S,MAT_INITIAL_MATRIX,PETSC_DEFAULT_REAL,R,ierr)
125: if (eta .ne. 1.0) then
126: alpha = 1.0/(eta*eta)
127: call MatScale(R,alpha,ierr)
128: endif
129: alpha = -tau
130: call MatAXPY(R,alpha,A,SAME_NONZERO_PATTERN,ierr)
131: call MatNorm(R,NORM_FROBENIUS,nrm,ierr)
132: if (nrm<100*PETSC_MACHINE_EPSILON) then
133: write (*,*) '||S*S-A||_F < 100*eps'
134: else
135: write (*,130) nrm
136: endif
137: 130 format ('||S*S-A||_F = ',F8.5)
139: ! *** Clean up
140: call MatDestroy(S,ierr)
141: call MatDestroy(R,ierr)
142: call MatDestroy(A,ierr)
143: call FNDestroy(fn,ierr)
144: call SlepcFinalize(ierr)
145: end
147: ! -----------------------------------------------------------------
149: subroutine FillUpMatrix(n,X)
150: PetscInt n,i,j
151: PetscScalar X(n,n)
153: do i=1,n
154: X(i,i) = 2.5
155: end do
156: do j=1,2
157: do i=1,n-j
158: X(i,i+j) = 1.d0
159: X(i+j,i) = 1.d0
160: end do
161: end do
162: return
163: end
165: !/*TEST
166: !
167: ! test:
168: ! suffix: 1
169: ! nsize: 1
170: ! args: -fn_scale .13,2 -n 19 -fn_method {{0 1 2 3}shared output}
171: ! filter: grep -v "computing matrix functions"
172: ! output_file: output/test7f_1.out
173: !
174: !TEST*/