FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
fstr_EIG_tridiag.f90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
8 use hecmw
9
10 implicit none
11
12 public
13
15 real(kind=kreal), pointer :: q(:) => null()
16 end type fstr_eigen_vec
17
19 real(kind=kreal), allocatable :: alpha(:)
20 real(kind=kreal), allocatable :: beta(:)
21 end type fstr_tri_diag
22
23contains
24
25 subroutine tridiag(hecMESH, hecMAT, fstrEIG, Q, Tri, iter, is_converge)
26 use hecmw
27 use m_fstr
29 implicit none
30 type(hecmwst_local_mesh) :: hecmesh
31 type(hecmwst_matrix) :: hecMAT
32 type(fstr_eigen) :: fstrEIG
33 type(fstr_tri_diag) :: Tri
34 type(fstr_eigen_vec), pointer :: Q(:)
35
36 integer(kind=kint), intent(in) :: iter
37 integer(kind=kint) :: N, NP, NDOF, NNDOF, NPNDOF
38 integer(kind=kint) :: i, j, k, in, jn, kn, nget
39 integer(kind=kint) :: iter2, ierr, maxiter
40 real(kind=kreal) :: resid, chk, sigma, tolerance, max_eigval
41 real(kind=kreal), allocatable :: alpha(:), beta(:), temp(:)
42 real(kind=kreal), allocatable :: l(:,:)
43
44 integer(kind=kint), allocatable :: iparm(:)
45 real(kind=kreal), pointer :: eigvec(:,:)
46 real(kind=kreal), pointer :: eigval(:)
47 logical :: is_converge
48
49 n = hecmat%N
50 np = hecmat%NP
51 ndof = hecmesh%n_dof
52 nndof = n *ndof
53 npndof = np*ndof
54 eigval => fstreig%eigval
55 eigvec => fstreig%eigvec
56 nget = fstreig%nget
57 maxiter = fstreig%maxiter
58 tolerance = fstreig%tolerance
59
60 allocate( iparm(maxiter) )
61 allocate( temp(maxiter) )
62 allocate( alpha(iter) )
63 allocate( beta(iter) )
64 allocate( l(iter, iter) )
65
66 do j=1, iter
67 do i = 1,iter
68 l(i,j) = 0.0d0
69 enddo
70 enddo
71
72 do i=1, iter
73 alpha(i) = tri%alpha(i)
74 l(i,i) = 1.0d0
75 enddo
76
77 beta(1) = 0.0d0
78 do i=2, iter
79 beta(i) = tri%beta(i)
80 enddo
81
82 call ql_decomposition(iter, iter, alpha, beta, l, ierr)
83
84 is_converge = .true.
85 chk = 0.0d0
86 max_eigval = maxval(alpha)
87 do i = 1, min(nget, iter)
88 resid = dabs(tri%beta(iter+1)*l(iter,i))/max_eigval
89 chk = max(chk, resid)
90 if(tolerance < resid) is_converge = .false.
91 enddo
92 if(myrank == 0) write(*,"(i8,1pe12.5)")iter, chk
93
94 if(iter < nget) is_converge = .false.
95 if(iter == maxiter-1) is_converge = .true.
96
97 if(is_converge)then
98 sigma = 0.0d0
99 if(fstreig%is_free) sigma = 0.1d0
100 do i = 1, iter
101 if(alpha(i) /= 0.0d0)then
102 eigval(i) = 1.0d0/alpha(i) + sigma
103 endif
104 enddo
105
106 call evsort(eigval, iparm, iter)
107
108 temp = eigval
109
110 eigvec = 0.0d0
111 do k=1, iter
112 in = iparm(k)
113 eigval(k) = temp(in)
114 do j=1, iter
115 do i=1, npndof
116 eigvec(i, k) = eigvec(i, k) + q(j)%q(i) * l(j, in)
117 enddo
118 enddo
119 enddo
120
121 do j=1, iter
122 chk = maxval(eigvec(:,j))
123 call hecmw_allreduce_r1(hecmesh, chk, hecmw_max)
124 if(chk /= 0.0d0)then
125 chk = 1.0d0/chk
126 do i = 1, nndof
127 eigvec(i,j) = eigvec(i,j) * chk
128 enddo
129 endif
130 enddo
131 endif
132
133 deallocate(iparm)
134 deallocate(temp)
135 deallocate(alpha)
136 deallocate(beta)
137 deallocate(l)
138 end subroutine tridiag
139
140 !======================================================================!
141 ! Description !
142 !======================================================================!
152 !
153 !on input
154 !
155 !nm must be set to the row dimension of two-dimensional
156 !array parameters as declared in the calling program
157 !dimension statement.
158 !
159 !is the order of the matrix.
160 !
161 !contains the diagonal elements of the input matrix.
162 !
163 !contains the subdiagonal elements of the input matrix
164 !in its last n-1 positions. e(1) is arbitrary.
165 !
166 !contains the transformation matrix produced in the
167 !reduction by tred2, if performed. if the eigenvectors
168 !of the tridiagonal matrix are desired, z must contain
169 !the identity matrix.
170 !
171 !on output
172 !
173 !d
174 !contains the eigenvalues in ascending order. if an
175 !error exit is made, the eigenvalues are correct but
176 !unordered for indices 1,2,...,ierror-1.
177 !-----------------------------------------------------
178 !GP
179 !du
180 !contains the unordered eigenvalues. if an
181 !error exit is made, the eigenvalues are correct and
182 !unordered for indices 1,2,...,ierror-1.
183 !-----------------------------------------------------
184 !e
185 !has been destroyed.
186 !
187 !z
188 !contains orthonormal eigenvectors of the symmetric
189 !tridiagonal (or full) matrix. if an error exit is made,
190 !z contains the eigenvectors associated with the stored
191 !eigenvalues.
192 !-----------------------------------------------------
193 !zu
194 !contains unordered eigenvectors of the symm. tridiag. matrix
195 !-----------------------------------------------------
196 !ierror is set to
197 ! zero for normal return,
198 ! j if the j-th eigenvalue has not been
199 ! determined after 30 iterations.
200 !
201 !calls a2b2 for dsqrt(a*a + b*b) .
202 !=======================================================================
203
204 !call QL_decomposition(iter, iter, alpha, beta, L, ierr)
205 subroutine ql_decomposition(nm, n, d, e, z, ierror)
206 use hecmw
207 implicit none
208 integer(kind=kint) :: i, j, k, l, m, n, ii, l1, l2, nm, mml, ierror
209 real(kind=kreal) :: d(n), e(n), z(nm, n)
210 real(kind=kreal) :: c, c2, c3, dl1, el1, f, g, h, p, r, s, s2, tst1, tst2
211
212 ierror = 0
213 if (n .eq. 1) go to 1001
214
215 do i = 2, n
216 e(i-1) = e(i)
217 enddo
218
219 f = 0.0d0
220 tst1 = 0.0d0
221 e(n) = 0.0d0
222
223 do 240 l = 1, n
224 j = 0
225 h = dabs(d(l)) + dabs(e(l))
226 if (tst1 .lt. h) tst1 = h
227 ! .......... look for small sub-diagonal element ..........
228 bb:do m = l, n
229 tst2 = tst1 + dabs(e(m))
230 if (tst2 .eq. tst1) exit bb
231 ! .......... e(n) is always zero, so there is no exit
232 ! through the bottom of the loop ..........
233 enddo bb
234
235 if (m .eq. l) go to 220
236
237 130 if (j .eq. 30) go to 1000
238 j = j + 1
239 ! .......... form shift ..........
240 l1 = l + 1
241 l2 = l1 + 1
242 g = d(l)
243 p = (d(l1) - g) / (2.0d0 * e(l))
244 r = a2b2(p,1.0d0)
245 d(l) = e(l) / (p + dsign(r,p))
246 d(l1) = e(l) * (p + dsign(r,p))
247 dl1 = d(l1)
248 h = g - d(l)
249 if (l2 .gt. n) go to 145
250
251 do i = l2, n
252 d(i) = d(i) - h
253 enddo
254
255 145 f = f + h
256 ! .......... ql transformation ..........
257 p = d(m)
258 c = 1.0d0
259 c2 = c
260 el1 = e(l1)
261 s = 0.0d0
262 s2 = 0.0d0
263 c3 = 0.0d0
264 mml = m - l
265 ! .......... for i=m-1 step -1 until l do -- ..........
266 do ii = 1, mml
267 c3 = c2
268 c2 = c
269 s2 = s
270 i = m - ii
271 g = c * e(i)
272 h = c * p
273 r = a2b2(p,e(i))
274 e(i+1) = s * r
275 s = e(i) / r
276 c = p / r
277 p = c * d(i) - s * g
278 d(i+1) = h + s * (c * g + s * d(i))
279 ! .......... form vector ..........
280 do k = 1, n
281 h = z(k,i+1)
282 z(k,i+1) = s * z(k,i) + c * h
283 z(k,i) = c * z(k,i) - s * h
284 enddo
285 enddo
286
287 p = -s * s2 * c3 * el1 * e(l) / dl1
288 e(l) = s * p
289 d(l) = c * p
290 tst2 = tst1 + dabs(e(l))
291 if (tst2 .gt. tst1) go to 130
292 220 d(l) = d(l) + f
293 240 continue
294
295 ! .......... order eigenvalues and eigenvectors ..........
296 !GP: Get unordered eigenvalues and eigenvectors----------------
297
298 do 300 ii = 2, n
299 i = ii - 1
300 k = i
301 p = d(i)
302
303 aa:do j = ii, n
304 if (d(j) .ge. p) exit aa
305 k = j
306 p = d(j)
307 enddo aa
308
309 if (k .eq. i) go to 300
310 d(k) = d(i)
311 d(i) = p
312
313 do j = 1, n
314 p = z(j,i)
315 z(j,i) = z(j,k)
316 z(j,k) = p
317 enddo
318
319 300 continue
320
321 go to 1001
322 ! .......... set error -- no convergence to an
323 ! eigenvalue after 30 iterations ..........
324 1000 ierror = l
325 1001 return
326 end subroutine ql_decomposition
327
328 function a2b2(a,b)
329 use hecmw
330 implicit none
331
332 real(kind=kreal) :: a2b2
333 real(kind=kreal) :: a, b
334 real(kind=kreal) :: p, q, r, s, t, u
335
336 p = dmax1(dabs(a), dabs(b))
337 if (p /= 0.0d0) then
338 r = (dmin1(dabs(a),dabs(b))/p) ** 2
339 do
340 t = 4.0d0 + r
341 if (t == 4.0d0) exit
342 s = r/t
343 u = 1.0d0 + 2.0d0*s
344 p = u * p
345 q = s/u
346 r = q * q * r
347 end do
348 end if
349 a2b2 = p
350 return
351
352 end function a2b2
353
354end module m_fstr_eig_tridiag
Definition: hecmw.f90:6
subroutine evsort(eig, new, neig)
Sort eigenvalues.
This module provides a subroutine to find the eigenvalues and eigenvectors of a symmetric tridiagonal...
subroutine tridiag(hecmesh, hecmat, fstreig, q, tri, iter, is_converge)
subroutine ql_decomposition(nm, n, d, e, z, ierror)
This subroutine has been adapted from the eispack routine tql2, which is a translation of the algol p...
real(kind=kreal) function a2b2(a, b)
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
integer(kind=kint) myrank
PARALLEL EXECUTION.
Definition: m_fstr.f90:80
Package of data used by Lanczos eigenvalue solver.
Definition: m_fstr.f90:562