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 (iv_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 if (trace_use) call da_trace_entry("da_minimise_cg")
45
46 write(unit=stdout,fmt='(A)') 'Minimize cost function using CG method'
47 if (calculate_cg_cost_fn) then
48 write(unit=stdout,fmt='(A)') &
49 'For this run cost function diagnostics will be written'
50 else
51 write(unit=stdout,fmt='(A)') &
52 'For this run cost function diagnostics will not be written'
53 end if
54 write(unit=stdout,fmt=*) ' '
55
56 ! Initialize temporary cv structures:
57 if ( .not. anal_type_randomcv ) xhat = 0.0
58 j_grad = 0.0
59 fhat = 0.0
60
61 call da_allocate_y(iv, jo_grad_y)
62
63
64 call da_calculate_j(it, 0, cv_size, be % cv % size_jb, be % cv % size_je, &
65 xbx, be, iv, xhat, cv, &
66 re, y, j, j_grad, &
67 grid, config_flags )
68
69 ghat = j_grad
70 phat = - ghat
71
72 rrmold = da_dot_cv(cv_size, cv_size_domain, ghat, ghat, grid, &
73 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
74
75 if (j%total == 0.0) go to 1004
76
77 if (it == 1) j_grad_norm_target = sqrt (rrmold)
78
79 write(unit=stdout,fmt='("Starting outer iteration : ",i3)') it
80 write(unit=stdout,fmt=11) j_grad_norm_target,eps(it)*j_grad_norm_target, &
81 j%total, sqrt(rrmold)
82 11 format('Original gradient is ',1PD15.8,/,&
83 'For this outer iteration gradient target is ',1PD15.8,/,&
84 'Starting cost function: ' ,1PD15.8,' gradient= ',1PD15.8)
85 write(unit=stdout,fmt='(A)') &
86 '----------------------------------------------------------'
87 if (calculate_cg_cost_fn) then
88 write(unit=stdout,fmt='(A)') &
89 'Iter Cost Function Gradient Step'
90 else
91 write(unit=stdout,fmt='(A)')'Iter Gradient Step'
92 end if
93
94 !-------------------------------------------------------------------------
95 ! [2.0] iteratively solve for minimum of cost function:
96 !-------------------------------------------------------------------------
97
98 do iter=1, ntmax
99 if (rrmold == 0.0) go to 1002
100 fhat = phat
101
102 call da_transform_vtoy(cv_size, be, grid%ep, fhat, iv, grid%vp, &
103 grid%vv, xbx, y, &
104 grid, config_flags )
105
106 call da_calculate_grady(iv, y , jo_grad_y)
107
108 call da_transform_vtoy_adj(iter, cv_size, be, grid%ep, fhat, iv, &
109 grid%vp, grid%vv, xbx, jo_grad_y, &
110 grid, config_flags)
111
112 je_start = be % cv % size_jb + 1
113 je_end = be % cv % size_jb + be % cv % size_je
114 cv_size_jb = be % cv % size_jb
115
116 fhat(1:cv_size_jb) = - fhat(1:cv_size_jb) + jb_factor*phat(1:cv_size_jb)
117 fhat(je_start:je_end) = - fhat(je_start:je_end) + &
118 je_factor*phat(je_start:je_end)
119 apdotp = da_dot_cv(cv_size, cv_size_domain, fhat, phat, grid, &
120 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
121
122 step = 0.0
123 if (apdotp .gt. 0.0) step = rrmold/apdotp
124 ghat = ghat + step * fhat
125 xhat = xhat + step * phat
126 if (calculate_cg_cost_fn) then
127 j_grad = 0.0
128 call da_calculate_j(it, iter, cv_size, be % cv % size_jb, &
129 be % cv % size_je, &
130 xbx, be, iv, xhat, cv, &
131 re, y, j, j_grad, &
132 grid, config_flags )
133 ob_grad = da_dot_cv(cv_size, cv_size_domain, j_grad, j_grad, &
134 grid, &
135 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
136 ob_grad = sqrt(ob_grad)
137 end if
138
139 rrmnew = da_dot_cv (cv_size, cv_size_domain, ghat, ghat, grid, &
140 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
141
142 rrmnew_norm = sqrt(rrmnew)
143
144 if (rrmnew_norm < eps(it) * j_grad_norm_target) go to 1002
145 ratio = 0.0
146 if (rrmold .gt. 0.0) ratio = rrmnew/rrmold
147
148 phat = - ghat + ratio * phat
149
150 rrmold=rrmnew
151 if (calculate_cg_cost_fn) then
152 write(unit=stdout,fmt=12)iter, j%total, ob_grad, step
153 else
154 write(unit=stdout,fmt=14)iter, rrmnew_norm , step
155 end if
156 12 format(i3,5x,1PD15.8,5x,1PD15.8,5x,1PD15.8)
157 14 format(i3,5x,1PD15.8,5x,1PD15.8)
158 end do
159
160 !-------------------------------------------------------------------------
161 ! End of the minimization of cost function
162 !-------------------------------------------------------------------------
163 iter = iter -1
164 go to 1003
165 1002 continue
166 if (calculate_cg_cost_fn) then
167 write(unit=stdout,fmt=12)iter, j%total, ob_grad, step
168 else
169 write(unit=stdout,fmt=14)iter, rrmnew_norm , step
170 end if
171 1003 continue
172 write(unit=stdout,fmt='(A)') &
173 '----------------------------------------------------------'
174 write(unit=stdout,fmt='(A)') " "
175 write(unit=stdout, &
176 fmt='("Inner iteration stopped after ",i4," iterations")') iter
177 write(unit=stdout,fmt='(A)') " "
178
179 if (calculate_cg_cost_fn) then
180 rrmnew_norm = ob_grad
181 else
182 call da_calculate_j(it, iter, cv_size, be % cv % size_jb, &
183 be % cv % size_je, &
184 xbx, be, iv, xhat, cv, &
185 re, y, j, j_grad, &
186 grid, config_flags )
187
188 rrmnew_norm = da_dot_cv(cv_size, cv_size_domain, j_grad, j_grad, &
189 grid, &
190 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
191 rrmnew_norm = sqrt(rrmnew_norm)
192 end if
193
194 write(unit=stdout,fmt=15) iter, j%total , rrmnew_norm
195 15 format('Final: ',I3,' iter, J=',1PD15.8,', g=',1PD15.8)
196 write(unit=stdout,fmt='(A)') &
197 '----------------------------------------------------------'
198 1004 continue
199
200 cv = cv + xhat
201
202 call da_deallocate_y(jo_grad_y)
203
204 if (trace_use) call da_trace_exit("da_minimise_cg")
205
206 end subroutine da_minimise_cg
207
208