Actual source code: test7f.F90
slepc-3.21.1 2024-04-26
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,A,ierr))
65: PetscCallA(PetscObjectSetName(A,'A',ierr))
66: PetscCallA(MatDenseGetArrayF90(A,aa,ierr))
67: call FillUpMatrix(n,aa)
68: PetscCallA(MatDenseRestoreArrayF90(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,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))
140: end
142: ! -----------------------------------------------------------------
144: subroutine FillUpMatrix(n,X)
145: PetscInt n,i,j
146: PetscScalar X(n,n)
148: do i=1,n
149: X(i,i) = 2.5
150: end do
151: do j=1,2
152: do i=1,n-j
153: X(i,i+j) = 1.d0
154: X(i+j,i) = 1.d0
155: end do
156: end do
158: end
160: !/*TEST
161: !
162: ! test:
163: ! suffix: 1
164: ! nsize: 1
165: ! args: -fn_scale .13,2 -n 19 -fn_method {{0 1 2 3}shared output}
166: ! filter: grep -v "computing matrix functions"
167: ! output_file: output/test7f_1.out
168: !
169: !TEST*/