EMPIRE DA  v1.9.1
Data assimilation codes using EMPIRE communication
 All Classes Files Functions Variables Pages
tests.f90
Go to the documentation of this file.
1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 !!! Time-stamp: <2016-08-15 19:32:47 pbrowne>
3 !!!
4 !!! Collection of subroutines to perform checks of user supplied routines
5 !!! Copyright (C) 2014 Philip A. Browne
6 !!!
7 !!! This program is free software: you can redistribute it and/or modify
8 !!! it under the terms of the GNU General Public License as published by
9 !!! the Free Software Foundation, either version 3 of the License, or
10 !!! (at your option) any later version.
11 !!!
12 !!! This program is distributed in the hope that it will be useful,
13 !!! but WITHOUT ANY WARRANTY; without even the implied warranty of
14 !!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 !!! GNU General Public License for more details.
16 !!!
17 !!! You should have received a copy of the GNU General Public License
18 !!! along with this program. If not, see <http://www.gnu.org/licenses/>.
19 !!!
20 !!! Email: p.browne @ reading.ac.uk
21 !!! Mail: School of Mathematical and Physical Sciences,
22 !!! University of Reading,
23 !!! Reading, UK
24 !!! RG6 6BB
25 !!!
26 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
27 
28 !!$ it turns out that these tests don't make any sense for H...
29 !!commenting out
30 !!$subroutine H_tests()
31 !!$ !> @brief These are some tests to check that the observation
32 !!$ !> operator is implemented correctly
33 !!$ use pf_control
34 !!$ use sizes
35 !!$ implicit none
36 !!$ integer, parameter :: rk = kind(1.0d0)
37 !!$ real(kind=rk), dimension(obs_dim) :: a,s
38 !!$ real(kind=rk), dimension(obs_dim,10) :: aa,ss
39 !!$ real(kind=rk), dimension(state_dim) :: j
40 !!$ real(kind=rk), dimension(state_dim,10) :: jj
41 !!$ real(kind=rk) :: dnrm2,r
42 !!$ integer :: pfcount,i,k
43 !!$
44 !!$ !let us save pf%count for later
45 !!$ pfcount = pf%count
46 !!$
47 !!$
48 !!$ write(6,*) 'TESTING H'
49 !!$ write(6,*) 'First doing H^T then doing H'
50 !!$ write(6,*) 'TESTING H WITH SINGLE RIGHT HAND SIDES'
51 !!$
52 !!$ !FIRST TESTS: HHT
53 !!$ !test one - zeros
54 !!$ a = 0.0d0
55 !!$ pf%count = 1
56 !!$ call HT(obs_dim,pf%count,a,j,1)
57 !!$ call H(obs_dim,pf%count,j,s,1)
58 !!$ !s should be a
59 !!$ r = dnrm2(obs_dim,s,1)
60 !!$ write(6,'(A)',advance='no') 'Test 1: H H^T(0) ... '
61 !!$ if(r .lt. 1.0d-15) then
62 !!$ !passed the test
63 !!$ write(6,'(A)',advance='yes') 'passed'
64 !!$ else
65 !!$ write(6,'(A,es24.17)',advance='yes') 'failed with r = ',r
66 !!$ end if
67 !!$
68 !!$ a = 1.0d0
69 !!$ call HT(obs_dim,pf%count,a,j,1)
70 !!$ call H(obs_dim,pf%count,j,s,1)
71 !!$ !s should be a
72 !!$ r = dnrm2(obs_dim,s-a,1)
73 !!$ write(6,'(A)',advance='no') 'Test 2: H H^T(1) ... '
74 !!$ if(r .lt. 1.0d-15) then
75 !!$ !passed the test
76 !!$ write(6,'(A)',advance='yes') 'passed'
77 !!$ else
78 !!$ write(6,'(A,es24.17)',advance='yes') 'failed with r = ',r
79 !!$ end if
80 !!$
81 !!$ a = -1.0d0
82 !!$ call HT(obs_dim,pf%count,a,j,1)
83 !!$ call H(obs_dim,pf%count,j,s,1)
84 !!$ !s should be a
85 !!$ r = dnrm2(obs_dim,s-a,1)
86 !!$ write(6,'(A)',advance='no') 'Test 3: H H^T(-1) ... '
87 !!$ if(r .lt. 1.0d-15) then
88 !!$ !passed the test
89 !!$ write(6,'(A)',advance='yes') 'passed'
90 !!$ else
91 !!$ write(6,'(A,es24.17)',advance='yes') 'failed with r = ',r
92 !!$ end if
93 !!$
94 !!$
95 !!$ call NormalRandomNumbers1D(0.0_rk,1.0_rk,obs_dim,a)
96 !!$ call HT(obs_dim,pf%count,a,j,1)
97 !!$ call H(obs_dim,pf%count,j,s,1)
98 !!$ !s should be a
99 !!$ r = dnrm2(obs_dim,s-a,1)
100 !!$ write(6,'(A)',advance='no') 'Test 4: H H^T( N(0,1) ) ... '
101 !!$ if(r .lt. 1.0d-15) then
102 !!$ !passed the test
103 !!$ write(6,'(A)',advance='yes') 'passed'
104 !!$ else
105 !!$ write(6,'(A,es24.17)',advance='yes') 'failed with r = ',r
106 !!$ end if
107 !!$
108 !!$ write(6,'(A)',advance='yes') 'TESTING H WITH MULTIPLE RIGHT HAND SIDES'
109 !!$ do i = 1,10
110 !!$ aa = 0.0_rk
111 !!$ pf%count = i
112 !!$ call HT(obs_dim,pf%count,aa(:,1:i),jj(:,1:i),i)
113 !!$ call H(obs_dim,pf%count,jj(:,1:i),ss(:,1:i),i)
114 !!$ write(6,'(A,i2,A,i2,A)',advance='no') 'Test 5 ',i,' RHS: H H^T(0) ... '
115 !!$ do k = 1,i
116 !!$ !s should be a
117 !!$ r = dnrm2(obs_dim,ss(:,k)-aa(:,k),1)
118 !!$! write(6,'(A,i2,A,i2,A)',advance='no') 'Test 5 ',i,',',k,': H H^T(0) ... '
119 !!$ if(r .lt. 1.0d-15) then
120 !!$ !passed the test
121 !!$ if(k .eq. 1 .and. i .eq. 1) then
122 !!$ write(6,'(A,i2)',advance='yes') 'passed: ',k
123 !!$ elseif(k .eq. 1) then
124 !!$ write(6,'(A,i2)',advance='no') 'passed: ',k
125 !!$ elseif(k .eq. i) then
126 !!$ write(6,'(A,i2)',advance='yes') ' ',k
127 !!$ else
128 !!$ write(6,'(A,i2)',advance='no') ' ',k
129 !!$ end if
130 !!$ else
131 !!$ !write(6,'(A)',advance='yes') ''
132 !!$ write(6,'(A,es9.2)',advance='yes') 'failed with r = ',r
133 !!$ end if
134 !!$ end do
135 !!$ end do
136 !!$
137 !!$
138 !!$ do i = 1,10
139 !!$ aa = 1.0_rk
140 !!$ pf%count = i
141 !!$ call HT(obs_dim,pf%count,aa(:,1:i),jj(:,1:i),i)
142 !!$ call H(obs_dim,pf%count,jj(:,1:i),ss(:,1:i),i)
143 !!$ write(6,'(A,i2,A,i2,A)',advance='no') 'Test 6 ',i,' RHS: H H^T(1) ... '
144 !!$ do k = 1,i
145 !!$ !s should be a
146 !!$ r = dnrm2(obs_dim,ss(:,k)-aa(:,k),1)
147 !!$! write(6,'(A,i2,A,i2,A)',advance='no') 'Test 5 ',i,',',k,': H H^T(0) ... '
148 !!$ if(r .lt. 1.0d-15) then
149 !!$ !passed the test
150 !!$ if(k .eq. 1 .and. i .eq. 1) then
151 !!$ write(6,'(A,i2)',advance='yes') 'passed: ',k
152 !!$ elseif(k .eq. 1) then
153 !!$ write(6,'(A,i2)',advance='no') 'passed: ',k
154 !!$ elseif(k .eq. i) then
155 !!$ write(6,'(A,i2)',advance='yes') ' ',k
156 !!$ else
157 !!$ write(6,'(A,i2)',advance='no') ' ',k
158 !!$ end if
159 !!$ else
160 !!$ write(6,'(A)',advance='yes') ''
161 !!$ write(6,'(A,es24.17)',advance='yes') 'failed with r = ',r
162 !!$ end if
163 !!$ end do
164 !!$ end do
165 !!$
166 !!$ do i = 1,10
167 !!$ aa = 0.0_rk
168 !!$ pf%count = i
169 !!$ call HT(obs_dim,pf%count,aa(:,1:i),jj(:,1:i),i)
170 !!$ call H(obs_dim,pf%count,jj(:,1:i),ss(:,1:i),i)
171 !!$ write(6,'(A,i2,A,i2,A)',advance='no') 'Test 7 ',i,' RHS: H H^T(-1) ... '
172 !!$ do k = 1,i
173 !!$ !s should be a
174 !!$ r = dnrm2(obs_dim,ss(:,k)-aa(:,k),1)
175 !!$! write(6,'(A,i2,A,i2,A)',advance='no') 'Test 5 ',i,',',k,': H H^T(0) ... '
176 !!$ if(r .lt. 1.0d-15) then
177 !!$ !passed the test
178 !!$ if(k .eq. 1 .and. i .eq. 1) then
179 !!$ write(6,'(A,i2)',advance='yes') 'passed: ',k
180 !!$ elseif(k .eq. 1) then
181 !!$ write(6,'(A,i2)',advance='no') 'passed: ',k
182 !!$ elseif(k .eq. i) then
183 !!$ write(6,'(A,i2)',advance='yes') ' ',k
184 !!$ else
185 !!$ write(6,'(A,i2)',advance='no') ' ',k
186 !!$ end if
187 !!$ else
188 !!$ write(6,'(A)',advance='yes') ''
189 !!$ write(6,'(A,es24.17)',advance='yes') 'failed with r = ',r
190 !!$ end if
191 !!$ end do
192 !!$ end do
193 !!$
194 !!$
195 !!$ call NormalRandomNumbers2D(0.0_rk,1.0_rk,obs_dim,10,aa)
196 !!$ do i = 1,10
197 !!$ pf%count = i
198 !!$ call HT(obs_dim,pf%count,aa(:,1:i),jj(:,1:i),i)
199 !!$ call H(obs_dim,pf%count,jj(:,1:i),ss(:,1:i),i)
200 !!$ write(6,'(A,i2,A,i2,A)',advance='no') 'Test 8 ',i,' RHS: H H^T( N&
201 !!$ &(0,1) ) ... '
202 !!$ do k = 1,i
203 !!$ !s should be a
204 !!$ r = dnrm2(obs_dim,ss(:,k)-aa(:,k),1)
205 !!$! write(6,'(A,i2,A,i2,A)',advance='no') 'Test 5 ',i,',',k,': H H^T(0) ... '
206 !!$ if(r .lt. 1.0d-15) then
207 !!$ !passed the test
208 !!$ if(k .eq. 1 .and. i .eq. 1) then
209 !!$ write(6,'(A,i2)',advance='yes') 'passed: ',k
210 !!$ elseif(k .eq. 1) then
211 !!$ write(6,'(A,i2)',advance='no') 'passed: ',k
212 !!$ elseif(k .eq. i) then
213 !!$ write(6,'(A,i2)',advance='yes') ' ',k
214 !!$ else
215 !!$ write(6,'(A,i2)',advance='no') ' ',k
216 !!$ end if
217 !!$ else
218 !!$ write(6,'(A)',advance='yes') ''
219 !!$ write(6,'(A,es9.2)',advance='yes') 'failed with r = ',r
220 !!$ end if
221 !!$ end do
222 !!$ end do
223 !!$
224 !!$
225 !!$
226 !!$ !put the count back to what it was before
227 !!$ pf%count = pfcount
228 !!$
229 !!$
230 !!$end subroutine H_tests
231 
232 
233 
234 
235 
236 
237 
238 
239 
240 
241 
242 
243 
244 
245 
246 
247 
248 
249 
250 
251 
252 
253 
254 
255 
256 
257 subroutine r_tests()
258 
260  use sizes
261  use pf_control
262  implicit none
263  integer, parameter :: rk = kind(1.0d0)
264  real(kind=rk), dimension(obs_dim) :: a,s,q,w,e
265  real(kind=rk), dimension(obs_dim,10) :: aa,ss,qq,ww,ee
266  real(kind=rk) :: dnrm2,rr
267  real(kind=rk), parameter :: pass = 1.0d-15
268  real(kind=rk), parameter :: warn = 1.0d-13
269  integer :: pfcount,i,k,l
270 
271  !let us save pf%count for later
272  pfcount = pf%count
273 
274 
275  write(6,*) 'TESTING R'
276  write(6,*) 'TESTING R WITH SINGLE RHS'
277 
278 
279  !FIRST TESTS: HHT
280  !test one - zeros
281  a = 0.0d0
282  pf%count = 1
283  call r(obs_dim,1,a,s,1)
284 
285  !s should be 0
286  rr = dnrm2(obs_dim,s,1)
287  write(6,'(A)',advance='no') 'Test 1: R(0) ... '
288  if(rr .lt. pass) then
289  !passed the test
290  write(6,'(A)',advance='yes') 'passed'
291  elseif(rr .lt. warn) then
292  write(6,'(A,es24.17)',advance='yes') 'passed with warning rr = ',rr
293  else
294  write(6,'(A,es24.17)',advance='yes') 'failed with rr = ',rr
295  end if
296 
297  a = 0.0d0
298  call rhalf(obs_dim,1,a,s,1)
299  !s should be a
300  rr = dnrm2(obs_dim,s,1)
301  write(6,'(A)',advance='no') 'Test 2: Rhalf(0) ... '
302  if(rr .lt. pass) then
303  !passed the test
304  write(6,'(A)',advance='yes') 'passed'
305  elseif(rr .lt. warn) then
306  write(6,'(A,es24.17)',advance='yes') 'passed with warning rr = ',rr
307  else
308  write(6,'(A,es24.17)',advance='yes') 'failed with rr = ',rr
309  end if
310 
311  a = 1.0d0
312  s = 1.0d0
313  call r(obs_dim,1,a,q,1)
314  call rhalf(obs_dim,1,s,w,1)
315 
316  if(all(w .eq. 0.0d0)) then
317  write(6,'(A)') 'SERIOUS ERROR: Rhalf*e = 0 i.e. Rhalf is the zero&
318  & matrix'
319  stop 'MEGA FAIL'
320  end if
321 
322  call rhalf(obs_dim,1,w,e,1)
323  !q should be e
324  rr = dnrm2(obs_dim,q-e,1)
325  write(6,'(A)',advance='no') 'Test 3: [RhalfRhalf(1)]-[R(1)] ... '
326  if(rr .lt. pass) then
327  !passed the test
328  write(6,'(A)',advance='yes') 'passed'
329  elseif(rr .lt. warn) then
330  write(6,'(A,es24.17)',advance='yes') 'passed with warning rr = ',rr
331  else
332  write(6,'(A,es24.17)',advance='yes') 'failed with rr = ',rr
333  end if
334 
335  a = -1.0d0
336  s = -1.0d0
337  call r(obs_dim,1,a,q,1)
338  call rhalf(obs_dim,1,s,w,1)
339  call rhalf(obs_dim,1,w,e,1)
340  !q should be e
341  rr = dnrm2(obs_dim,q-e,1)
342  write(6,'(A)',advance='no') 'Test 4: [RhalfRhalf(-1)]-[R(-1)] ... '
343  if(rr .lt. pass) then
344  !passed the test
345  write(6,'(A)',advance='yes') 'passed'
346  elseif(rr .lt. warn) then
347  write(6,'(A,es24.17)',advance='yes') 'passed with warning rr = ',rr
348  else
349  write(6,'(A,es24.17)',advance='yes') 'failed with rr = ',rr
350  end if
351 
352  call normalrandomnumbers1d(0.0d0,1.0d0,obs_dim,a)
353  s = a
354  call r(obs_dim,1,a,q,1)
355  call rhalf(obs_dim,1,s,w,1)
356  call rhalf(obs_dim,1,w,e,1)
357  !q should be e
358  rr = dnrm2(obs_dim,q-e,1)
359  write(6,'(A)',advance='no') 'Test 5: [RhalfRhalf( N(0,1) )]-[R( N(0,1) )] ... '
360  if(rr .lt. pass) then
361  !passed the test
362  write(6,'(A)',advance='yes') 'passed'
363  elseif(rr .lt. warn) then
364  write(6,'(A,es24.17)',advance='yes') 'passed with warning rr = ',rr
365  else
366  write(6,'(A,es24.17)',advance='yes') 'failed with rr = ',rr
367  end if
368 
369  write(6,'(A)',advance='yes') 'TESTING R WITH MULTIPLE RIGHT HAND SIDES'
370 
371 
372  do l = 6,9
373 
374 
375  do i = 1,10
376  if( l .eq. 6) then
377  aa = 0.0_rk
378  ss = 0.0_rk
379  write(6,'(A,i2,A,i2,A)',advance='no') 'Test 6 ',i,' RHS: R(0)-&
380  &Rhalf(Rhalf(0)) ... '
381  elseif( l .eq. 7) then
382  aa = 1.0_rk
383  ss = 1.0_rk
384  write(6,'(A,i2,A,i2,A)',advance='no') 'Test 7 ',i,' RHS: R(1)-&
385  &Rhalf(Rhalf(1)) ... '
386  elseif( l .eq. 8) then
387  aa = -1.0_rk
388  ss = -1.0_rk
389  write(6,'(A,i2,A,i2,A)',advance='no') 'Test 8 ',i,' RHS: R(-1)-Rha&
390  &lf(Rhalf(-1)) ... '
391 
392  elseif( l .eq. 9) then
393  call normalrandomnumbers2d(0.0d0,1.0d0,obs_dim,i,aa(:,1:i))
394  ss = aa
395  write(6,'(A,i2,A,i2,A)',advance='no') 'Test 9 ',i,' RHS: R( N(&
396  &0,1) )-&
397  &Rhalf(Rhalf( N(0,1) )) ... '
398  end if
399  call r(obs_dim,i,aa(:,1:i),qq(:,1:i),1)
400  call rhalf(obs_dim,i,ss(:,1:i),ww(:,1:i),1)
401  call rhalf(obs_dim,i,ww(:,1:i),ee(:,1:i),1)
402 
403  do k = 1,i
404  !qq should be ee
405  rr = dnrm2(obs_dim,qq(:,k)-ee(:,k),1)
406 ! write(6,'(A,i2,A,i2,A)',advance='no') 'Test 5 ',i,',',k,': H H^T(0) ... '
407  if(rr .lt. pass) then
408  !passed the test
409  if(k .eq. 1 .and. i .eq. 1) then
410  write(6,'(A,i2)',advance='yes') 'passed: ',k
411  elseif(k .eq. 1) then
412  write(6,'(A,i2)',advance='no') 'passed: ',k
413  elseif(k .eq. i) then
414  write(6,'(A,i2)',advance='yes') ' ',k
415  else
416  write(6,'(A,i2)',advance='no') ' ',k
417  end if
418  elseif(rr .lt. warn) then
419 
420  if(k .eq. 1 .and. i .eq. 1) then
421  write(6,'(A,i2,A,es9.2)',advance='yes') 'passed: ',k,' warn ',rr
422  elseif(k .eq. 1) then
423  write(6,'(A,i2,A,es9.2)',advance='no') 'passed: ',k,' warn ',rr
424  elseif(k .eq. i) then
425  write(6,'(A,i2,A,es9.2)',advance='yes') ' ',k,' warn ',rr
426  else
427  write(6,'(A,i2,A,es9.2)',advance='no') ' ',k,' warn ',rr
428  end if
429  else
430  !write(6,'(A)',advance='yes') ''
431  write(6,'(i2,A,es9.2)',advance='yes') k,' failed with rr = ',rr
432  end if
433  end do
434  end do
435 
436  end do
437 
438  write(6,*) 'TESTING SOLVE R WITH SINGLE RHS'
439 
440  a = 0.0d0
441  pf%count = 1
442  call solve_r(obs_dim,pf%count,a,s,1)
443 
444  !s should be 0
445  rr = dnrm2(obs_dim,s,1)
446  write(6,'(A)',advance='no') 'Test 10: R^(-1)(0) ... '
447  if(rr .lt. pass) then
448  !passed the test
449  write(6,'(A)',advance='yes') 'passed'
450  elseif(rr .lt. warn) then
451  write(6,'(A,es24.17)',advance='yes') 'passed with warning rr = ',rr
452  else
453  write(6,'(A,es24.17)',advance='yes') 'failed with rr = ',rr
454  end if
455 
456  do l = 11,13
457  if( l .eq. 11) then
458  a = 1.0d0
459  write(6,'(A)',advance='no') 'Test 11: R^(-1)[R(1)] ... '
460  elseif (l .eq. 12) then
461  a = -1.0d0
462  write(6,'(A)',advance='no') 'Test 12: R^(-1)[R(-1)] ... '
463  elseif (l .eq. 13) then
464  call normalrandomnumbers1d(0.0d0,1.0d0,obs_dim,a)
465  write(6,'(A)',advance='no') 'Test 13: R^(-1)[R( N(0,1) )] ... '
466  end if
467 
468 
469  call r(obs_dim,1,a,q,1)
470  call solve_r(obs_dim,pf%count,q,w,1)
471  !w should be a
472  rr = dnrm2(obs_dim,w-a,1)
473 
474 
475  if(rr .lt. pass) then
476  !passed the test
477  write(6,'(A)',advance='yes') 'passed'
478  elseif(rr .lt. warn) then
479  write(6,'(A,es24.17)',advance='yes') 'passed with warning rr = ',rr
480  else
481  write(6,'(A,es24.17)',advance='yes') 'failed with rr = ',rr
482  end if
483 
484  end do
485 
486  write(6,*) 'TESTING SOLVE R WITH MULTIPLE RHS'
487 
488  do l = 14,16
489 
490  do i = 1,10
491  if( l .eq. 14) then
492  aa = 1.0d0
493  write(6,'(A)',advance='no') 'Test 14: R^(-1)[R(1)] ... '
494  elseif (l .eq. 15) then
495  aa = -1.0d0
496  write(6,'(A)',advance='no') 'Test 15: R^(-1)[R(-1)] ... '
497  elseif (l .eq. 16) then
498  call normalrandomnumbers2d(0.0d0,1.0d0,obs_dim,10,aa)
499  write(6,'(A)',advance='no') 'Test 16: R^(-1)[R( N(0,1) )] ... '
500  end if
501 
502  pf%count = i
503  call r(obs_dim,i,aa(:,1:i),qq(:,1:i),1)
504  call solve_r(obs_dim,pf%count,qq(:,1:i),ww(:,1:i),1)
505  !w should be a
506  do k = 1,i
507  rr = dnrm2(obs_dim,ww(:,k)-aa(:,k),1)
508 
509  if(rr .lt. pass) then
510  !passed the test
511  if(k .eq. 1 .and. i .eq. 1) then
512  write(6,'(A,i2)',advance='yes') 'passed: ',k
513  elseif(k .eq. 1) then
514  write(6,'(A,i2)',advance='no') 'passed: ',k
515  elseif(k .eq. i) then
516  write(6,'(A,i2)',advance='yes') ' ',k
517  else
518  write(6,'(A,i2)',advance='no') ' ',k
519  end if
520  elseif(rr .lt. warn) then
521 
522  if(k .eq. 1 .and. i .eq. 1) then
523  write(6,'(A,i2,A,es9.2)',advance='yes') 'passed: ',k,' warn ',rr
524  elseif(k .eq. 1) then
525  write(6,'(A,i2,A,es9.2)',advance='no') 'passed: ',k,' warn ',rr
526  elseif(k .eq. i) then
527  write(6,'(A,i2,A,es9.2)',advance='yes') ' ',k,' warn ',rr
528  else
529  write(6,'(A,i2,A,es9.2)',advance='no') ' ',k,' warn ',rr
530  end if
531  else
532  !write(6,'(A)',advance='yes') ''
533  write(6,'(i2,A,es9.2)',advance='yes') k,' failed with rr = ',rr
534  end if
535 
536 
537  end do
538 
539  end do
540 
541  end do
542 
543 
544  write(6,*) 'TESTING SOLVE Rhalf WITH SINGLE RHS'
545 
546  a = 0.0d0
547  pf%count = 1
548  call solve_rhalf(obs_dim,pf%count,a,s,1)
549 
550  !s should be 0
551  rr = dnrm2(obs_dim,s,1)
552  write(6,'(A)',advance='no') 'Test 17: R^(-1/2)(0) ... '
553  if(rr .lt. pass) then
554  !passed the test
555  write(6,'(A)',advance='yes') 'passed'
556  elseif(rr .lt. warn) then
557  write(6,'(A,es24.17)',advance='yes') 'passed with warning rr = ',rr
558  else
559  write(6,'(A,es24.17)',advance='yes') 'failed with rr = ',rr
560  end if
561 
562  do l = 18,20
563  if( l .eq. 18) then
564  a = 1.0d0
565  write(6,'(A)',advance='no') 'Test 18: R^{-1/2}[R^{1/2}(1)] ... '
566  elseif (l .eq. 19) then
567  a = -1.0d0
568  write(6,'(A)',advance='no') 'Test 19: R^{-1/2}[R^{1/2}(-1)] ... '
569  elseif (l .eq. 20) then
570  call normalrandomnumbers1d(0.0d0,1.0d0,obs_dim,a)
571  write(6,'(A)',advance='no') 'Test 20: R^{-1/2}[R^{1/2}( N(0,1) )] ... '
572  end if
573 
574 
575  call rhalf(obs_dim,1,a,q,1)
576  call solve_rhalf(obs_dim,pf%count,q,w,1)
577  !w should be a
578  rr = dnrm2(obs_dim,w-a,1)
579 
580 
581  if(rr .lt. pass) then
582  !passed the test
583  write(6,'(A)',advance='yes') 'passed'
584  elseif(rr .lt. warn) then
585  write(6,'(A,es24.17)',advance='yes') 'passed with warning rr = ',rr
586  else
587  write(6,'(A,es24.17)',advance='yes') 'failed with rr = ',rr
588  end if
589 
590  end do
591 
592  write(6,*) 'TESTING SOLVE Rhalf WITH MULTIPLE RHS'
593 
594  do l = 21,23
595 
596  do i = 1,10
597  if( l .eq. 21) then
598  aa = 1.0d0
599  write(6,'(A)',advance='no') 'Test 21: R^{-1/2}[R^{1/2}(1)] ... '
600  elseif (l .eq. 22) then
601  aa = -1.0d0
602  write(6,'(A)',advance='no') 'Test 22: R^{-1/2}[R^{1/2}(-1)] ... '
603  elseif (l .eq. 23) then
604  call normalrandomnumbers2d(0.0d0,1.0d0,obs_dim,10,aa)
605  write(6,'(A)',advance='no') 'Test 23: R^{-1/2}[R^{1/2}( N(0,1) )] ... '
606  end if
607 
608  pf%count = i
609  call rhalf(obs_dim,i,aa(:,1:i),qq(:,1:i),1)
610  call solve_rhalf(obs_dim,pf%count,qq(:,1:i),ww(:,1:i),1)
611  !w should be a
612  do k = 1,i
613  rr = dnrm2(obs_dim,ww(:,k)-aa(:,k),1)
614 
615  if(rr .lt. pass) then
616  !passed the test
617  if(k .eq. 1 .and. i .eq. 1) then
618  write(6,'(A,i2)',advance='yes') 'passed: ',k
619  elseif(k .eq. 1) then
620  write(6,'(A,i2)',advance='no') 'passed: ',k
621  elseif(k .eq. i) then
622  write(6,'(A,i2)',advance='yes') ' ',k
623  else
624  write(6,'(A,i2)',advance='no') ' ',k
625  end if
626  elseif(rr .lt. warn) then
627 
628  if(k .eq. 1 .and. i .eq. 1) then
629  write(6,'(A,i2,A,es9.2)',advance='yes') 'passed: ',k,' warn ',rr
630  elseif(k .eq. 1) then
631  write(6,'(A,i2,A,es9.2)',advance='no') 'passed: ',k,' warn ',rr
632  elseif(k .eq. i) then
633  write(6,'(A,i2,A,es9.2)',advance='yes') ' ',k,' warn ',rr
634  else
635  write(6,'(A,i2,A,es9.2)',advance='no') ' ',k,' warn ',rr
636  end if
637  else
638  !write(6,'(A)',advance='yes') ''
639  write(6,'(i2,A,es9.2)',advance='yes') k,' failed with rr = ',rr
640  end if
641 
642 
643  end do
644 
645  end do
646 
647  end do
648 
649 
650 
651 
652 
653  pf%count=pfcount
654 
655 end subroutine r_tests
656 
657 
658 
659 
660 
661 
662 
663 
664 
665 
666 
667 
668 
669 
670 
671 
672 
673 
674 
675 subroutine q_tests()
676 
678  use sizes
679  use pf_control
680  implicit none
681  integer, parameter :: rk = kind(1.0d0)
682  real(kind=rk), dimension(state_dim) :: a,s,qqq,w,e
683  real(kind=rk), dimension(state_dim,10) :: aa,ss,qq,ww,ee
684  real(kind=rk) :: dnrm2,rr
685  real(kind=rk), parameter :: pass = 1.0d-15
686  real(kind=rk), parameter :: warn = 1.0d-13
687  integer :: pfcount,i,k,l
688 
689  !let us save pf%count for later
690  pfcount = pf%count
691 
692 
693  write(6,*) 'TESTING Q'
694  write(6,*) 'TESTING Q WITH SINGLE RHS'
695 
696 
697  !FIRST TESTS: HHT
698  !test one - zeros
699  a = 0.0d0
700  pf%count = 1
701  call q(1,a,s)
702 
703  !s should be 0
704  rr = dnrm2(state_dim,s,1)
705  write(6,'(A)',advance='no') 'Test 1: Q(0) ... '
706  if(rr .lt. pass) then
707  !passed the test
708  write(6,'(A)',advance='yes') 'passed'
709  elseif(rr .lt. warn) then
710  write(6,'(A,es24.17)',advance='yes') 'passed with warning rr = ',rr
711  else
712  write(6,'(A,es24.17)',advance='yes') 'failed with rr = ',rr
713  end if
714 
715  a = 0.0d0
716  call qhalf(1,a,s)
717  !s should be a
718  rr = dnrm2(state_dim,s,1)
719  write(6,'(A)',advance='no') 'Test 2: Qhalf(0) ... '
720  if(rr .lt. pass) then
721  !passed the test
722  write(6,'(A)',advance='yes') 'passed'
723  elseif(rr .lt. warn) then
724  write(6,'(A,es24.17)',advance='yes') 'passed with warning rr = ',rr
725  else
726  write(6,'(A,es24.17)',advance='yes') 'failed with rr = ',rr
727  end if
728 
729  a = 1.0d0
730  s = 1.0d0
731  call q(1,a,qqq)
732  call qhalf(1,s,w)
733 
734  if(all(w .eq. 0.0d0)) then
735  write(6,'(A)') 'SERIOUS ERROR: Qhalf*e = 0 i.e. Qhalf is the zero&
736  & matrix'
737  stop 'MEGA FAIL'
738  end if
739 
740  call qhalf(1,w,e)
741  !q should be e
742  rr = dnrm2(state_dim,qqq-e,1)
743  write(6,'(A)',advance='no') 'Test 3: [QhalfQhalf(1)]-[Q(1)] ... '
744  if(rr .lt. pass) then
745  !passed the test
746  write(6,'(A)',advance='yes') 'passed'
747  elseif(rr .lt. warn) then
748  write(6,'(A,es24.17)',advance='yes') 'passed with warning rr = ',rr
749  else
750  write(6,'(A,es24.17)',advance='yes') 'failed with rr = ',rr
751  end if
752 
753  a = -1.0d0
754  s = -1.0d0
755  call q(1,a,qqq)
756  call qhalf(1,s,w)
757  call qhalf(1,w,e)
758  !q should be e
759  rr = dnrm2(state_dim,qqq-e,1)
760  write(6,'(A)',advance='no') 'Test 4: [QhalfQhalf(-1)]-[R(-1)] ... '
761  if(rr .lt. pass) then
762  !passed the test
763  write(6,'(A)',advance='yes') 'passed'
764  elseif(rr .lt. warn) then
765  write(6,'(A,es24.17)',advance='yes') 'passed with warning rr = ',rr
766  else
767  write(6,'(A,es24.17)',advance='yes') 'failed with rr = ',rr
768  end if
769 
770  call normalrandomnumbers1d(0.0d0,1.0d0,state_dim,a)
771  s = a
772  call q(1,a,qqq)
773  call qhalf(1,s,w)
774  call qhalf(1,w,e)
775  !q should be e
776  rr = dnrm2(state_dim,qqq-e,1)
777  write(6,'(A)',advance='no') 'Test 5: [QhalfQhalf( N(0,1) )]-[Q( N(0,1) )] ... '
778  if(rr .lt. pass) then
779  !passed the test
780  write(6,'(A)',advance='yes') 'passed'
781  elseif(rr .lt. warn) then
782  write(6,'(A,es24.17)',advance='yes') 'passed with warning rr = ',rr
783  else
784  write(6,'(A,es24.17)',advance='yes') 'failed with rr = ',rr
785  end if
786 
787  write(6,'(A)',advance='yes') 'TESTING Q WITH MULTIPLE RIGHT HAND SIDES'
788 
789 
790  do l = 6,8
791 
792 
793  do i = 1,10
794  if( l .eq. 6) then
795  aa = 1.0_rk
796  ss = 1.0_rk
797  write(6,'(A,i2,A,i2,A)',advance='no') 'Test 6 ',i,' RHS: Q(1)-&
798  &Qhalf(Qhalf(1)) ... '
799  elseif( l .eq. 7) then
800  aa = -1.0_rk
801  ss = -1.0_rk
802  write(6,'(A,i2,A,i2,A)',advance='no') 'Test 7 ',i,' RHS: Q(-1)-Qha&
803  &lf(Qhalf(-1)) ... '
804 
805  elseif( l .eq. 8) then
806  call normalrandomnumbers2d(0.0d0,1.0d0,state_dim,i,aa(:,1:i))
807  ss = aa
808  write(6,'(A,i2,A,i2,A)',advance='no') 'Test 8 ',i,' RHS: Q( N(&
809  &0,1) )-Qhalf(Qhalf( N(0,1) )) ... '
810  end if
811  call q(i,aa(:,1:i),qq(:,1:i))
812  call qhalf(i,ss(:,1:i),ww(:,1:i))
813  call qhalf(i,ww(:,1:i),ee(:,1:i))
814 
815  do k = 1,i
816  !qq should be ee
817  rr = dnrm2(state_dim,qq(:,k)-ee(:,k),1)
818 ! write(6,'(A,i2,A,i2,A)',advance='no') 'Test 5 ',i,',',k,': H H^T(0) ... '
819  if(rr .lt. pass) then
820  !passed the test
821  if(k .eq. 1 .and. i .eq. 1) then
822  write(6,'(A,i2)',advance='yes') 'passed: ',k
823  elseif(k .eq. 1) then
824  write(6,'(A,i2)',advance='no') 'passed: ',k
825  elseif(k .eq. i) then
826  write(6,'(A,i2)',advance='yes') ' ',k
827  else
828  write(6,'(A,i2)',advance='no') ' ',k
829  end if
830  elseif(rr .lt. warn) then
831 
832  if(k .eq. 1 .and. i .eq. 1) then
833  write(6,'(A,i2,A,es9.2)',advance='yes') 'passed: ',k,' warn ',rr
834  elseif(k .eq. 1) then
835  write(6,'(A,i2,A,es9.2)',advance='no') 'passed: ',k,' warn ',rr
836  elseif(k .eq. i) then
837  write(6,'(A,i2,A,es9.2)',advance='yes') ' ',k,' warn ',rr
838  else
839  write(6,'(A,i2,A,es9.2)',advance='no') ' ',k,' warn ',rr
840  end if
841  else
842  !write(6,'(A)',advance='yes') ''
843  write(6,'(i2,A,es9.2)',advance='yes') k,' failed with rr = ',rr
844  end if
845  end do
846  end do
847 
848  end do
849 
850  pf%count = pfcount
851 
852 end subroutine q_tests
853 
854 
855 
856 
857 
858 
859 
860 
861 
862 
863 
864 
865 
866 
867 
868 
869 
870 
871 
872 
873 
874 
875 
876 
877 
878 
879 
880 
881 subroutine hqhtr_tests()
882 
887  use pf_control
888  use sizes
889  implicit none
890  integer, parameter :: rk = kind(1.0d0)
891  real(kind=rk), dimension(obs_dim) :: a,s,qq,hqha,hqha_r
892  real(kind=rk), dimension(state_dim) :: qha,ha
893  real(kind=rk) :: dnrm2,rr
894  real(kind=rk), parameter :: pass = 1.0d-15
895  real(kind=rk), parameter :: warn = 1.0d-13
896 
897  integer :: pfcount,l
898 
899  !let us save pf%count for later
900  pfcount = pf%count
901 
902 
903 write(6,*) 'TESTING (HQH^T+R)^(-1)'
904 
905  !FIRST TESTS: HHT
906  !test one - zeros
907  a = 0.0d0
908  pf%count = 1
909  call solve_hqht_plus_r(obs_dim,a,s,1)
910  !s should be a
911  rr = dnrm2(obs_dim,s,1)
912  write(6,'(A)',advance='no') 'Test 1: (HQH^T+R)^(-1)[0] ... '
913  if(rr .lt. pass) then
914  !passed the test
915  write(6,'(A)',advance='yes') 'passed'
916  elseif( rr .lt. warn) then
917  write(6,'(A,es24.17)',advance='yes') 'passed with warning rr = ',rr
918  else
919  write(6,'(A,es24.17)',advance='yes') 'failed with rr = ',rr
920  end if
921 
922 
923  do l = 2,10
924  if( l .eq. 2) then
925  a = 1.0d0
926  write(6,'(A)',advance='no') 'Test 2: (HQH^T+R)^(-1)[(HQH^T+R)(1)] ... '
927  elseif (l .eq. 3) then
928  a = -1.0d0
929  write(6,'(A)',advance='no') 'Test 3: (HQH^T+R)^(-1)[(HQH^T+R)(1)] ... '
930  elseif (l .ge. 4) then
931  call normalrandomnumbers1d(0.0d0,1.0d0,obs_dim,a)
932  write(6,'(A,i2,A)',advance='no') 'Test ',l,': (HQH^T+R)^(-1)[(&
933  &HQH^T+R)( N(0,1) )] ... '
934  end if
935 
936 
937  call r(obs_dim,1,a,qq,1)
938  call ht(obs_dim,1,a,ha,1)
939  call q(1,ha,qha)
940  call h(obs_dim,1,qha,hqha,1)
941  hqha_r = hqha+qq
942  call solve_hqht_plus_r(obs_dim,hqha_r,s,1)
943  !w should be a
944  rr = dnrm2(obs_dim,s-a,1)
945 
946 
947  if(rr .lt. pass) then
948  !passed the test
949  write(6,'(A)',advance='yes') 'passed'
950  elseif(rr .lt. warn) then
951  write(6,'(A,es9.2)',advance='yes') 'warning rr = ',rr
952  else
953  write(6,'(A,es9.2)',advance='yes') 'failed rr = ',rr
954  end if
955 
956  end do
957 
958 
959 
960 
961 
962 end subroutine hqhtr_tests
963 
964 
965 subroutine b_tests()
966 
968  use sizes
969  use pf_control
970  implicit none
971  integer, parameter :: rk = kind(1.0d0)
972  real(kind=rk), dimension(state_dim) :: a,s,q,qt,w
973  real(kind=rk), dimension(state_dim,10) :: aa,qq,qqt,ww
974  real(kind=rk) :: dnrm2,rr
975  real(kind=rk), parameter :: pass = 1.0d-15
976  real(kind=rk), parameter :: warn = 1.0d-13
977  integer :: pfcount,i,k,l
978 
979  !let us save pf%count for later
980  pfcount = pf%count
981 
982 
983  write(6,*) 'TESTING B'
984  write(6,*) 'TESTING B WITH SINGLE RHS'
985 
986 
987  !FIRST TESTS: HHT
988  !test one - zeros
989  a = 0.0d0
990  pf%count = 1
991  call bhalf(1,a,s)
992 
993  !s should be 0
994  rr = dnrm2(state_dim,s,1)
995  write(6,'(A)',advance='no') 'Test 1: Bhalf(0) ... '
996  if(rr .lt. pass) then
997  !passed the test
998  write(6,'(A)',advance='yes') 'passed'
999  elseif(rr .lt. warn) then
1000  write(6,'(A,es24.17)',advance='yes') 'passed with warning rr = ',rr
1001  else
1002  write(6,'(A,es24.17)',advance='yes') 'failed with rr = ',rr
1003  end if
1004 
1005 
1006  a = 1.0d0
1007  s = 1.0d0
1008  call bhalf(1,s,w)
1009 
1010  if(all(w .eq. 0.0d0)) then
1011  write(6,'(A)') 'SERIOUS ERROR: Bhalf*e = 0 i.e. Bhalf is the zero&
1012  & matrix'
1013  stop 'MEGA FAIL'
1014  end if
1015 
1016 
1017 
1018  write(6,*) 'TESTING SOLVE B WITH SINGLE RHS'
1019 
1020  a = 0.0d0
1021  pf%count = 1
1022  call solve_b(pf%count,a,s)
1023 
1024  !s should be 0
1025  rr = dnrm2(state_dim,s,1)
1026  write(6,'(A)',advance='no') 'Test 10: B^(-1)(0) ... '
1027  if(rr .lt. pass) then
1028  !passed the test
1029  write(6,'(A)',advance='yes') 'passed'
1030  elseif(rr .lt. warn) then
1031  write(6,'(A,es24.17)',advance='yes') 'passed with warning rr = ',rr
1032  else
1033  write(6,'(A,es24.17)',advance='yes') 'failed with rr = ',rr
1034  end if
1035 
1036  do l = 11,13
1037  if( l .eq. 11) then
1038  a = 1.0d0
1039  write(6,'(A)',advance='no') 'Test 11: B^(-1)[BhalfBhalf(1)] ... '
1040  elseif (l .eq. 12) then
1041  a = -1.0d0
1042  write(6,'(A)',advance='no') 'Test 12: B^(-1)[BhalfBhalf(-1)] ... '
1043  elseif (l .eq. 13) then
1044  call normalrandomnumbers1d(0.0d0,1.0d0,obs_dim,a)
1045  write(6,'(A)',advance='no') 'Test 13: B^(-1)[BhalfBhalf( N(0,1) )] ... '
1046  end if
1047 
1048 
1049  call bhalf(1,a,qt)
1050  call bhalf(1,qt,q)
1051  call solve_b(pf%count,q,w)
1052  !w should be a
1053  rr = dnrm2(state_dim,w-a,1)
1054 
1055 
1056  if(rr .lt. pass) then
1057  !passed the test
1058  write(6,'(A)',advance='yes') 'passed'
1059  elseif(rr .lt. warn) then
1060  write(6,'(A,es24.17)',advance='yes') 'passed with warning rr = ',rr
1061  else
1062  write(6,'(A,es24.17)',advance='yes') 'failed with rr = ',rr
1063  end if
1064 
1065  end do
1066 
1067  write(6,*) 'TESTING SOLVE B WITH MULTIPLE RHS'
1068 
1069  do l = 14,16
1070 
1071  do i = 1,10
1072  if( l .eq. 14) then
1073  aa = 1.0d0
1074  write(6,'(A)',advance='no') 'Test 14: B^(-1)[BhalfBhalf(1)] ... '
1075  elseif (l .eq. 15) then
1076  aa = -1.0d0
1077  write(6,'(A)',advance='no') 'Test 15: B^(-1)[BhalfBhalf(-1)] ... '
1078  elseif (l .eq. 16) then
1079  call normalrandomnumbers2d(0.0d0,1.0d0,state_dim,10,aa)
1080  write(6,'(A)',advance='no') 'Test 16: B^(-1)[BhalfBhalf( N(0,1) )] ... '
1081  end if
1082 
1083  pf%count = i
1084  call bhalf(i,aa(:,1:i),qqt(:,1:i))
1085  call bhalf(i,qqt(:,1:i),qq(:,1:i))
1086  call solve_b(pf%count,qq(:,1:i),ww(:,1:i))
1087  !w should be a
1088  do k = 1,i
1089  rr = dnrm2(state_dim,ww(:,k)-aa(:,k),1)
1090 
1091  if(rr .lt. pass) then
1092  !passed the test
1093  if(k .eq. 1 .and. i .eq. 1) then
1094  write(6,'(A,i2)',advance='yes') 'passed: ',k
1095  elseif(k .eq. 1) then
1096  write(6,'(A,i2)',advance='no') 'passed: ',k
1097  elseif(k .eq. i) then
1098  write(6,'(A,i2)',advance='yes') ' ',k
1099  else
1100  write(6,'(A,i2)',advance='no') ' ',k
1101  end if
1102  elseif(rr .lt. warn) then
1103 
1104  if(k .eq. 1 .and. i .eq. 1) then
1105  write(6,'(A,i2,A,es9.2)',advance='yes') 'passed: ',k,' warn ',rr
1106  elseif(k .eq. 1) then
1107  write(6,'(A,i2,A,es9.2)',advance='no') 'passed: ',k,' warn ',rr
1108  elseif(k .eq. i) then
1109  write(6,'(A,i2,A,es9.2)',advance='yes') ' ',k,' warn ',rr
1110  else
1111  write(6,'(A,i2,A,es9.2)',advance='no') ' ',k,' warn ',rr
1112  end if
1113  else
1114  !write(6,'(A)',advance='yes') ''
1115  write(6,'(i2,A,es9.2)',advance='yes') k,' failed with rr = ',rr
1116  end if
1117 
1118 
1119  end do
1120 
1121  end do
1122 
1123 end do
1124 
1125 
1126  pf%count = pfcount
1127 
1128 end subroutine b_tests
normalrandomnumbers2d
subroutine normalrandomnumbers2d(mean, stdev, n, k, phi)
generate two dimensional Normal random numbers
Definition: gen_rand.f90:109
b_tests
subroutine b_tests()
Definition: tests.f90:965
bhalf
subroutine bhalf(nrhs, x, bx)
subroutine to take a full state vector x and return in state space.
Definition: model_specific.f90:281
solve_r
subroutine solve_r(obsDim, nrhs, y, v, t)
subroutine to take an observation vector y and return v in observation space.
Definition: model_specific.f90:68
r_tests
subroutine r_tests()
Definition: tests.f90:257
sizes
Module that stores the dimension of observation and state spaces.
Definition: sizes.f90:29
normalrandomnumbers1d
subroutine normalrandomnumbers1d(mean, stdev, n, phi)
generate one dimension of Normal random numbers
Definition: gen_rand.f90:74
h
subroutine h(obsDim, nrhs, x, hx, t)
subroutine to take a full state vector x and return H(x) in observation space.
Definition: model_specific.f90:221
solve_b
subroutine solve_b(nrhs, x, v)
subroutine to take a state vector x and return v in state space.
Definition: model_specific.f90:301
q_tests
subroutine q_tests()
Definition: tests.f90:675
solve_rhalf
subroutine solve_rhalf(obsdim, nrhs, y, v, t)
subroutine to take an observation vector y and return v in observation space.
Definition: model_specific.f90:89
hqhtr_tests
subroutine hqhtr_tests()
Definition: tests.f90:881
ht
subroutine ht(obsDim, nrhs, y, x, t)
subroutine to take an observation vector y and return x in full state space.
Definition: model_specific.f90:243
rhalf
subroutine rhalf(obsDim, nrhs, y, Ry, t)
subroutine to take an observation vector x and return Rx in observation space.
Definition: model_specific.f90:198
pf_control
module pf_control holds all the information to control the the main program
Definition: pf_control.f90:29
solve_hqht_plus_r
subroutine solve_hqht_plus_r(obsdim, y, v, t)
subroutine to take an observation vector y and return v in observation space.
Definition: model_specific.f90:111
qhalf
subroutine qhalf(nrhs, x, Qx)
subroutine to take a full state vector x and return in state space.
Definition: model_specific.f90:156
q
subroutine q(nrhs, x, Qx)
subroutine to take a full state vector x and return Qx in state space.
Definition: model_specific.f90:131
r
subroutine r(obsDim, nrhs, y, Ry, t)
subroutine to take an observation vector x and return Rx in observation space.
Definition: model_specific.f90:176