Neko-TOP
A portable framework for high-order spectral element flow toplogy optimization.
Loading...
Searching...
No Matches
mma_vector.f90
Go to the documentation of this file.
1submodule(mma) mma_vector
2
3contains
4 module subroutine mma_gensub_vector(this, iter, x, df0dx, fval, dfdx)
5 ! ----------------------------------------------------- !
6 ! Generate the approximation sub problem by computing !
7 ! the lower and upper asymtotes and the other necessary !
8 ! parameters (alpha, beta, p0j, q0j, pij, qij, ...). !
9 ! ----------------------------------------------------- !
10 class(mma_t), intent(inout) :: this
11 type(vector_t), intent(in) :: x
12 type(vector_t), intent(in) :: df0dx
13 type(vector_t), intent(in) :: fval
14 type(matrix_t), intent(in) :: dfdx
15 integer, intent(in) :: iter
16 integer :: i, j, ierr
17 type(vector_t) :: globaltmp_m
18
19 call globaltmp_m%init(this%m)
20
21 if (iter .lt. 3) then
22 this%low = x - this%asyinit * (this%xmax - this%xmin)
23 this%upp = x + this%asyinit * (this%xmax - this%xmin)
24 else
25 !Move asymptotes low and upp
26 ! Todo: Port to vectorized operations
27 do j = 1, this%n
28 if ((x%x(j) - this%xold1%x(j))*(this%xold1%x(j) - this%xold2%x(j)) &
29 .lt. 0) then
30 this%low%x(j) = x%x(j) - &
31 this%asydecr * (this%xold1%x(j) - this%low%x(j))
32 this%upp%x(j) = x%x(j) + &
33 this%asydecr * (this%upp%x(j) - this%xold1%x(j))
34
35 else if ((x%x(j) - this%xold1%x(j))* &
36 (this%xold1%x(j) - this%xold2%x(j)) .gt. 0) then
37 this%low%x(j) = x%x(j) - &
38 this%asyincr * (this%xold1%x(j) - this%low%x(j))
39 this%upp%x(j) = x%x(j) + &
40 this%asyincr * (this%upp%x(j) - this%xold1%x(j))
41 else
42 this%low%x(j) = x%x(j) - (this%xold1%x(j) - this%low%x(j))
43 this%upp%x(j) = x%x(j) + (this%upp%x(j) - this%xold1%x(j))
44 end if
45
46 ! setting a minimum and maximum for the low and upp
47 ! asymptotes (eq3.9)
48 this%low%x(j) = max(this%low%x(j), &
49 x%x(j) - 10*(this%xmax%x(j) - this%xmin%x(j)))
50 this%low%x(j) = min(this%low%x(j), &
51 x%x(j) - 0.01*(this%xmax%x(j) - this%xmin%x(j)))
52
53 this%upp%x(j) = min(this%upp%x(j), &
54 x%x(j) + 10*(this%xmax%x(j) - this%xmin%x(j)))
55 this%upp%x(j) = max(this%upp%x(j), &
56 x%x(j) + 0.01*(this%xmax%x(j) - this%xmin%x(j)))
57 end do
58 end if
59 ! we can move alpha and beta out of the following loop if needed as:
60 ! this%alpha = max(this%xmin, this%low + &
61 ! 0.1*(this%x- this%low), this- 0.5*(this%xmax - this%xmin))
62 ! this%beta = min(this%xmax, this%upp - &
63 ! 0.1*(this%upp - this), this+ 0.5*(this%xmax - this%xmin))
64 ! Todo: Port to vectorized operations
65 do j = 1, this%n
66 ! set the the bounds and coefficients for the approximation
67 ! the move bounds (alpha and beta )are slightly more restrictive
68 ! than low and upp. This is done based on eq(3.6)--eq(3.10).
69 ! also check
70 ! https://comsolyar.com/wp-content/uploads/2020/03/gcmma.pdf
71 ! eq (2.8) and (2.9)
72 this%alpha%x(j) = max(this%xmin%x(j), this%low%x(j) + &
73 0.1*(x%x(j)- this%low%x(j)), &
74 x%x(j) - 0.5*(this%xmax%x(j) - this%xmin%x(j)))
75 this%beta%x(j) = min(this%xmax%x(j), this%upp%x(j) - &
76 0.1*(this%upp%x(j) - x%x(j)), &
77 x%x(j) + 0.5*(this%xmax%x(j) - this%xmin%x(j)))
78
79 !Calculate p0j, q0j, pij, qij
80 !where j = 1,2,...,n and i = 1,2,...,m (eq(2.3)-eq(2.5))
81 this%p0j%x(j) = (this%upp%x(j) - x%x(j))**2 * &
82 (1.001_rp*max(df0dx%x(j), 0.0_rp) + &
83 0.001_rp*max(-df0dx%x(j), 0.0_rp) + &
84 (0.00001_rp/(max(0.00001_rp, &
85 (this%xmax%x(j) - this%xmin%x(j))))))
86
87 this%q0j%x(j) = (x%x(j) - this%low%x(j))**2 * &
88 (0.001_rp*max(df0dx%x(j),0.0_rp) + &
89 1.001_rp*max(-df0dx%x(j),0.0_rp) + &
90 (0.00001_rp/(max(0.00001_rp, &
91 (this%xmax%x(j) - this%xmin%x(j))))))
92
93 do i = 1, this%m
94 this%pij%x(i,j) = (this%upp%x(j) - x%x(j))**2 * &
95 (1.001_rp*max(dfdx%x(i, j), 0.0_rp) + &
96 0.001_rp*max(-dfdx%x(i, j), 0.0_rp) + &
97 (0.00001_rp/(max(0.00001_rp, &
98 (this%xmax%x(j) - this%xmin%x(j))))))
99 this%qij%x(i,j) = (x%x(j) - this%low%x(j))**2 * &
100 (0.001_rp*max(dfdx%x(i, j), 0.0_rp) + &
101 1.001_rp*max(-dfdx%x(i, j), 0.0_rp) + &
102 (0.00001_rp/(max(0.00001_rp, &
103 (this%xmax%x(j) - this%xmin%x(j))))))
104 end do
105 end do
106
107 !computing bi as defined in page 5
108 ! Todo: Port to vectorized operations
109 this%bi = 0.0_rp
110 do i = 1, this%m
111 !MPI: here this%n is the global n
112 do j = 1, this%n
113 this%bi%x(i) = this%bi%x(i) + &
114 this%pij%x(i,j) / (this%upp%x(j) - x%x(j)) + &
115 this%qij%x(i,j) / (x%x(j) - this%low%x(j))
116 end do
117 end do
118
119
120 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
121 !!!!Showing that for double precision, bi will be different when!!!!!!!!
122 !!!!!!!!!!!computed in parallel compare to sequential!!!!!!!!!!!!!!!!!!!
123 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
124 ! this%bi= 0_rp
125 ! longbi = 0.0
126 ! do i = 1, this%m
127 ! !MPI: here this%n is the global n
128 ! do j = 1, this%n
129 ! this%bi(i) = this%bi(i) + &
130 ! this%pij%x(i,j)/ (this%upp%x(j) - x(j)) + &
131 ! this%qij%x(i,j)/(x(j) - this%low%x(j))
132 ! longbi(i) = longbi(i) + &
133 ! this%pij%x(i,j)/ (this%upp%x(j) - x(j)) + &
134 ! this%qij%x(i,j)/(x(j) - this%low%x(j))
135 ! end do
136 ! end do
137 ! print *, "bi = ", this%bi%x, "this%n = ", this%n
138 ! print *, "longbi = ", longbi
139 ! ierr = 2160
140 ! longbi = 0.0
141 ! this%bi= 0_rp
142 ! do i = 1, this%m
143 ! do j = 1, ierr
144 ! this%bi(i) = this%bi(i) + &
145 ! this%pij%x(i,j)/ (this%upp%x(j) - x(j)) + &
146 ! this%qij%x(i,j)/(x(j) - this%low%x(j))
147 ! longbi(i) = longbi(i) + &
148 ! this%pij%x(i,j)/ (this%upp%x(j) - x(j)) + &
149 ! this%qij%x(i,j)/(x(j) - this%low%x(j))
150 ! end do
151 ! end do
152 ! print *, "bi = ", this%bi%x, "first batch(1-ierr)"
153 ! print *, "longbi = ", longbi, "first batch(1-ierr)"
154 ! longbiglobal = longbi
155 ! longbi = 0.0
156 ! globaltmp_m= this%bi
157 ! this%bi= 0_rp
158 ! do i = 1, this%m
159 ! do j = ierr+1, this%n
160 ! this%bi(i) = this%bi(i) + &
161 ! this%pij%x(i,j)/ (this%upp%x(j) - x(j)) + &
162 ! this%qij%x(i,j)/(x(j) - this%low%x(j))
163 ! longbi(i) = longbi(i) + &
164 ! this%pij%x(i,j)/ (this%upp%x(j) - x(j)) + &
165 ! this%qij%x(i,j)/(x(j) - this%low%x(j))
166 ! end do
167 ! end do
168 ! print *, "bi = ", this%bi%x, "second batch(ierr+1:end)"
169 ! print *, "longbi = ", longbi, "second batch(ierr+1:end)"
170 ! print *, "bi = ", this%bi+globaltmp_m, "first + second"
171 ! print *, "longbi = ", longbi+longbiglobal, "first + second"
172 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
173 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
174
175
176 globaltmp_m = 0.0_rp
177 call mpi_allreduce(this%bi, globaltmp_m, this%m, &
178 mpi_real_precision, mpi_sum, neko_comm, ierr)
179 this%bi = globaltmp_m - fval
180
181 end subroutine mma_gensub_vector
182
183 subroutine mma_subsolve_dpip_vector(this, designx)
184 ! ------------------------------------------------------- !
185 ! Dual-primal interior point method using Newton's step !
186 ! to solve MMA sub problem. !
187 ! A Backtracking Line Search approach is used to compute !
188 ! the step size; starting with the full Newton's step !
189 ! (delta = 1) and deviding by 2 until we have a step size !
190 ! that leads to a feasible point while ensuring a !
191 ! decrease in the residue. !
192 ! ------------------------------------------------------- !
193 class(mma_t), intent(inout) :: this
194 type(vector_t), intent(inout) :: designx
195 !Note that there is a local dummy "x" in this subroutine, thus, we call
196 !the current design "designx" instead of just "x"
197 integer :: i, j, k, iter, ggdumiter, itto, ierr
198 real(kind=rp) :: epsi, residumax, residunorm, &
199 z, zeta, rez, rezeta, &
200 delz, dz, dzeta, &
201 steg, dummy_one, zold, zetaold, newresidu
202 type(vector_t) :: y, lambda, s, mu, &
203 rey, relambda, remu, res, &
204 dely, dellambda, &
205 dy, dlambda, ds, dmu, &
206 yold, lambdaold, sold, muold, &
207 globaltmp_m
208 type(vector_t) :: x, xsi, eta, &
209 rex, rexsi, reeta, &
210 delx, diagx, dx, dxsi, deta, &
211 xold, xsiold, etaold
212 type(vector_t) :: residu
213 type(vector_t) :: residu_small
214 type(vector_t) :: xx, dxx
215
216 type(vector_t) :: bb
217 type(matrix_t) :: GG
218 type(matrix_t) :: AA
219 type(matrix_t) :: globaltmp_mm
220
221 ! using DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO ) in lapack to solve
222 ! the linear system which needs the following parameters
223 integer :: info
224 integer, dimension(this%m+1) :: ipiv
225
226 real(kind=rp) :: re_xstuff_squ_global
227
228 integer :: nglobal
229
230 call y%init(this%m)
231 call lambda%init(this%m)
232 call s%init(this%m)
233 call mu%init(this%m)
234 call rey%init(this%m)
235 call relambda%init(this%m)
236 call remu%init(this%m)
237 call res%init(this%m)
238 call dely%init(this%m)
239 call dellambda%init(this%m)
240 call dy%init(this%m)
241 call dlambda%init(this%m)
242 call ds%init(this%m)
243 call dmu%init(this%m)
244 call yold%init(this%m)
245 call lambdaold%init(this%m)
246 call sold%init(this%m)
247 call muold%init(this%m)
248 call globaltmp_m%init(this%m)
249
250 call x%init(this%n)
251 call xsi%init(this%n)
252 call eta%init(this%n)
253 call rex%init(this%n)
254 call rexsi%init(this%n)
255 call reeta%init(this%n)
256 call delx%init(this%n)
257 call diagx%init(this%n)
258 call dx%init(this%n)
259 call dxsi%init(this%n)
260 call deta%init(this%n)
261 call xold%init(this%n)
262 call xsiold%init(this%n)
263 call etaold%init(this%n)
264
265 call residu%init(3*this%n+4*this%m+2)
266 call residu_small%init(4*this%m+2)
267 call xx%init(2*this%n+4*this%m+2)
268 call dxx%init(2*this%n+4*this%m+2)
269
270 call bb%init(this%m+1)
271
272 call gg%init(this%m, this%n)
273 call aa%init(this%m+1, this%m+1)
274 call globaltmp_mm%init(this%m, this%m)
275
276 ! intial value for the parameters in the subsolve based on
277 ! page 15 of "https://people.kth.se/~krille/mmagcmma.pdf"
278 dummy_one = 1
279 epsi = 1 !100
280 x = 0.5_rp * (this%alpha + this%beta)
281 y = 1.0_rp
282 z = 1.0_rp
283 zeta = 1.0_rp
284 lambda = 1.0_rp
285 s = 1.0_rp
286 xsi = max(1.0_rp, 1.0_rp / (x - this%alpha))
287 eta = max(1.0_rp, 1.0_rp / (this%beta - x))
288 mu = max(1.0_rp, 0.5_rp * this%c)
289
290 do while (epsi .gt. 0.9*this%epsimin)
291 ! calculating residuals based on
292 ! "https://people.kth.se/~krille/mmagcmma.pdf" for the variables
293 ! x, y, z, lambda residuals based on eq(5.9a)-(5.9d), respectively.
294 rex = ((this%p0j + matmul(transpose(this%pij), &
295 lambda))/(this%upp - x)**2 - &
296 (this%q0j + matmul(transpose(this%qij), &
297 lambda))/(x - this%low)**2 ) - &
298 xsi + eta
299
300 call mpi_allreduce(this%n, nglobal, 1, &
301 mpi_integer, mpi_sum, neko_comm, ierr)
302
303 !!!! computing without matmul and transpose
304 ! rex = 0.0_rp
305 ! do j = 1, this%n
306 ! do i = 1, this%m
307 ! rex(j) = rex(j) + this%pij%x(i,j) * &
308 ! lambda(i)/(this%upp%x(j) - x(j))**2 - &
309 ! this%qij%x(i,j) * lambda(i)/(x(j) - this%low%x(j))**2
310 ! end do
311 ! rex(j) = rex(j) + this%p0j%x(j)/(this%upp%x(j) - x(j))**2 &
312 ! - this%q0j%x(j)/(x(j) - this%low%x(j))**2 &
313 ! - xsi(j) + eta(j)
314 ! end do
315
316
317 rey = this%c + this%d*y - lambda - mu
318 rez = this%a0 - zeta - dot_product(lambda, this%a)
319
320 ! relambda = matmul(this%pij,1.0/(this%upp - x)) + &
321 ! matmul(this%qij, 1.0/(x - this%low)) - &
322 ! this%a*z - y + s - this%bi
323 relambda = 0.0_rp
324 do i = 1, this%m
325 do j = 1, this%n !this n is global
326 ! Accumulate sums for relambda (the term gi(x))
327 relambda(i) = relambda(i) + &
328 this%pij%x(i,j)/(this%upp%x(j) - x(j)) &
329 + this%qij%x(i,j)/(x(j) - this%low%x(j))
330 end do
331 end do
332
333
334 globaltmp_m = 0.0_rp
335 call mpi_allreduce(relambda, globaltmp_m, this%m, &
336 mpi_real_precision, mpi_sum, neko_comm, ierr)
337 relambda = globaltmp_m - this%a*z - y + s - this%bi
338
339
340 rexsi = xsi*(x - this%alpha) - epsi
341 reeta = eta*(this%beta - x) - epsi
342 remu = mu*y - epsi
343 rezeta = zeta*z - epsi
344 res = lambda*s - epsi
345
346 residu = [rex, rey, rez, relambda, rexsi, reeta, remu, rezeta, res]
347 residumax = 0_rp
348
349 call mpi_allreduce(maxval(abs(residu)), residumax, 1, &
350 mpi_real_precision, mpi_max, neko_comm, ierr)
351
352 re_xstuff_squ_global = 0_rp
353 call mpi_allreduce(norm2(rex)**2+norm2(rexsi)**2+norm2(reeta)**2,&
354 re_xstuff_squ_global, 1, mpi_real_precision, mpi_sum,&
355 neko_comm, ierr)
356 residu_small = [rey, rez, relambda, &
357 remu, rezeta, res]
358 residunorm = sqrt(norm2(residu_small)**2 + re_xstuff_squ_global)
359
360
361 do iter = 1, this%max_iter !ittt
362 if (iter .gt. (this%max_iter -2)) then
363 ! print *, "The mma inner loop seems not to converge"
364 ! print *, "residumax = ", residumax, "for epsi = ", epsi, &
365 ! ", ittt = ", iter, "out of ", this%max_iter
366 end if
367 !Check the condition
368 if (residumax .lt. epsi) exit
369
370 delx = 0.0_rp
371 do j = 1, this%n
372 do i = 1, this%m
373 delx(j) = delx(j) + this%pij%x(i,j) * &
374 lambda(i)/(this%upp%x(j) - x(j))**2 &
375 - this%qij%x(i,j) * lambda(i)/(x(j) - this%low%x(j))**2
376 end do
377 delx(j) = delx(j) + this%p0j%x(j)/(this%upp%x(j) - x(j))**2 &
378 - this%q0j%x(j)/(x(j) - this%low%x(j))**2 &
379 - epsi/(x(j) - this%alpha%x(j)) &
380 + epsi/(this%beta%x(j) - x(j))
381 end do
382 dely = this%c+ this%d*y - lambda - epsi/y
383 delz = this%a0 - dot_product(lambda, this%a) - epsi/z
384
385 dellambda = 0.0_rp
386 do i = 1, this%m
387 do j = 1, this%n !this n is global
388 ! Accumulate sums for dellambda (the term gi(x))
389 dellambda(i) = dellambda(i) + &
390 this%pij%x(i,j)/(this%upp%x(j) - x(j)) &
391 + this%qij%x(i,j)/(x(j) - this%low%x(j))
392 end do
393 end do
394
395 globaltmp_m = 0.0_rp
396 call mpi_allreduce(dellambda, globaltmp_m, this%m, &
397 mpi_real_precision, mpi_sum, neko_comm, ierr)
398
399 dellambda = globaltmp_m - this%a*z - y - this%bi + epsi / lambda
400
401 ! delx = ((this%p0j + matmul(transpose(this%pij), &
402 ! lambda))/(this%upp - x)**2 - &
403 ! (this%q0j + matmul(transpose(this%qij), &
404 ! lambda))/(x - this%low)**2 ) - &
405 ! epsi/(x - this%alpha) + epsi/(this%beta - x)
406
407 ! dely = this%c + this%d*y - lambda - epsi/y
408 ! delz = this%a0 - dot_product(lambda, this%a) - epsi/z
409 ! dellambda = matmul(this%pij,1.0/(this%upp - x))+&
410 ! matmul(this%qij, 1.0/(x - this%low)) - &
411 ! this%a*z - y - this%bi + epsi/lambda
412
413 do ggdumiter = 1, this%m
414 gg(ggdumiter, :) = this%pij%x(ggdumiter,:)/ &
415 (this%upp - x)**2 - &
416 this%qij%x(ggdumiter,:)/(x - this%low)**2
417 end do
418
419 diagx = ((this%p0j + matmul(transpose(this%pij), &
420 lambda))/(this%upp - x)**3 + &
421 (this%q0j + matmul(transpose(this%qij), &
422 lambda))/(x - this%low)**3 )
423 diagx = 2*diagx + xsi/(x - this%alpha) + &
424 eta/(this%beta- x)
425
426
427 !Here we only consider the case m<n in the matlab code
428 !assembling the right hand side matrix based on eq(5.20)
429 ! bb = [dellambda + dely/(this%d + &
430 ! (mu/y)) - matmul(GG,delx/diagx), delz ]
431 !!!!!!!!!!!!!!for MPI computation of bb!!!!!!!!!!!!!!!!!!!!!!!!!
432 bb = 0.0_rp
433 do i = 1, this%m
434 do j = 1, this%n ! this n is global
435 bb(i) = bb(i) + gg(i, j) * (delx(j) / diagx(j))
436 end do
437 end do
438 globaltmp_m = 0.0_rp
439 call mpi_allreduce(bb(1:this%m), globaltmp_m, this%m, &
440 mpi_real_precision, mpi_sum, neko_comm, ierr)
441 bb(1:this%m) = globaltmp_m
442
443 bb(1:this%m) = dellambda + dely/(this%d+ (mu/y)) - bb(1:this%m)
444 bb(this%m +1) = delz
445 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
446 ! !assembling the coefficients matrix AA based on eq(5.20)
447 ! AA(1:this%m,1:this%m) = &
448 ! matmul(matmul(GG,mma_diag(1/diagx)), transpose(GG))
449 ! !update diag(AA)
450 ! AA(1:this%m,1:this%m) = AA(1:this%m,1:this%m) + &
451 ! mma_diag(s/lambda + 1.0/(this%d + (mu/y)))
452
453 aa = 0.0_rp
454 !Direct computation of the matrix multiplication
455 !(for better performance)
456 do i = 1, this%m
457 do j = 1, this%m
458 ! Compute the (i, j) element of AA
459 do k = 1, this%n !this n is global
460 aa(i, j) = aa(i, j) + gg(i, k) * &
461 (1.0_rp / diagx(k)) * gg(j, k)
462 end do
463 end do
464 end do
465
466 globaltmp_mm = 0.0_rp
467 call mpi_allreduce(aa(1:this%m, 1:this%m), globaltmp_mm, &
468 this%m*this%m, mpi_real_precision, mpi_sum, neko_comm, ierr)
469 aa(1:this%m,1:this%m) = globaltmp_mm
470 do i = 1, this%m
471 !update the diag AA
472 aa(i, i) = aa(i, i) + (s(i) / lambda(i) + &
473 1.0_rp / (this%d%x(i) + mu(i) / y(i)))
474 end do
475
476 aa(1:this%m, this%m+1) = this%a
477 aa(this%m+1, 1:this%m) = this%a
478 aa(this%m+1, this%m+1) = -zeta/z
479
480
481
482
483 call dgesv(this%m+1, 1, aa, this%m+1, ipiv, bb, this%m+1, info)
484 ! if info! = 0 then DGESV is failed.
485 if (info .ne. 0) then
486 write(stderr, *) "DGESV failed to solve the linear system in MMA."
487 write(stderr, *) "Please check mma_subsolve_dpip in mma.f90"
488 error stop
489 end if
490
491 dlambda = bb(1:this%m)
492 dz = bb(this%m + 1)
493 ! based on eq(5.19)
494 dx = -delx / diagx - matmul(transpose(gg), dlambda)/diagx
495
496 dy = (-dely + dlambda) / (this%d + (mu / y))
497 dxsi = -xsi + (epsi - dx*xsi) / (x - this%alpha)
498 deta = -eta + (epsi + dx*eta) / (this%beta - x)
499 dmu = -mu + (epsi - mu*dy) / y
500 dzeta = -zeta + (epsi-zeta*dz)/z
501 ds = -s + (epsi - dlambda*s) / lambda
502
503 !2*this%n+4*this%m+2
504 dxx = [dy, dz, dlambda, dxsi, deta, dmu, dzeta, ds]
505 xx = [y, z, lambda, xsi, eta, mu, zeta, s]
506 steg = maxval([dummy_one, -1.01*dxx/xx, -1.01*dx/ &
507 (x - this%alpha), 1.01 * dx / (this%beta - x)])
508 steg = 1.0/steg
509
510 call mpi_allreduce(steg, steg, 1, &
511 mpi_real_precision, mpi_min, neko_comm, ierr)
512
513 xold = x
514 yold = y
515 zold = z
516 lambdaold = lambda
517 xsiold = xsi
518 etaold = eta
519 muold = mu
520 zetaold = zeta
521 sold = s
522
523 !The innermost loop to determine the suitable step length
524 !using the Backtracking Line Search approach
525 newresidu = 2*residunorm
526 itto = 0
527 do while ((newresidu .gt. residunorm) .and. (itto .lt. 50))
528 itto = itto + 1
529 !update the variables
530 x = xold + steg*dx
531 y = yold + steg*dy
532 z = zold + steg*dz
533 lambda = lambdaold + steg*dlambda
534 xsi = xsiold + steg*dxsi
535 eta = etaold + steg*deta
536 mu = muold + steg*dmu
537 zeta = zetaold + steg*dzeta
538 s = sold + steg*ds
539 !recompute the newresidu to see if this stepsize improves
540 !the residue
541 rex = ((this%p0j + matmul(transpose(this%pij), &
542 lambda))/(this%upp - x)**2 - &
543 (this%q0j + matmul(transpose(this%qij), &
544 lambda))/(x - this%low)**2 ) - &
545 xsi + eta
546 rey = this%c + this%d*y - lambda - mu
547 rez = this%a0 - zeta - dot_product(lambda, this%a)
548 ! relambda = matmul(this%pij,1.0/&
549 ! (this%upp - x)) + matmul(this%qij, &
550 ! 1.0/(x - this%low)) - this%a*z - &
551 ! y + s - this%bi
552 relambda = 0.0_rp
553 do i = 1, this%m
554 do j = 1, this%n !this n is global
555 ! Accumulate sums for relambda (the term gi(x))
556 relambda(i) = relambda(i) + &
557 this%pij%x(i,j)/(this%upp%x(j) - x(j)) &
558 + this%qij%x(i,j)/(x(j) - this%low%x(j))
559 end do
560 end do
561 globaltmp_m = 0.0_rp
562 call mpi_allreduce(relambda, globaltmp_m, this%m, &
563 mpi_real_precision, mpi_sum, neko_comm, ierr)
564 relambda = globaltmp_m
565
566
567 relambda = relambda - this%a*z - y + s - this%bi
568
569 rexsi = xsi*(x - this%alpha) - epsi
570 reeta = eta*(this%beta - x) - epsi
571 remu = mu*y - epsi
572 rezeta = zeta*z - epsi
573 res = lambda*s - epsi
574
575 residu = [rex, rey, rez, relambda, &
576 rexsi, reeta, remu, rezeta, res]
577
578 re_xstuff_squ_global = 0_rp
579 call mpi_allreduce(norm2(rex)**2 + &
580 norm2(rexsi)**2+norm2(reeta)**2, re_xstuff_squ_global, &
581 1, mpi_real_precision, mpi_sum, neko_comm, ierr)
582 residu_small = [rey, rez, relambda, &
583 remu, rezeta, res]
584 newresidu = sqrt(norm2(residu_small)**2 + &
585 re_xstuff_squ_global)
586
587 steg = steg/2
588 end do
589
590 residunorm = newresidu
591 residumax = 0_rp
592 call mpi_allreduce(maxval(abs(residu)), residumax, 1, &
593 mpi_real_precision, mpi_max, neko_comm, ierr)
594
595 !correct the step size for the extra devision by 2 in the final
596 !loop
597 steg = 2*steg
598
599 ! print *,"Processor ",pe_rank, "iter = ", iter, "epsi = ", epsi, &
600 ! "steg = ", steg, "residunorm = ",residunorm, &
601 ! "residumax = ",residumax
602 end do
603 epsi = 0.1*epsi
604
605 end do
606
607 ! Save the new design
608 this%xold2 = this%xold1
609 this%xold1 = designx
610 designx = x
611
612 !update the parameters of the MMA object nesessary to compute KKT residu
613 this%y = y
614 this%z = z
615 this%lambda = lambda
616 this%zeta = zeta
617 this%xsi = xsi
618 this%eta = eta
619 this%mu = mu
620 this%s = s
621
622 end subroutine mma_subsolve_dpip_vector
623
624 subroutine mma_kkt_vector(this, x, df0dx, fval, dfdx)
625 ! ----------------------------------------------------- !
626 ! Compute the KKT condition right hand side for a given !
627 ! design x and set the max and norm values of the !
628 ! residue of KKT system to this%residumax and !
629 ! this%residunorm. !
630 ! !
631 ! The left hand sides of the KKT conditions are computed!
632 ! for the following nonlinear programming problem: !
633 ! Minimize f_0(x) + a_0*z + !
634 ! sum( c_i*y_i + 0.5*d_i*(y_i)^2 )!
635 ! subject to f_i(x) - a_i*z - y_i <= 0, i = 1,...,m !
636 ! xmax_j <= x_j <= xmin_j, j = 1,...,n !
637 ! z >= 0, y_i >= 0, i = 1,...,m !
638 ! !
639 ! !
640 ! Note that before calling this function, the function !
641 ! values (f0val, fval, dfdx, ...) should be updated !
642 ! using the new x values. !
643 ! ----------------------------------------------------- !
644 class(mma_t), intent(inout) :: this
645 type(vector_t), intent(in) :: x
646 type(vector_t), intent(in) :: fval
647 type(vector_t), intent(in) :: df0dx
648 type(matrix_t), intent(in) :: dfdx
649
650
651 real(kind=rp) :: rez, rezeta
652 type(vector_t) :: rey, relambda, remu, res
653 type(vector_t) :: rex, rexsi, reeta
654 type(vector_t) :: residu
655
656 type(vector_t) :: residu_small
657 integer :: ierr
658 real(kind=rp) :: re_xstuff_squ_global
659
660 call rey%init(this%m)
661 call relambda%init(this%m)
662 call remu%init(this%m)
663 call res%init(this%m)
664
665 call rex%init(this%n)
666 call rexsi%init(this%n)
667 call reeta%init(this%n)
668
669 call residu%init(3*this%n+4*this%m+2)
670 call residu_small%init(4*this%m+2)
671
672
673 rex = df0dx+ matmul(transpose(dfdx), this%lambda) - this%xsi + &
674 this%eta
675 rey = this%c + this%d*this%y - this%lambda - this%mu
676 rez = this%a0 - this%zeta - dot_product(this%lambda, this%a)
677
678 relambda = fval- this%a*this%z - this%y + this%s
679 rexsi = this%xsi*(x - this%xmin)
680 reeta = this%eta*(this%xmax - x)
681 remu = this%mu*this%y
682 rezeta = this%zeta*this%z
683 res = this%lambda*this%s
684
685 residu = [rex, rey, rez, relambda, rexsi, reeta, remu, rezeta, res]
686
687 call mpi_allreduce(maxval(abs(residu)), this%residumax, 1, &
688 mpi_real_precision, mpi_max, neko_comm, ierr)
689
690 call mpi_allreduce(norm2(rex)**2 + norm2(rexsi)**2 + norm2(reeta)**2, &
691 re_xstuff_squ_global, 1, mpi_real_precision, mpi_sum, neko_comm, ierr)
692
693 residu_small = [rey, rez, relambda, remu, rezeta, res]
694
695 this%residunorm = sqrt(norm2(residu_small)**2 + re_xstuff_squ_global)
696
697 end subroutine mma_kkt_vector
698end submodule mma_vector
mma
Definition mma.f90:34