1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
5module m_fstr_EIG_lanczos_util
6
7contains
8
9  !> Initialize Lanczos iterations
10  subroutine lanczos_set_initial_value(hecMESH, hecMAT, fstrEIG, eigvec, p, q, beta)
11    use m_fstr
12    use hecmw_util
13    implicit none
14    type(hecmwST_local_mesh) :: hecMESH
15    type(hecmwST_matrix)     :: hecMAT
16    type(fstr_eigen)         :: fstrEIG
17    integer(kind=kint) :: N, NP, NDOF, NDOF2, NNDOF, NPNDOF
18    integer(kind=kint) :: i, j
19    real(kind=kreal)   :: eigvec(:, :), p(:), beta, chk, sigma
20    real(kind=kreal), allocatable :: temp(:)
21    real(kind=kreal), pointer     :: q(:), mass(:), filter(:)
22    logical :: is_free
23
24    N      = hecMAT%N
25    NP     = hecMAT%NP
26    NDOF   = hecMESH%n_dof
27    NDOF2  = NDOF*NDOF
28    NNDOF  = N *NDOF
29    NPNDOF = NP*NDOF
30
31    sigma  = 0.1d0
32    mass   => fstrEIG%mass
33    filter => fstrEIG%filter
34
35    allocate(temp(NNDOF))
36    temp = 0.0d0
37
38    !> shifting
39    if(fstrEIG%is_free)then
40      do i = 1, NP
41        do j = 1, NDOF
42          hecMAT%D(NDOF2*(i-1) + (NDOF+1)*(j-1) + 1) = hecMAT%D(NDOF2*(i-1) + (NDOF+1)*(j-1) + 1) + sigma * mass(NDOF*(i-1) + j)
43        enddo
44      enddo
45    endif
46
47    call URAND1(NNDOF, temp, hecMESH%my_rank)
48
49    do i = 1, NNDOF
50      temp(i) = temp(i) * filter(i)
51    enddo
52
53    !> M-orthogonalization
54    do i = 1, NNDOF
55      eigvec(i,1) = mass(i) * temp(i)
56    enddo
57
58    chk = 0.0d0
59    do i = 1, NNDOF
60      chk = chk + temp(i) * eigvec(i,1)
61    enddo
62    call hecmw_allreduce_R1(hecMESH, chk, hecmw_sum)
63    beta = dsqrt(chk)
64
65    if(beta == 0.0d0)then
66      call hecmw_finalize()
67      stop "Self-orthogonal"
68    endif
69
70    chk = 1.0d0/beta
71    do i = 1, NNDOF
72      q(i) = temp(i) * chk
73    enddo
74
75    do i = 1, NNDOF
76      p(i) = mass(i) * q(i)
77    enddo
78  end subroutine lanczos_set_initial_value
79
80  !> Sort eigenvalues
81  subroutine EVSORT(EIG, NEW, NEIG)
82    use hecmw
83    implicit none
84    integer(kind=kint) :: i, j, n, ip, minloc, NEIG, IBAF, NEW(NEIG)
85    real(kind=kreal) :: EMIN, EIG(NEIG)
86
87    do i = 1, NEIG
88      NEW(i) = i
89    enddo
90
91    n = NEIG-1
92    do i = 1, n
93      minloc = i
94      EMIN = dabs(EIG(NEW(I)))
95      IP = I+1
96      do J=IP,NEIG
97        if(dabs(EIG(NEW(J))).LT.EMIN)then
98          minloc = J
99          EMIN = dabs(EIG(NEW(J)))
100        endif
101      enddo
102      IBAF = NEW(I)
103      NEW(I) = NEW(minloc)
104      NEW(minloc) = IBAF
105    enddo
106  end subroutine EVSORT
107
108  subroutine URAND1(N, X, SHIFT)
109    use hecmw
110    implicit none
111    real(kind=kreal) :: X(N), INVM
112    integer(kind=kint), parameter :: MM = 1664501
113    integer(kind=kint), parameter :: LAMBDA = 1229
114    integer(kind=kint), parameter :: MU = 351750
115    integer(kind=kint) :: i, N, IR, SHIFT
116
117    IR = 0
118    INVM = 1.0D0 / MM
119    do I = 1, SHIFT
120      IR = mod( LAMBDA * IR + MU, MM)
121    enddo
122    do I = SHIFT+1, SHIFT+N
123      IR = mod( LAMBDA * IR + MU, MM)
124      X(I-SHIFT) = INVM * IR
125    enddo
126  end subroutine URAND1
127
128end module m_fstr_EIG_lanczos_util
129
130