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