da_minimise_cg.inc
References to this file elsewhere.
1 subroutine da_minimise_cg(grid, config_flags, &
2 it, cv_size, xbx, be, iv, &
3 j_grad_norm_target, xhat, cv, &
4 re, y, j)
5
6 !-------------------------------------------------------------------------
7 ! Purpose: Main Conjugate Gradient minimisation routine
8 !
9 ! Here
10 ! cv is updated in outer-loop.
11 ! xhat is the control variable in inner-loop.
12 !-------------------------------------------------------------------------
13
14 implicit none
15
16 integer, intent(in) :: it ! external iteration.
17 integer, intent(in) :: cv_size ! Total cv size
18 type (xbx_type),intent(in) :: xbx ! Header & non-gridded vars.
19 type (be_type), intent(in) :: be ! background error structure.
20 type (ob_type), intent(inout) :: iv ! ob. increment vector.
21 real, intent(inout) :: j_grad_norm_target ! Target norm.
22 real, intent(out) :: xhat(1:cv_size) ! control variable (local).
23 real, intent(inout) :: cv(1:cv_size) ! control variable (local).
24 type (y_type), intent(inout) :: re ! residual (o-a) structure.
25 type (y_type), intent(inout) :: y ! y = H(x_inc) structure.
26
27 type (j_type), intent(out) :: j ! cost function
28
29 type(domain), intent(inout) :: grid
30 type(grid_config_rec_type), intent(inout) :: config_flags
31
32 integer :: iter
33 integer :: je_start, je_end ! Start/end indices of Je.
34 integer :: cv_size_jb ! end indices of Jb.
35 real :: j_grad(1:cv_size) ! grad_v j (local-grid)
36 real :: fhat(1:cv_size) ! cv copy.
37 real :: ghat(1:cv_size) ! cv copy.
38 real :: phat(1:cv_size) ! cv copy.
39 real :: apdotp,step,rrmold,rrmnew,ratio
40 real :: ob_grad,rrmnew_norm
41
42 type (y_type) :: jo_grad_y ! Grad_y(Jo)
43
44 integer :: cg_jcdf
45
46 if (trace_use) call da_trace_entry("da_minimise_cg")
47
48 write(unit=stdout,fmt='(A)') 'Minimize cost function using CG method'
49 if (calculate_cg_cost_fn) then
50 write(unit=stdout,fmt='(A)') &
51 'For this run cost function diagnostics will be written'
52 else
53 write(unit=stdout,fmt='(A)') &
54 'For this run cost function diagnostics will not be written'
55 end if
56 write(unit=stdout,fmt=*) ' '
57
58 ! Initialize temporary cv structures:
59 if ( .not. anal_type_randomcv ) xhat = 0.0
60 j_grad = 0.0
61 fhat = 0.0
62
63 call da_allocate_y(iv, jo_grad_y)
64
65
66 call da_calculate_j(it, 0, cv_size, be % cv % size_jb, be % cv % size_je, &
67 xbx, be, iv, xhat, cv, &
68 re, y, j, j_grad, &
69 grid, config_flags )
70
71 ghat = j_grad
72 phat = - ghat
73
74 rrmold = da_dot_cv(cv_size, cv_size_domain, ghat, ghat, grid%xp, &
75 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
76
77 if (j%total == 0.0) go to 1004
78
79 if (it == 1) j_grad_norm_target = sqrt (rrmold)
80
81 write(unit=stdout,fmt='("Starting outer iteration : ",i3)') it
82 write(unit=stdout,fmt=11) j_grad_norm_target,eps(it)*j_grad_norm_target, &
83 j%total, sqrt(rrmold)
84 11 format('Original gradient is ',1PD15.8,/,&
85 'For this outer iteration gradient target is ',1PD15.8,/,&
86 'Starting cost function: ' ,1PD15.8,' gradient= ',1PD15.8)
87 write(unit=stdout,fmt='(A)') &
88 '----------------------------------------------------------'
89 if (calculate_cg_cost_fn) then
90 write(unit=stdout,fmt='(1x,/,A)') &
91 'Iter Cost Function Gradient Step'
92 else
93 write(unit=stdout,fmt='(1x,/,A)')'Iter Gradient Step'
94 end if
95 write(unit=stdout,fmt=*) ' '
96
97 !-------------------------------------------------------------------------
98 ! [2.0] iteratively solve for minimum of cost function:
99 !-------------------------------------------------------------------------
100
101 do iter=1, ntmax
102 if (rrmold == 0.0) go to 1002
103 fhat = phat
104 cg_jcdf = 1
105
106 call da_transform_vtoy(cv_size, be, grid%ep, fhat, iv, grid%vp, &
107 grid%vv, grid%xa, grid%xb, xbx, grid%xp, y, &
108 grid, config_flags )
109
110 call da_calculate_grady(iv, y , jo_grad_y)
111
112 call da_transform_vtoy_adj(iter, cv_size, be, grid%ep, fhat, iv, &
113 grid%vp, grid%vv, grid%xa, grid%xb, xbx, grid%xp, jo_grad_y, &
114 grid, config_flags, cg_jcdf)
115
116 je_start = be % cv % size_jb + 1
117 je_end = be % cv % size_jb + be % cv % size_je
118 cv_size_jb = be % cv % size_jb
119
120 fhat(1:cv_size_jb) = - fhat(1:cv_size_jb) + jb_factor*phat(1:cv_size_jb)
121 fhat(je_start:je_end) = - fhat(je_start:je_end) + &
122 je_factor*phat(je_start:je_end)
123 apdotp = da_dot_cv(cv_size, cv_size_domain, fhat, phat, grid%xp, &
124 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
125
126 step = 0.0
127 if (apdotp .gt. 0.0) step = rrmold/apdotp
128 ghat = ghat + step * fhat
129 xhat = xhat + step * phat
130 if (calculate_cg_cost_fn) then
131 j_grad = 0.
132 call da_calculate_j(it, iter, cv_size, be % cv % size_jb, &
133 be % cv % size_je, &
134 xbx, be, iv, xhat, cv, &
135 re, y, j, j_grad, &
136 grid, config_flags )
137 ob_grad = da_dot_cv(cv_size, cv_size_domain, j_grad, j_grad, &
138 grid%xp, &
139 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
140 ob_grad = sqrt(ob_grad)
141 end if
142
143 rrmnew = da_dot_cv (cv_size, cv_size_domain, ghat, ghat, grid%xp, &
144 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
145
146 rrmnew_norm = sqrt(rrmnew)
147
148 if (rrmnew_norm < eps(it) * j_grad_norm_target) go to 1002
149 ratio = 0.0
150 if (rrmold .gt. 0.0) ratio = rrmnew/rrmold
151
152 phat = - ghat + ratio * phat
153
154 rrmold=rrmnew
155 if (calculate_cg_cost_fn) then
156 write(unit=stdout,fmt=12)iter, j%total, ob_grad, step
157 else
158 write(unit=stdout,fmt=14)iter, rrmnew_norm , step
159 end if
160 12 format(i3,5x,1PD15.8,5x,1PD15.8,5x,1PD15.8)
161 14 format(i3,5x,1PD15.8,5x,1PD15.8)
162 end do
163
164 !-------------------------------------------------------------------------
165 ! End of the minimization of cost function
166 !-------------------------------------------------------------------------
167 iter = iter -1
168 go to 1003
169 1002 continue
170 if (calculate_cg_cost_fn) then
171 write(unit=stdout,fmt=12)iter, j%total, ob_grad, step
172 else
173 write(unit=stdout,fmt=14)iter, rrmnew_norm , step
174 end if
175 1003 continue
176 write(unit=stdout,fmt='(A)') &
177 '----------------------------------------------------------'
178 write(unit=stdout, &
179 fmt='("Inner iteration stopped after ",i4," iterations")') iter
180
181 if (calculate_cg_cost_fn) then
182 rrmnew_norm = ob_grad
183 else
184 call da_calculate_j(it, iter, cv_size, be % cv % size_jb, &
185 be % cv % size_je, &
186 xbx, be, iv, xhat, cv, &
187 re, y, j, j_grad, &
188 grid, config_flags )
189
190 rrmnew_norm = da_dot_cv(cv_size, cv_size_domain, j_grad, j_grad, &
191 grid%xp, &
192 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
193 rrmnew_norm = sqrt(rrmnew_norm)
194 end if
195
196 write(unit=stdout,fmt=15) iter, j%total , rrmnew_norm
197 15 format('Final: ',I3,' iter, J=',1PD15.8,', g=',1PD15.8)
198 write(unit=stdout,fmt='(A)') &
199 '----------------------------------------------------------'
200 1004 continue
201
202 cv = cv + xhat
203
204 call da_deallocate_y(jo_grad_y)
205
206 if (trace_use) call da_trace_exit("da_minimise_cg")
207
208 end subroutine da_minimise_cg
209
210