GCC Code Coverage Report

Directory: ./
File: src/eps/tutorials/ex6f.F90
Date: 2025-10-03 04:28:47
Exec Total Coverage
Lines: 96 101 95.0%
Functions: 5 5 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 8 program main
22 #include <slepc/finclude/slepceps.h>
23 4 use slepceps
24 implicit none
25
26 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
27 ! Declarations
28 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
29 !
30 ! Variables:
31 ! A operator matrix
32 ! eps eigenproblem solver context
33
34 Mat A
35 EPS eps
36 EPSType tname
37 PetscReal tol
38 PetscInt N, m
39 PetscInt nev, maxit, its
40 PetscMPIInt sz, rank
41 PetscErrorCode ierr
42 PetscBool flg, terse
43
44 ! This is the routine to use for matrix-free approach
45 !
46 external MatMult_Ising, MatMultTranspose_Ising
47
48 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
49 ! Beginning of program
50 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51
52
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER,ierr))
53 #if defined(PETSC_USE_COMPLEX)
54 SETERRA(PETSC_COMM_SELF,PETSC_ERR_SUP,'This example requires real numbers')
55 #endif
56
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,sz,ierr))
57
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
58
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 if (sz .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only!'); endif
59 4 m = 30
60
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr))
61 4 N = 2*m
62
63
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (rank .eq. 0) then
64 4 write(*,'(/A,I6,A/)') 'Ising Model Eigenproblem, m=',m,', (N=2*m)'
65 endif
66
67 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68 ! Register the matrix-vector subroutine for the operator that defines
69 ! the eigensystem, Ax=kx
70 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71
72
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))
73
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(MatShellSetOperation(A,MATOP_MULT,MatMult_Ising,ierr))
74
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(MatShellSetOperation(A,MATOP_MULT_TRANSPOSE,MatMultTranspose_Ising,ierr))
75
76 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
77 ! Create the eigensolver and display info
78 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79
80 ! ** Create eigensolver context
81
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSCreate(PETSC_COMM_WORLD,eps,ierr))
82
83 ! ** Set operators. In this case, it is a standard eigenvalue problem
84
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSSetOperators(eps,A,PETSC_NULL_MAT,ierr))
85
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSSetProblemType(eps,EPS_NHEP,ierr))
86
87 ! ** Set solver parameters at runtime
88
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSSetFromOptions(eps,ierr))
89
90 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91 ! Solve the eigensystem
92 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93
94
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSSolve(eps,ierr))
95
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSGetIterationNumber(eps,its,ierr))
96
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (rank .eq. 0) then
97 4 write(*,'(A,I4)') ' Number of iterations of the method: ',its
98 endif
99
100 ! ** Optional: Get some information from the solver and display it
101
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSGetType(eps,tname,ierr))
102
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (rank .eq. 0) then
103 4 write(*,'(A,A)') ' Solution method: ', tname
104 endif
105
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSGetDimensions(eps,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
106
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (rank .eq. 0) then
107 4 write(*,'(A,I2)') ' Number of requested eigenvalues:',nev
108 endif
109
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSGetTolerances(eps,tol,maxit,ierr))
110
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (rank .eq. 0) then
111 4 write(*,'(A,1PE11.4,A,I6)') ' Stopping condition: tol=',tol,', maxit=', maxit
112 endif
113
114 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115 ! Display solution and clean up
116 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117
118 ! ** show detailed info unless -terse option is given by user
119
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr))
120
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (terse) then
121
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr))
122 else
123 ✗ PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr))
124 ✗ PetscCallA(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD,ierr))
125 ✗ PetscCallA(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD,ierr))
126 ✗ PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr))
127 endif
128
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(EPSDestroy(eps,ierr))
129
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 PetscCallA(MatDestroy(A,ierr))
130
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 4 PetscCallA(SlepcFinalize(ierr))
131 ✗ end
132
133 ! -------------------------------------------------------------------
134 !
135 ! MatMult_Ising - user provided matrix-vector multiply
136 !
137 ! Input Parameters:
138 ! A - matrix
139 ! x - input vector
140 !
141 ! Output Parameter:
142 ! y - output vector
143 !
144 524 subroutine MatMult_Ising(A,x,y,ierr)
145 #include <petsc/finclude/petscmat.h>
146 use petscmat
147 implicit none
148
149 Mat A
150 Vec x,y
151 PetscInt trans,one,N
152 524 PetscScalar, pointer :: xx(:),yy(:)
153 PetscErrorCode ierr
154
155 ! The actual routine for the matrix-vector product
156 external mvmisg
157
158
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 524 PetscCall(MatGetSize(A,N,PETSC_NULL_INTEGER,ierr))
159
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 524 PetscCall(VecGetArrayRead(x,xx,ierr))
160
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 524 PetscCall(VecGetArray(y,yy,ierr))
161
162 524 trans = 0
163 524 one = 1
164
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)
165
166
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 524 PetscCall(VecRestoreArrayRead(x,xx,ierr))
167
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 524 PetscCall(VecRestoreArray(y,yy,ierr))
168
169 524 end
170
171 ! -------------------------------------------------------------------
172 !
173 ! MatMultTranspose_Ising - user provided transpose matrix-vector multiply
174 !
175 ! Input Parameters:
176 ! A - matrix
177 ! x - input vector
178 !
179 ! Output Parameter:
180 ! y - output vector
181 !
182 344 subroutine MatMultTranspose_Ising(A,x,y,ierr)
183 #include <petsc/finclude/petscmat.h>
184 use petscmat
185 implicit none
186
187 Mat A
188 Vec x,y
189 PetscInt trans,one,N
190 344 PetscScalar, pointer :: xx(:),yy(:)
191 PetscErrorCode ierr
192
193 ! The actual routine for the transpose matrix-vector product
194 external mvmisg
195
196
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 344 PetscCall(MatGetSize(A,N,PETSC_NULL_INTEGER,ierr))
197
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 344 PetscCall(VecGetArrayRead(x,xx,ierr))
198
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 344 PetscCall(VecGetArray(y,yy,ierr))
199
200 344 trans = 1
201 344 one = 1
202
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)
203
204
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 344 PetscCall(VecRestoreArrayRead(x,xx,ierr))
205
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 344 PetscCall(VecRestoreArray(y,yy,ierr))
206
207 344 end
208
209 ! -------------------------------------------------------------------
210 ! The actual routine for the matrix-vector product
211 ! See https://math.nist.gov/MatrixMarket/data/NEP/mvmisg/mvmisg.html
212
213 868 SUBROUTINE MVMISG(TRANS, N, M, X, LDX, Y, LDY)
214 #include <petsc/finclude/petscsys.h>
215 use petscsys
216 ! ..
217 ! .. Scalar Arguments ..
218 PetscInt LDY, LDX, M, N, TRANS
219 ! ..
220 ! .. Array Arguments ..
221 PetscScalar Y(LDY, *), X(LDX, *)
222 ! ..
223 !
224 ! Purpose
225 ! =======
226 !
227 ! Compute
228 !
229 ! Y(:,1:M) = op(A)*X(:,1:M)
230 !
231 ! where op(A) is A or A' (the transpose of A). The A is the Ising
232 ! matrix.
233 !
234 ! Arguments
235 ! =========
236 !
237 ! TRANS (input) INTEGER
238 ! If TRANS = 0, compute Y(:,1:M) = A*X(:,1:M)
239 ! If TRANS = 1, compute Y(:,1:M) = A'*X(:,1:M)
240 !
241 ! N (input) INTEGER
242 ! The order of the matrix A. N has to be an even number.
243 !
244 ! M (input) INTEGER
245 ! The number of columns of X to multiply.
246 !
247 ! X (input) DOUBLE PRECISION array, dimension (LDX, M)
248 ! X contains the matrix (vectors) X.
249 !
250 ! LDX (input) INTEGER
251 ! The leading dimension of array X, LDX >= max(1, N)
252 !
253 ! Y (output) DOUBLE PRECISION array, dimension (LDX, M)
254 ! contains the product of the matrix op(A) with X.
255 !
256 ! LDY (input) INTEGER
257 ! The leading dimension of array Y, LDY >= max(1, N)
258 !
259 ! ===================================================================
260 !
261 !
262
263 ! .. Local Variables ..
264 PetscInt I, K
265 PetscReal ALPHA, BETA
266 PetscReal COSA, COSB, SINA, SINB
267 PetscScalar TEMP, TEMP1
268 !
269 ! .. Intrinsic functions ..
270 INTRINSIC COS, SIN
271 !
272 868 ALPHA = PETSC_PI/4
273 868 BETA = PETSC_PI/4
274 868 COSA = COS(ALPHA)
275 868 SINA = SIN(ALPHA)
276 868 COSB = COS(BETA)
277 868 SINB = SIN(BETA)
278 !
279
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 868 IF (TRANS.EQ.0) THEN
280 !
281 ! Compute Y(:,1:M) = A*X(:,1:M)
282
283
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 1048 DO 30 K = 1, M
284 !
285 524 Y(1, K) = COSB*X(1, K) - SINB*X(N, K)
286
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 524 DO 10 I = 2, N-1, 2
287 15196 Y(I, K) = COSB*X(I, K) + SINB*X(I+1, K)
288
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)
289 524 10 CONTINUE
290 524 Y(N, K) = SINB*X(1, K) + COSB*X(N, K)
291 !
292
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 524 DO 20 I = 1, N, 2
293 15720 TEMP = COSA*Y(I, K) + SINA*Y(I+1, K)
294 15720 Y(I+1, K) = -SINA*Y(I, K) + COSA*Y(I+1, K)
295
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 15720 Y(I, K) = TEMP
296 524 20 CONTINUE
297 !
298 524 30 CONTINUE
299 !
300
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 344 ELSE IF (TRANS.EQ.1) THEN
301 !
302 ! Compute Y(:1:M) = A'*X(:,1:M)
303 !
304
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 688 DO 60 K = 1, M
305 !
306
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 344 DO 40 I = 1, N, 2
307 10320 Y(I, K) = COSA*X(I, K) - SINA*X(I+1, K)
308
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)
309 344 40 CONTINUE
310 344 TEMP = COSB*Y(1,K) + SINB*Y(N,K)
311
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 344 DO 50 I = 2, N-1, 2
312 9976 TEMP1 = COSB*Y(I, K) - SINB*Y(I+1, K)
313 9976 Y(I+1, K) = SINB*Y(I, K) + COSB*Y(I+1, K)
314
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 9976 Y(I, K) = TEMP1
315 344 50 CONTINUE
316 344 Y(N, K) = -SINB*Y(1, K) + COSB*Y(N, K)
317 344 Y(1, K) = TEMP
318 !
319 344 60 CONTINUE
320 !
321 END IF
322 !
323 868 RETURN
324 !
325 ! END OF MVMISG
326 END
327
328 !/*TEST
329 !
330 ! testset:
331 ! args: -eps_max_it 1000 -eps_ncv 12 -eps_tol 1e-5 -eps_nev 4 -eps_largest_imaginary -terse
332 ! requires: !complex
333 ! output_file: output/ex6f_1.out
334 ! filter: grep -v iterations | sed -e 's/-0.00000/0.00000/g'
335 ! test:
336 ! suffix: 1
337 ! test:
338 ! suffix: 1_ts
339 ! args: -eps_two_sided
340 ! requires: !single
341 !
342 !TEST*/
343