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