Line data Source code
1 :
2 : #if defined HAVE_CONFIG_H
3 : #include "config.h"
4 : #endif
5 : #include "abi_common.h"
6 :
7 : module m_mathfuncs
8 :
9 : use defs_basis, only: dp, PI, std_out
10 : use m_errors
11 : use m_abicore
12 : use m_random_xoroshiro128plus
13 : implicit none
14 :
15 :
16 : ! integer/real
17 : ! change a diagonal (a array) to a matrix
18 : ! or get the diagonal of a 2D matrix
19 : ! similar to matlab diag function
20 : interface diag
21 : procedure diag_mat_int
22 : procedure diag_mat_real
23 : procedure diag_array_int
24 : procedure diag_array_real
25 : end interface diag
26 :
27 : contains
28 :
29 : ! vector cross production
30 865728 : function cross(a, b) result(c)
31 : real(dp), intent(in) :: a(3), b(3)
32 : real(dp) :: c(3)
33 865728 : c(1) = a(2)*b(3) - a(3)*b(2)
34 865728 : c(2) = a(3)*b(1) - a(1)*b(3)
35 865728 : c(3) = a(1)*b(2) - a(2)*b(1)
36 865728 : return
37 : end function cross
38 :
39 :
40 : ! defines outer product of two vectors.
41 0 : function outer_product(a,b) ! result (c)
42 :
43 : real(dp), intent(in) :: a(:), b(:)
44 : real(dp) :: outer_product(size(b, dim=1), size(a, dim=1))
45 : integer:: i, j
46 0 : do i=1, size(a, dim=1)
47 0 : do j=1, size(b, dim=1)
48 0 : outer_product(j,i) = a(i)*b(j)
49 : enddo
50 : enddo
51 0 : return
52 : end function outer_product
53 :
54 : !----------------------------------------------------------------------
55 : !> @brief compare two arrays a<b
56 : !>
57 : !> @param[in] a: integer array
58 : !> @param[in] b: integer array
59 : !> @param[in] N: length of a and b
60 : !> @return bool.
61 : !----------------------------------------------------------------------
62 29097478 : function array_lessthan(a, b, N) result (y)
63 : integer, intent(in) :: a(:), b(:), N
64 : logical :: y
65 : integer :: i
66 29097478 : y=.False.
67 66325592 : do i =1, N
68 66325592 : if (a(i).lt.b(i)) then
69 : y=.True.
70 : exit
71 48697614 : elseif (a(i).gt.b(i)) then
72 : y=.False.
73 : exit
74 : end if
75 : end do
76 29097478 : end function array_lessthan
77 :
78 :
79 : !----------------------------------------------------------------------
80 : !> @brief compare two arrays a>b
81 : !>
82 : !> @param[in] a: integer array
83 : !> @param[in] b: integer array
84 : !> @param[in] N: length of a and b
85 : !> @return bool.
86 : !----------------------------------------------------------------------
87 1278480 : function array_morethan(a, b, N) result (y)
88 : integer, intent(in) :: a(:), b(:), N
89 : logical :: y
90 : integer :: i
91 1278480 : y=.False.
92 1303790 : do i =1, N
93 1303790 : if (a(i).gt.b(i)) then
94 : y=.True.
95 : exit
96 1254295 : elseif (a(i).lt.b(i)) then
97 : y=.False.
98 : exit
99 : end if
100 : end do
101 1278480 : end function array_morethan
102 :
103 : !----------------------------------------------------------------------
104 : !> @brief compare two arrays a <= b
105 : !>
106 : !> @param[in] a: integer array
107 : !> @param[in] b: integer array
108 : !> @param[in] N: length of a and b
109 : !> @return bool.
110 : !----------------------------------------------------------------------
111 15325708 : function array_le(a, b, N) result (y)
112 : integer, intent(in) :: a(:), b(:), N
113 : logical :: y
114 : integer :: i
115 15325708 : y=.True.
116 28412720 : do i =1, N
117 28412720 : if (a(i).gt.b(i)) then
118 : y=.False.
119 : exit
120 20742931 : elseif (a(i).lt.b(i)) then
121 : y=.True.
122 : exit
123 : end if
124 : end do
125 15325708 : end function array_le
126 :
127 :
128 :
129 :
130 : !----------------------------------------------------------------------
131 : !> @brief find an integer from a array
132 : !>
133 : !> @param[in] a: the array to find from
134 : !> @param[in] x: the value to be find
135 : !> @param[out] ix: the index found. 0 if not found.
136 : !----------------------------------------------------------------------
137 0 : function find_int(a, x) result(ix)
138 : integer, intent(in):: a(:), x
139 : integer :: ix, i
140 0 : ix=0
141 0 : do i=1, size(a)
142 0 : if( a(i)==x ) then
143 0 : ix=i
144 : endif
145 : end do
146 0 : end function find_int
147 :
148 : !----------------------------------------------------------------------
149 : !> @brief find an integer from a SORTED array using binary search
150 : !>
151 : !> @param[in] a: the array to find from
152 : !> @param[in] x: the value to be find
153 : !> @param[out] ix: the index found. 0 if not found.
154 : !----------------------------------------------------------------------
155 0 : function binsearch_left_integer(a, x) result(ix)
156 : integer, intent(in):: a(:), x
157 : integer :: n,ix, ub, lb
158 : integer , save :: i=1
159 0 : n=size(a)
160 0 : if (i<0 .or. i>n) i=(size(a)/2+1)
161 0 : if (a(i)==x) then
162 : ix=i
163 : else
164 : ub=n
165 : lb=1
166 0 : do while (lb<ub)
167 0 : i=floor((lb+ub)/2.0)
168 0 : if (a(i)< x) then
169 0 : lb=i+1
170 : else
171 : ub=i
172 : end if
173 : end do
174 0 : if(a(lb)==x) then
175 0 : ix=lb
176 0 : i=ix
177 : else
178 : ix=0
179 : end if
180 : endif
181 0 : end function binsearch_left_integer
182 :
183 : !-------------------------------------------------------------------!
184 : !Binaray search in a interger list.
185 : ! Once it find one match, the index is returned
186 : ! Input:
187 : ! a: the list.
188 : ! x: the element to search for
189 : ! Output:
190 : ! ix: the index of x. If x is not in a, ix =0.
191 : !-------------------------------------------------------------------!
192 3617120 : function binsearch_left_integerlist(a, x) result(ix)
193 : integer, intent(in):: a(:,:), x(:)
194 : integer :: n,ix, ub, lb, nx
195 : integer , save :: i=1
196 3617120 : nx=size(x)
197 3617120 : n=size(a, dim=2)
198 3617120 : if (i<0 .or. i>n) i=(size(a, dim=2)/2+1)
199 10187864 : if (all(a(:,i)==x(:))) then
200 : ix=i
201 : else
202 : ub=n
203 : lb=1
204 32714449 : do while (lb<ub)
205 29097478 : i=floor((lb+ub)/2.0)
206 32714449 : if (array_lessthan(a(:, i), x, nx)) then
207 14026128 : lb=i+1
208 : else
209 : ub=i
210 : end if
211 : end do
212 14467884 : if(all(a(:, lb)==x)) then
213 3616971 : ix=lb
214 3616971 : i=ix
215 : else
216 : ix=0
217 : end if
218 : endif
219 3617120 : end function binsearch_left_integerlist
220 :
221 0 : subroutine set_random_seed(seed)
222 :
223 : integer , intent(in) :: seed(:)
224 : write(std_out,*) "Warning! Currently I'm not sure about how this function &
225 0 : &(set_random_seed,which calls RANDOM_SEED) works. Do test it!"
226 0 : call RANDOM_SEED(put=seed(:))
227 0 : end subroutine set_random_seed
228 :
229 : !----------------------------------------------------------------------
230 : !> @brief get the diagonal of a matrix
231 : !>
232 : !> @param[in] mat: matrix
233 : !> @param[out] ret: the diagonal, a 1d array
234 : !----------------------------------------------------------------------
235 0 : pure function diag_mat_real(mat) result (ret)
236 : real(dp), intent(in) :: mat(:, :)
237 : real(dp):: ret(size(mat, dim=1))
238 : integer :: n, i
239 0 : n=size(mat, dim=1)
240 0 : do i=1, n
241 0 : ret(i)=mat(i,i)
242 : end do
243 0 : end function diag_mat_real
244 :
245 : !----------------------------------------------------------------------
246 : !> @brief get the diagonal of a integer matrix
247 : !>
248 : !> @param[in] mat: matrix
249 : !> @param[out] ret: the diagonal, a 1d array
250 : !----------------------------------------------------------------------
251 0 : pure function diag_mat_int(mat) result (ret)
252 : integer, intent(in) :: mat(:, :)
253 : integer:: ret(size(mat, dim=1))
254 : integer :: n, i
255 0 : n=size(mat, dim=1)
256 0 : do i=1, n
257 0 : ret(i)=mat(i,i)
258 : end do
259 0 : end function diag_mat_int
260 :
261 : !----------------------------------------------------------------------
262 : !> @brief build a matrix from the diagonal (int)
263 : !>
264 : !> @param[in] a: the diagonal array
265 : !> @param[out] the matrix
266 : !----------------------------------------------------------------------
267 0 : pure function diag_array_int(a) result (ret)
268 : integer, intent(in) :: a(:)
269 : integer:: ret(size(a), size(a))
270 : integer :: i
271 0 : ret(:,:)=0
272 0 : do i=1, size(a)
273 0 : ret(i, i)=a(i)
274 : end do
275 0 : end function diag_array_int
276 :
277 : !----------------------------------------------------------------------
278 : !> @brief build a matrix from the diagonal (real)
279 : !>
280 : !> @param[in] a: the diagonal array
281 : !> @param[out] the matrix
282 : !----------------------------------------------------------------------
283 0 : pure function diag_array_real(a) result (ret)
284 : real(dp), intent(in) :: a(:)
285 : real(dp):: ret(size(a), size(a))
286 : integer :: i
287 0 : ret(:,:)=0.0_dp
288 0 : do i=1, size(a)
289 0 : ret(i, i)=a(i)
290 : end do
291 0 : end function diag_array_real
292 :
293 :
294 :
295 : ! Random number generator; Normal (Gaussian) dist.
296 : ! This should NOT be used in practice.
297 : ! There is no standart builtin random number in fortran compiler,
298 : ! which could be a very bad one.
299 : ! also the algorithm used here is not efficient.
300 : ! Only for test purpose
301 0 : subroutine rand_normal_builtin(a)
302 :
303 : real(dp), intent(out)::a(:,:)
304 0 : real(dp), allocatable :: b(:,:)
305 0 : ABI_MALLOC(b, (size(a,dim=1), size(a, dim=2)))
306 0 : call random_number(a)
307 0 : b(:,:) = sqrt(-2*dlog(1.0-a(:,:)))
308 0 : call random_number(a)
309 0 : a(:,:)=b(:,:)*cos(PI*a(:,:))
310 0 : ABI_FREE(b)
311 0 : end subroutine rand_normal_builtin
312 :
313 : !-------------------------------------------------------------------!
314 : ! Shortcut for Hermitian matrix eigen value and eigen vectors.
315 : ! Input:
316 : ! evecs
317 : ! Output:
318 : ! evals: eigen values
319 : ! evecs: eigen vectors
320 : !-------------------------------------------------------------------!
321 0 : subroutine eigensh(evals, evecs)
322 : real(dp), intent(inout) :: evals(:)
323 : complex(dp), intent(inout) :: evecs(:,:)
324 0 : complex(dp), allocatable :: work(:)
325 : integer :: info, lwork
326 0 : real(dp) :: Rwork(3*size(evecs,1)-2)
327 : integer :: ndim
328 : external ZHEEV
329 0 : ndim = size(evecs, 1)
330 0 : lwork = -1
331 0 : ABI_MALLOC(work , (1))
332 0 : call ZHEEV('V', 'U', ndim, evecs, ndim, evals, work, lwork, rwork, info)
333 0 : lwork= INT(work(1))
334 0 : ABI_SFREE(work)
335 :
336 0 : ABI_MALLOC(work , (lwork))
337 0 : call ZHEEV('V', 'U', ndim, evecs, ndim, evals, work, lwork, rwork, info)
338 0 : ABI_SFREE(work)
339 0 : IF( INFO.gt.0 ) THEN
340 0 : ABI_ERROR('The zheev algorithm failed to compute eigenvalues.')
341 0 : STOP
342 : END IF
343 0 : end subroutine eigensh
344 :
345 : !-----------------------------------------------------------------------
346 : !> @brief rotate vector around axis by angle
347 : !> Using the quaternion rotation algorithm
348 : !> \vec{v}_new = \vec{v} + 2 \vec{r} \cross (\vec{r}\cross\vec{v}
349 : !> + w \vec{v})
350 : !> see https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation
351 : !> @param [in] angle : angle
352 : !> @param [in] axis : axis vector. The norm of axis does not matter.
353 : !> @param [in] vec : vector to roate
354 : !> @param [out] vec2 : result vector
355 : !-----------------------------------------------------------------------
356 432 : function rotate_by_angle_around_axis(angle, axis, vec) result(vec2)
357 : real(dp), intent(in) :: angle, axis(3), vec(3)
358 : real(dp) :: vec2(3)
359 : real(dp) :: half_angle, r(3), w, norm
360 432 : half_angle=angle/2.0_dp
361 1728 : norm=norm2(axis)
362 1728 : r(:)=axis(:)/norm * sin(half_angle)
363 432 : w=cos(half_angle)
364 : ! (w, r) is the quaternion
365 3456 : vec2(:) = vec(:) + 2.0 * cross(r, (cross(r, vec) + w*vec))
366 432 : end function rotate_by_angle_around_axis
367 :
368 : end module m_mathfuncs
|
