Neko-TOP
A portable framework for high-order spectral element flow toplogy optimization.
Loading...
Searching...
No Matches
mma_cpu.f90
1
33
34submodule(mma) mma_cpu
35 use lapack_interfaces, only: dgesv
36 use mpi_f08, only: mpi_in_place, mpi_max, mpi_min
37 use comm, only: neko_comm, pe_rank, mpi_real_precision
38 use math, only: neko_eps
39 use profiler, only: profiler_start_region, profiler_end_region
40
41 implicit none
42
43contains
44
45 module subroutine mma_update_cpu(this, iter, x, df0dx, fval, dfdx)
46 ! ----------------------------------------------------- !
47 ! Update the design variable x by solving the convex !
48 ! approximation of the problem. !
49 ! !
50 ! This subroutine is called in each iteration of the !
51 ! optimization loop !
52 ! ----------------------------------------------------- !
53 class(mma_t), intent(inout) :: this
54 integer, intent(in) :: iter
55 real(kind=rp), dimension(this%n), intent(inout) :: x
56 real(kind=rp), dimension(this%n), intent(in) :: df0dx
57 real(kind=rp), dimension(this%m), intent(in) :: fval
58 real(kind=rp), dimension(this%m, this%n), intent(in) :: dfdx
59
60 if (.not. this%is_initialized) then
61 call neko_error("The MMA object is not initialized.")
62 end if
63
64 call profiler_start_region("MMA update")
65
66 ! generate a convex approximation of the problem
67 call profiler_start_region("MMA gensub")
68 call mma_gensub_cpu(this, iter, x, df0dx, fval, dfdx)
69 call profiler_end_region("MMA gensub")
70
71 !solve the approximation problem using interior point method
72 call profiler_start_region("MMA subsolve")
73 if (this%subsolver .eq. "dip") then
74 call mma_subsolve_dip_cpu(this, x)
75 else if (this%subsolver .eq. "pdip") then
76 call mma_subsolve_pdip_cpu(this, x)
77 else
78 call neko_error("Unrecognized subsolver for MMA in mma_cpu.")
79 end if
80 call profiler_end_region("MMA subsolve")
81
82 call profiler_end_region("MMA update")
83
84 this%is_updated = .true.
85 end subroutine mma_update_cpu
86
88 module subroutine mma_kkt_cpu(this, x, df0dx, fval, dfdx)
89 ! ----------------------------------------------------- !
90 ! Compute the KKT condition right hand side for a given !
91 ! designx x and set the max and norm values of the !
92 ! residue of KKT system to this%residumax and !
93 ! this%residunorm. !
94 ! !
95 ! The left hand sides of the KKT conditions are computed!
96 ! for the following nonlinear programming problem: !
97 ! Minimize f_0(x) + a_0*z + !
98 ! sum(c_i*y_i + 0.5*d_i*(y_i)^2)!
99 ! subject to f_i(x) - a_i*z - y_i <= 0, i = 1,...,m !
100 ! xmax_j <= x_j <= xmin_j, j = 1,...,n !
101 ! z >= 0, y_i >= 0, i = 1,...,m !
102 ! !
103 ! !
104 ! Note that before calling this function, the function !
105 ! values (f0val, fval, dfdx, ...) should be updated !
106 ! using the new x values. !
107 ! ----------------------------------------------------- !
108 class(mma_t), intent(inout) :: this
109 real(kind=rp), dimension(this%n), intent(in) :: x
110 real(kind=rp), dimension(this%n), intent(in) :: df0dx
111 real(kind=rp), dimension(this%m), intent(in) :: fval
112 real(kind=rp), dimension(this%m, this%n), intent(in) :: dfdx
113
114 call profiler_start_region("MMA KKT computation")
115
116 if (this%subsolver .eq. "dip") then
117 call mma_dip_kkt_cpu(this, x, df0dx, fval, dfdx)
118 else if (this%subsolver .eq. "pdip") then
119 call mma_pdip_kkt_cpu(this, x, df0dx, fval, dfdx)
120 else
121 call neko_error("Unrecognized subsolver for MMA in mma_cpu.")
122 end if
123
124 call profiler_end_region("MMA KKT computation")
125 end subroutine mma_kkt_cpu
126
128 ! point method (pdip) subsolve of MMA algorithm.
129 module subroutine mma_pdip_kkt_cpu(this, x, df0dx, fval, dfdx)
130 ! ----------------------------------------------------- !
131 ! Compute the KKT condition right hand side for a given !
132 ! designx x and set the max and norm values of the !
133 ! residue of KKT system to this%residumax and !
134 ! this%residunorm. !
135 ! !
136 ! The left hand sides of the KKT conditions are computed!
137 ! for the following nonlinear programming problem: !
138 ! Minimize f_0(x) + a_0*z + !
139 ! sum( c_i*y_i + 0.5*d_i*(y_i)^2 )!
140 ! subject to f_i(x) - a_i*z - y_i <= 0, i = 1,...,m !
141 ! xmax_j <= x_j <= xmin_j, j = 1,...,n !
142 ! z >= 0, y_i >= 0, i = 1,...,m !
143 ! !
144 ! !
145 ! Note that before calling this function, the function !
146 ! values (f0val, fval, dfdx, ...) should be updated !
147 ! using the new x values. !
148 ! ----------------------------------------------------- !
149 class(mma_t), intent(inout) :: this
150 real(kind=rp), dimension(this%n), intent(in) :: x
151 real(kind=rp), dimension(this%n), intent(in) :: df0dx
152 real(kind=rp), dimension(this%m), intent(in) :: fval
153 real(kind=rp), dimension(this%m, this%n), intent(in) :: dfdx
154
155 real(kind=rp) :: rez, rezeta
156 real(kind=rp), dimension(this%m) :: rey, relambda, remu, res
157 real(kind=rp), dimension(this%n) :: rex, rexsi, reeta
158 real(kind=rp), dimension(3*this%n+4*this%m+2) :: residual
159
160 real(kind=rp), dimension(4*this%m+2) :: residual_small
161 integer :: ierr
162 real(kind=rp) :: re_sq_norm
163
164 rex = df0dx + matmul(transpose(dfdx), this%lambda%x) &
165 - this%xsi%x + this%eta%x
166 rey = this%c%x + this%d%x*this%y%x - this%lambda%x - this%mu%x
167 rez = this%a0 - this%zeta - dot_product(this%lambda%x, this%a%x)
168
169 relambda = fval - this%a%x * this%z - this%y%x + this%s%x
170 rexsi = this%xsi%x * (x - this%xmin%x)
171 reeta = this%eta%x * (this%xmax%x - x)
172 remu = this%mu%x * this%y%x
173 rezeta = this%zeta * this%z
174 res = this%lambda%x * this%s%x
175
176 residual = [rex, rey, rez, relambda, rexsi, reeta, remu, rezeta, res]
177 residual_small = [rey, rez, relambda, remu, rezeta, res]
178
179 this%residumax = maxval(abs(residual))
180 re_sq_norm = norm2(rex)**2 + norm2(rexsi)**2 + norm2(reeta)**2
181
182 call mpi_allreduce(mpi_in_place, this%residumax, 1, &
183 mpi_real_precision, mpi_max, neko_comm, ierr)
184
185 call mpi_allreduce(mpi_in_place, re_sq_norm, 1, &
186 mpi_real_precision, mpi_sum, neko_comm, ierr)
187
188 this%residunorm = sqrt(norm2(residual_small)**2 + re_sq_norm)
189 end subroutine mma_pdip_kkt_cpu
190
192 ! point method (dip) subsolve of MMA algorithm.
193 module subroutine mma_dip_kkt_cpu(this, x, df0dx, fval, dfdx)
194 ! ----------------------------------------------------- !
195 ! Compute the KKT condition right hand side for a given !
196 ! designx x and set the max and norm values of the !
197 ! residue of KKT system to this%residumax and !
198 ! this%residunorm. !
199 ! !
200 ! The left hand sides of the KKT conditions are computed!
201 ! for the following nonlinear programming problem: !
202 ! Minimize f_0(x) + a_0*z + !
203 ! sum( c_i*y_i + 0.5*d_i*(y_i)^2 )!
204 ! subject to f_i(x) - a_i*z - y_i <= 0, i = 1,...,m !
205 ! xmax_j <= x_j <= xmin_j, j = 1,...,n !
206 ! z >= 0, y_i >= 0, i = 1,...,m !
207 ! !
208 ! !
209 ! Note that before calling this function, the function !
210 ! values (f0val, fval, dfdx, ...) should be updated !
211 ! using the new x values. !
212 ! ----------------------------------------------------- !
213 class(mma_t), intent(inout) :: this
214 real(kind=rp), dimension(this%n), intent(in) :: x
215 real(kind=rp), dimension(this%n), intent(in) :: df0dx
216 real(kind=rp), dimension(this%m), intent(in) :: fval
217 real(kind=rp), dimension(this%m, this%n), intent(in) :: dfdx
218
219 real(kind=rp), dimension(this%m) :: relambda, remu
220 real(kind=rp), dimension(2*this%m) :: residual
221
222
223 relambda = fval - this%a%x * this%z - this%y%x + this%mu%x
224 ! Compute residual for mu (eta in the paper)
225 remu = this%lambda%x * this%mu%x
226
227 residual = abs([relambda, remu])
228 this%residumax = maxval(residual)
229 this%residunorm = norm2(residual)
230
231 end subroutine mma_dip_kkt_cpu
232
233 !============================================================================!
234 ! private internal subroutines
235
237 subroutine mma_gensub_cpu(this, iter, x, df0dx, fval, dfdx)
238 ! ----------------------------------------------------- !
239 ! Generate the approximation sub problem by computing !
240 ! the lower and upper asymtotes and the other necessary !
241 ! parameters (alpha, beta, p0j, q0j, pij, qij, ...). !
242 ! ----------------------------------------------------- !
243 class(mma_t), intent(inout) :: this
244 real(kind=rp), dimension(this%n), intent(in) :: x
245 real(kind=rp), dimension(this%n), intent(in) :: df0dx
246 real(kind=rp), dimension(this%m), intent(in) :: fval
247 real(kind=rp), dimension(this%m, this%n), intent(in) :: dfdx
248 integer, intent(in) :: iter
249 integer :: i, j, ierr
250 real(kind=rp), dimension(this%n) :: xmin_eff, xmax_eff
251 real(kind=rp), dimension(this%n) :: x_diff
252 real(kind=rp) :: asy_factor
253
254 xmin_eff = this%xmin%x
255 xmax_eff = this%xmax%x
256 if (this%move_limit .gt. 0.0_rp) then
257 xmin_eff = max(xmin_eff, x - this%move_limit)
258 xmax_eff = min(xmax_eff, x + this%move_limit)
259 end if
260
261 x_diff = max(xmax_eff - xmin_eff, 1.0e-5_rp)
262
263 ! ------------------------------------------------------------------------ !
264 ! Setup the current asymptotes
265 associate(low => this%low%x, upp => this%upp%x, &
266 x_1 => this%xold1%x, x_2 => this%xold2%x)
267
268 if (iter .lt. 3) then
269 ! Initialize the lower and upper asymptotes
270 low = x - this%asyinit * x_diff
271 upp = x + this%asyinit * x_diff
272 else
273 do j = 1, this%n
274 if ((x(j) - x_1(j)) * (x_1(j) - x_2(j)) .lt. 0.0_rp) then
275 asy_factor = this%asydecr
276 else if ((x(j) - x_1(j)) * (x_1(j) - x_2(j)) .gt. 0.0_rp) then
277 asy_factor = this%asyincr
278 else
279 asy_factor = 1.0_rp
280 end if
281
282 low(j) = x(j) - asy_factor * (x_1(j) - low(j))
283 upp(j) = x(j) + asy_factor * (upp(j) - x_1(j))
284 end do
285
286
287 ! Setting a minimum and maximum for the low and upp
288 ! asymptotes (eq3.9)
289 low = max(low, x - 10.0_rp * x_diff)
290 low = min(low, x - 0.01_rp * x_diff)
291
292 upp = min(upp, x + 10.0_rp * x_diff)
293 upp = max(upp, x + 0.01_rp * x_diff)
294 end if
295
296 end associate
297
298
299 ! ------------------------------------------------------------------------ !
300 ! Set the the bounds and coefficients for the approximation
301 ! the move bounds (alpha and beta) are slightly more restrictive
302 ! than low and upp. This is done based on eq(3.6)--eq(3.10).
303 ! also check
304 ! https://comsolyar.com/wp-content/uploads/2020/03/gcmma.pdf
305 ! eq (2.8) and (2.9)
306
307 associate(alpha => this%alpha%x, beta => this%beta%x, &
308 xmin => xmin_eff, xmax => xmax_eff, &
309 low => this%low%x, upp => this%upp%x, x => x)
310
311 alpha = max(xmin, low + 0.1_rp*(x - low), x - 0.5_rp*x_diff)
312 beta = min(xmax, upp - 0.1_rp*(upp - x), x + 0.5_rp*x_diff)
313 end associate
314
315 ! ------------------------------------------------------------------------ !
316 ! Calculate p0j, q0j, pij, qij
317 ! where j = 1,2,...,n and i = 1,2,...,m (eq(2.3)-eq(2.5))
318
319 associate(p0j => this%p0j%x, q0j => this%q0j%x, &
320 pij => this%pij%x, qij => this%qij%x, &
321 low => this%low%x, upp => this%upp%x)
322
323 p0j = ( &
324 1.001_rp * max(df0dx, 0.0_rp) &
325 + 0.001_rp * max(-df0dx, 0.0_rp) &
326 + 0.00001_rp / max(x_diff, 0.00001_rp) &
327 ) * (upp - x)**2
328
329 q0j = ( &
330 0.001_rp * max(df0dx, 0.0_rp) &
331 + 1.001_rp * max(-df0dx, 0.0_rp) &
332 + 0.00001_rp / max(x_diff, 0.00001_rp)&
333 ) * (x - low)**2
334
335 do j = 1, this%n
336 do i = 1, this%m
337 pij(i, j) = ( &
338 1.001_rp * max(dfdx(i, j), 0.0_rp) &
339 + 0.001_rp * max(-dfdx(i, j), 0.0_rp) &
340 + 0.00001_rp / max(x_diff(j), 0.00001_rp) &
341 ) * (upp(j) - x(j))**2
342
343 qij(i, j) = ( &
344 0.001_rp * max(dfdx(i, j), 0.0_rp) &
345 + 1.001_rp * max(-dfdx(i, j), 0.0_rp) &
346 + 0.00001_rp / max(x_diff(j), 0.00001_rp) &
347 ) * (x(j) - low(j))**2
348 end do
349 end do
350
351 end associate
352
353 ! ------------------------------------------------------------------------ !
354 ! Computing bi as defined in page 5
355
356 associate(bi => this%bi%x, &
357 pij => this%pij%x, qij => this%qij%x, &
358 low => this%low%x, upp => this%upp%x)
359
360 bi = 0.0_rp
361 do i = 1, this%m
362 do j = 1, this%n
363 bi(i) = bi(i) &
364 + pij(i, j) / (upp(j) - x(j)) &
365 + qij(i, j) / (x(j) - low(j))
366 end do
367 end do
368
369 call mpi_allreduce(mpi_in_place, bi, this%m, &
370 mpi_real_precision, mpi_sum, neko_comm, ierr)
371 bi = bi - fval
372
373 end associate
374 end subroutine mma_gensub_cpu
375
378 subroutine mma_subsolve_pdip_cpu(this, designx)
379 ! ------------------------------------------------------- !
380 ! Dual-primal interior point method using Newton's step !
381 ! to solve MMA sub problem. !
382 ! A Backtracking Line Search approach is used to compute !
383 ! the step size; starting with the full Newton's step !
384 ! (delta = 1) and dividing by 2 until we have a step size !
385 ! that leads to a feasible point while ensuring a !
386 ! decrease in the residue. !
387 ! ------------------------------------------------------- !
388
389 class(mma_t), intent(inout) :: this
390 real(kind=rp), dimension(this%n), intent(inout) :: designx
391 ! Note that there is a local dummy "x" in this subroutine, thus, we call
392 ! the current design "designx" instead of just "x"
393 integer :: i, j, k, iter, itto, ierr
394 real(kind=rp) :: epsi, residual_max, residual_norm, &
395 z, zeta, rez, rezeta, &
396 delz, dz, dzeta, &
397 steg, zold, zetaold, new_residual
398 real(kind=rp), dimension(this%m) :: y, lambda, s, mu, &
399 rey, relambda, remu, res, &
400 dely, dellambda, &
401 dy, dlambda, ds, dmu, &
402 yold, lambdaold, sold, muold
403 real(kind=rp), dimension(this%n) :: x, xsi, eta, &
404 rex, rexsi, reeta, &
405 delx, diagx, dx, dxsi, deta, &
406 xold, xsiold, etaold
407 real(kind=rp), dimension(4*this%m + 2) :: residual_small
408 real(kind=rp), dimension(3*this%n + 4*this%m + 2) :: residual
409 real(kind=rp), dimension(2*this%n + 4*this%m + 2) :: xx, dxx
410
411 real(kind=rp), dimension(this%m, this%n) :: gg
412 real(kind=rp), dimension(this%m+1) :: bb
413 real(kind=rp), dimension(this%m+1, this%m+1) :: aa
414 real(kind=rp), dimension(this%m * this%m) :: aa_buffer
415
416 ! using DGESV in lapack to solve
417 ! the linear system which needs the following parameters
418 integer :: info
419 integer, dimension(this%m+1) :: ipiv
420
421 ! Parameters for global communication
422 real(kind=rp) :: re_sq_norm
423 real(kind=rp) :: minimal_epsilon
424
425 ! ------------------------------------------------------------------------ !
426 ! initial value for the parameters in the subsolve based on
427 ! page 15 of "https://people.kth.se/~krille/mmagcmma.pdf"
428
429 epsi = 1.0_rp !100
430 x = 0.5_rp * (this%alpha%x + this%beta%x)
431 y = 1.0_rp
432 z = 1.0_rp
433 zeta = 1.0_rp
434 lambda = 1.0_rp
435 s = 1.0_rp
436 xsi = max(1.0_rp, 1.0_rp / (x - this%alpha%x))
437 eta = max(1.0_rp, 1.0_rp / (this%beta%x - x))
438 mu = max(1.0_rp, 0.5_rp * this%c%x)
439
440 ! ------------------------------------------------------------------------ !
441 ! Computing the minimal epsilon and choose the most conservative one
442
443 minimal_epsilon = max(0.9_rp * this%epsimin, 1.0e-12_rp)
444 call mpi_allreduce(mpi_in_place, minimal_epsilon, 1, &
445 mpi_real_precision, mpi_min, neko_comm, ierr)
446
447 ! ------------------------------------------------------------------------ !
448 ! The main loop of the dual-primal interior point method.
449
450 do while (epsi .gt. minimal_epsilon)
451
452 ! --------------------------------------------------------------------- !
453 ! Calculating residuals based on
454 ! "https://people.kth.se/~krille/mmagcmma.pdf" for the variables
455 ! x, y, z, lambda residuals based on eq(5.9a)-(5.9d), respectively.
456
457 associate(p0j => this%p0j%x, q0j => this%q0j%x, &
458 pij => this%pij%x, qij => this%qij%x, &
459 low => this%low%x, upp => this%upp%x, &
460 alpha => this%alpha%x, beta => this%beta%x, &
461 c => this%c%x, d => this%d%x, &
462 a0 => this%a0, a => this%a%x)
463
464 rex = (p0j + matmul(transpose(pij), lambda)) / (upp - x)**2 &
465 - (q0j + matmul(transpose(qij), lambda)) / (x - low)**2 &
466 - xsi + eta
467
468 rey = c + d * y - lambda - mu
469 rez = a0 - zeta - dot_product(lambda, a)
470
471 relambda = 0.0_rp
472 do i = 1, this%m
473 do j = 1, this%n
474 ! Accumulate sums for relambda (the term gi(x))
475 relambda(i) = relambda(i) &
476 + pij(i, j) / (upp(j) - x(j)) &
477 + qij(i, j) / (x(j) - low(j))
478 end do
479 end do
480
481 end associate
482
483 ! --------------------------------------------------------------------- !
484 ! Computing the norm of the residuals
485
486 ! Complete the computations of lambda residuals
487 call mpi_allreduce(mpi_in_place, relambda, this%m, &
488 mpi_real_precision, mpi_sum, neko_comm, ierr)
489 relambda = relambda - this%a%x*z - y + s - this%bi%x
490
491 rexsi = xsi * (x - this%alpha%x) - epsi
492 reeta = eta * (this%beta%x - x) - epsi
493 remu = mu * y - epsi
494 rezeta = zeta * z - epsi
495 res = lambda * s - epsi
496
497 ! Setup vectors of residuals and their norms
498 residual = [rex, rey, rez, relambda, rexsi, reeta, remu, rezeta, res]
499 residual_small = [rey, rez, relambda, remu, rezeta, res]
500
501 residual_max = maxval(abs(residual))
502 re_sq_norm = norm2(rex)**2 + norm2(rexsi)**2 + norm2(reeta)**2
503
504 call mpi_allreduce(mpi_in_place, residual_max, 1, &
505 mpi_real_precision, mpi_max, neko_comm, ierr)
506
507 call mpi_allreduce(mpi_in_place, re_sq_norm, &
508 1, mpi_real_precision, mpi_sum, neko_comm, ierr)
509
510 residual_norm = sqrt(norm2(residual_small)**2 + re_sq_norm)
511
512 ! --------------------------------------------------------------------- !
513 ! Internal loop
514
515 do iter = 1, this%max_iter
516
517 !Check the condition
518 if (residual_max .lt. epsi) exit
519
520 delx = 0.0_rp
521 do j = 1, this%n
522 do i = 1, this%m
523 delx(j) = delx(j) &
524 + this%pij%x(i,j) * lambda(i) / (this%upp%x(j) - x(j))**2 &
525 - this%qij%x(i,j) * lambda(i) / (x(j) - this%low%x(j))**2
526 end do
527 end do
528
529 delx = delx &
530 + this%p0j%x / (this%upp%x - x)**2 &
531 - this%q0j%x / (x - this%low%x)**2 &
532 - epsi / (x - this%alpha%x) &
533 + epsi / (this%beta%x - x)
534
535 dely = this%c%x + this%d%x * y - lambda - epsi / y
536 delz = this%a0 - dot_product(lambda, this%a%x) - epsi / z
537
538 ! Accumulate sums for dellambda (the term gi(x))
539 dellambda = 0.0_rp
540 do i = 1, this%m
541 do j = 1, this%n
542 dellambda(i) = dellambda(i) &
543 + this%pij%x(i, j) / (this%upp%x(j) - x(j)) &
544 + this%qij%x(i, j) / (x(j) - this%low%x(j))
545 end do
546 end do
547
548 call mpi_allreduce(mpi_in_place, dellambda, this%m, &
549 mpi_real_precision, mpi_sum, neko_comm, ierr)
550
551 dellambda = dellambda - this%a%x*z - y - this%bi%x + epsi / lambda
552
553 do i = 1, this%m
554 gg(i,:) = this%pij%x(i,:) / (this%upp%x - x)**2 &
555 - this%qij%x(i,:) / (x - this%low%x)**2
556 end do
557
558 diagx = &
559 (this%p0j%x + matmul(transpose(this%pij%x), lambda)) &
560 / (this%upp%x - x)**3 &
561 + (this%q0j%x + matmul(transpose(this%qij%x), lambda)) &
562 / (x - this%low%x)**3
563
564 diagx = 2.0_rp * diagx &
565 + xsi / (x - this%alpha%x) &
566 + eta / (this%beta%x - x)
567
568
569 !Here we only consider the case m<n in the matlab code
570 !assembling the right hand side matrix based on eq(5.20)
571 ! bb = [dellambda + dely/(this%d%x + &
572 ! (mu/y)) - matmul(GG,delx/diagx), delz ]
573
574 !--------------------------------------------------------------------!
575 ! for MPI computation of bb
576
577 bb = 0.0_rp
578 do i = 1, this%m
579 do j = 1, this%n
580 bb(i) = bb(i) + gg(i, j) * (delx(j) / diagx(j))
581 end do
582 end do
583
584 call mpi_allreduce(mpi_in_place, bb, this%m, &
585 mpi_real_precision, mpi_sum, neko_comm, ierr)
586
587 bb(1:this%m) = dellambda + dely / (this%d%x + mu / y) - bb(1:this%m)
588 bb(this%m + 1) = delz
589
590 !--------------------------------------------------------------------!
591 ! assembling the coefficients matrix AA based on eq(5.20)
592 ! AA(1:this%m,1:this%m) = &
593 ! matmul(matmul(GG,mma_diag(1/diagx)), transpose(GG))
594 ! !update diag(AA)
595 ! AA(1:this%m,1:this%m) = AA(1:this%m,1:this%m) + &
596 ! mma_diag(s/lambda + 1.0/(this%d%x + (mu/y)))
597
598 aa = 0.0_rp
599 ! Direct computation of the matrix multiplication
600 ! (for better performance)
601 do i = 1, this%m
602 do j = 1, this%m
603 ! Compute the (i, j) element of AA
604 do k = 1, this%n !this n is global
605 aa(i, j) = aa(i, j) &
606 + gg(i, k) * (1.0_rp / diagx(k)) * gg(j, k)
607 end do
608 end do
609 end do
610
611 aa_buffer = reshape(aa(1:this%m, 1:this%m), [this%m * this%m])
612
613 call mpi_allreduce(mpi_in_place, aa_buffer, &
614 this%m*this%m, mpi_real_precision, mpi_sum, neko_comm, ierr)
615
616 aa(1:this%m, 1:this%m) = reshape(aa_buffer, [this%m, this%m])
617
618 do i = 1, this%m
619 ! update the diag AA
620 aa(i, i) = aa(i, i) &
621 + s(i) / lambda(i) &
622 + 1.0_rp / (this%d%x(i) + mu(i) / y(i))
623 end do
624
625 aa(1:this%m, this%m+1) = this%a%x
626 aa(this%m+1, 1:this%m) = this%a%x
627 aa(this%m+1, this%m+1) = - zeta/z
628
629 call dgesv(this%m + 1, 1, aa, this%m + 1, ipiv, bb, this%m + 1, info)
630
631 if (info .ne. 0) then
632 call neko_error("DGESV failed to solve the linear system in " // &
633 "mma_subsolve_pdip.")
634 end if
635
636
637 dlambda = bb(1:this%m)
638 dz = bb(this%m + 1)
639
640 ! based on eq(5.19)
641 dx = - delx / diagx - matmul(transpose(gg), dlambda) / diagx
642 dy = (-dely + dlambda) / (this%d%x + mu / y)
643
644 dxsi = -xsi + (epsi - dx * xsi) / (x - this%alpha%x)
645 deta = -eta + (epsi + dx * eta) / (this%beta%x - x)
646 dmu = -mu + (epsi - mu * dy) / y
647 dzeta = -zeta + (epsi - zeta * dz) / z
648 ds = -s + (epsi - dlambda * s) / lambda
649
650 dxx = [dy, dz, dlambda, dxsi, deta, dmu, dzeta, ds]
651 xx = [y, z, lambda, xsi, eta, mu, zeta, s]
652
653 steg = 1.0_rp / maxval([ &
654 1.0_rp, &
655 -1.01_rp * dxx / xx, &
656 -1.01_rp * dx / (x - this%alpha%x), &
657 1.01_rp * dx / (this%beta%x - x) &
658 ])
659
660 ! Save the old values
661 xold = x
662 yold = y
663 zold = z
664 lambdaold = lambda
665 xsiold = xsi
666 etaold = eta
667 muold = mu
668 zetaold = zeta
669 sold = s
670
671 new_residual = 2.0_rp * residual_norm
672
673 ! Share the new_residual and steg values
674 call mpi_allreduce(mpi_in_place, steg, 1, &
675 mpi_real_precision, mpi_min, neko_comm, ierr)
676 call mpi_allreduce(mpi_in_place, new_residual, 1, &
677 mpi_real_precision, mpi_min, neko_comm, ierr)
678
679 ! The innermost loop to determine the suitable step length
680 ! using the Backtracking Line Search approach
681 itto = 0
682 do while ((new_residual .gt. residual_norm) .and. (itto .lt. 50))
683 itto = itto + 1
684
685 ! update the variables
686 x = xold + steg*dx
687 y = yold + steg*dy
688 z = zold + steg*dz
689
690 lambda = lambdaold + steg*dlambda
691 xsi = xsiold + steg*dxsi
692 eta = etaold + steg*deta
693 mu = muold + steg*dmu
694 zeta = zetaold + steg*dzeta
695 s = sold + steg*ds
696
697 ! Recompute the new_residual to see if this stepsize improves
698 ! the residue
699 rex = (this%p0j%x + matmul(transpose(this%pij%x), lambda)) &
700 / (this%upp%x - x)**2 &
701 - (this%q0j%x + matmul(transpose(this%qij%x), lambda)) &
702 / (x - this%low%x)**2 &
703 - xsi + eta
704
705 rey = this%c%x + this%d%x*y - lambda - mu
706 rez = this%a0 - zeta - dot_product(lambda, this%a%x)
707
708 ! Accumulate sums for relambda (the term gi(x))
709 relambda = 0.0_rp
710 do i = 1, this%m
711 do j = 1, this%n
712 relambda(i) = relambda(i) &
713 + this%pij%x(i, j) / (this%upp%x(j) - x(j)) &
714 + this%qij%x(i, j) / (x(j) - this%low%x(j))
715 end do
716 end do
717
718 call mpi_allreduce(mpi_in_place, relambda, this%m, &
719 mpi_real_precision, mpi_sum, neko_comm, ierr)
720
721 relambda = relambda - this%a%x*z - y + s - this%bi%x
722
723 rexsi = xsi * (x - this%alpha%x) - epsi
724 reeta = eta * (this%beta%x - x) - epsi
725 remu = mu * y - epsi
726 rezeta = zeta * z - epsi
727 res = lambda * s - epsi
728
729 ! Compute squared norms for the residuals
730 re_sq_norm = norm2(rex)**2 + norm2(rexsi)**2 + norm2(reeta)**2
731 call mpi_allreduce(mpi_in_place, re_sq_norm, &
732 1, mpi_real_precision, mpi_sum, neko_comm, ierr)
733
734 residual_small = [rey, rez, relambda, remu, rezeta, res]
735 new_residual = sqrt(norm2(residual_small)**2 + re_sq_norm)
736 call mpi_allreduce(mpi_in_place, new_residual, 1, &
737 mpi_real_precision, mpi_min, neko_comm, ierr)
738
739 steg = steg / 2.0_rp
740 end do
741 steg = 2.0_rp * steg ! Correction for the final division by 2
742
743 residual = [rex, rey, rez, relambda, rexsi, reeta, remu, rezeta, res]
744
745 ! Update the maximum and norm of the residuals
746 residual_norm = new_residual
747 residual_max = maxval(abs(residual))
748 call mpi_allreduce(mpi_in_place, residual_max, 1, &
749 mpi_real_precision, mpi_max, neko_comm, ierr)
750 end do
751
752 epsi = 0.1_rp * epsi
753 end do
754
755 ! Save the new designx
756 this%xold2%x = this%xold1%x
757 this%xold1%x = designx
758 designx = x
759
760 !update the parameters of the MMA object nesessary to compute KKT residual
761 this%y%x = y
762 this%z = z
763 this%lambda%x = lambda
764 this%zeta = zeta
765 this%xsi%x = xsi
766 this%eta%x = eta
767 this%mu%x = mu
768 this%s%x = s
769
770 end subroutine mma_subsolve_pdip_cpu
771
774 subroutine mma_subsolve_dip_cpu(this, designx)
775 ! ------------------------------------------------------------------------ !
776 ! -------------------------------Dual Solver------------------------------ !
777 ! ------------------------------------------------------------------------ !
778 ! This implementation is based on: !
779 ! https://doi.org/10.1007/s00158-012-0869-2 !
780 ! Definition of the Lagrangian function: !
781 ! (Note that the equation is slightly different with d(i)=1 and a quadratic!
782 ! term for z. This is done to ensure that we have quadratic terms for both !
783 ! y and z.) !
784 ! !
785 ! L(x, y, z, λ) = !
786 ! sum_{j=1}^{n} [ (p_{0j} + sum_{i=1}^{m} λ_i * p_{ij}) / (u_j - x_j)!
787 ! + (q_{0j} + sum_{i=1}^{m} λ_i * q_{ij}) / (x_j - l_j) ]!
788 ! - sum_{i=1}^{m} λ_i * b_i !
789 ! + sum_{i=1}^{m} [ (c_i - λ_i) * y_i + 0.5 * y_i^2 ] !
790 ! + (a_0 - sum_{i=1}^{m} λ_i * a_i) * z + 0.5 * z^2 !
791 ! !
792 ! Breakdown of terms: !
793 ! - Terms related to x: L_x (the first three lines of L(x, y, z, λ)) !
794 ! - Terms related to y: L_y (the fourth line of L(x, y, z, λ)) !
795 ! - Terms related to z: L_z (the last line of L(x, y, z, λ)) !
796 ! !
797 ! Optimization problem if λ is given: !
798 ! !
799 ! Minimize L(x, y, z, λ) !
800 ! subject to: α_j ≤ x_j ≤ β_j, z ≥ 0, and y_i ≥ 0 for all i, j. !
801 ! !
802 ! Since the problem is separable: !
803 ! Ψ(λ) = !
804 ! sum_{j=1}^{n} min_xj {L_x(x_j, λ) | α_j ≤ x_j ≤ β_j} !
805 ! + min_z {L_z(z, λ) | z ≥ 0} !
806 ! + sum_{i=1}^{m} min_yi {L_y(y_i, λ) | y_i ≥ 0} !
807 ! !
808 ! Maximize Ψ(λ) subject to λ_i ≥ 0 for i = 1, ..., m. !
809 ! !
810 ! ------------------------------------------------------------------------ !
811
812 class(mma_t), intent(inout) :: this
813 real(kind=rp), dimension(this%n), intent(inout) :: designx
814 ! Note that there is a local dummy "x" in this subroutine, thus, we call
815 ! the current design "designx" instead of just "x"
816 integer :: i, j, k, iter, ierr
817 real(kind=rp) :: epsi, residual_max, z, steg
818 real(kind=rp), dimension(this%m) :: y, lambda, mu, &
819 relambda, remu, dlambda, dmu, gradlambda
820 real(kind=rp), dimension(this%n) :: x, pjlambda, qjlambda
821
822 ! To compute the Hessian based on eq(13)
823 ! https://doi.org/10.1007/s00158-012-0869-2
824
825 ! inverse of a diag matrix:
826 real(kind=rp), dimension(this%n) :: ljjxinv ! [∇_x^2 Ljj]−1
827 real(kind=rp), dimension(this%m,this%n) :: hijx ! ∇_x hij
828 real(kind=rp), dimension(this%m,this%m) :: hess
829 real(kind=rp) :: hesstrace
830
831
832 ! using DGESV in lapack to solve
833 ! the linear system which needs the following parameters
834 integer :: info
835 integer, dimension(this%m+1) :: ipiv
836
837 ! Parameters for global communication
838 real(kind=rp) :: minimal_epsilon
839
840 ! ------------------------------------------------------------------------ !
841 ! initial value for the parameters in the subsolve based on
842 ! page 15 of "https://people.kth.se/~krille/mmagcmma.pdf"
843
844 epsi = 1.0_rp !100
845 ! x = 0.5_rp * (this%alpha%x + this%beta%x)
846 y = 1.0_rp
847 z = 0.0_rp
848 lambda = max(1.0_rp, 0.5_rp * this%c%x)
849 mu = 1.0_rp !this parameter is eta in Niel's paper
850 ! note that mu in the paper translates to epsi in the code following the
851 ! same style as the Cpp code by Neils
852
853 ! ------------------------------------------------------------------------ !
854 ! Computing the minimal epsilon and choose the most conservative one
855
856 minimal_epsilon = max(0.9_rp * this%epsimin, 1.0e-12_rp)
857 call mpi_allreduce(mpi_in_place, minimal_epsilon, 1, &
858 mpi_real_precision, mpi_min, neko_comm, ierr)
859
860 ! ------------------------------------------------------------------------ !
861 ! The main loop of the dual-primal interior point method.
862
863 do while (epsi .gt. minimal_epsilon)
864
865 ! --------------------------------------------------------------------- !
866 ! Calculating residuals based on
867 ! "https://people.kth.se/~krille/mmagcmma.pdf" for the variables
868 ! x, y, z, lambda residuals based on eq(5.9a)-(5.9d), respectively.
869
870 associate(p0j => this%p0j%x, q0j => this%q0j%x, &
871 pij => this%pij%x, qij => this%qij%x, &
872 low => this%low%x, upp => this%upp%x, &
873 alpha => this%alpha%x, beta => this%beta%x, &
874 c => this%c%x, d => this%d%x, &
875 a0 => this%a0, a => this%a%x, &
876 bi => this%bi%x)
877 ! minimize(L_x, L_y, L_z) and compute x(λ), y(λ), z(λ) for
878 ! the initial value of λ
879
880 ! Comput the value of y that minimizes L_y for the current λ
881 ! minimize (sum_{i=1}^{m} [ (c_i - λ_i) * y_i + 0.5 * y_i^2 ])
882 ! dL_y/dy =0 => y= (λ_i - c_i), ensure y>=0
883 y = max(0.0_rp, lambda - c)
884
885 ! Comput the value of z that minimizes L_z for the current λ
886 ! minimize ((a_0 - sum_{i=1}^{m} λ_i * a_i) * z + 0.5 * z^2)
887 ! ensure z>=0
888 z = max(0.0_rp,dot_product(lambda, a) - a0)
889
890 ! Comput the value of x that minimizes L_x for the current λ
891 ! minimize( sum_{j=1}^{n} [ (p_{0j} + sum_{i=1}^{m} ����_i *
892 ! p_{ij}) / (u_j - x_j) + (q_{0j} + sum_{i=1}^{m} λ_i * q_{ij}) /
893 ! (x_j - l_j) ] - sum_{i=1}^{m} λ_i * b_i)
894 pjlambda = (p0j + matmul(transpose(pij), lambda))
895 qjlambda = (q0j + matmul(transpose(qij), lambda))
896 x = (sqrt(pjlambda) * low + sqrt(qjlambda) * upp) / &
897 (sqrt(pjlambda) + sqrt(qjlambda))
898
899 ! Ensure that x is feasible (alpha<=x<=beta)
900 x = merge(alpha, x, x .lt. alpha)
901 x = merge(beta, x, x .gt. beta)
902
903 ! Compute the residual for the lambda and mu using eq(9) and eq(15)
904 relambda = matmul(pij, 1.0_rp / (upp - x)) + &
905 matmul(qij, 1.0_rp / (x - low))
906
907 ! Global comminucation for relambda values
908 call mpi_allreduce(mpi_in_place, relambda, this%m, &
909 mpi_real_precision, mpi_sum, neko_comm, ierr)
910 relambda = relambda - bi - y - a * z + mu
911
912 ! Compute residual for mu (eta in the paper)
913 remu = mu * lambda - epsi
914
915 residual_max = maxval(abs([relambda, remu]))
916 call mpi_allreduce(mpi_in_place, residual_max, 1, &
917 mpi_real_precision, mpi_max, neko_comm, ierr)
918
919 ! ------------------------------------------------------------------- !
920 ! Internal loop
921 do iter = 1, this%max_iter
922
923 !Check the condition
924 if (residual_max .lt. epsi) exit
925
926 ! Compute dL(x, y, z, λ)/dλ for the updated x(λ), y(λ), z(λ)
927
928 gradlambda = matmul(pij, 1.0_rp / (upp - x)) + &
929 matmul(qij, 1.0_rp/(x - low))
930
931 ! Global comminucation for gradlambda values
932 call mpi_allreduce(mpi_in_place, gradlambda, this%m, &
933 mpi_real_precision, mpi_sum, neko_comm, ierr)
934 gradlambda = gradlambda - bi - y - a * z
935
936 ! Update gradlambda as the right hand side for Newton's method(eq10)
937 gradlambda = - gradlambda - epsi / lambda
938
939 ! Computing the Hessian as in equation (13) in
940 !! https://doi.org/10.1007/s00158-012-0869-2
941
942 !--------------contributions of x terms to Hess--------------------!
943 ljjxinv= - 1.0_rp / ( (2.0_rp * pjlambda/(upp - x)**3) + &
944 (2.0_rp * qjlambda/(x - low)**3))
945
946
947 ! Remove the sensitivity for the active primal constraints
948 ljjxinv = merge(0.0_rp, ljjxinv, x - alpha < neko_eps)
949 ljjxinv = merge(0.0_rp, ljjxinv, beta - x < neko_eps)
950
951 do i = 1, this%m
952 hijx(i,:) = pij(i,:) / (upp - x)**2 &
953 - qij(i,:) / (x - low)**2
954 end do
955
956 hess = 0.0_rp
957 ! Direct computation of the matrix multiplication
958 ! (for better performance)
959 do i = 1, this%m
960 do j = 1, this%m
961 ! Compute the (i, j) element of AA
962 do k = 1, this%n !this n is global
963 hess(i, j) = hess(i, j) &
964 + hijx(i, k) * (ljjxinv(k)) * hijx(j, k)
965 end do
966 end do
967 end do
968
969 call mpi_allreduce(mpi_in_place, hess, &
970 this%m * this%m, mpi_real_precision, mpi_sum, neko_comm, ierr)
971
972 !---------------contributions of z terms to Hess-------------------!
973 ! Only for inactive constraint, we consider contributions to Hess
974 ! based on the cpp code by Niels.
975 if (dot_product(lambda, a) .gt. 0.0_rp) then
976 do i = 1, this%m
977 do j = 1, this%m
978 hess(i, j) = hess(i, j) - a(i) * a(j)
979 end do
980 end do
981 end if
982
983 !---------------contributions of y terms to Hess-------------------!
984 ! Only for inactive constraint, we consider contributions to Hess.
985 do i = 1, this%m
986 if (y(i) .gt. 0.0_rp) then
987 hess(i, i) = hess(i, i) - 1.0_rp
988 end if
989 ! Based on eq(10), note the term (-\Omega \Lambda)
990 hess(i, i) = hess(i, i) - mu(i) / lambda(i)
991 end do
992
993 ! Improve the robustness by stablizing the Hess using
994 ! Levenberg-Marquardt algorithm (heuristically)
995 hesstrace = 0.0_rp
996 do i=1, this%m
997 hesstrace = hesstrace + hess(i, i)
998 end do
999 do i=1, this%m
1000 hess(i,i) = hess(i, i) - &
1001 max(-1.0e-4_rp*hesstrace/this%m, 1.0e-7_rp)
1002 end do
1003
1004 call dgesv(this%m , 1, hess, this%m , ipiv, &
1005 gradlambda, this%m, info)
1006
1007 if (info .ne. 0) then
1008 call neko_error("DGESV failed to solve the linear system in " // &
1009 "mma_subsolve_dip.")
1010 end if
1011 dlambda = gradlambda
1012
1013 ! based on eq(11) for delta eta
1014 dmu = -mu + epsi / lambda - dlambda * mu / lambda
1015
1016 ! Compute the stepsize and update lambda and mu (eta in the paper)
1017
1018 steg = 1.005_rp
1019 do i = 1, this%m
1020 steg = merge(-1.01_rp * dlambda(i) / lambda(i), steg, &
1021 steg < -1.01_rp * dlambda(i) / lambda(i))
1022 steg = merge(-1.01_rp * dmu(i) / mu(i), &
1023 steg, steg < -1.01_rp * dmu(i) / mu(i))
1024 end do
1025
1026 steg = 1.0_rp / steg
1027
1028 lambda = lambda + steg*dlambda
1029 mu = mu + steg*dmu
1030
1031 ! minimize(L_x, L_y, L_z) and compute x(λ), y(λ), z(λ) for
1032 ! the updated values of λ
1033
1034 ! Comput the value of y that minimizes L_y for the current λ
1035 ! minimize (sum_{i=1}^{m} [ (c_i - λ_i) * y_i + 0.5 * y_i^2 ])
1036 ! dL_y/dy =0 => y= (λ_i - c_i), ensure y>=0
1037 y = max(0.0_rp, lambda - c)
1038
1039
1040 ! Comput the value of z that minimizes L_z for the current λ
1041 ! minimize ((a_0 - sum_{i=1}^{m} λ_i * a_i) * z + 0.5 * z^2)
1042 ! ensure z>=0
1043 z = max(0.0_rp,dot_product(lambda, a) - a0)
1044
1045 ! Comput the value of x that minimizes L_x for the current λ
1046 ! minimize( sum_{j=1}^{n} [ (p_{0j} + sum_{i=1}^{m} λ_i *
1047 ! p_{ij}) / (u_j - x_j) + (q_{0j} + sum_{i=1}^{m} λ_i * q_{ij}) /
1048 ! (x_j - l_j) ] - sum_{i=1}^{m} λ_i * b_i)
1049 pjlambda = (p0j + matmul(transpose(pij), lambda))
1050 qjlambda = (q0j + matmul(transpose(qij), lambda))
1051 x = (sqrt(pjlambda) * low + sqrt(qjlambda) * upp) / &
1052 (sqrt(pjlambda) + sqrt(qjlambda))
1053
1054 ! Ensure that x is feasible (alpha<=x<=beta)
1055 x = merge(alpha, x, x .lt. alpha)
1056 x = merge(beta, x, x .gt. beta)
1057
1058 ! Compute the residual for the lambda and mu using eq(9) and eq(15)
1059 relambda = matmul(pij, 1.0_rp / (upp - x)) + &
1060 matmul(qij, 1.0_rp / (x - low))
1061 ! Global comminucation for relambda values
1062 call mpi_allreduce(mpi_in_place, relambda, this%m, &
1063 mpi_real_precision, mpi_sum, neko_comm, ierr)
1064 relambda = relambda - bi - y - a * z + mu
1065
1066 ! Compute residual for mu (eta in the paper)
1067 remu = mu * lambda - epsi
1068
1069 residual_max = maxval(abs([relambda, remu]))
1070 call mpi_allreduce(mpi_in_place, residual_max, 1, &
1071 mpi_real_precision, mpi_max, neko_comm, ierr)
1072 end do
1073 end associate
1074
1075 epsi = 0.1_rp * epsi
1076 end do
1077
1078 ! Save the new designx
1079 this%xold2%x = this%xold1%x
1080 this%xold1%x = designx
1081 designx = x
1082
1083 !update the parameters of the MMA object nesessary to compute KKT residual
1084
1085 this%y%x = y
1086 this%z = z
1087 this%lambda%x = lambda
1088 this%mu%x = mu
1089 end subroutine mma_subsolve_dip_cpu
1090
1091end submodule mma_cpu
lapack_interfaces::dgesv
Definition lapack_interfaces.f90:44
mma
MMA module.
Definition mma.f90:69