Ticket #349: ropp_pp_utils.f90_SSY

File ropp_pp_utils.f90_SSY, 13.0 KB (added by Ian Culverwell, 12 years ago)

ropp_pp_utils.f90_SSY

Line 
1! $Id: ropp_pp_utils.f90 3491 2013-02-06 12:43:43Z idculv $
2
3MODULE ropp_pp_utils
4
5!****m* Modules/ropp_pp_utils *
6!
7! NAME
8! ropp_pp_utils - Utility routines for ionospheric correction
9!
10! SYNOPSIS
11! use ropp_pp_utils
12!
13! DESCRIPTION
14! This module provides signal and matrix processing routines used by
15! the ROPP pre-processor ionospheric correction package.
16!
17! AUTHOR
18! M Gorbunov, Russian Academy of Sciences, Russia.
19! Any comments on this software should be given via the ROM SAF
20! Helpdesk at http://www.romsaf.org
21!
22! COPYRIGHT
23! Copyright (c) 1998-2010 Michael Gorbunov <michael.gorbunov@zmaw.de>
24! For further details please refer to the file COPYRIGHT
25! which you should have received as part of this distribution.
26!
27!****
28
29 USE typesizes, ONLY: wp => EightByteReal
30
31 INTERFACE ropp_pp_matmul
32 MODULE PROCEDURE ropp_pp_matmul_arr
33 MODULE PROCEDURE ropp_pp_matmul_vecl
34 MODULE PROCEDURE ropp_pp_matmul_vecr
35 END INTERFACE
36
37CONTAINS
38
39!****s* PPUtils/ropp_pp_matmul *
40!
41! NAME
42! ropp_pp_matmul - Matrix multiplication
43!
44! SYNOPSIS
45! C = ropp_pp_matmul(A, B)
46!
47! DESCRIPTION
48! Replicates intrinsic 'MATMUL' function
49!
50!****
51!
52SUBROUTINE ropp_pp_matmul_sub(A, B, C)
53
54 IMPLICIT NONE
55
56 REAL(wp), DIMENSION(:,:), INTENT(in) :: A
57 REAL(wp), DIMENSION(:,:), INTENT(in) :: B
58 REAL(wp), DIMENSION(:,:), INTENT(out) :: C
59 INTEGER :: i, j, k
60
61 DO i=1,SIZE(A,1)
62 DO j=1, SIZE(B,2)
63 C(i,j) = 0.0_wp
64 DO k=1, SIZE(B,1)
65 C(i,j) = C(i,j) + A(i,k)*B(k,j)
66 ENDDO
67 ENDDO
68 ENDDO
69
70END SUBROUTINE ropp_pp_matmul_sub
71
72FUNCTION ropp_pp_matmul_arr(A, B) RESULT (C)
73
74 IMPLICIT NONE
75
76 REAL(wp), DIMENSION(:,:), INTENT(in) :: A
77 REAL(wp), DIMENSION(:,:), INTENT(in) :: B
78 REAL(wp), DIMENSION(SIZE(A,1),SIZE(B,2)) :: C
79 INTEGER :: i, j, k
80
81 DO i=1,SIZE(A,1)
82 DO j=1, SIZE(B,2)
83 C(i,j) = 0.0_wp
84 DO k=1, SIZE(B,1)
85 C(i,j) = C(i,j) + A(i,k)*B(k,j)
86 ENDDO
87 ENDDO
88 ENDDO
89
90END FUNCTION ropp_pp_matmul_arr
91
92FUNCTION ropp_pp_matmul_vecr(A, B) RESULT (C)
93
94 IMPLICIT NONE
95
96 REAL(wp), DIMENSION(:,:), INTENT(in) :: A
97 REAL(wp), DIMENSION(:), INTENT(in) :: B
98 REAL(wp), DIMENSION(SIZE(A,1)) :: C
99 INTEGER :: j, k
100
101 DO j=1, SIZE(A,1)
102 C(j) = 0.0_wp
103 DO k=1, SIZE(B)
104 C(j) = C(j) + A(j,k)*B(k)
105 ENDDO
106 ENDDO
107
108END FUNCTION ropp_pp_matmul_vecr
109
110FUNCTION ropp_pp_matmul_vecl(A, B) RESULT (C)
111
112 IMPLICIT NONE
113
114 REAL(wp), DIMENSION(:), INTENT(in) :: A
115 REAL(wp), DIMENSION(:,:), INTENT(in) :: B
116 REAL(wp), DIMENSION(SIZE(B,2)) :: C
117 INTEGER :: j, k
118
119 DO j=1, SIZE(B,2)
120 C(j) = 0.0_wp
121 DO k=1, SIZE(B,1)
122 C(j) = C(j) + A(k)*B(k,j)
123 ENDDO
124 ENDDO
125
126END FUNCTION ropp_pp_matmul_vecl
127
128
129!****s* PPUtils/ropp_pp_invert_matrix *
130!
131! NAME
132! ropp_pp_invert_matrix - Invert a matrix A(dimY, dimX)
133!
134! SYNOPSIS
135! call ropp_pp_invert_matrix(A, B)
136!
137! DESCRIPTION
138! Gauss elimination
139!
140!****
141!
142SUBROUTINE ropp_pp_invert_matrix(A, B)
143
144 IMPLICIT NONE
145
146 ! 2.1 Declarations
147
148 REAL(wp), DIMENSION(:,:), INTENT(inout) :: A ! Matrix to invert
149 REAL(wp), DIMENSION(:,:), INTENT(out) :: B ! Inverted matrix
150
151 INTEGER :: N ! Matrix dimension
152 INTEGER :: i, k ! Matrix indeces
153 INTEGER :: m(1) ! Search index
154 REAL(wp) :: Alpha ! Row combination coefficien
155 REAL(wp) :: Det ! Matrix determinant
156 REAL(wp) :: Amin, Amax ! Max and min diagonal elements
157 REAL(wp), DIMENSION(:), ALLOCATABLE :: F ! Exchange work space
158
159 ! 2.2 Check array shape
160
161 IF (SIZE(A,1) /= SIZE(A,2)) THEN
162 RETURN
163 ENDIF
164 IF (ANY(SHAPE(B) /= SHAPE(B))) THEN
165 RETURN
166 ENDIF
167
168 ALLOCATE(F(SIZE(A,1)))
169
170 ! 2.3 Initialization
171
172 N = SIZE(A,1)
173 B(:,:) = 0.0_wp
174
175 DO i=1,N
176 B(i,i) = 1.0_wp
177 ENDDO
178
179 ! 2.4 Elimination of sub-diagonal elements
180
181 Lower: DO k=1,N-1
182
183 IF (A(k,k) == 0.0_wp) THEN
184 m(:) = MAXLOC(A(k+1:N,k), A(k+1:N,k) /= 0.0_wp)
185 IF (m(1) == 0) THEN
186 EXIT Lower
187 END IF
188 i = k + m(1)
189 F(:) = A(k,:)
190 A(k,:) = A(i,:)
191 A(i,:) = F(:)
192 F(:) = B(k,:)
193 B(k,:) = B(i,:)
194 B(i,:) = F(:)
195 END IF
196
197 DO i=k+1,N
198 Alpha = A(i,k)/A(k,k)
199 A(i,:) = A(i,:) - Alpha*A(k,:)
200 B(i,:) = B(i,:) - Alpha*B(k,:)
201 END DO
202
203 END DO Lower
204
205 ! 2.5 Checking for degenerated matrix
206
207 Det = 1.0_wp
208 Amin = ABS (A(1,1))
209 Amax = ABS (A(1,1))
210 Diagonal: DO i=1,N
211 Amin = MIN(ABS (A(i,i)), Amin)
212 Amax = MAX(ABS (A(i,i)), Amax)
213 Det = Det*A(i,i)
214 END DO Diagonal
215
216 IF (Det == 0.0_wp) THEN
217 RETURN
218 END IF
219
220 ! 2.6 Elimination of super-diagonal elements
221
222 Upper: DO k=N,2,-1
223 DO i=k-1,1,-1
224 Alpha = A(i,k)/A(k,k)
225 A(i,:) = A(i,:) - Alpha*A(k,:)
226 B(i,:) = B(i,:) - Alpha*B(k,:)
227 END DO
228 END DO Upper
229
230 DO i=1,N
231 B(i,:) = B(i,:)/A(i,i)
232 END DO
233
234 DEALLOCATE(F)
235
236END SUBROUTINE ropp_pp_invert_matrix
237
238!****s* PPUtils/ropp_pp_quasi_invert *
239!
240! NAME
241! ropp_pp_quasi_invert - Quasi-inverse of a matrix K(dimY, dimX)
242!
243! SYNOPSIS
244! call ropp_pp_quasi_invert(K, Q)
245!
246! DESCRIPTION
247! 1. dimY >= dimX:
248! x = Qy - vector minimizing ||Kx - y||
249! QK = E; Q is left inverse operator.
250! Q = (K^T K)^-1 K^T
251! 2. dimY <= dimX:
252! x = Qy - solution minimizing ||x||
253! KQ = E; Q is right inverse operator.
254! Q = K^T (KK^T)^-1
255!
256!****
257!
258 SUBROUTINE ropp_pp_quasi_invert(K, Q)
259
260 IMPLICIT NONE
261
262 ! 3.1 Declarations
263
264 REAL(wp), DIMENSION(:,:), INTENT(in) :: K ! Matrix to quasi-invert
265 REAL(wp), DIMENSION(:,:), INTENT(out) :: Q ! Quasi-inverse
266
267 INTEGER :: DimX, DimY ! Dimension of x and y spaces
268 INTEGER :: N ! Work matrix dimension
269 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: W ! Work arrays for inversion
270 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: WI ! Work arrays for inversion
271 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: KT
272
273 ! 3.2 Initialization
274
275 DimX = SIZE(K,2)
276 DimY = SIZE(K,1)
277
278 IF (ANY(SHAPE(Q) /= (/ DimX, DimY /))) THEN
279 RETURN
280 END IF
281
282 N = MIN(DimX, DimY)
283 ALLOCATE (W(N,N))
284 ALLOCATE (WI(N,N))
285 ALLOCATE (KT(DimX,DimY))
286
287 KT = TRANSPOSE(K)
288
289 ! 3.3 Quasi-Inversion
290
291 IF (DimY >= DimX) THEN
292
293 ! 3.3.1 Left inversion W=K^TK, WI=(K^TK)^-1, Q=(K^TK)^-1 K^T
294
295 W = ropp_pp_matmul(KT, K)
296 CALL ropp_pp_invert_matrix(W, WI)
297! Q = ropp_pp_matmul(WI, KT)
298 CALL ropp_pp_matmul_sub(WI, KT, Q)
299
300 ELSE
301
302 ! 3.3.2. Right inversion W=KK^T, WI=(KK^T)^-1, Q=K^T(KK^T)^-1
303
304 W = ropp_pp_matmul(K, KT)
305 CALL ropp_pp_invert_matrix(W, WI)
306! Q = ropp_pp_matmul(KT, WI)
307 CALL ropp_pp_matmul_sub(KT, WI, Q)
308
309 END IF
310
311 ! 3.4 Clean up
312
313 DEALLOCATE (W)
314 DEALLOCATE (WI)
315 DEALLOCATE (KT)
316
317 END SUBROUTINE ropp_pp_quasi_invert
318
319!****s* PPUtils/ropp_pp_polynomial *
320!
321! NAME
322! ropp_pp_polynomial - Calculation of a polynomial and its derivative
323!
324! SYNOPSIS
325! call ropp_pp_polynomial(c, x, P, DP)
326!
327! DESCRIPTION
328! Horner scheme
329!
330!****
331!
332 SUBROUTINE ropp_pp_polynomial(c, x, P, DP)
333
334 USE typesizes, ONLY: wp => EightByteReal
335
336 ! 4.1 Declarations
337
338 IMPLICIT NONE
339
340 REAL(wp), DIMENSION(0:), INTENT(in) :: c ! polynomial coefficients
341 REAL(wp), INTENT(in) :: x ! polynomial argument
342 REAL(wp), INTENT(out) :: P ! polynomial value
343 REAL(wp), OPTIONAL, INTENT(out) :: DP ! polynomial derivative
344
345 INTEGER :: n, i
346
347 ! 4.2 Calculate polynomial value
348
349 n = UBOUND(c,1)
350
351 P = c(n)
352 DO i = n-1, 0, -1
353 P = c(i) + x*P
354 ENDDO
355
356 ! 4.3 Calculate polynomial derivative
357
358 IF ( PRESENT(DP) ) THEN
359 DP = c(n)*REAL(n,wp)
360 DO i = n-1, 1, -1
361 DP = c(i)*REAL(i,wp) + x*DP
362 ENDDO
363 ENDIF
364
365 END SUBROUTINE ropp_pp_polynomial
366
367!****s* PPUtils/ropp_pp_init_polynomial *
368!
369! NAME
370! ropp_pp_init_polynomial - Generate matrix of basic polynomials
371! for regression
372!
373! SYNOPSIS
374! call ropp_pp_init_polynomial(x, K)
375!
376! DESCRIPTION
377! Compute basic polynomials as
378! K_ij = f_j(x_i) = (x_i)^j for j=0..UBound(K,2)
379!
380!****
381!
382 SUBROUTINE ropp_pp_init_polynomial(x, K)
383
384 USE typesizes, ONLY: wp => EightByteReal
385
386 ! 5.1 Declarations
387
388 IMPLICIT NONE
389
390 REAL(wp), DIMENSION(1:), INTENT(in) :: x ! grid of argument x
391 REAL(wp), DIMENSION(1:,0:), INTENT(out) :: K ! matrix of polynomials
392 INTEGER :: i,j
393
394 ! 5.2 Generate polynomials
395
396 K(:,0) = 1.0_wp
397 DO j=1,UBOUND(K,2)
398 DO i=1,SIZE(x)
399 IF (x(i) > 0) THEN
400 K(i, j) = x(i)**REAL(j,wp)
401 ELSE
402 K(i, j) = x(i)**INT(REAL(j,wp))
403 ENDIF
404 ENDDO
405 ENDDO
406
407 END SUBROUTINE ropp_pp_init_polynomial
408
409!****s* PPUtils/ropp_pp_regression *
410!
411! NAME
412! ropp_pp_regression - Linear regression
413!
414! SYNOPSIS
415! call ropp_pp_regression(K, y, a)
416!
417! DESCRIPTION
418! Quasi-inversion of matrix of basic functions:
419! || Sum_j a_j f_j(x_i) - y_i || = min
420! Solution is a = Q y where Q is left inverse of K_ij = f_j(x_i)
421!
422!****
423!
424 SUBROUTINE ropp_pp_regression(K, y, a)
425
426 USE typesizes, ONLY: wp => EightByteReal
427
428 ! 6.1 Declarations
429
430 IMPLICIT NONE
431
432 REAL(wp), DIMENSION(:,:), INTENT(in) :: K ! matrix of functions
433 REAL(wp), DIMENSION(:), INTENT(in) :: y ! regression data
434 REAL(wp), DIMENSION(:), INTENT(out) :: a ! regression coefficients
435
436 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: Q ! quasi-inverse matrix
437
438 ALLOCATE(Q(SIZE(K,2), SIZE(K,1)))
439
440 ! 6.2 Invert matrix of functions
441
442 CALL ropp_pp_quasi_invert(K, Q)
443
444 ! 6.3 Solve to find coefficients
445
446
447 a(:) = ropp_pp_matmul(Q, y(:))
448
449 DEALLOCATE(Q)
450
451 END SUBROUTINE ropp_pp_regression
452
453!****s* PPUtils/ropp_pp_residual_regression *
454!
455! NAME
456! ropp_pp_residual_regression - linear regression on residual
457!
458! SYNOPSIS
459! call ropp_pp_residual_regression(K, x, y, a)
460!
461! DESCRIPTION
462! Updates the regression coefficients using
463! the residual from a previous regression estimate
464!
465!****
466!
467
468 SUBROUTINE ropp_pp_residual_regression(K, x, y, a)
469
470 USE typesizes, ONLY: wp => EightByteReal
471
472 ! 7.1 Declarations
473
474 IMPLICIT NONE
475
476 REAL(wp), DIMENSION(:,:),INTENT(in) :: K ! matrix of functions
477 REAL(wp), DIMENSION(:), INTENT(in) :: x ! independent variable
478 REAL(wp), DIMENSION(:), INTENT(in) :: y ! regression data
479 REAL(wp), DIMENSION(:), INTENT(inout):: a ! regression coefficients
480
481 REAL(wp), DIMENSION(:), ALLOCATABLE :: y_a ! polynomial data
482 REAL(wp), DIMENSION(:), ALLOCATABLE :: y_res ! residual data
483 REAL(wp), DIMENSION(:), ALLOCATABLE :: a_res ! residual coefficients
484 INTEGER :: i
485
486 ALLOCATE(y_a(SIZE(y)))
487 ALLOCATE(y_res(SIZE(y)))
488 ALLOCATE(a_res(SIZE(a)))
489
490 ! 7.2 Calculate polynomial data from input regression coefficients
491
492 DO i=1,SIZE(y)
493 CALL ropp_pp_polynomial(a, x(i), y_a(i))
494 END DO
495
496 ! 7.3 Perform regression on the residuals
497
498 y_res(:) = y(:) - y_a(:)
499
500 CALL ropp_pp_regression(K, y_res, a_res)
501
502 ! 7.4 Update the regression coefficients
503
504 a(:) = a(:) + a_res(:)
505
506 DEALLOCATE(y_a)
507 DEALLOCATE(y_res)
508 DEALLOCATE(a_res)
509
510 END SUBROUTINE ropp_pp_residual_regression
511
512!****s* PPUtils/ropp_pp_nearest_power2 *
513!
514! NAME
515! ropp_pp_nearest_power2 - Find power of 2 nearest input integer
516!
517! SYNOPSIS
518! n2 = ropp_pp_nearest_power2(n)
519!
520!****
521
522 FUNCTION ropp_pp_nearest_power2(N) RESULT(nearest_power2)
523
524 IMPLICIT NONE
525 INTEGER, INTENT(in) :: N ! given integer number
526 INTEGER :: nearest_power2 ! power of 2 right-nearest to N
527 INTEGER :: i
528
529 i = 1
530
531 DO
532 IF (i >= N) THEN
533 Nearest_Power2 = i
534 RETURN
535 ENDIF
536 i = 2*i
537 ENDDO
538
539 END FUNCTION ropp_pp_nearest_power2
540
541
542!****s* PPUtils/ropp_pp_isnan *
543!
544! NAME
545! ropp_pp_isnan - Test if variable value is NaN
546!
547! SYNOPSIS
548! true_or_false = ropp_pp_isnan(x)
549!
550!****
551
552 FUNCTION ropp_pp_isnan(x) RESULT(it_is)
553
554 USE typesizes, ONLY: wp => EightByteReal
555 IMPLICIT NONE
556
557 REAL(wp), INTENT(in) :: x
558 LOGICAL :: it_is
559
560 it_is = .FALSE.
561
562 IF ( x /= x ) it_is = .TRUE.
563 IF ( x + 1.0_wp == x ) it_is = .TRUE.
564 IF ((x > 0) .EQV. (x <= 0)) it_is=.true.
565
566 END FUNCTION ropp_pp_isnan
567
568END MODULE ropp_pp_utils