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: #include <slepc/finclude/slepcfn.h>
17: program test7f
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 == 0) then
47: write (*, '(/a,i3,a)') 'Matrix square root, n =', n, ' (Fortran)'
48: end if
50: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51: ! Create FN object and matrix
52: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: PetscCallA(FNCreate(PETSC_COMM_WORLD, fn, ierr))
55: PetscCallA(FNSetType(fn, FNSQRT, ierr))
56: tau = 0.15
57: eta = 1.0
58: PetscCallA(FNSetScale(fn, tau, eta, ierr))
59: PetscCallA(FNSetFromOptions(fn, ierr))
60: PetscCallA(FNGetScale(fn, tau, eta, ierr))
61: PetscCallA(FNView(fn, PETSC_NULL_VIEWER, ierr))
63: PetscCallA(MatCreateSeqDense(PETSC_COMM_SELF, n, n, PETSC_NULL_SCALAR_ARRAY, A, ierr))
64: PetscCallA(PetscObjectSetName(A, 'A', ierr))
65: PetscCallA(MatDenseGetArray(A, aa, ierr))
66: call FillUpMatrix(n, aa)
67: PetscCallA(MatDenseRestoreArray(A, aa, ierr))
68: PetscCallA(MatSetOption(A, MAT_HERMITIAN, PETSC_TRUE, ierr))
70: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71: ! Scalar evaluation
72: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74: x = 2.2
75: PetscCallA(FNEvaluateFunction(fn, x, y, ierr))
76: PetscCallA(FNEvaluateDerivative(fn, x, yp, ierr))
78: if (rank == 0) then
79: re = PetscRealPart(y)
80: im = PetscImaginaryPart(y)
81: if (abs(im) < 1.d-10) then
82: write (*, '(a3,f3.1,a,f8.5)') 'f(', PetscRealPart(x), ') = ', re
83: else
84: write (*, '(a3,f3.1,a,f10.5,sp,f9.5,a)') 'f(', PetscRealPart(x), ') = ', re, im, 'i'
85: end if
86: re = PetscRealPart(yp)
87: im = PetscImaginaryPart(yp)
88: if (abs(im) < 1.d-10) then
89: write (*, '(a3,f3.1,a,f8.5)') 'f''(', PetscRealPart(x), ') = ', re
90: else
91: write (*, '(a3,f3.1,a,f8.5,sp,f8.5,a)') 'f''(', PetscRealPart(x), ') = ', re, im, 'i'
92: end if
93: end if
95: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: ! Compute matrix square root
97: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99: PetscCallA(MatCreateSeqDense(PETSC_COMM_SELF, n, n, PETSC_NULL_SCALAR_ARRAY, S, ierr))
100: PetscCallA(PetscObjectSetName(S, 'S', ierr))
101: if (inplace) then
102: PetscCallA(MatCopy(A, S, SAME_NONZERO_PATTERN, ierr))
103: PetscCallA(MatSetOption(S, MAT_HERMITIAN, PETSC_TRUE, ierr))
104: PetscCallA(FNEvaluateFunctionMat(fn, S, PETSC_NULL_MAT, ierr))
105: else
106: PetscCallA(FNEvaluateFunctionMat(fn, A, S, ierr))
107: end if
108: if (verbose) then
109: if (rank == 0) write (*, *) 'Matrix A - - - - - - - -'
110: PetscCallA(MatView(A, PETSC_NULL_VIEWER, ierr))
111: if (rank == 0) write (*, *) 'Computed sqrtm(A) - - - - - - - -'
112: PetscCallA(MatView(S, PETSC_NULL_VIEWER, ierr))
113: end if
115: ! *** check error ||S*S-A||_F
116: PetscCallA(MatMatMult(S, S, MAT_INITIAL_MATRIX, PETSC_DEFAULT_REAL, R, ierr))
117: if (eta /= 1.0) then
118: alpha = 1.0/(eta*eta)
119: PetscCallA(MatScale(R, alpha, ierr))
120: end if
121: alpha = -tau
122: PetscCallA(MatAXPY(R, alpha, A, SAME_NONZERO_PATTERN, ierr))
123: PetscCallA(MatNorm(R, NORM_FROBENIUS, nrm, ierr))
124: if (nrm < 100*PETSC_MACHINE_EPSILON) then
125: write (*, *) '||S*S-A||_F < 100*eps'
126: else
127: write (*, '(a,f8.5)') '||S*S-A||_F = ', nrm
128: end if
130: ! *** Clean up
131: PetscCallA(MatDestroy(S, ierr))
132: PetscCallA(MatDestroy(R, ierr))
133: PetscCallA(MatDestroy(A, ierr))
134: PetscCallA(FNDestroy(fn, ierr))
135: PetscCallA(SlepcFinalize(ierr))
137: contains
139: subroutine FillUpMatrix(n, X)
140: PetscInt :: n, i, j
141: PetscScalar :: X(n, n)
143: do i = 1, n
144: X(i, i) = 2.5
145: end do
146: do j = 1, 2
147: do i = 1, n - j
148: X(i, i + j) = 1.d0
149: X(i + j, i) = 1.d0
150: end do
151: end do
153: end
155: end program test7f
157: !/*TEST
158: !
159: ! test:
160: ! suffix: 1
161: ! nsize: 1
162: ! args: -fn_scale .13,2 -n 19 -fn_method {{0 1 2 3}} -inplace {{0 1}}
163: ! filter: grep -v "computing matrix functions"
164: ! output_file: output/test7f_1.out
165: !
166: !TEST*/