FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
fstr_EIG_lanczos.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 !-------------------------------------------------------------------------------
7 contains
8 
10  subroutine fstr_solve_lanczos(hecMESH, hecMAT, fstrSOLID, fstrEIG)
11  use m_fstr
12  use hecmw_util
13  use m_eigen_lib
16 
17  implicit none
18 
19  type(hecmwst_local_mesh) :: hecMESH
20  type(hecmwst_matrix) :: hecMAT
21  type(fstr_solid) :: fstrSOLID
22  type(fstr_eigen) :: fstrEIG
23  type(fstr_tri_diag) :: Tri
24  type(fstr_eigen_vec), pointer :: Q(:)
25  integer(kind=kint) :: N, NP, NDOF, NNDOF, NPNDOF
26  integer(kind=kint) :: iter, maxiter, nget, ierr
27  integer(kind=kint) :: i, j, k, in, jn, kn, ik, it
28  integer(kind=kint) :: ig, ig0, is0, ie0, its0, ite0
29  real(kind=kreal) :: t1, t2, tolerance
30  real(kind=kreal) :: alpha, beta, beta0
31  real(kind=kreal), allocatable :: s(:), t(:), p(:)
32  logical :: is_converge
33 
34  n = hecmat%N
35  np = hecmat%NP
36  ndof = hecmesh%n_dof
37  nndof = n *ndof
38  npndof = np*ndof
39 
40  allocate(fstreig%filter(npndof))
41  fstreig%filter = 1.0d0
42 
43  jn = 0
44  do ig0 = 1, fstrsolid%BOUNDARY_ngrp_tot
45  ig = fstrsolid%BOUNDARY_ngrp_ID(ig0)
46  is0 = hecmesh%node_group%grp_index(ig-1) + 1
47  ie0 = hecmesh%node_group%grp_index(ig )
48  it = fstrsolid%BOUNDARY_ngrp_type(ig0)
49  its0 = (it - mod(it,10))/10
50  ite0 = mod(it,10)
51 
52  do ik = is0, ie0
53  in = hecmesh%node_group%grp_item(ik)
54  if(ndof < ite0) ite0 = ndof
55  do i = its0, ite0
56  jn = jn + 1
57  fstreig%filter((in-1)*ndof+i) = 0.0d0
58  enddo
59  enddo
60  enddo
61 
62  call hecmw_allreduce_i1(hecmesh, jn, hecmw_sum)
63  if(jn == 0)then
64  fstreig%is_free = .true.
65  if(myrank == 0)then
66  write(*,*) '** free modal analysis: shift factor = 0.1'
67  endif
68  endif
69 
70  call hecmw_update_m_r(hecmesh, fstreig%filter, np, ndof)
71 
72  in = 0
73  do i = 1, nndof
74  if(fstreig%filter(i) == 1.0d0) in = in + 1
75  enddo
76  call hecmw_allreduce_i1(hecmesh, in, hecmw_sum)
77 
78  fstreig%maxiter = fstreig%maxiter + 1
79  if(in < fstreig%maxiter)then
80  if(myrank == 0)then
81  write(imsg,*) '** changed maxiter to system matrix size.'
82  endif
83  fstreig%maxiter = in
84  endif
85 
86  if(in < fstreig%nget)then
87  fstreig%nget = in
88  endif
89 
90  maxiter = fstreig%maxiter
91 
92  allocate(q(0:maxiter))
93  allocate(q(0)%q(npndof))
94  allocate(q(1)%q(npndof))
95  allocate(fstreig%eigval(maxiter))
96  allocate(fstreig%eigvec(npndof, maxiter))
97  allocate(tri%alpha(maxiter))
98  allocate(tri%beta(maxiter))
99  allocate(t(npndof))
100  allocate(s(npndof))
101  allocate(p(npndof))
102 
103  fstreig%eigval = 0.0d0
104  fstreig%eigvec = 0.0d0
105  t = 0.0d0
106  p = 0.0d0
107  s = 0.0d0
108  q(0)%q = 0.0d0
109  q(1)%q = 0.0d0
110  tri%alpha = 0.0d0
111  tri%beta = 0.0d0
112  hecmat%X = 0.0d0
113 
114  call lanczos_set_initial_value(hecmesh, hecmat, fstreig, fstreig%eigvec, p, q(1)%q, tri%beta(1))
115 
116  hecmat%Iarray(98) = 1 !Assmebly complete
117  hecmat%Iarray(97) = 1 !Need numerical factorization
118 
119  if(myrank == 0)then
120  write(imsg,*)
121  write(imsg,*) ' ***** STAGE Begin Lanczos loop **'
122  endif
123 
124  do iter = 1, maxiter-1
126  do i = 1, npndof
127  hecmat%B(i) = p(i)
128  enddo
129 
130  call solve_lineq(hecmesh, hecmat)
131 
132  allocate(q(iter+1)%q(npndof))
133 
134  do i = 1, npndof
135  t(i) = hecmat%X(i) * fstreig%filter(i)
136  enddo
137 
140  do i = 1, npndof
141  t(i) = t(i) - tri%beta(iter) * q(iter-1)%q(i)
142  enddo
143 
144  alpha = 0.0d0
145  do i = 1, nndof
146  alpha = alpha + p(i) * t(i)
147  enddo
148  call hecmw_allreduce_r1(hecmesh, alpha, hecmw_sum)
149  tri%alpha(iter) = alpha
150 
152  do i = 1, npndof
153  t(i) = t(i) - tri%alpha(iter) * q(iter)%q(i)
154  enddo
155 
157  s = 0.0d0
158 
159  do i = 1, npndof
160  s(i) = fstreig%mass(i) * t(i)
161  enddo
162 
163  do j = 0, iter
164  t1 = 0.0d0
165  do i = 1, nndof
166  t1 = t1 + q(j)%q(i) * s(i)
167  enddo
168  call hecmw_allreduce_r1(hecmesh, t1, hecmw_sum)
169  do i = 1, npndof
170  t(i) = t(i) - t1 * q(j)%q(i)
171  enddo
172  enddo
173 
175  do i = 1, npndof
176  s(i) = fstreig%mass(i) * t(i)
177  enddo
178 
179  beta = 0.0d0
180  do i = 1, nndof
181  beta = beta + s(i) * t(i)
182  enddo
183  call hecmw_allreduce_r1(hecmesh, beta, hecmw_sum)
184  tri%beta(iter+1) = dsqrt(beta)
185 
188  beta = 1.0d0/tri%beta(iter+1)
189  do i = 1, npndof
190  p(i) = s(i) * beta
191  q(iter+1)%q(i) = t(i) * beta
192  enddo
193 
194  fstreig%iter = iter
195  if(iter == 1) beta0 = tri%beta(iter+1)
196 
197  call tridiag(hecmesh, hecmat, fstreig, q, tri, iter, is_converge)
198 
199  if(is_converge) exit
200  enddo
201 
202  do i = 0, iter
203  if(associated(q(i)%q)) deallocate(q(i)%q)
204  enddo
205  deallocate(tri%alpha)
206  deallocate(tri%beta)
207  deallocate(t)
208  deallocate(s)
209  deallocate(p)
210 
211  t2 = hecmw_wtime()
212 
213  if(myrank == 0)then
214  write(imsg,*)
215  write(imsg,*) ' * STAGE Output and postprocessing **'
216  write(idbg,'(a,f10.2)') 'Lanczos loop (sec) :', t2 - t1
217  endif
218 
219  end subroutine fstr_solve_lanczos
220 
221 end module m_fstr_eig_lanczos
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=kint), parameter hecmw_sum
integer(kind=4), parameter kreal
real(kind=kreal) function hecmw_wtime()
This modules just summarizes all modules used in eigen analysis.
Definition: eigen_LIB.f90:6
subroutine lanczos_set_initial_value(hecMESH, hecMAT, fstrEIG, eigvec, p, q, beta)
Initialize Lanczos iterations.
Lanczos iteration calculation.
subroutine fstr_solve_lanczos(hecMESH, hecMAT, fstrSOLID, fstrEIG)
SOLVE EIGENVALUE PROBLEM.
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)
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
integer(kind=kint), parameter imsg
Definition: m_fstr.f90:94
integer(kind=kint), parameter idbg
Definition: m_fstr.f90:95
Package of data used by Lanczos eigenvalue solver.
Definition: m_fstr.f90:562