Branch data 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
|
