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