Actual source code: ex6f.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> ./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>
23: module ex6fmodule
24: use slepceps
25: implicit none
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: subroutine mvmisg(trans, n, m, x, ldx, y, ldy)
58: use petscsys
59: implicit none
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
67: alpha = PETSC_PI/4
68: beta = PETSC_PI/4
69: cosa = cos(alpha)
70: sina = sin(alpha)
71: cosb = cos(beta)
72: sinb = sin(beta)
74: if (trans == 0) then
76: ! ** Compute y(:,1:m) = A*x(:,1:m)
78: do k = 1, m
79: y(1, k) = cosb*x(1, k) - sinb*x(n, k)
80: do i = 2, n - 1, 2
81: y(i, k) = cosb*x(i, k) + sinb*x(i + 1, k)
82: y(i + 1, k) = -sinb*x(i, k) + cosb*x(i + 1, k)
83: end do
84: y(n, k) = sinb*x(1, k) + cosb*x(n, k)
85: do i = 1, n, 2
86: temp = cosa*y(i, k) + sina*y(i + 1, k)
87: y(i + 1, k) = -sina*y(i, k) + cosa*y(i + 1, k)
88: y(i, k) = temp
89: end do
90: end do
92: else if (trans == 1) then
94: ! ** Compute y(:1:m) = A'*x(:,1:m)
96: do k = 1, m
97: do i = 1, n, 2
98: y(i, k) = cosa*x(i, k) - sina*x(i + 1, k)
99: y(i + 1, k) = sina*x(i, k) + cosa*x(i + 1, k)
100: end do
101: temp = cosb*y(1, k) + sinb*y(n, k)
102: do i = 2, n - 1, 2
103: temp1 = cosb*y(i, k) - sinb*y(i + 1, k)
104: y(i + 1, k) = sinb*y(i, k) + cosb*y(i + 1, k)
105: y(i, k) = temp1
106: end do
107: y(n, k) = -sinb*y(1, k) + cosb*y(n, k)
108: y(1, k) = temp
109: end do
111: end if
112: end subroutine
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: subroutine MatMult_Ising(A, x, y, ierr)
125: use petscmat
126: implicit none
128: Mat :: A
129: Vec :: x, y
130: PetscInt :: trans, N
131: PetscScalar, pointer :: xx(:), yy(:)
132: PetscErrorCode, intent(out) :: ierr
134: PetscCall(MatGetSize(A, N, PETSC_NULL_INTEGER, ierr))
135: PetscCall(VecGetArrayRead(x, xx, ierr))
136: PetscCall(VecGetArray(y, yy, ierr))
138: trans = 0
139: call mvmisg(trans, N, 1_PETSC_INT_KIND, xx, N, yy, N)
141: PetscCall(VecRestoreArrayRead(x, xx, ierr))
142: PetscCall(VecRestoreArray(y, yy, ierr))
143: end subroutine
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: subroutine MatMultTranspose_Ising(A, x, y, ierr)
156: use petscmat
157: implicit none
159: Mat :: A
160: Vec :: x, y
161: PetscInt :: trans, N
162: PetscScalar, pointer :: xx(:), yy(:)
163: PetscErrorCode, intent(out) :: ierr
165: PetscCall(MatGetSize(A, N, PETSC_NULL_INTEGER, ierr))
166: PetscCall(VecGetArrayRead(x, xx, ierr))
167: PetscCall(VecGetArray(y, yy, ierr))
169: trans = 1
170: call mvmisg(trans, N, 1_PETSC_INT_KIND, xx, N, yy, N)
172: PetscCall(VecRestoreArrayRead(x, xx, ierr))
173: PetscCall(VecRestoreArray(y, yy, ierr))
174: end subroutine
176: end module ex6fmodule
178: program ex6f
179: use slepceps
180: use ex6fmodule
181: implicit none
183: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184: ! Declarations
185: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
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
197: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198: ! Beginning of program
199: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
201: 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: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, sz, ierr))
206: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
207: PetscCheckA(sz == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, 'This is a uniprocessor example only!')
208: m = 30
209: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-m', m, flg, ierr))
210: N = 2*m
212: if (rank == 0) then
213: write (*, '(/a,i6,a/)') 'Ising Model Eigenproblem, m=', m, ', (N=2*m)'
214: end if
216: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
217: ! Register the matrix-vector subroutine for the operator that defines
218: ! the eigensystem, Ax=kx
219: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
221: PetscCallA(MatCreateShell(PETSC_COMM_WORLD, N, N, N, N, PETSC_NULL_INTEGER, A, ierr))
222: PetscCallA(MatShellSetOperation(A, MATOP_MULT, MatMult_Ising, ierr))
223: PetscCallA(MatShellSetOperation(A, MATOP_MULT_TRANSPOSE, MatMultTranspose_Ising, ierr))
225: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
226: ! Create the eigensolver and display info
227: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
229: ! ** Create eigensolver context
230: PetscCallA(EPSCreate(PETSC_COMM_WORLD, eps, ierr))
232: ! ** Set operators. In this case, it is a standard eigenvalue problem
233: PetscCallA(EPSSetOperators(eps, A, PETSC_NULL_MAT, ierr))
234: PetscCallA(EPSSetProblemType(eps, EPS_NHEP, ierr))
236: ! ** Set solver parameters at runtime
237: PetscCallA(EPSSetFromOptions(eps, ierr))
239: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
240: ! Solve the eigensystem
241: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
243: PetscCallA(EPSSolve(eps, ierr))
244: PetscCallA(EPSGetIterationNumber(eps, its, ierr))
245: if (rank == 0) then
246: write (*, '(a,i4)') ' Number of iterations of the method: ', its
247: end if
249: ! ** Optional: Get some information from the solver and display it
250: PetscCallA(EPSGetType(eps, tname, ierr))
251: if (rank == 0) then
252: write (*, '(a,a)') ' Solution method: ', tname
253: end if
254: PetscCallA(EPSGetDimensions(eps, nev, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
255: if (rank == 0) then
256: write (*, '(a,i2)') ' Number of requested eigenvalues:', nev
257: end if
258: PetscCallA(EPSGetTolerances(eps, tol, maxit, ierr))
259: if (rank == 0) then
260: write (*, '(a,1pe11.4,a,i6)') ' Stopping condition: tol=', tol, ', maxit=', maxit
261: end if
263: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
264: ! Display solution and clean up
265: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
267: ! ** show detailed info unless -terse option is given by user
268: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-terse', terse, ierr))
269: if (terse) then
270: 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: PetscCallA(EPSDestroy(eps, ierr))
278: PetscCallA(MatDestroy(A, ierr))
279: PetscCallA(SlepcFinalize(ierr))
280: end program ex6f
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*/