A Modular Arbitrary-Order Ocean-Atmosphere Model -- Stochastic implementation
util.f90
Go to the documentation of this file.
1 
2 ! util.f90
3 !
4 !> Utility module
5 !
6 !> @copyright
7 !> 2018 Lesley De Cruz & Jonathan Demaeyer.
8 !> See LICENSE.txt for license information.
9 !
10 !---------------------------------------------------------------------------!
11 
12 MODULE util
13  IMPLICIT NONE
14 
15  PRIVATE
16 
18  PUBLIC :: mat_trace, mat_contract, choldc
20  PUBLIC :: vector_outer
21 
22 CONTAINS
23 
24  ! SUBROUTINE scalar_allocate(x)
25  ! INTEGER :: AllocStat
26  ! IF (.NOT. ALLOCATED(x)) THEN
27  ! ALLOCATE(x, STAT=AllocStat)
28 
29  !> Convert an integer to string.
30  CHARACTER(len=20) FUNCTION str(k)
31  INTEGER, INTENT(IN) :: k
32  WRITE (str, *) k
33  str = adjustl(str)
34  END FUNCTION str
35 
36  !> Convert a real to string with a given format
37  CHARACTER(len=40) FUNCTION rstr(x,fm)
38  REAL(KIND=8), INTENT(IN) :: x
39  CHARACTER(len=20), INTENT(IN) :: fm
40  WRITE (rstr, trim(adjustl(fm))) x
41  rstr = adjustl(rstr)
42  END FUNCTION rstr
43 
44  !> Determine if a character is in a string and where
45  !> @remark: return positions in a vector if found and 0 vector if not found
46  FUNCTION isin(c,s)
47  CHARACTER, INTENT(IN) :: c
48  CHARACTER, DIMENSION(:), INTENT(IN) :: s
49  INTEGER, DIMENSION(size(s)) :: isin
50  INTEGER :: i,j
51 
52  isin=0
53  j=0
54  DO i=size(s),1,-1
55  IF (c==s(i)) THEN
56  j=j+1
57  isin(j)=i
58  END IF
59  END DO
60  END FUNCTION isin
61 
62  !> Random generator initialization routine
63  SUBROUTINE init_random_seed()
64  USE iso_fortran_env, only: int64
65  USE ifport !, only: getpid
66  IMPLICIT NONE
67  INTEGER, ALLOCATABLE :: seed_loc(:)
68  INTEGER :: i, n, un, istat, dt(8), pid
69  INTEGER(int64) :: t
70 
71  CALL random_seed(size = n)
72  ALLOCATE(seed_loc(n))
73  ! First try IF the OS provides a random number generator
74  OPEN(newunit=un, file="/dev/urandom", access="stream", &
75  form="unformatted", action="read", status="old", iostat=istat)
76  IF (istat == 0) THEN
77  READ(un) seed_loc
78  CLOSE(un)
79  ELSE
80  ! Fallback to XOR:ing the current time and pid. The PID is
81  ! useful in case one launches multiple instances of the same
82  ! program in parallel.
83  CALL system_clock(t)
84  IF (t == 0) THEN
85  CALL date_and_time(values=dt)
86  t = (dt(1) - 1970) * 365_int64 * 24 * 60 * 60 * 1000 &
87  + dt(2) * 31_int64 * 24 * 60 * 60 * 1000 &
88  + dt(3) * 24_int64 * 60 * 60 * 1000 &
89  + dt(5) * 60 * 60 * 1000 &
90  + dt(6) * 60 * 1000 + dt(7) * 1000 &
91  + dt(8)
92  END IF
93  pid = getpid()
94  t = ieor(t, int(pid, kind(t)))
95  DO i = 1, n
96  seed_loc(i) = lcg(t)
97  END DO
98  END IF
99  CALL random_seed(put=seed_loc)
100  contains
101  ! This simple PRNG might not be good enough for real work, but is
102  ! sufficient for seeding a better PRNG.
103  FUNCTION lcg(s)
104  integer :: lcg
105  integer(int64) :: s
106  IF (s == 0) THEN
107  s = 104729
108  ELSE
109  s = mod(s, 4294967296_int64)
110  END IF
111  s = mod(s * 279470273_int64, 4294967291_int64)
112  lcg = int(mod(s, int(huge(0), int64)), kind(0))
113  END FUNCTION lcg
114  END SUBROUTINE init_random_seed
115 
116  !> Simple card player sorting function
117  SUBROUTINE piksrt(k,arr,par)
118  INTEGER, INTENT(IN) :: k
119  INTEGER, DIMENSION(k), INTENT(INOUT) :: arr
120  INTEGER, INTENT(OUT) :: par
121  INTEGER :: i,j,a,b
122 
123  par=1
124 
125  DO j=2,k
126  a=arr(j)
127  DO i=j-1,1,-1
128  if(arr(i).le.a) EXIT
129  arr(i+1)=arr(i)
130  par=-par
131  END DO
132  arr(i+1)=a
133  ENDDO
134  RETURN
135  END SUBROUTINE piksrt
136 
137  !> Initialize a square matrix A as a unit matrix
138  SUBROUTINE init_one(A)
139  REAL(KIND=8), DIMENSION(:,:),INTENT(INOUT) :: A
140  INTEGER :: i,n
141  n=size(a,1)
142  a=0.0d0
143  DO i=1,n
144  a(i,i)=1.0d0
145  END DO
146 
147  END SUBROUTINE init_one
148 
149  FUNCTION mat_trace(A)
150  REAL(KIND=8), DIMENSION(:,:) :: A
151  REAL(KIND=8) :: mat_trace
152  INTEGER :: i,n
153  n=size(a,1)
154  mat_trace=0.d0
155  DO i=1,n
156  mat_trace=mat_trace+a(i,i)
157  END DO
158  RETURN
159  END FUNCTION mat_trace
160 
161  FUNCTION mat_contract(A,B)
162  REAL(KIND=8), DIMENSION(:,:) :: A,B
163  REAL(KIND=8) :: mat_contract
164  INTEGER :: i,j,n
165  n=size(a,1)
166  mat_contract=0.d0
167  DO i=1,n
168  DO j=1,n
169  mat_contract=mat_contract+a(i,j)*b(i,j)
170  END DO
171  ENDDO
172  RETURN
173  END FUNCTION mat_contract
174 
175  SUBROUTINE choldc(a,p)
176  REAL(KIND=8), DIMENSION(:,:) :: a
177  REAL(KIND=8), DIMENSION(:) :: p
178  INTEGER :: n
179  INTEGER :: i,j,k
180  REAL(KIND=8) :: sum
181  n=size(a,1)
182  DO i=1,n
183  DO j=i,n
184  sum=a(i,j)
185  DO k=i-1,1,-1
186  sum=sum-a(i,k)*a(j,k)
187  END DO
188  IF (i.eq.j) THEN
189  IF (sum.le.0.) stop 'choldc failed'
190  p(i)=sqrt(sum)
191  ELSE
192  a(j,i)=sum/p(i)
193  ENDIF
194  END DO
195  END DO
196  RETURN
197  END SUBROUTINE choldc
198 
199  SUBROUTINE printmat(A) ! to be moved to util
200  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: A
201  INTEGER :: i
202 
203  DO i=1,SIZE(a,1)
204  print*, a(i,:)
205  END DO
206  END SUBROUTINE printmat
207 
208  SUBROUTINE cprintmat(A) ! to be moved to util
209  COMPLEX(KIND=16), DIMENSION(:,:), INTENT(IN) :: A
210  INTEGER :: i
211 
212  DO i=1,SIZE(a,1)
213  print*, a(i,:)
214  END DO
215  END SUBROUTINE cprintmat
216 
217  FUNCTION invmat(A) RESULT(Ainv)
218  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: A
219  REAL(KIND=8), DIMENSION(SIZE(A,1),SIZE(A,2)) :: Ainv
220 
221  REAL(KIND=8), DIMENSION(SIZE(A,1)) :: work ! work array for LAPACK
222  INTEGER, DIMENSION(SIZE(A,1)) :: ipiv ! pivot indices
223  INTEGER :: n, info
224 
225  ! Store A in Ainv to prevent it from being overwritten by LAPACK
226  ainv = a
227  n = size(a,1)
228 
229  ! DGETRF computes an LU factorization of a general M-by-N matrix A
230  ! using partial pivoting with row interchanges.
231  CALL dgetrf(n, n, ainv, n, ipiv, info)
232 
233  IF (info /= 0) THEN
234  stop 'Matrix is numerically singular!'
235  ENDIF
236 
237  ! DGETRI computes the inverse of a matrix using the LU factorization
238  ! computed by DGETRF.
239  CALL dgetri(n, ainv, n, ipiv, work, n, info)
240 
241  IF (info /= 0) THEN
242  stop 'Matrix inversion failed!'
243  ENDIF
244  END FUNCTION invmat
245 
246  SUBROUTINE triu(A,T)
247  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: A
248  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: T
249  INTEGER i,j
250  t=0.d0
251  DO i=1,SIZE(a,1)
252  DO j=i,SIZE(a,1)
253  t(i,j)=a(i,j)
254  END DO
255  END DO
256  END SUBROUTINE triu
257 
258  SUBROUTINE diag(A,d)
259  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: A
260  REAL(KIND=8), DIMENSION(:), INTENT(OUT) :: d
261  INTEGER :: i
262 
263  DO i=1,SIZE(a,1)
264  d(i)=a(i,i)
265  END DO
266  END SUBROUTINE diag
267 
268  SUBROUTINE cdiag(A,d)
269  COMPLEX(KIND=16), DIMENSION(:,:), INTENT(IN) :: A
270  COMPLEX(KIND=16), DIMENSION(:), INTENT(OUT) :: d
271  INTEGER :: i
272 
273  DO i=1,SIZE(a,1)
274  d(i)=a(i,i)
275  END DO
276  END SUBROUTINE cdiag
277 
278 
279  FUNCTION floordiv(i,j)
280  INTEGER :: i,j,floordiv
281  floordiv=int(floor(real(i)/real(j)))
282  RETURN
283  END FUNCTION floordiv
284 
285  SUBROUTINE reduce(A,Ared,n,ind,rind)
286  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: A
287  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: Ared
288  INTEGER, INTENT(OUT) :: n
289  INTEGER, DIMENSION(:), INTENT(OUT) :: ind,rind
290  LOGICAL, DIMENSION(SIZE(A,1)) :: sel
291  INTEGER :: i,j
292 
293  ind=0
294  rind=0
295  sel=.false.
296  n=0
297  DO i=1,SIZE(a,1)
298  IF (any(a(i,:)/=0)) THEN
299  n=n+1
300  sel(i)=.true.
301  ind(n)=i
302  rind(i)=n
303  ENDIF
304  END DO
305  ared=0.d0
306  DO i=1,SIZE(a,1)
307  DO j=1,SIZE(a,1)
308  IF (sel(i).and.sel(j)) ared(rind(i),rind(j))=a(i,j)
309  ENDDO
310  ENDDO
311  END SUBROUTINE reduce
312 
313  SUBROUTINE ireduce(A,Ared,n,ind,rind)
314  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: A
315  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: Ared
316  INTEGER, INTENT(IN) :: n
317  INTEGER, DIMENSION(:), INTENT(IN) :: ind,rind
318  INTEGER :: i,j
319  a=0.d0
320  DO i=1,n
321  DO j=1,n
322  a(ind(i),ind(j))=ared(i,j)
323  END DO
324  END DO
325  END SUBROUTINE ireduce
326 
327  SUBROUTINE vector_outer(u,v,A)
328  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: u,v
329  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: A
330  INTEGER :: i,j
331 
332  a=0.d0
333  DO i=1,SIZE(u)
334  DO j=1,SIZE(v)
335  a(i,j)=u(i)*v(j)
336  ENDDO
337  ENDDO
338  END SUBROUTINE vector_outer
339 
340 END MODULE util
subroutine, public choldc(a, p)
Definition: util.f90:176
Utility module.
Definition: util.f90:12
subroutine, public init_random_seed()
Random generator initialization routine.
Definition: util.f90:64
integer function, public floordiv(i, j)
Definition: util.f90:280
subroutine, public cdiag(A, d)
Definition: util.f90:269
subroutine, public piksrt(k, arr, par)
Simple card player sorting function.
Definition: util.f90:118
character(len=40) function, public rstr(x, fm)
Convert a real to string with a given format.
Definition: util.f90:38
integer function lcg(s)
Definition: util.f90:104
subroutine, public printmat(A)
Definition: util.f90:200
subroutine, public diag(A, d)
Definition: util.f90:259
subroutine, public ireduce(A, Ared, n, ind, rind)
Definition: util.f90:314
real(kind=8) function, public mat_trace(A)
Definition: util.f90:150
character(len=20) function, public str(k)
Convert an integer to string.
Definition: util.f90:31
subroutine, public vector_outer(u, v, A)
Definition: util.f90:328
real(kind=8) function, public mat_contract(A, B)
Definition: util.f90:162
integer function, dimension(size(s)), public isin(c, s)
Determine if a character is in a string and where.
Definition: util.f90:47
subroutine, public reduce(A, Ared, n, ind, rind)
Definition: util.f90:286
subroutine, public triu(A, T)
Definition: util.f90:247
subroutine, public init_one(A)
Initialize a square matrix A as a unit matrix.
Definition: util.f90:139
real(kind=8) function, dimension(size(a, 1), size(a, 2)), public invmat(A)
Definition: util.f90:218
subroutine, public cprintmat(A)
Definition: util.f90:209