GCC Code Coverage Report

Directory: ./
File: src/eps/tutorials/ex6f.F90
Date: 2025-12-10 04:20:18
Exec Total Coverage
Lines: 90 95 94.7%
Functions: 2 2 100.0%
Branches: 60 116 51.7%
Line Branch Exec Source
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> ./ex6f [-help] [-m <m>] [all SLEPc options]
11 !
12 ! Description: Eigensystem from the Ising model for ferromagnetic materials.
13 ! Information about the model can be found at the following
14 ! site https://math.nist.gov/MatrixMarket/data/NEP
15 !
16 ! The command line options are:
17 ! -m <m>, where <m> is the number of 2x2 blocks, i.e. matrix size N=2*m
18 !
19 ! ----------------------------------------------------------------------
20 !
21 #include <slepc/finclude/slepceps.h>
22
23 module ex6fmodule
24 use slepceps
25 implicit none
26
27 contains
28 ! -------------------------------------------------------------------
29 ! The actual routine for the matrix-vector product
30 ! See https://math.nist.gov/MatrixMarket/data/NEP/mvmisg/mvmisg.html
31 !
32 ! Computes Y(:,1:M) = op(A)*X(:,1:M)
33 ! where op(A) is A or A' and A is the Ising matrix.
34 !
35 ! trans (input) integer
36 ! If trans = 0, compute y(:,1:M) = A*x(:,1:M)
37 ! If trans = 1, compute y(:,1:M) = A'*x(:,1:M)
38 !
39 ! n (input) integer
40 ! The order of the matrix, it has to be an even number.
41 !
42 ! m (input) integeR
43 ! The number of columns of x to multiply.
44 !
45 ! x (input) double precision array, dimension (ldx, m)
46 ! x contains the matrix (vectors) x.
47 !
48 ! ldx (input) integer
49 ! The leading dimension of array x, ldx >= max(1, n)
50 !
51 ! y (output) double precision array, dimension (ldx, m)
52 ! contains the product of the matrix op(A) with x.
53 !
54 ! ldy (input) integer
55 ! The leading dimension of array y, ldy >= max(1, n)
56 !
57 868 subroutine mvmisg(trans, n, m, x, ldx, y, ldy)
58 use petscsys
59 implicit none
60
61 PetscInt :: ldy, ldx, m, n, trans
62 PetscScalar :: y(ldy, *), x(ldx, *)
63 PetscInt :: i, k
64 PetscReal :: alpha, beta, cosa, cosb, sina, sinb
65 PetscScalar :: temp, temp1
66
67 868 alpha = PETSC_PI/4
68 868 beta = PETSC_PI/4
69 868 cosa = cos(alpha)
70 868 sina = sin(alpha)
71 868 cosb = cos(beta)
72 868 sinb = sin(beta)
73
74
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 868 if (trans == 0) then
75
76 ! ** Compute y(:,1:m) = A*x(:,1:m)
77
78
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 1048 do k = 1, m
79 524 y(1, k) = cosb*x(1, k) - sinb*x(n, k)
80
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 524 do i = 2, n - 1, 2
81 15196 y(i, k) = cosb*x(i, k) + sinb*x(i + 1, k)
82
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 15196 y(i + 1, k) = -sinb*x(i, k) + cosb*x(i + 1, k)
83 end do
84 524 y(n, k) = sinb*x(1, k) + cosb*x(n, k)
85
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 1048 do i = 1, n, 2
86 15720 temp = cosa*y(i, k) + sina*y(i + 1, k)
87 15720 y(i + 1, k) = -sina*y(i, k) + cosa*y(i + 1, k)
88
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 15720 y(i, k) = temp
89 end do
90 end do
91
92
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 344 else if (trans == 1) then
93
94 ! ** Compute y(:1:m) = A'*x(:,1:m)
95
96
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 688 do k = 1, m
97
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 344 do i = 1, n, 2
98 10320 y(i, k) = cosa*x(i, k) - sina*x(i + 1, k)
99
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 10320 y(i + 1, k) = sina*x(i, k) + cosa*x(i + 1, k)
100 end do
101 344 temp = cosb*y(1, k) + sinb*y(n, k)
102
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 344 do i = 2, n - 1, 2
103 9976 temp1 = cosb*y(i, k) - sinb*y(i + 1, k)
104 9976 y(i + 1, k) = sinb*y(i, k) + cosb*y(i + 1, k)
105
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 9976 y(i, k) = temp1
106 end do
107 344 y(n, k) = -sinb*y(1, k) + cosb*y(n, k)
108 688 y(1, k) = temp
109 end do
110
111 end if
112 868 end subroutine
113
114 ! -------------------------------------------------------------------
115 ! MatMult_Ising - user provided matrix-vector multiply
116 !
117 ! Input Parameters:
118 ! A - matrix
119 ! x - input vector
120 !
121 ! Output Parameter:
122 ! y - output vector
123 !
124 524 subroutine MatMult_Ising(A, x, y, ierr)
125 use petscmat
126 implicit none
127
128 Mat :: A
129 Vec :: x, y
130 PetscInt :: trans, one, N
131 524 PetscScalar, pointer :: xx(:), yy(:)
132 PetscErrorCode :: ierr
133
134
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 524 PetscCall(MatGetSize(A, N, PETSC_NULL_INTEGER, ierr))
135
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 524 PetscCall(VecGetArrayRead(x, xx, ierr))
136
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 524 PetscCall(VecGetArray(y, yy, ierr))
137
138 524 trans = 0
139 524 one = 1
140
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
 524 call mvmisg(trans, N, one, xx, N, yy, N)
141
142
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 524 PetscCall(VecRestoreArrayRead(x, xx, ierr))
143
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 524 PetscCall(VecRestoreArray(y, yy, ierr))
144 524 end subroutine
145
146 ! -------------------------------------------------------------------
147 ! MatMultTranspose_Ising - user provided transpose matrix-vector multiply
148 !
149 ! Input Parameters:
150 ! A - matrix
151 ! x - input vector
152 !
153 ! Output Parameter:
154 ! y - output vector
155 !
156 344 subroutine MatMultTranspose_Ising(A, x, y, ierr)
157 use petscmat
158 implicit none
159
160 Mat :: A
161 Vec :: x, y
162 PetscInt :: trans, one, N
163 344 PetscScalar, pointer :: xx(:), yy(:)
164 PetscErrorCode :: ierr
165
166
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 344 PetscCall(MatGetSize(A, N, PETSC_NULL_INTEGER, ierr))
167
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 344 PetscCall(VecGetArrayRead(x, xx, ierr))
168
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 344 PetscCall(VecGetArray(y, yy, ierr))
169
170 344 trans = 1
171 344 one = 1
172
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
 344 call mvmisg(trans, N, one, xx, N, yy, N)
173
174
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 344 PetscCall(VecRestoreArrayRead(x, xx, ierr))
175
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 344 PetscCall(VecRestoreArray(y, yy, ierr))
176 344 end subroutine
177
178 end module ex6fmodule
179
180 8 program ex6f
181 4 use slepceps
182 use ex6fmodule
183 implicit none
184
185 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186 ! Declarations
187 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188
189 Mat :: A ! operator matrix
190 EPS :: eps ! aeigenproblem solver context
191 EPSType :: tname
192 PetscReal :: tol
193 PetscInt :: N, m
194 PetscInt :: nev, maxit, its
195 PetscMPIInt :: sz, rank
196 PetscErrorCode :: ierr
197 PetscBool :: flg, terse
198
199 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
200 ! Beginning of program
201 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202
203
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
204 #if defined(PETSC_USE_COMPLEX)
205 SETERRA(PETSC_COMM_SELF, PETSC_ERR_SUP, 'This example requires real numbers')
206 #endif
207
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, sz, ierr))
208
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
209
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCheckA(sz == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, 'This is a uniprocessor example only!')
210 4 m = 30
211
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-m', m, flg, ierr))
212 4 N = 2*m
213
214
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (rank == 0) then
215 4 write (*, '(/a,i6,a/)') 'Ising Model Eigenproblem, m=', m, ', (N=2*m)'
216 end if
217
218 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
219 ! Register the matrix-vector subroutine for the operator that defines
220 ! the eigensystem, Ax=kx
221 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
222
223
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(MatCreateShell(PETSC_COMM_WORLD, N, N, N, N, PETSC_NULL_INTEGER, A, ierr))
224
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(MatShellSetOperation(A, MATOP_MULT, MatMult_Ising, ierr))
225
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(MatShellSetOperation(A, MATOP_MULT_TRANSPOSE, MatMultTranspose_Ising, ierr))
226
227 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
228 ! Create the eigensolver and display info
229 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
230
231 ! ** Create eigensolver context
232
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSCreate(PETSC_COMM_WORLD, eps, ierr))
233
234 ! ** Set operators. In this case, it is a standard eigenvalue problem
235
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSSetOperators(eps, A, PETSC_NULL_MAT, ierr))
236
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSSetProblemType(eps, EPS_NHEP, ierr))
237
238 ! ** Set solver parameters at runtime
239
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSSetFromOptions(eps, ierr))
240
241 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
242 ! Solve the eigensystem
243 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
244
245
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSSolve(eps, ierr))
246
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSGetIterationNumber(eps, its, ierr))
247
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (rank == 0) then
248 4 write (*, '(a,i4)') ' Number of iterations of the method: ', its
249 end if
250
251 ! ** Optional: Get some information from the solver and display it
252
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSGetType(eps, tname, ierr))
253
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (rank == 0) then
254 4 write (*, '(a,a)') ' Solution method: ', tname
255 end if
256
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSGetDimensions(eps, nev, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
257
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (rank == 0) then
258 4 write (*, '(a,i2)') ' Number of requested eigenvalues:', nev
259 end if
260
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSGetTolerances(eps, tol, maxit, ierr))
261
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (rank == 0) then
262 4 write (*, '(a,1pe11.4,a,i6)') ' Stopping condition: tol=', tol, ', maxit=', maxit
263 end if
264
265 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
266 ! Display solution and clean up
267 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
268
269 ! ** show detailed info unless -terse option is given by user
270
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-terse', terse, ierr))
271
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (terse) then
272
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSErrorView(eps, EPS_ERROR_RELATIVE, PETSC_NULL_VIEWER, ierr))
273 else
274 ✗ PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL, ierr))
275 ✗ PetscCallA(EPSConvergedReasonView(eps, PETSC_VIEWER_STDOUT_WORLD, ierr))
276 ✗ PetscCallA(EPSErrorView(eps, EPS_ERROR_RELATIVE, PETSC_VIEWER_STDOUT_WORLD, ierr))
277 ✗ PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD, ierr))
278 end if
279
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSDestroy(eps, ierr))
280
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(MatDestroy(A, ierr))
281
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 4 PetscCallA(SlepcFinalize(ierr))
282 ✗ end program ex6f
283
284 !/*TEST
285 !
286 ! testset:
287 ! args: -eps_max_it 1000 -eps_ncv 12 -eps_tol 1e-5 -eps_nev 4 -eps_largest_imaginary -terse
288 ! requires: !complex
289 ! output_file: output/ex6f_1.out
290 ! filter: grep -v iterations | sed -e 's/-0.00000/0.00000/g'
291 ! test:
292 ! suffix: 1
293 ! test:
294 ! suffix: 1_ts
295 ! args: -eps_two_sided
296 ! requires: !single
297 !
298 !TEST*/
299