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