GCC Code Coverage Report

Directory: ./
File: src/eps/tutorials/ex6f.F90
Date: 2026-02-22 03:58:10
Exec Total Coverage
Lines: 88 93 94.6%
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 2620 subroutine MatMult_Ising(A, x, y, ierr)
125 use petscmat
126 implicit none
127
128 Mat :: A
129 Vec :: x, y
130 PetscInt :: trans, N
131 524 PetscScalar, pointer :: xx(:), yy(:)
132 PetscErrorCode, intent(out) :: 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
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, 1_PETSC_INT_KIND, xx, N, yy, N)
140
141
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 524 PetscCall(VecRestoreArrayRead(x, xx, ierr))
142
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 524 PetscCall(VecRestoreArray(y, yy, ierr))
143 524 end subroutine
144
145 ! -------------------------------------------------------------------
146 ! MatMultTranspose_Ising - user provided transpose matrix-vector multiply
147 !
148 ! Input Parameters:
149 ! A - matrix
150 ! x - input vector
151 !
152 ! Output Parameter:
153 ! y - output vector
154 !
155 1720 subroutine MatMultTranspose_Ising(A, x, y, ierr)
156 use petscmat
157 implicit none
158
159 Mat :: A
160 Vec :: x, y
161 PetscInt :: trans, N
162 344 PetscScalar, pointer :: xx(:), yy(:)
163 PetscErrorCode, intent(out) :: ierr
164
165
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 344 PetscCall(MatGetSize(A, N, PETSC_NULL_INTEGER, ierr))
166
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 344 PetscCall(VecGetArrayRead(x, xx, ierr))
167
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 344 PetscCall(VecGetArray(y, yy, ierr))
168
169 344 trans = 1
170
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, 1_PETSC_INT_KIND, xx, N, yy, N)
171
172
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 344 PetscCall(VecRestoreArrayRead(x, xx, ierr))
173
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 344 PetscCall(VecRestoreArray(y, yy, ierr))
174 344 end subroutine
175
176 end module ex6fmodule
177
178 52 program ex6f
179 4 use slepceps
180 use ex6fmodule
181 implicit none
182
183 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184 ! Declarations
185 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186
187 Mat :: A ! operator matrix
188 EPS :: eps ! aeigenproblem solver context
189 EPSType :: tname
190 PetscReal :: tol
191 PetscInt :: N, m
192 PetscInt :: nev, maxit, its
193 PetscMPIInt :: sz, rank
194 PetscErrorCode :: ierr
195 PetscBool :: flg, terse
196
197 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198 ! Beginning of program
199 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
200
201
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
202 #if defined(PETSC_USE_COMPLEX)
203 SETERRA(PETSC_COMM_SELF, PETSC_ERR_SUP, 'This example requires real numbers')
204 #endif
205
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, sz, ierr))
206
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
207
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!')
208 4 m = 30
209
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-m', m, flg, ierr))
210 4 N = 2*m
211
212
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (rank == 0) then
213 4 write (*, '(/a,i6,a/)') 'Ising Model Eigenproblem, m=', m, ', (N=2*m)'
214 end if
215
216 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
217 ! Register the matrix-vector subroutine for the operator that defines
218 ! the eigensystem, Ax=kx
219 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
220
221
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))
222
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(MatShellSetOperation(A, MATOP_MULT, MatMult_Ising, ierr))
223
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(MatShellSetOperation(A, MATOP_MULT_TRANSPOSE, MatMultTranspose_Ising, ierr))
224
225 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
226 ! Create the eigensolver and display info
227 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
228
229 ! ** Create eigensolver context
230
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSCreate(PETSC_COMM_WORLD, eps, ierr))
231
232 ! ** Set operators. In this case, it is a standard eigenvalue problem
233
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSSetOperators(eps, A, PETSC_NULL_MAT, ierr))
234
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSSetProblemType(eps, EPS_NHEP, ierr))
235
236 ! ** Set solver parameters at runtime
237
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSSetFromOptions(eps, ierr))
238
239 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
240 ! Solve the eigensystem
241 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
242
243
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSSolve(eps, ierr))
244
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSGetIterationNumber(eps, its, ierr))
245
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (rank == 0) then
246 4 write (*, '(a,i4)') ' Number of iterations of the method: ', its
247 end if
248
249 ! ** Optional: Get some information from the solver and display it
250
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSGetType(eps, tname, ierr))
251
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (rank == 0) then
252 4 write (*, '(a,a)') ' Solution method: ', tname
253 end if
254
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSGetDimensions(eps, nev, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
255
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (rank == 0) then
256 4 write (*, '(a,i2)') ' Number of requested eigenvalues:', nev
257 end if
258
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSGetTolerances(eps, tol, maxit, ierr))
259
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (rank == 0) then
260 4 write (*, '(a,1pe11.4,a,i6)') ' Stopping condition: tol=', tol, ', maxit=', maxit
261 end if
262
263 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
264 ! Display solution and clean up
265 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
266
267 ! ** show detailed info unless -terse option is given by user
268
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-terse', terse, ierr))
269
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (terse) then
270
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSErrorView(eps, EPS_ERROR_RELATIVE, PETSC_NULL_VIEWER, ierr))
271 else
272 ✗ PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL, ierr))
273 ✗ PetscCallA(EPSConvergedReasonView(eps, PETSC_VIEWER_STDOUT_WORLD, ierr))
274 ✗ PetscCallA(EPSErrorView(eps, EPS_ERROR_RELATIVE, PETSC_VIEWER_STDOUT_WORLD, ierr))
275 ✗ PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD, ierr))
276 end if
277
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSDestroy(eps, ierr))
278
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(MatDestroy(A, ierr))
279
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 4 PetscCallA(SlepcFinalize(ierr))
280 ✗ end program ex6f
281
282 !/*TEST
283 !
284 ! testset:
285 ! args: -eps_max_it 1000 -eps_ncv 12 -eps_tol 1e-5 -eps_nev 4 -eps_largest_imaginary -terse
286 ! requires: !complex
287 ! output_file: output/ex6f_1.out
288 ! filter: grep -v iterations | sed -e 's/-0.00000/0.00000/g'
289 ! test:
290 ! suffix: 1
291 ! test:
292 ! suffix: 1_ts
293 ! args: -eps_two_sided
294 ! requires: !single
295 !
296 !TEST*/
297