MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
mma.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12#include "mma.hpp"
13
14#include "vector.hpp"
15#include "../general/communication.hpp"
16#include "../general/error.hpp"
17
18#include <fstream>
19#include <math.h>
20
21namespace
22{
23// check if C++ 14 or beyond
24#if __cplusplus >= 201402L
25inline ::std::unique_ptr<::mfem::real_t[]> allocArray(int size) { return ::std::make_unique<::mfem::real_t[]>(size); }
26#else
27inline ::std::unique_ptr<::mfem::real_t[]> allocArray(int size) { return ::std::unique_ptr<::mfem::real_t[]>(new ::mfem::real_t[size]); }
28#endif
29}
30
31#ifdef MFEM_USE_LAPACK
32extern "C" void dgesv_(int* nLAP, int* nrhs, double* AA, int* lda,
33 int* ipiv,
34 double* bb, int* ldb, int* info);
35
36extern "C" void sgesv_(int* nLAP, int* nrhs, float* AA, int* lda,
37 int* ipiv,
38 float* bb, int* ldb, int* info);
39#endif
40
41namespace mfem
42{
43
44void solveLU(int nCon, real_t* AA1, real_t* bb1)
45{
46 // Solve linear system with LU decomposition ifndef LAPACK
47 int nLAP = nCon + 1;
48
49 // Convert AA1 to matrix A and bb1 to vector B
50 ::std::unique_ptr<::std::unique_ptr<real_t[]>[]> A(
51 new ::std::unique_ptr<real_t[]>[nLAP]);
52 for (int i = 0; i < nLAP; ++i)
53 {
54 A[i] = allocArray(nLAP);
55 }
56 ::std::unique_ptr<real_t[]> B = allocArray(nLAP);
57 for (int i = 0; i < nLAP; ++i)
58 {
59 for (int j = 0; j < nLAP; ++j)
60 {
61 A[i][j] = AA1[j * nLAP + i];
62 }
63 B[i] = bb1[i];
64 }
65
66 // Perform LU decomposition
67 ::std::unique_ptr<::std::unique_ptr<real_t[]>[]> L(
68 new ::std::unique_ptr<real_t[]>[nLAP]),
69 U(new ::std::unique_ptr<real_t[]>[nLAP]);
70 for (int i = 0; i < nLAP; ++i)
71 {
72 L[i] = allocArray(nLAP);
73 U[i] = allocArray(nLAP);
74 for (int j = 0; j < nLAP; ++j)
75 {
76 L[i][j] = 0.0;
77 U[i][j] = 0.0;
78 }
79 }
80
81 for (int i = 0; i < nLAP; ++i)
82 {
83 for (int k = i; k < nLAP; ++k)
84 {
85 real_t sum = 0.0;
86 for (int j = 0; j < i; ++j)
87 {
88 sum += (L[i][j] * U[j][k]);
89 }
90 U[i][k] = A[i][k] - sum;
91 }
92 for (int k = i; k < nLAP; ++k)
93 {
94 if (i == k)
95 {
96 L[i][i] = 1.0;
97 }
98 else
99 {
100 real_t sum = 0.0;
101 for (int j = 0; j < i; ++j)
102 {
103 sum += (L[k][j] * U[j][i]);
104 }
105 L[k][i] = (A[k][i] - sum) / U[i][i];
106 }
107 }
108 }
109
110 // Check for singular matrix
111 for (int i = 0; i < nLAP; ++i)
112 {
113 if (U[i][i] == 0.0)
114 {
115 MFEM_ABORT("Error: matrix in MMA LU Solve is singular.");
116 }
117 }
118
119 // Forward substitution to solve L * Y = B
120 ::std::unique_ptr<real_t[]> Y=allocArray(nLAP);
121 for (int i = 0; i < nLAP; ++i)
122 {
123 real_t sum = 0.0;
124 for (int j = 0; j < i; ++j)
125 {
126 sum += L[i][j] * Y[j];
127 }
128 Y[i] = (B[i] - sum) / L[i][i];
129 }
130
131 // Backward substitution to solve U * X = Y
132 ::std::unique_ptr<real_t[]> X=allocArray(nLAP);
133 for (int i = nLAP - 1; i >= 0; --i)
134 {
135 real_t sum = 0.0;
136 for (int j = i + 1; j < nLAP; ++j)
137 {
138 sum += U[i][j] * X[j];
139 }
140 X[i] = (Y[i] - sum) / U[i][i];
141 }
142
143
144 // Copy results back to bb1
145 for (int i = 0; i < (nCon + 1); i++)
146 {
147 bb1[i] = X[i];
148 }
149}
150
151void MMA::MMASubSvanberg::AllocSubData(int nvar, int ncon)
152{
153 epsi = 1.0;
154 ittt = itto = itera = 0;
155 raa0 = 0.00001;
156 move = 0.5;
157 albefa = 0.1;
158 xmamieps = 1e-5;
159 ux1 = allocArray(nvar); // ini
160 xl1 = allocArray(nvar); // ini
161 plam = allocArray(nvar); // ini
162 qlam = allocArray(nvar); // ini
163 gvec = allocArray(ncon); // ini
164 residu = allocArray(3 * nvar + 4 * ncon + 2); // ini
165 GG = allocArray(nvar * ncon); // ini
166 delx = allocArray(nvar); // ini
167 dely = allocArray(ncon); // ini
168 dellam = allocArray(ncon); // ini
169 dellamyi = allocArray(ncon);
170 diagx = allocArray(nvar); // ini
171 diagy = allocArray(ncon); // ini
172 diaglamyi = allocArray(ncon); // ini
173 bb = allocArray(nvar + 1);
174 bb1 = allocArray(ncon + 1); // ini
175 Alam = allocArray(ncon * ncon); // ini
176 AA = allocArray((nvar + 1) * (nvar + 1));
177 AA1 = allocArray((ncon + 1) * (ncon + 1)); // ini
178 dlam = allocArray(ncon); // ini
179 dx = allocArray(nvar); // ini
180 dy = allocArray(ncon); // ini
181 dxsi = allocArray(nvar); // ini
182 deta = allocArray(nvar); // ini
183 dmu = allocArray(ncon); // ini
184 Axx = allocArray(nvar * ncon); // ini
185 axz = allocArray(nvar); // ini
186 ds = allocArray(ncon); // ini
187 xx = allocArray(4 * ncon + 2 * nvar + 2); // ini
188 dxx = allocArray(4 * ncon + 2 * nvar + 2); // ini
189 stepxx = allocArray(4 * ncon + 2 * nvar + 2); // ini
190 sum = 0;
191 sum1 = allocArray(nvar);
192 stepalfa = allocArray(nvar); // ini
193 stepbeta = allocArray(nvar); // ini
194 xold = allocArray(nvar); // ini
195 yold = allocArray(ncon); // ini
196 lamold = allocArray(ncon); // ini
197 xsiold = allocArray(nvar); // ini
198 etaold = allocArray(nvar); // ini
199 muold = allocArray(ncon); // ini
200 sold = allocArray(ncon); // ini
201 q0 = allocArray(nvar); // ini
202 p0 = allocArray(nvar); // ini
203 P = allocArray(ncon * nvar); // ini
204 Q = allocArray(ncon * nvar); // ini
205 alfa = allocArray(nvar); // ini
206 beta = allocArray(nvar); // ini
207 xmami = allocArray(nvar);
208 b = allocArray(ncon); // ini
209
210 b_local = allocArray(ncon);
211 gvec_local = allocArray(ncon);
212 Alam_local = allocArray(ncon * ncon);
213 sum_local = allocArray(ncon);
214 sum_global = allocArray(ncon);
215
216
217 for (int i=0; i<(3 * nvar + 4 * ncon + 2); i++)
218 {
219 residu[i]=0.0;
220 }
221}
222
223void MMA::MMASubSvanberg::Update(const real_t* dfdx,
224 const real_t* gx,
225 const real_t* dgdx,
226 const real_t* xmin,
227 const real_t* xmax,
228 const real_t* xval)
229{
230 MMA& mma = this->mma_ref;
231
232 int rank = 0;
233#ifdef MFEM_USE_MPI
234 MPI_Comm_rank(mma.comm, &rank);
235#endif
236
237 int ncon = mma.nCon;
238 int nvar = mma.nVar;
239
240 real_t zero = 0.0;
241
242 ittt = 0;
243 itto = 0;
244 epsi = 1.0;
245 itera = 0;
246 mma.z = 1.0;
247 mma.zet = 1.0;
248
249 for (int i = 0; i < ncon; i++)
250 {
251 b[i] = 0.0;
252 b_local[i] = 0.0;
253 }
254
255 for (int i = 0; i < nvar; i++)
256 {
257 // Calculation of bounds alfa and beta according to:
258 // alfa = max{xmin, low + 0.1(xval-low), xval-0.5(xmax-xmin)}
259 // beta = min{xmax, upp - 0.1(upp-xval), xval+0.5(xmax-xmin)}
260
261 alfa[i] = std::max(std::max(mma.low[i] + albefa * (xval[i] - mma.low[i]),
262 xval[i] - move * (xmax[i] - xmin[i])), xmin[i]);
263 beta[i] = std::min(std::min(mma.upp[i] - albefa * (mma.upp[i] - xval[i]),
264 xval[i] + move * (xmax[i] - xmin[i])), xmax[i]);
265 xmami[i] = std::max(xmax[i] - xmin[i], xmamieps);
266
267 // Calculations of p0, q0, P, Q, and b
268 ux1[i] = mma.upp[i] - xval[i];
269 if (std::fabs(ux1[i]) <= mma.machineEpsilon)
270 {
271 ux1[i] = mma.machineEpsilon;
272 }
273 xl1[i] = xval[i] - mma.low[i];
274 if (std::fabs(xl1[i]) <= mma.machineEpsilon)
275 {
276 xl1[i] = mma.machineEpsilon;
277 }
278 p0[i] = ( std::max(dfdx[i], zero) + 0.001 * (std::max(dfdx[i],
279 zero) + std::max(-dfdx[i], zero)) + raa0 / xmami[i]) * ux1[i] * ux1[i];
280 q0[i] = ( std::max(-dfdx[i], zero) + 0.001 * (std::max(dfdx[i],
281 zero) + std::max(-dfdx[i], zero)) + raa0 / xmami[i]) * xl1[i] * xl1[i];
282 }
283
284 // P = max(dgdx,0)
285 // Q = max(-dgdx,0)
286 // P = P + 0.001(P+Q) + raa0/xmami
287 // Q = Q + 0.001(P+Q) + raa0/xmami
288 for (int i = 0; i < ncon; i++)
289 {
290 for (int j = 0; j < nvar; j++)
291 {
292 // P = P * spdiags(ux2,0,n,n)
293 // Q = Q * spdiags(xl2,0,n,n)
294 P[i * nvar + j] = (std::max(dgdx[i * nvar + j],
295 zero) + 0.001 * (std::max(dgdx[i * nvar + j],
296 zero) + std::max(-1*dgdx[i * nvar + j],
297 zero)) + raa0 / xmami[j]) * ux1[j] * ux1[j];
298 Q[i * nvar + j] = (std::max(-1*dgdx[i * nvar + j],
299 zero) + 0.001 * (std::max(dgdx[i * nvar + j],
300 zero) + std::max(-1*dgdx[i * nvar + j],
301 zero)) + raa0 / xmami[j]) * xl1[j] * xl1[j];
302 // b = P/ux1 + Q/xl1 - gx
303 b_local[i] = b_local[i] + P[i * nvar + j] / ux1[j] + Q[i * nvar + j] / xl1[j];
304 }
305 }
306
307 std::copy(b_local.get(), b_local.get() + ncon, b.get());
308
309#ifdef MFEM_USE_MPI
310 MPI_Allreduce(b_local.get(), b.get(), ncon, MPITypeMap<real_t>::mpi_type,
311 MPI_SUM,
312 mma.comm);
313#endif
314
315 for (int i = 0; i < ncon; i++)
316 {
317 b[i] = b[i] - gx[i];
318 }
319
320
321 for (int i = 0; i < nvar; i++)
322 {
323 mma.x[i] = 0.5 * (alfa[i] + beta[i]);
324 mma.xsi[i] = 1.0/(mma.x[i] - alfa[i]);
325 mma.xsi[i] = std::max(mma.xsi[i], static_cast<real_t>(1.0));
326 mma.eta[i] = 1.0/(beta[i] - mma.x[i]);
327 mma.eta[i] = std::max(mma.eta[i], static_cast<real_t>(1.0));
328 ux1[i] = 0.0;
329 xl1[i] = 0.0;
330 }
331
332 for (int i = 0; i < ncon; i++)
333 {
334 mma.y[i] = 1.0;
335 mma.lam[i] = 1.0;
336 mma.mu[i] = std::max(1.0, 0.5 * mma.c[i]);
337 mma.s[i] = 1.0;
338 }
339
340 while (epsi > mma.epsimin)
341 {
342 residu[nvar + ncon] = mma.a0 - mma.zet; // rez
343 for (int i = 0; i < nvar; i++)
344 {
345 ux1[i] = mma.upp[i] - mma.x[i];
346 if (std::fabs(ux1[i]) < mma.machineEpsilon)
347 {
348 ux1[i] = mma.machineEpsilon;
349 }
350
351 xl1[i] = mma.x[i] - mma.low[i];
352 if (std::fabs(xl1[i]) < mma.machineEpsilon)
353 {
354 xl1[i] = mma.machineEpsilon;
355 }
356
357 // plam = P' * lam, qlam = Q' * lam
358 plam[i] = p0[i];
359 qlam[i] = q0[i];
360 for (int j = 0; j < ncon; j++)
361 {
362 plam[i] += P[j * nvar + i] * mma.lam[j];
363 qlam[i] += Q[j * nvar + i] * mma.lam[j];
364 residu[nvar + ncon] -= mma.a[j] * mma.lam[j]; // rez
365 }
366 residu[i] = plam[i] / (ux1[i] * ux1[i]) - qlam[i] / (xl1[i] * xl1[i]) -
367 mma.xsi[i] + mma.eta[i]; // rex
368 // residu[nvar + ncon] -= mma.a[i] * mma.lam[i]; // rez
369 residu[nvar + ncon + 1 + ncon + i] = mma.xsi[i] * (mma.x[i] - alfa[i]) -
370 epsi; // rexsi
371 if (std::fabs(mma.x[i]-alfa[i]) < mma.machineEpsilon)
372 {
373 residu[nvar + ncon + 1 + ncon + i] = mma.xsi[i] * mma.machineEpsilon - epsi;
374 }
375 residu[nvar + ncon + 1 + ncon + nvar + i] = mma.eta[i] *
376 (beta[i] - mma.x[i]) - epsi; // reeta
377 if (std::fabs(beta[i] - mma.x[i]) < mma.machineEpsilon)
378 {
379 residu[nvar + ncon + 1 + ncon + nvar + i] = mma.eta[i] * mma.machineEpsilon -
380 epsi;
381 }
382 }
383 for (int i = 0; i < ncon; i++)
384 {
385 gvec_local[i] = 0.0;
386 // gvec = P/ux + Q/xl
387 for (int j = 0; j < nvar; j++)
388 {
389 gvec_local[i] = gvec_local[i] + P[i * nvar + j] / ux1[j] + Q[i * nvar + j] /
390 xl1[j];
391 }
392 }
393
394 std::copy(gvec_local.get(), gvec_local.get() + ncon, gvec.get());
395
396#ifdef MFEM_USE_MPI
397 MPI_Allreduce(gvec_local.get(), gvec.get(), ncon, MPITypeMap<real_t>::mpi_type,
398 MPI_SUM,
399 mma.comm);
400#endif
401
402 if ( rank == 0)
403 {
404 for (int i = 0; i < ncon; i++)
405 {
406 residu[nvar + i] = mma.c[i] + mma.d[i] * mma.y[i] - mma.mu[i] -
407 mma.lam[i]; // rey
408 residu[nvar + ncon + 1 + i] = gvec[i] - mma.a[i] * mma.z - mma.y[i] +
409 mma.s[i] - b[i]; // relam
410 residu[nvar + ncon + 1 + ncon + 2 * nvar + i] = mma.mu[i] * mma.y[i] -
411 epsi; // remu
412 residu[nvar + ncon + 1 + 2 * nvar + 2 * ncon + 1 + i] = mma.lam[i] * mma.s[i]
413 - epsi; // res
414 }
415 residu[nvar + ncon + 1 + 2 * nvar + 2 * ncon] = mma.zet * mma.z - epsi;
416 }
417
418 // Get vector product and maximum absolute value
419 residunorm = 0.0;
420 residumax = 0.0;
421 for (int i = 0; i < (3 * nvar + 4 * ncon + 2); i++)
422 {
423 residunorm += residu[i] * residu[i];
424 residumax = std::max(residumax, std::abs(residu[i]));
425 }
426
427 global_norm = residunorm;
428 global_max = residumax;
429
430#ifdef MFEM_USE_MPI
431 MPI_Allreduce(&residunorm, &global_norm, 1
432 , MPITypeMap<real_t>::mpi_type, MPI_SUM, mma.comm);
433 MPI_Allreduce(&residumax, &global_max, 1
434 , MPITypeMap<real_t>::mpi_type, MPI_MAX, mma.comm);
435#endif
436 // Norm of the residual
437 residunorm = std::sqrt(global_norm);
438 residumax = global_max;
439
440 ittt = 0;
441
442 while (residumax > 0.9 * epsi && ittt < 200)
443 {
444 ittt++;
445 for (int i = 0; i < nvar; i++)
446 {
447 ux1[i] = mma.upp[i] - mma.x[i];
448 if (std::fabs(ux1[i]) < mma.machineEpsilon)
449 {
450 ux1[i] = mma.machineEpsilon;
451 }
452
453 xl1[i] = mma.x[i] - mma.low[i];
454 if (std::fabs(xl1[i]) <= mma.machineEpsilon)
455 {
456 xl1[i] = mma.machineEpsilon;
457 }
458 // plam = P' * lam, qlam = Q' * lam
459 plam[i] = p0[i];
460 qlam[i] = q0[i];
461 for (int j = 0; j < ncon; j++)
462 {
463 plam[i] += P[j * nvar + i] * mma.lam[j];
464 qlam[i] += Q[j * nvar + i] * mma.lam[j];
465 }
466 // NaN-Avoidance
467 if (std::fabs(mma.x[i] - alfa[i]) < mma.machineEpsilon)
468 {
469 if (std::fabs(beta[i] - mma.x[i]) < mma.machineEpsilon)
470 {
471 delx[i] = plam[i] / (ux1[i] * ux1[i]) - qlam[i] / (xl1[i] * xl1[i]);
472 diagx[i] = 2 * (plam[i] / (ux1[i] * ux1[i] * ux1[i]) + qlam[i] /
473 (xl1[i] * xl1[i] * xl1[i])) + mma.xsi[i] / mma.machineEpsilon + mma.eta[i] /
474 mma.machineEpsilon;
475 }
476 else
477 {
478 delx[i] = plam[i] / (ux1[i] * ux1[i]) - qlam[i] / (xl1[i] * xl1[i]) - epsi /
479 mma.machineEpsilon + epsi / (beta[i] - mma.x[i]);
480 diagx[i] = 2 * (plam[i] / (ux1[i] * ux1[i] * ux1[i]) + qlam[i] /
481 (xl1[i] * xl1[i] * xl1[i])) + mma.xsi[i] / (mma.x[i] - alfa[i]) +
482 mma.eta[i] / (beta[i] - mma.x[i]);
483 }
484 }
485 else if (std::fabs(beta[i] - mma.x[i]) < mma.machineEpsilon)
486 {
487 delx[i] = plam[i] / (ux1[i] * ux1[i]) - qlam[i] / (xl1[i] * xl1[i]) - epsi /
488 (mma.x[i] - alfa[i]) + epsi / mma.machineEpsilon;
489 diagx[i] = 2 * (plam[i] / (ux1[i] * ux1[i] * ux1[i]) + qlam[i] /
490 (xl1[i] * xl1[i] * xl1[i])) + mma.xsi[i] / (mma.x[i] - alfa[i]) +
491 mma.eta[i] / mma.machineEpsilon;
492 }
493 else
494 {
495 delx[i] = plam[i] / (ux1[i] * ux1[i]) - qlam[i] / (xl1[i] * xl1[i]) - epsi /
496 (mma.x[i] - alfa[i]) + epsi / (beta[i] - mma.x[i]);
497 diagx[i] = 2 * (plam[i] / (ux1[i] * ux1[i] * ux1[i]) + qlam[i] /
498 (xl1[i] * xl1[i] * xl1[i])) + mma.xsi[i] / (mma.x[i] - alfa[i]) +
499 mma.eta[i] / (beta[i] - mma.x[i]);
500 }
501 }
502
503 for (int i = 0; i < ncon; i++)
504 {
505 gvec_local[i] = 0.0;
506 // gvec = P/ux + Q/xl
507 for (int j = 0; j < nvar; j++)
508 {
509 gvec_local[i] = gvec_local[i] + P[i * nvar + j] / ux1[j] + Q[i * nvar + j] /
510 xl1[j];
511 GG[i * nvar + j] = P[i * nvar + j] / (ux1[j] * ux1[j]) - Q[i * nvar + j] /
512 (xl1[j] * xl1[j]);
513 }
514 }
515
516 std::copy(gvec_local.get(), gvec_local.get() + ncon, gvec.get());
517#ifdef MFEM_USE_MPI
518 MPI_Allreduce(gvec_local.get(), gvec.get(), ncon,
519 MPITypeMap<real_t>::mpi_type, MPI_SUM, mma.comm);
520#endif
521
522 delz = mma.a0 - epsi / mma.z;
523 for (int i = 0; i < ncon; i++)
524 {
525 dely[i] = mma.c[i] + mma.d[i] * mma.y[i] - mma.lam[i] - epsi / mma.y[i];
526 delz -= mma.a[i] * mma.lam[i];
527 dellam[i] = gvec[i] - mma.a[i] * mma.z - mma.y[i] - b[i] + epsi /
528 mma.lam[i];
529 diagy[i] = mma.d[i] + mma.mu[i] / mma.y[i];
530 diaglamyi[i] = mma.s[i] / mma.lam[i] + 1.0 / diagy[i];
531 }
532
533 if (ncon < nVar_global)
534 {
535 // bb1 = dellam + dely./diagy - GG*(delx./diagx);
536 // bb1 = [bb1; delz];
537 for (int j = 0; j < ncon; j++)
538 {
539 sum_local[j] = 0.0;
540 for (int i = 0; i < nvar; i++)
541 {
542 sum_local[j] = sum_local[j] + GG[j * nvar + i] * (delx[i] / diagx[i]);
543 }
544 }
545
546 std::copy(sum_local.get(), sum_local.get() + ncon, sum_global.get());
547
548#ifdef MFEM_USE_MPI
549 MPI_Allreduce(sum_local.get(), sum_global.get(), ncon,
550 MPITypeMap<real_t>::mpi_type, MPI_SUM, mma.comm);
551#endif
552
553 for (int j = 0; j < ncon; j++)
554 {
555 bb1[j] = - sum_global[j] + dellam[j] + dely[j] / diagy[j];
556 }
557 bb1[ncon] = delz;
558
559 // Alam = spdiags(diaglamyi,0,m,m) + GG*spdiags(diagxinv,0,n,n)*GG';
560 for (int i = 0; i < ncon; i++)
561 {
562 // Axx = GG*spdiags(diagxinv,0,n,n);
563 for (int k = 0; k < nvar; k++)
564 {
565 Axx[i * nvar + k] = GG[k * ncon + i] / diagx[k];
566 }
567 }
568 // Alam = spdiags(diaglamyi,0,m,m) + Axx*GG';
569 for (int i = 0; i < ncon; i++)
570 {
571 for (int j = 0; j < ncon; j++)
572 {
573 Alam_local[i * ncon + j] = 0.0;
574 for (int k = 0; k < nvar; k++)
575 {
576 Alam_local[i * ncon + j] += Axx[i * nvar + k] * GG[j * nvar + k];
577 }
578 }
579 }
580
581 std::copy(Alam_local.get(), Alam_local.get() + ncon * ncon, Alam.get());
582#ifdef MFEM_USE_MPI
583 MPI_Reduce(Alam_local.get(), Alam.get(), ncon * ncon,
584 MPITypeMap<real_t>::mpi_type, MPI_SUM, 0, mma.comm);
585#endif
586
587 if (0 == rank)
588 {
589 for (int i = 0; i < ncon; i++)
590 {
591 for (int j = 0; j < ncon; j++)
592 {
593 if (i == j)
594 {
595 Alam[i * ncon + j] += diaglamyi[i];
596 }
597 }
598 }
599 // AA1 = [Alam a ]
600 // [ a' -zet/z]
601 for (int i = 0; i < ncon; i++)
602 {
603 for (int j = 0; j < ncon; j++)
604 {
605 AA1[i * (ncon + 1) + j] = Alam[i * ncon + j];
606 }
607 AA1[i * (ncon + 1) + ncon] = mma.a[i];
608 }
609 for (int i = 0; i < ncon; i++)
610 {
611 AA1[ncon * (ncon + 1) + i] = mma.a[i];
612 }
613 AA1[(ncon + 1) * (ncon + 1) - 1] = -mma.zet / mma.z;
614
615#ifdef MFEM_USE_LAPACK
616 // bb1 = AA1\bb1 --> solve linear system of equations using LAPACK
617 int info;
618 int nLAP = ncon + 1;
619 int nrhs = 1;
620 int lda = nLAP;
621 int ldb = nLAP;
622 int* ipiv = new int[nLAP];
623#if defined(MFEM_USE_DOUBLE)
624 dgesv_(&nLAP, &nrhs, AA1.get(), &lda, ipiv, bb1.get(), &ldb, &info);
625#elif defined(MFEM_USE_SINGLE)
626 sgesv_(&nLAP, &nrhs, AA1.get(), &lda, ipiv, bb1.get(), &ldb, &info);
627#else
628#error "Only single and double precision are supported!"
629#endif
630 if (info == 0)
631 {
632 delete[] ipiv;
633 }
634 else if (info > 0)
635 {
636 MFEM_ABORT("MMA: matrix is singular.");
637 }
638 else
639 {
640 MFEM_ABORT("MMA: Argument " << info <<
641 " in linear system solve has illegal value.");
642 }
643#else
644 solveLU(ncon, AA1.get(), bb1.get());
645#endif
646 }
647#ifdef MFEM_USE_MPI
648 MPI_Bcast(bb1.get(), ncon + 1,
649 MPITypeMap<real_t>::mpi_type, 0, mma.comm);
650#endif
651 // Reassign results
652 for (int i = 0; i < ncon; i++)
653 {
654 dlam[i] = bb1[i];
655 }
656 dz = bb1[ncon];
657
658 // dx = -(GG'*dlam)./diagx - delx./diagx;
659 for (int i = 0; i < nvar; i++)
660 {
661 sum = 0.0;
662 for (int j = 0; j < ncon; j++)
663 {
664 sum = sum + GG[j * nvar + i] * dlam[j];
665 }
666 dx[i] = -sum / diagx[i] - delx[i] / diagx[i];
667 }
668 }
669 else
670 {
671 MFEM_ABORT("MMA: Optimization problem case which has more constraints than design variables is not implemented!");
672 }
673
674 for (int i = 0; i < ncon; i++)
675 {
676 dy[i] = -dely[i] / diagy[i] + dlam[i] / diagy[i];
677 dmu[i] = -mma.mu[i] + epsi / mma.y[i] - (mma.mu[i] * dy[i]) / mma.y[i];
678 ds[i] = -mma.s[i] + epsi / mma.lam[i] - (mma.s[i] * dlam[i]) / mma.lam[i];
679 // xx = [y z lam xsi eta mu zet s]
680 // dxx = [dy dz dlam dxsi deta dmu dzet ds]
681 xx[i] = mma.y[i];
682 xx[ncon + 1 + i] = mma.lam[i];
683 xx[2 * ncon + 1 + 2 * nvar + i] = mma.mu[i];
684 xx[3 * ncon + 2 * nvar + 2 + i] = mma.s[i];
685
686 dxx[i] = dy[i];
687 dxx[ncon + 1 + i] = dlam[i];
688 dxx[2 * ncon + 1 + 2 * nvar + i] = dmu[i];
689 dxx[3 * ncon + 2 * nvar + 2 + i] = ds[i];
690 }
691 xx[ncon] = mma.z;
692 xx[3 * ncon + 2 * nvar + 1] = mma.zet;
693 dxx[ncon] = dz;
694 for (int i = 0; i < nvar; i++)
695 {
696 // NaN-Avoidance
697 if (std::fabs(mma.x[i] - alfa[i]) < mma.machineEpsilon)
698 {
699 if (std::fabs(beta[i] - mma.x[i]) < mma.machineEpsilon)
700 {
701 dxsi[i] = -mma.xsi[i] + epsi / mma.machineEpsilon - (mma.xsi[i] * dx[i]) /
702 mma.machineEpsilon;
703 deta[i] = -mma.eta[i] + epsi / mma.machineEpsilon + (mma.eta[i] * dx[i]) /
704 mma.machineEpsilon;
705 }
706 else
707 {
708 dxsi[i] = -mma.xsi[i] + epsi / mma.machineEpsilon - (mma.xsi[i] * dx[i]) /
709 mma.machineEpsilon;
710 deta[i] = -mma.eta[i] + epsi / (beta[i] - mma.x[i]) +
711 (mma.eta[i] * dx[i]) / (beta[i] - mma.x[i]);
712 }
713 }
714 else if (std::fabs(beta[i] - mma.x[i]) < mma.machineEpsilon)
715 {
716 dxsi[i] = -mma.xsi[i] + epsi / (mma.x[i] - alfa[i]) -
717 (mma.xsi[i] * dx[i]) / (mma.x[i] - alfa[i]);
718 deta[i] = -mma.eta[i] + epsi / mma.machineEpsilon + (mma.eta[i] * dx[i]) /
719 mma.machineEpsilon;
720 }
721 else
722 {
723 dxsi[i] = -mma.xsi[i] + epsi / (mma.x[i] - alfa[i]) -
724 (mma.xsi[i] * dx[i]) / (mma.x[i] - alfa[i]);
725 deta[i] = -mma.eta[i] + epsi / (beta[i] - mma.x[i]) +
726 (mma.eta[i] * dx[i]) / (beta[i] - mma.x[i]);
727 }
728 xx[ncon + 1 + ncon + i] = mma.xsi[i];
729 xx[ncon + 1 + ncon + nvar + i] = mma.eta[i];
730 dxx[ncon + 1 + ncon + i] = dxsi[i];
731 dxx[ncon + 1 + ncon + nvar + i] = deta[i];
732 }
733 dzet = -mma.zet + epsi / mma.z - mma.zet * dz / mma.z;
734 dxx[3 * ncon + 2 * nvar + 1] = dzet;
735
736 stmxx = 0.0;
737 for (int i = 0; i < (4 * ncon + 2 * nvar + 2); i++)
738 {
739 stepxx[i] = -1.01*dxx[i] / xx[i];
740 stmxx = std::max(stepxx[i], stmxx);
741 }
742 stmxx_global = stmxx;
743#ifdef MFEM_USE_MPI
744 MPI_Allreduce(&stmxx, &stmxx_global, 1,
745 MPITypeMap<real_t>::mpi_type, MPI_MAX, mma.comm);
746#endif
747
748 stmalfa = 0.0;
749 stmbeta = 0.0;
750 for (int i = 0; i < nvar; i++)
751 {
752 // NaN-Avoidance
753 if (std::fabs(mma.x[i] - alfa[i]) < mma.machineEpsilon)
754 {
755 stepalfa[i] = -1.01*dx[i] / mma.machineEpsilon;
756 }
757 else
758 {
759 stepalfa[i] = -1.01*dx[i] / (mma.x[i] - alfa[i]);
760 }
761 if (std::fabs(beta[i] - mma.x[i]) < mma.machineEpsilon)
762 {
763 stepbeta[i] = 1.01*dx[i] / mma.machineEpsilon;
764 }
765 else
766 {
767 stepbeta[i] = 1.01*dx[i] / (beta[i] - mma.x[i]);
768 }
769 stmalfa = std::max(stepalfa[i], stmalfa);
770 stmbeta = std::max(stepbeta[i], stmbeta);
771 }
772 stmalfa_global = stmalfa;
773 stmbeta_global = stmbeta;
774#ifdef MFEM_USE_MPI
775 MPI_Allreduce(&stmalfa, &stmalfa_global, 1,
776 MPITypeMap<real_t>::mpi_type, MPI_MAX, mma.comm);
777 MPI_Allreduce(&stmbeta, &stmbeta_global, 1,
778 MPITypeMap<real_t>::mpi_type, MPI_MAX, mma.comm);
779#endif
780 stminv = std::max(std::max(std::max(stmalfa_global, stmbeta_global),
781 stmxx_global), static_cast<real_t>(1.0));
782 steg = 1.0 / stminv;
783
784 for (int i = 0; i < nvar; i++)
785 {
786 xold[i] = mma.x[i];
787 xsiold[i] = mma.xsi[i];
788 etaold[i] = mma.eta[i];
789 }
790 for (int i = 0; i < ncon; i++)
791 {
792 yold[i] = mma.y[i];
793 lamold[i] = mma.lam[i];
794 muold[i] = mma.mu[i];
795 sold[i] = mma.s[i];
796 }
797 zold = mma.z;
798 zetold = mma.zet;
799
800 itto = 0;
801 resinew = 2.0 * residunorm;
802 while (resinew > residunorm && itto < 50)
803 {
804 itto++;
805
806 for (int i = 0; i < ncon; ++i)
807 {
808 mma.y[i] = yold[i] + steg * dy[i];
809 if (std::fabs(mma.y[i])< mma.machineEpsilon)
810 {
811 mma.y[i] = mma.machineEpsilon;
812 }
813
814 mma.lam[i] = lamold[i] + steg * dlam[i];
815 if (std::fabs(mma.lam[i])< mma.machineEpsilon )
816 {
817 mma.lam[i] = mma.machineEpsilon;
818 }
819 mma.mu[i] = muold[i] + steg * dmu[i];
820 mma.s[i] = sold[i] + steg * ds[i];
821 }
822
823 residu[nvar + ncon] = mma.a0 - mma.zet; // rez
824 for (int i = 0; i < nvar; ++i)
825 {
826 mma.x[i] = xold[i] + steg * dx[i];
827 mma.xsi[i] = xsiold[i] + steg * dxsi[i];
828 mma.eta[i] = etaold[i] + steg * deta[i];
829
830 ux1[i] = mma.upp[i] - mma.x[i];
831 if (std::fabs(ux1[i]) < mma.machineEpsilon)
832 {
833 ux1[i] = mma.machineEpsilon;
834 }
835 xl1[i] = mma.x[i] - mma.low[i];
836 if (std::fabs(xl1[i]) < mma.machineEpsilon )
837 {
838 xl1[i] = mma.machineEpsilon;
839 }
840 // plam & qlam
841 plam[i] = p0[i];
842 qlam[i] = q0[i];
843 for (int j = 0; j < ncon; j++)
844 {
845 plam[i] += P[j * nvar + i] * mma.lam[j];
846 qlam[i] += Q[j * nvar + i] * mma.lam[j];
847 residu[nvar + ncon] -= mma.a[j] * mma.lam[j]; // rez
848 }
849
850 // Assembly starts here
851
852 residu[i] = plam[i] / (ux1[i] * ux1[i]) - qlam[i] / (xl1[i] * xl1[i]) -
853 mma.xsi[i] + mma.eta[i]; // rex
854 // residu[nvar + ncon] -= mma.a[i] * mma.lam[i]; // rez
855 residu[nvar + ncon + 1 + ncon + i] = mma.xsi[i] * (mma.x[i] - alfa[i]) -
856 epsi; // rexsi
857 if (std::fabs(mma.x[i] - alfa[i]) < mma.machineEpsilon)
858 {
859 residu[nvar + ncon + 1 + ncon + i] = mma.xsi[i] * mma.machineEpsilon - epsi;
860 }
861 residu[nvar + ncon + 1 + ncon + nvar + i] = mma.eta[i] *
862 (beta[i] - mma.x[i]) - epsi; // reeta
863 if (std::fabs(beta[i] - mma.x[i]) < mma.machineEpsilon)
864 {
865 residu[nvar + ncon + 1 + ncon + nvar + i] = mma.eta[i] * mma.machineEpsilon -
866 epsi;
867 }
868 }
869 mma.z = zold + steg * dz;
870 if (std::fabs(mma.z) < mma.machineEpsilon)
871 {
872 mma.z = mma.machineEpsilon;
873 }
874 mma.zet = zetold + steg * dzet;
875
876 // gvec = P/ux + Q/xl
877 for (int i = 0; i < ncon; i++)
878 {
879 gvec_local[i] = 0.0;
880 for (int j = 0; j < nvar; j++)
881 {
882 gvec_local[i] = gvec_local[i] + P[i * nvar + j] / ux1[j] + Q[i * nvar + j] /
883 xl1[j];
884 }
885 }
886 std::copy(gvec_local.get(), gvec_local.get() + ncon, gvec.get());
887
888#ifdef MFEM_USE_MPI
889 MPI_Allreduce(gvec_local.get(), gvec.get(), ncon,
890 MPITypeMap<real_t>::mpi_type, MPI_SUM, mma.comm);
891#endif
892 if (rank == 0)
893 {
894 for (int i = 0; i < ncon; i++)
895 {
896 residu[nvar + i] = mma.c[i] + mma.d[i] * mma.y[i]
897 - mma.mu[i] - mma.lam[i]; // rey
898 residu[nvar + ncon + 1 + i] = gvec[i] - mma.a[i] * mma.z
899 - mma.y[i] + mma.s[i] - b[i];
900 // relam
901 residu[nvar + ncon + 1 + ncon + 2 * nvar + i] = mma.mu[i]
902 * mma.y[i] -
903 epsi; // remu
904 residu[nvar + ncon + 1 + 2 * nvar + 2 * ncon + 1 + i] =
905 mma.lam[i] * mma.s[i] - epsi; // res
906 }
907 residu[nvar + ncon + 1 + 2 * nvar + 2 * ncon] =
908 mma.zet * mma.z - epsi; // rezet
909 }
910
911 // Get vector product and maximum absolute value
912 resinew = 0.0;
913 for (int i = 0; i < (3 * nvar + 4 * ncon + 2); i++)
914 {
915 resinew = resinew + residu[i] * residu[i];
916 }
917
918 global_norm = resinew;
919#ifdef MFEM_USE_MPI
920 MPI_Allreduce(&resinew, &global_norm, 1,
921 MPITypeMap<real_t>::mpi_type, MPI_SUM, mma.comm);
922#endif
923
924 // Norm of the residual
925 resinew = std::sqrt(global_norm);
926
927 steg = steg / 2.0;
928 }
929
930 residunorm = resinew;
931 residumax = 0.0;
932 for (int i = 0; i < (3 * nvar + 4 * ncon + 2); i++)
933 {
934 residumax = std::max(residumax, std::abs(residu[i]));
935 }
936 global_max = residumax;
937#ifdef MFEM_USE_MPI
938 MPI_Allreduce(&residumax, &global_max, 1,
939 MPITypeMap<real_t>::mpi_type, MPI_MAX, mma.comm);
940#endif
941 residumax = global_max;
942 steg = steg * 2.0;
943
944 }
945 if (ittt > 198 && mma.print_level>=2)
946 {
947 out << "Warning: Max number of iterations reached in MMA subsolve.\n";
948 }
949 epsi = 0.1 * epsi;
950 }
951
952 // returns x, y, z, lam, xsi, eta, mu, zet, s
953}
954
955void MMA::InitData(real_t *xval)
956{
957 for (int i = 0; i < nVar; i++)
958 {
959 x[i]=xval[i];
960 xo1[i] = 0.0;
961 xo2[i] = 0.0;
962 }
963
964
965 for (int i = 0; i < nCon; i++)
966 {
967 a[i] = 0.0;
968 c[i] = 1000.0;
969 d[i] = 1.0;
970 }
971 a0 = 1.0;
972}
973
974/// Serial MMA
975MMA::MMA(int nVar, int nCon, real_t *xval, int iter)
976{
977#ifdef MFEM_USE_MPI
978 comm=MPI_COMM_SELF;
979#endif
980
981 AllocData(nVar,nCon);
982 InitData(xval);
983 // allocate the serial subproblem
984#if __cplusplus >= 201402L
985 mSubProblem = ::std::make_unique<MMA::MMASubSvanberg>(*this, nVar, nCon);
986#else
987 mSubProblem.reset(new MMA::MMASubSvanberg(*this, nVar, nCon));
988#endif
989}
990
991MMA::MMA(const int nVar, int nCon, Vector &xval, int iter) : MMA(nVar, nCon,
992 xval.GetData(), iter)
993{}
994
995#ifdef MFEM_USE_MPI
996MMA::MMA(MPI_Comm comm_, int nVar, int nCon, real_t *xval, int iter)
997{
998 int rank = 0;
999 MPI_Comm_rank(comm_, &rank);
1000
1001 // create new communicator
1002 int colour;
1003
1004 if ( 0 != nVar)
1005 {
1006 colour = 0;
1007 }
1008 else
1009 {
1010 colour = MPI_UNDEFINED;
1011 }
1012
1013 // Split the global communicator
1014 MPI_Comm_split(comm_, colour, rank, &comm);
1015
1016 AllocData(nVar,nCon);
1017 InitData(xval);
1018 // allocate the serial subproblem
1019 mSubProblem.reset(new MMA::MMASubSvanberg(*this, nVar, nCon));
1020}
1021
1022MMA::MMA(MPI_Comm comm_, const int & nVar, const int & nCon,
1023 const Vector & xval, int iter) : MMA(comm_, nVar, nCon, xval.GetData(), iter)
1024{}
1025#endif
1026
1027
1028MMA::~MMA()
1029{
1030}
1031
1032void MMA::AllocData(int nVariables,int nConstr)
1033{
1034 // accessed by the subproblems
1035 nVar = nVariables;
1036 nCon = nConstr;
1037
1038 x= allocArray(nVar); // ini
1039 xo1 = allocArray(nVar); // ini
1040 xo2 = allocArray(nVar); // ini
1041
1042 y = allocArray(nCon); // ini
1043 c = allocArray(nCon); // ini
1044 d = allocArray(nCon); // ini
1045 a = allocArray(nCon); // ini
1046
1047 lam = allocArray(nCon); // ini
1048
1049 xsi = allocArray(nVar); // ini
1050 eta = allocArray(nVar); // ini
1051
1052 mu = allocArray(nCon); // ini
1053 s = allocArray(nCon); // ini
1054
1055 z = zet = 1.0;
1056 kktnorm = 10;
1057 machineEpsilon = 1e-10;
1058
1059
1060 // accessed by MMA
1061 epsimin = 1e-7;
1062 asyinit = 0.5;
1063 asyincr = 1.1;
1064 asydecr = 0.7;
1065 low = allocArray(nVar); // ini
1066 upp = allocArray(nVar); // ini
1067 factor = allocArray(nVar); // ini
1068 lowmin = lowmax = uppmin = uppmax = zz = 0.0;
1069
1070}
1071
1072void MMA::Update( const Vector& dfdx,
1073 const Vector& gx, const Vector& dgdx,
1074 const Vector& xmin, const Vector& xmax,
1075 Vector& xval)
1076{
1077 this->Update(dfdx.GetData(),
1078 gx.GetData(),dgdx.GetData(),
1079 xmin.GetData(), xmax.GetData(),
1080 xval.GetData());
1081}
1082
1083void MMA::Update( const Vector& dfdx,
1084 const Vector& xmin, const Vector& xmax,
1085 Vector& xval)
1086{
1087 MFEM_ASSERT(0 == nCon,
1088 "MMA nCon != 0. Provide constraint values and gradients");
1089
1090 this->Update(dfdx.GetData(),
1091 nullptr,nullptr,
1092 xmin.GetData(), xmax.GetData(),
1093 xval.GetData());
1094}
1095
1096void MMA::Update(const real_t* dfdx,
1097 const real_t* gx,const real_t* dgdx,
1098 const real_t* xmin, const real_t* xmax,
1099 real_t* xval)
1100{
1101 // Calculation of the asymptotes low and upp
1102 if (iter < 3)
1103 {
1104 for (int i = 0; i < nVar; i++)
1105 {
1106 low[i] = xval[i] - asyinit * (xmax[i] - xmin[i]);
1107 upp[i] = xval[i] + asyinit * (xmax[i] - xmin[i]);
1108 }
1109 }
1110 else
1111 {
1112 for (int i = 0; i < nVar; i++)
1113 {
1114 // Determine sign
1115 zz = (xval[i] - xo1[i]) * (xo1[i] - xo2[i]);
1116 if ( zz > 0.0)
1117 {
1118 factor[i] = asyincr;
1119 }
1120 else if ( zz < 0.0)
1121 {
1122 factor[i] = asydecr;
1123 }
1124 else
1125 {
1126 factor[i] = 1.0;
1127 }
1128
1129
1130 // Find new asymptote
1131 low[i] = xval[i] - factor[i] * (xo1[i] - low[i]);
1132 upp[i] = xval[i] + factor[i] * (upp[i] - xo1[i]);
1133
1134 lowmin = xval[i] - 10.0 * (xmax[i] - xmin[i]);
1135 lowmax = xval[i] - 0.01 * (xmax[i] - xmin[i]);
1136 uppmin = xval[i] + 0.01 * (xmax[i] - xmin[i]);
1137 uppmax = xval[i] + 10.0 * (xmax[i] - xmin[i]);
1138
1139 low[i] = std::max(low[i], lowmin);
1140 low[i] = std::min(low[i], lowmax);
1141 upp[i] = std::max(upp[i], uppmin);
1142 upp[i] = std::min(upp[i], uppmax);
1143 }
1144 }
1145
1146 mSubProblem->Update(dfdx,gx,dgdx,xmin,xmax,xval);
1147 // Update design variables
1148 for (int i = 0; i < nVar; i++)
1149 {
1150 xo2[i] = xo1[i];
1151 xo1[i] = xval[i];
1152 xval[i] = x[i];
1153 }
1154
1155 iter++;
1156}
1157
1158}
mfem::MMA
MMA (Method of Moving Asymptotes) solves an optimization problem of the form:
Definition mma.hpp:53
mfem::MMA::mu
::std::unique_ptr< real_t[]> mu
Definition mma.hpp:112
mfem::MMA::b
::std::unique_ptr< real_t[]> b
Definition mma.hpp:99
mfem::MMA::eta
::std::unique_ptr< real_t[]> eta
Definition mma.hpp:112
mfem::MMA::d
::std::unique_ptr< real_t[]> d
Definition mma.hpp:99
mfem::MMA::a0
real_t a0
Definition mma.hpp:100
mfem::MMA::~MMA
~MMA()
Destructor.
Definition mma.cpp:1028
mfem::MMA::machineEpsilon
real_t machineEpsilon
Definition mma.hpp:100
mfem::MMA::zet
real_t zet
Definition mma.hpp:101
mfem::MMA::c
::std::unique_ptr< real_t[]> c
Definition mma.hpp:99
mfem::MMA::Update
void Update(const Vector &dfdx, const Vector &gx, const Vector &dgdx, const Vector &xmin, const Vector &xmax, Vector &xval)
Definition mma.cpp:1072
mfem::MMA::MMA
MMA(int nVar, int nCon, real_t *xval, int iterationNumber=0)
Serial MMA.
Definition mma.cpp:975
mfem::MMA::upp
::std::unique_ptr< real_t[]> upp
Definition mma.hpp:111
mfem::MMA::a
::std::unique_ptr< real_t[]> a
Definition mma.hpp:99
mfem::MMA::x
::std::unique_ptr< real_t[]> x
Definition mma.hpp:112
mfem::MMA::xsi
::std::unique_ptr< real_t[]> xsi
Definition mma.hpp:112
mfem::MMA::lam
::std::unique_ptr< real_t[]> lam
Definition mma.hpp:112
mfem::MMA::low
::std::unique_ptr< real_t[]> low
Definition mma.hpp:111
mfem::MMA::z
real_t z
Definition mma.hpp:101
mfem::MMA::y
::std::unique_ptr< real_t[]> y
Definition mma.hpp:112
mfem::MMA::nCon
int nCon
Definition mma.hpp:102
mfem::MMA::s
::std::unique_ptr< real_t[]> s
Definition mma.hpp:112
mfem::MMA::nVar
int nVar
Definition mma.hpp:102
mfem::MMA::iter
int iter
Definition mma.hpp:105
mfem::MMA::epsimin
real_t epsimin
Definition mma.hpp:100
mfem::Vector
Vector data type.
Definition vector.hpp:82
mfem::Vector::GetData
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
Definition vector.hpp:235
beta
Vector beta
Definition convection-diffusion.cpp:82
sgesv_
void sgesv_(int *nLAP, int *nrhs, float *AA, int *lda, int *ipiv, float *bb, int *ldb, int *info)
dgesv_
void dgesv_(int *nLAP, int *nrhs, double *AA, int *lda, int *ipiv, double *bb, int *ldb, int *info)
mfem::out
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition globals.hpp:66
mfem::real_t
float real_t
Definition config.hpp:43
mfem::solveLU
void solveLU(int nCon, real_t *AA1, real_t *bb1)
Definition mma.cpp:44
Generated on Wed Apr 9 2025 16:53:05 for MFEM by doxygen 1.11.0