61 if (dot_prod_type == 0)
75 int print_level_ = print_lvl;
78 if (dot_prod_type != 0)
81 MPI_Comm_rank(comm, &rank);
100 if (dot_prod_type != 0)
103 MPI_Comm_rank(comm, &rank);
106 derived_print_level = -1;
120 if (comm != MPI_COMM_NULL)
122 MPI_Comm_rank(comm, &rank);
126 switch (print_level_)
143 MFEM_WARNING(
"Unknown print level " << print_level_ <<
144 ". Defaulting to level 0.");
160 else if (print_options_.
summary)
192 const Vector& x,
bool final)
const
196 if (it == 0 && !
final)
209 ess_tdof_list(nullptr),
218 Solver(
a.FESpace()->GetTrueVSize()),
221 ess_tdof_list(&ess_tdofs),
226 a.AssembleDiagonal(diag);
240 ess_tdof_list(&ess_tdofs),
273 MFEM_VERIFY(
height ==
width,
"not a square matrix!");
275 ess_tdof_list =
nullptr;
288 auto D = diag.
Read();
289 auto DI = dinv.
Write();
290 const bool use_abs_diag_ = use_abs_diag;
295 MFEM_ABORT_KERNEL(
"Zero diagonal entry in OperatorJacobiSmoother");
297 if (!use_abs_diag_) { DI[i] =
delta / D[i]; }
298 else { DI[i] =
delta / std::abs(D[i]); }
300 if (ess_tdof_list && ess_tdof_list->
Size() > 0)
302 auto I = ess_tdof_list->
Read();
314 MFEM_VERIFY(x.
Size() ==
Width(),
"invalid input vector");
315 MFEM_VERIFY(y.
Size() ==
Height(),
"invalid output vector");
319 MFEM_VERIFY(oper,
"iterative_mode == true requires the forward operator");
320 oper->
Mult(y, residual);
329 auto DI = dinv.
Read();
330 auto R = residual.
Read();
334 Y[i] += DI[i] * R[i];
341 int order_,
real_t max_eig_estimate_)
345 max_eig_estimate(max_eig_estimate_),
350 ess_tdof_list(ess_tdofs),
352 oper(&oper_) {
Setup(); }
358 int order_, MPI_Comm comm,
359 int power_iterations,
367 int power_iterations,
377 ess_tdof_list(ess_tdofs),
390 MFEM_VERIFY(power_seed != 0,
"invalid seed!");
402 int order_,
real_t max_eig_estimate_)
409 int order_, MPI_Comm comm,
int power_iterations,
real_t power_tolerance)
411 power_iterations, power_tolerance) { }
416 int order_,
int power_iterations,
real_t power_tolerance)
426 auto D = diag.
Read();
427 auto X = dinv.
Write();
428 mfem::forall(N, [=] MFEM_HOST_DEVICE (
int i) { X[i] = 1.0 / D[i]; });
429 auto I = ess_tdof_list.
Read();
438 real_t upper_bound = 1.2 * max_eig_estimate;
439 real_t lower_bound = 0.3 * max_eig_estimate;
440 real_t theta = 0.5 * (upper_bound + lower_bound);
447 coeffs[0] = 1.0 / theta;
452 real_t tmp_0 = 1.0/(pow(
delta, 2) - 2*pow(theta, 2));
453 coeffs[0] = -4*theta*tmp_0;
460 real_t tmp_1 = pow(theta, 2);
461 real_t tmp_2 = 1.0/(-4*pow(theta, 3) + theta*tmp_0);
462 coeffs[0] = tmp_2*(tmp_0 - 12*tmp_1);
463 coeffs[1] = 12/(tmp_0 - 4*tmp_1);
464 coeffs[2] = -4*tmp_2;
470 real_t tmp_1 = pow(theta, 2);
472 real_t tmp_3 = 1.0/(pow(
delta, 4) + 8*pow(theta, 4) - tmp_1*tmp_2);
473 coeffs[0] = tmp_3*(32*pow(theta, 3) - 16*theta*tmp_0);
474 coeffs[1] = tmp_3*(-48*tmp_1 + tmp_2);
475 coeffs[2] = 32*theta*tmp_3;
476 coeffs[3] = -8*tmp_3;
482 real_t tmp_1 = pow(theta, 4);
483 real_t tmp_2 = pow(theta, 2);
487 real_t tmp_6 = 1.0/(16*pow(theta, 5) - pow(theta, 3)*tmp_5 + theta*tmp_0);
489 real_t tmp_8 = 1.0/(tmp_0 + 16*tmp_1 - tmp_2*tmp_5);
490 coeffs[0] = tmp_6*(tmp_0 + 80*tmp_1 - tmp_2*tmp_4);
491 coeffs[1] = tmp_8*(tmp_4 - tmp_7);
492 coeffs[2] = tmp_6*(-tmp_5 + tmp_7);
493 coeffs[3] = -80*tmp_8;
494 coeffs[4] = 16*tmp_6;
498 MFEM_ABORT(
"Chebyshev smoother not implemented for order = " << order);
506 MFEM_ABORT(
"Chebyshev smoother not implemented for iterative mode");
511 MFEM_ABORT(
"Chebyshev smoother requires operator");
521 for (
int k = 0; k < order; ++k)
526 oper->
Mult(residual, helperVector);
527 residual = helperVector;
532 auto Dinv = dinv.
Read();
534 mfem::forall(n, [=] MFEM_HOST_DEVICE (
int i) { R[i] *= Dinv[i]; });
538 auto C = coeffs.
Read();
539 mfem::forall(n, [=] MFEM_HOST_DEVICE (
int i) { Y[i] += C[k] * R[i]; });
554 const bool zero_b = (
b.Size() == 0);
558 "empty 'b' can be used only in iterative mode!");
569 if (!zero_b) { x +=
z; }
596 real_t r0, nom, nom0, nomold = 1, cf;
612 nom0 = nom = sqrt(
Dot(
z,
z));
616 nom0 = nom = sqrt(
Dot(
r,
r));
622 mfem::out <<
" Iteration : " << setw(3) << right << 0 <<
" ||Br|| = "
642 if (!zero_b) { x +=
z; }
647 if (!zero_b) { x +=
r; }
657 nom = sqrt(
Dot(
z,
z));
661 nom = sqrt(
Dot(
r,
r));
682 mfem::out <<
" Iteration : " << setw(3) << right << (i-1)
683 <<
" ||Br|| = " << setw(11) << left << nom
684 <<
"\tConv. rate: " << cf <<
'\n';
692 const auto rf = pow (nom/nom0, 1.0/
final_iter);
694 <<
"Conv. rate: " << cf <<
'\n'
695 <<
"Average reduction factor: "<< rf <<
'\n';
699 mfem::out <<
"SLI: No convergence!" <<
'\n';
706 int print_iter,
int max_num_iter,
721 int print_iter,
int max_num_iter,
777 nom0 = nom =
Dot(
d,
r);
779 MFEM_VERIFY(
IsFinite(nom),
"nom = " << nom);
782 mfem::out <<
" Iteration : " << setw(3) << 0 <<
" (B r, r) = "
790 mfem::out <<
"PCG: The preconditioner is not positive definite. (Br, r) = "
802 if (
Monitor(0, nom,
r, x) || nom <= r0)
814 MFEM_VERIFY(
IsFinite(den),
"den = " << den);
819 mfem::out <<
"PCG: The operator is not positive definite. (Ad, d) = "
851 MFEM_VERIFY(
IsFinite(betanom),
"betanom = " << betanom);
856 mfem::out <<
"PCG: The preconditioner is not positive definite. (Br, r) = "
866 mfem::out <<
" Iteration : " << setw(3) << i <<
" (B r, r) = "
867 << betanom << std::endl;
870 if (
Monitor(i, betanom,
r, x) || betanom <= r0)
893 MFEM_VERIFY(
IsFinite(den),
"den = " << den);
898 mfem::out <<
"PCG: The operator is not positive definite. (Ad, d) = "
921 const auto arf = pow (betanom/nom0, 0.5/
final_iter);
922 mfem::out <<
"Average reduction factor = " << arf <<
'\n';
926 mfem::out <<
"PCG: No convergence!" <<
'\n';
935 int print_iter,
int max_num_iter,
950 int print_iter,
int max_num_iter,
974 else if (fabs(dy) > fabs(dx))
977 sn = 1.0 / sqrt( 1.0 + temp*temp );
983 cs = 1.0 / sqrt( 1.0 + temp*temp );
990 real_t temp = cs * dx + sn * dy;
991 dy = -sn * dx + cs * dy;
1001 for (
int i = k; i >= 0; i--)
1004 for (
int j = i - 1; j >= 0; j--)
1006 y(j) -= h(j,i) * y(i);
1010 for (
int j = 0; j <= k; j++)
1084 <<
" Iteration : " << setw(3) << 0
1085 <<
" ||B r|| = " <<
beta
1098 v[0]->Set(1.0/
beta, r);
1099 s = 0.0; s(0) =
beta;
1101 for (i = 0; i <
m && j <=
max_iter; i++, j++)
1113 for (k = 0; k <= i; k++)
1115 H(k,i) =
Dot(w, *v[k]);
1116 w.
Add(-H(k,i), *v[k]);
1120 MFEM_VERIFY(
IsFinite(H(i+1,i)),
"Norm(w) = " << H(i+1,i));
1121 if (v[i+1] == NULL) { v[i+1] =
new Vector(n); }
1122 v[i+1]->Set(1.0/H(i+1,i), w);
1124 for (k = 0; k < i; k++)
1133 const real_t resid = fabs(s(i+1));
1134 MFEM_VERIFY(
IsFinite(resid),
"resid = " << resid);
1147 mfem::out <<
" Pass : " << setw(2) << (j-1)/
m+1
1148 <<
" Iteration : " << setw(3) << j
1149 <<
" ||B r|| = " << resid <<
'\n';
1188 mfem::out <<
" Pass : " << setw(2) << (j-1)/
m+1
1198 mfem::out <<
"GMRES: No convergence!\n";
1203 for (i = 0; i < v.
Size(); i++)
1252 <<
" Iteration : " << setw(3) << 0
1253 <<
" || r || = " <<
beta
1259 for (i= 0; i<=
m; i++)
1275 s = 0.0; s(0) =
beta;
1277 for (i = 0; i <
m && j <=
max_iter; i++, j++)
1297 for (k = 0; k <= i; k++)
1299 H(k,i) =
Dot( r, *v[k]);
1300 r.
Add(-H(k,i), (*v[k]));
1306 v[i+1] =
new Vector(
b.Size());
1310 v[i+1] ->
Add (1.0/H(i+1,i), r);
1312 for (k = 0; k < i; k++)
1321 const real_t resid = fabs(s(i+1));
1322 MFEM_VERIFY(
IsFinite(resid),
"resid = " << resid);
1326 mfem::out <<
" Pass : " << setw(2) << (j-1)/
m+1
1327 <<
" Iteration : " << setw(3) << j
1328 <<
" || r || = " << resid << endl;
1343 for (i= 0; i<=
m; i++)
1345 if (v[i]) {
delete v[i]; }
1346 if (z[i]) {
delete z[i]; }
1349 Monitor(j, resid, r, x,
true);
1374 for (i = 0; i <=
m; i++)
1376 if (v[i]) {
delete v[i]; }
1377 if (z[i]) {
delete z[i]; }
1386 mfem::out <<
" Pass : " << setw(2) << (j-1)/
m+1
1387 <<
" Iteration : " << setw(3) << j-1
1396 mfem::out <<
"FGMRES: No convergence!\n";
1404 int &max_iter,
int m,
real_t &tol,
real_t atol,
int printit)
1423 int print_iter,
int max_num_iter,
int m,
real_t rtol,
real_t atol)
1425 GMRES(A, x,
b, B, max_num_iter, m, rtol, atol, print_iter);
1481 MFEM_VERIFY(
IsFinite(resid),
"resid = " << resid);
1484 mfem::out <<
" Iteration : " << setw(3) << 0
1490 if (
Monitor(0, resid,
r, x) || resid <= tol_goal)
1507 mfem::out <<
" Iteration : " << setw(3) << i
1508 <<
" ||r|| = " << resid <<
'\n';
1520 mfem::out <<
"BiCGStab: No convergence!\n";
1548 MFEM_VERIFY(
IsFinite(resid),
"resid = " << resid);
1549 if (
Monitor(i, resid,
r, x) || resid < tol_goal)
1554 mfem::out <<
" Iteration : " << setw(3) << i
1555 <<
" ||s|| = " << resid <<
'\n';
1570 mfem::out <<
" Iteration : " << setw(3) << i
1571 <<
" ||s|| = " << resid;
1589 MFEM_VERIFY(
IsFinite(resid),
"resid = " << resid);
1592 mfem::out <<
" ||r|| = " << resid <<
'\n';
1594 if (
Monitor(i, resid,
r, x) || resid < tol_goal)
1601 mfem::out <<
" Iteration : " << setw(3) << i
1602 <<
" ||r|| = " << resid <<
'\n';
1619 mfem::out <<
" Iteration : " << setw(3) << i
1620 <<
" ||r|| = " << resid <<
'\n';
1628 mfem::out <<
"BiCGStab: No convergence!\n";
1643 <<
" ||r|| = " << resid <<
'\n';
1651 mfem::out <<
"BiCGStab: No convergence!\n";
1667 bicgstab.
Mult(
b, x);
1674 int print_iter,
int max_num_iter,
real_t rtol,
real_t atol)
1676 BiCGSTAB(A, x,
b, B, max_num_iter, rtol, atol, print_iter);
1716 real_t beta, eta, gamma0, gamma1, sigma0, sigma1;
1738 MFEM_VERIFY(
IsFinite(eta),
"eta = " << eta);
1739 gamma0 = gamma1 = 1.;
1740 sigma0 = sigma1 = 0.;
1744 if (
Monitor(0, eta, *z, x) || eta <= norm_goal)
1752 mfem::out <<
"MINRES: iteration " << setw(3) << 0 <<
": ||r||_B = "
1774 rho2 = sigma1*
alpha + gamma0*gamma1*
beta;
1789 w0.
Set(1./rho1, *z);
1793 add(1./rho1, *z, -rho2/rho1,
w1,
w0);
1797 add(-rho3/rho1,
w0, -rho2/rho1,
w1,
w0);
1798 w0.
Add(1./rho1, *z);
1802 gamma1 =
delta/rho1;
1804 x.
Add(gamma1*eta,
w0);
1810 MFEM_VERIFY(
IsFinite(eta),
"eta = " << eta);
1812 if (fabs(eta) <= norm_goal)
1819 mfem::out <<
"MINRES: iteration " << setw(3) << it <<
": ||r||_B = "
1820 << fabs(eta) <<
'\n';
1822 if (
Monitor(it, fabs(eta), *z, x)) {
goto loop_end; }
1840 mfem::out <<
"MINRES: iteration " << setw(3) << it <<
": ||r||_B = "
1841 << fabs(eta) <<
'\n';
1866 mfem::out <<
"MINRES: No convergence!\n";
1903 MFEM_VERIFY(
height ==
width,
"square Operator is required.");
1914 MFEM_VERIFY(
oper != NULL,
"the Operator is not set (use SetOperator).");
1915 MFEM_VERIFY(
prec != NULL,
"the Solver is not set (use SetSolver).");
1918 real_t norm0, norm, norm_goal;
1919 const bool have_b = (
b.Size() ==
Height());
1937 mfem::out <<
"Newton iteration " << setw(2) << 0
1938 <<
" : ||r|| = " << norm <<
"...\n";
1945 for (it = 0;
true; it++)
1947 MFEM_VERIFY(
IsFinite(norm),
"norm = " << norm);
1950 mfem::out <<
"Newton iteration " << setw(2) << it
1951 <<
" : ||r|| = " << norm;
1954 mfem::out <<
", ||r||/||r_0|| = " << norm/norm0;
1959 if (
Monitor(it, norm,
r, x) || norm <= norm_goal)
1992 add(x, -c_scale,
c, x);
2012 <<
", ||r||/||r_0|| = " <<
final_norm/norm0 <<
'\n';
2016 mfem::out <<
"Newton: No convergence!\n";
2031 this->
alpha = alpha_;
2032 this->
gamma = gamma_;
2037 const real_t fnorm)
const
2045 real_t sg_threshold = 0.1;
2065 MFEM_ABORT(
"Unknown adaptive linear solver rtol version");
2070 if (sg_eta > sg_threshold) { eta = std::max(eta, sg_eta); }
2074 iterative_solver->SetRelTol(eta);
2078 mfem::out <<
"Eisenstat-Walker rtol = " << eta <<
"\n";
2085 const real_t fnorm)
const
2105 MFEM_VERIFY(
oper != NULL,
"the Operator is not set (use SetOperator).");
2127 alpha.UseDevice(
true);
2129 int last_saved_id = -1;
2132 real_t norm0, norm, norm_goal;
2133 const bool have_b = (
b.Size() ==
Height());
2144 if (have_b) {
r -=
b; }
2151 mfem::out <<
"LBFGS iteration " << setw(2) << 0
2152 <<
" : ||r|| = " << norm <<
"...\n";
2155 for (it = 0;
true; it++)
2157 MFEM_VERIFY(
IsFinite(norm),
"norm = " << norm);
2161 <<
" : ||r|| = " << norm;
2164 mfem::out <<
", ||r||/||r_0|| = " << norm/norm0;
2169 if (norm <= norm_goal)
2188 add(x, -c_scale,
c, x);
2200 sk =
c; sk *= -c_scale;
2204 last_saved_id = (last_saved_id ==
m-1) ? 0 : last_saved_id+1;
2209 for (
int i = last_saved_id; i > -1; i--)
2217 for (
int i =
m-1; i > last_saved_id; i--)
2228 for (
int i = last_saved_id+1; i <
m ; i++)
2234 for (
int i = 0; i < last_saved_id+1 ; i++)
2251 <<
", ||r||/||r_0|| = " <<
final_norm/norm0 <<
'\n';
2255 mfem::out <<
"LBFGS: No convergence!\n";
2261 int m_max,
int m_min,
int m_step,
real_t cf,
2269 Vector s(m+1), cs(m+1), sn(m+1);
2281 real_t normb = w.Norml2();
2294 resid =
beta / normb;
2296 if (resid * resid <= tol)
2298 tol = resid * resid;
2306 <<
" Iteration : " << setw(3) << 0
2307 <<
" (r, r) = " <<
beta*
beta <<
'\n';
2310 tol *= (normb*normb);
2311 tol = (atol > tol) ? atol : tol;
2315 for (i= 0; i<=m; i++)
2323 while (j <= max_iter)
2327 s = 0.0; s(0) =
beta;
2331 for (i = 0; i < m && j <= max_iter; i++)
2336 for (k = 0; k <= i; k++)
2338 H(k,i) = w * (*v[k]);
2339 w.
Add(-H(k,i), (*v[k]));
2342 H(i+1,i) = w.Norml2();
2344 v[i+1] ->
Add (1.0/H(i+1,i), w);
2346 for (k = 0; k < i; k++)
2355 resid = fabs(s(i+1));
2359 <<
" Iteration : " << setw(3) << i+1
2360 <<
" (r, r) = " << resid*resid <<
'\n';
2363 if ( resid*resid < tol)
2366 tol = resid * resid;
2368 for (i= 0; i<=m; i++)
2387 if ( resid*resid < tol)
2389 tol = resid * resid;
2391 for (i= 0; i<=m; i++)
2400 if (m - m_step >= m_min)
2413 tol = resid * resid;
2414 for (i= 0; i<=m; i++)
2424 : C(C_), D(D_), c_e(NULL), d_lo(NULL), d_hi(NULL), x_lo(NULL), x_hi(NULL),
2433 MFEM_VERIFY(
C,
"The C operator is unspecified -- can't set constraints.");
2434 MFEM_VERIFY(c.
Size() ==
C->
Height(),
"Wrong size of the constraint.");
2442 MFEM_VERIFY(
D,
"The D operator is unspecified -- can't set constraints.");
2444 "Wrong size of the constraint.");
2452 "Wrong size of the constraint.");
2469 MFEM_WARNING(
"Objective functional is ignored as SLBQP always minimizes"
2470 "the l2 norm of (x - x_target).");
2472 MFEM_VERIFY(
prob.GetC(),
"Linear constraint is not set.");
2473 MFEM_VERIFY(
prob.GetC()->Height() == 1,
"Solver expects scalar constraint.");
2494 mfem::out <<
"SLBQP iteration " << it <<
": residual = " << r
2495 <<
", lambda = " << l <<
'\n';
2525 mfem::out <<
"SLBQP bracketing phase" <<
'\n';
2543 llow = l; rlow = r; l = l + dl;
2546 r =
solve(l,xt,x,nclip);
2549 while ((r < 0) && (nclip <
max_iter))
2553 if (s < smin) { s = smin; }
2558 r =
solve(l,xt,x,nclip);
2566 lupp = l; rupp = r; l = l - dl;
2569 r =
solve(l,xt,x,nclip);
2572 while ((r > 0) && (nclip <
max_iter))
2576 if (s < smin) { s = smin; }
2581 r =
solve(l,xt,x,nclip);
2594 mfem::out <<
"SLBQP secant phase" <<
'\n';
2597 s = 1.0 - rlow/rupp; dl = dl/s; l = lupp - dl;
2600 r =
solve(l,xt,x,nclip);
2603 while ( (fabs(r) > tol) && (nclip <
max_iter) )
2609 lupp = l; rupp = r; s = 1.0 - rlow/rupp;
2610 dl = (lupp - llow)/s; l = lupp - dl;
2615 if (s < smin) { s = smin; }
2617 lnew = 0.75*llow + 0.25*l;
2618 if (lnew < l-dl) { lnew = l-dl; }
2619 lupp = l; rupp = r; l = lnew;
2620 s = (lupp - llow)/(lupp - l);
2628 llow = l; rlow = r; s = 1.0 - rlow/rupp;
2629 dl = (lupp - llow)/s; l = lupp - dl;
2634 if (s < smin) { s = smin; }
2636 lnew = 0.75*lupp + 0.25*l;
2637 if (lnew < l+dl) { lnew = l+dl; }
2638 llow = l; rlow = r; l = lnew;
2639 s = (lupp - llow)/(lupp - l);
2644 r =
solve(l,xt,x,nclip);
2660 <<
" lambda = " << l <<
'\n'
2665 mfem::out <<
"SLBQP: No convergence!" <<
'\n';
2671 const std::vector<real_t> &w;
2672 std::vector<size_t> c;
2673 std::vector<int> loc;
2675 WeightMinHeap(
const std::vector<real_t> &w_) : w(w_)
2677 c.reserve(w.size());
2678 loc.resize(w.size());
2679 for (
size_t i=0; i<w.size(); ++i) { push(i); }
2682 size_t percolate_up(
size_t pos,
real_t val)
2684 for (; pos > 0 && w[c[(pos-1)/2]] > val; pos = (pos-1)/2)
2686 c[pos] = c[(pos-1)/2];
2687 loc[c[(pos-1)/2]] =
static_cast<int>(pos);
2692 size_t percolate_down(
size_t pos,
real_t val)
2694 while (2*pos+1 < c.size())
2696 size_t left = 2*pos+1;
2697 size_t right = left+1;
2699 if (right < c.size() && w[c[right]] < w[c[left]]) { tgt = right; }
2700 else { tgt = left; }
2701 if (w[c[tgt]] < val)
2704 loc[c[tgt]] =
static_cast<int>(pos);
2719 size_t pos = c.size()-1;
2720 pos = percolate_up(pos, val);
2722 loc[i] =
static_cast<int>(pos);
2728 size_t j = c.back();
2732 if (c.empty()) {
return static_cast<int>(i); }
2735 pos = percolate_down(pos, val);
2737 loc[j] =
static_cast<int>(pos);
2738 return static_cast<int>(i);
2741 void update(
size_t i)
2743 size_t pos = loc[i];
2745 pos = percolate_up(pos, val);
2746 pos = percolate_down(pos, val);
2748 loc[i] =
static_cast<int>(pos);
2751 bool picked(
size_t i)
2767 for (
int i=0; i<n; ++i)
2769 for (
int j=I[i]; j<I[i+1]; ++j)
2771 V[j] = abs(V[j]/D[i]);
2775 std::vector<real_t> w(n, 0.0);
2776 for (
int k=0; k<n; ++k)
2779 for (
int ii=I[k]; ii<I[k+1]; ++ii)
2784 for (
int kk=I[i]; kk<I[i+1]; ++kk)
2792 for (
int jj=I[k]; jj<I[k+1]; ++jj)
2795 if (j == k) {
continue; }
2797 bool ij_exists =
false;
2798 for (
int jj2=I[i]; jj2<I[i+1]; ++jj2)
2806 if (!ij_exists) { w[k] += pow(C_ik*C_kj,2); }
2812 WeightMinHeap w_heap(w);
2816 for (
int ii=0; ii<n; ++ii)
2818 int pi = w_heap.pop();
2821 for (
int kk=I[pi]; kk<I[pi+1]; ++kk)
2824 if (w_heap.picked(k)) {
continue; }
2828 for (
int ii2=I[k]; ii2<I[k+1]; ++ii2)
2831 if (w_heap.picked(i)) {
continue; }
2834 for (
int kk2=I[i]; kk2<I[i+1]; ++kk2)
2842 for (
int jj=I[k]; jj<I[k+1]; ++jj)
2845 if (j == k || w_heap.picked(j)) {
continue; }
2847 bool ij_exists =
false;
2848 for (
int jj2=I[i]; jj2<I[i+1]; ++jj2)
2856 if (!ij_exists) { w[k] += pow(C_ik*C_kj,2); }
2869 block_size(block_size_),
2871 reordering(reordering_)
2878 :
BlockILU(block_size_, reordering_, k_fill_)
2900 MFEM_ABORT(
"BlockILU must be created with a SparseMatrix or HypreParMatrix");
2905 MFEM_VERIFY(A->
Finalized(),
"Matrix must be finalized.");
2906 CreateBlockPattern(*A);
2910void BlockILU::CreateBlockPattern(
const SparseMatrix &A)
2912 MFEM_VERIFY(k_fill == 0,
"Only block ILU(0) is currently supported.");
2913 if (A.
Height() % block_size != 0)
2915 MFEM_ABORT(
"BlockILU: block size must evenly divide the matrix size");
2923 int nblockrows = nrows / block_size;
2925 std::vector<std::set<int>> unique_block_cols(nblockrows);
2927 for (
int iblock = 0; iblock < nblockrows; ++iblock)
2929 for (
int bi = 0; bi < block_size; ++bi)
2931 int i = iblock * block_size + bi;
2932 for (
int k = I[i]; k < I[i + 1]; ++k)
2934 unique_block_cols[iblock].insert(J[k] / block_size);
2937 nnz +=
static_cast<int>(unique_block_cols[iblock].size());
2942 SparseMatrix C(nblockrows, nblockrows);
2943 for (
int iblock = 0; iblock < nblockrows; ++iblock)
2945 for (
int jblock : unique_block_cols[iblock])
2947 for (
int bi = 0; bi < block_size; ++bi)
2949 int i = iblock * block_size + bi;
2950 for (
int k = I[i]; k < I[i + 1]; ++k)
2953 if (j >= jblock * block_size && j < (jblock + 1) * block_size)
2955 C.Add(iblock, jblock, V[k]*V[k]);
2962 real_t *CV = C.GetData();
2963 for (
int i=0; i<C.NumNonZeroElems(); ++i)
2965 CV[i] = sqrt(CV[i]);
2974 MFEM_ABORT(
"BlockILU: unknown reordering")
2981 for (
int i=0; i<nblockrows; ++i)
2989 for (
int i=0; i<nblockrows; ++i)
2995 std::vector<std::vector<int>> unique_block_cols_perminv(nblockrows);
2996 for (
int i=0; i<nblockrows; ++i)
2998 std::vector<int> &cols = unique_block_cols_perminv[i];
2999 for (
int j : unique_block_cols[P[i]])
3001 cols.push_back(Pinv[j]);
3003 std::sort(cols.begin(), cols.end());
3010 AB.
SetSize(block_size, block_size, nnz);
3011 DB.
SetSize(block_size, block_size, nblockrows);
3014 ipiv.
SetSize(block_size*nblockrows);
3017 for (
int iblock = 0; iblock < nblockrows; ++iblock)
3019 int iblock_perm = P[iblock];
3020 for (
int jblock : unique_block_cols_perminv[iblock])
3022 int jblock_perm = P[jblock];
3023 if (iblock == jblock)
3025 ID[iblock] = counter;
3027 JB[counter] = jblock;
3028 for (
int bi = 0; bi < block_size; ++bi)
3030 int i = iblock_perm*block_size + bi;
3031 for (
int k = I[i]; k < I[i + 1]; ++k)
3034 if (j >= jblock_perm*block_size && j < (jblock_perm + 1)*block_size)
3036 int bj = j - jblock_perm*block_size;
3038 AB(bi, bj, counter) = val;
3040 if (iblock == jblock)
3042 DB(bi, bj, iblock) = val;
3049 IB[iblock + 1] = counter;
3053void BlockILU::Factorize()
3055 int nblockrows =
Height()/block_size;
3058 for (
int i=0; i<nblockrows; ++i)
3060 LUFactors factorization(DB.
GetData(i), &ipiv[i*block_size]);
3061 factorization.Factor(block_size);
3067 DenseMatrix A_ik, A_ij, A_kj;
3069 for (
int i=1; i<nblockrows; ++i)
3072 for (
int kk=IB[i]; kk<IB[i+1]; ++kk)
3076 if (k == i) {
break; }
3079 MFEM_ABORT(
"Matrix must be sorted with nonzero diagonal");
3081 LUFactors A_kk_inv(DB.
GetData(k), &ipiv[k*block_size]);
3082 A_ik.UseExternalData(&AB(0,0,kk), block_size, block_size);
3084 A_kk_inv.RightSolve(block_size, block_size, A_ik.GetData());
3086 for (
int jj=kk+1; jj<IB[i+1]; ++jj)
3089 if (j <= k) {
continue; }
3090 A_ij.UseExternalData(&AB(0,0,jj), block_size, block_size);
3091 for (
int ll=IB[k]; ll<IB[k+1]; ++ll)
3096 A_kj.UseExternalData(&AB(0,0,ll), block_size, block_size);
3103 LUFactors factorization(DB.
GetData(i), &ipiv[i*block_size]);
3104 factorization.Factor(block_size);
3116 MFEM_VERIFY(
height > 0,
"BlockILU(0) preconditioner is not constructed");
3117 int nblockrows =
Height()/block_size;
3126 for (
int i=0; i<nblockrows; ++i)
3129 for (
int ib=0; ib<block_size; ++ib)
3131 yi[ib] =
b[ib + P[i]*block_size];
3133 for (
int k=IB[i]; k<ID[i]; ++k)
3143 for (
int i=nblockrows-1; i >= 0; --i)
3146 for (
int ib=0; ib<block_size; ++ib)
3148 xi[ib] = y[ib + i*block_size];
3150 for (
int k=ID[i]+1; k<IB[i+1]; ++k)
3158 LUFactors A_ii_inv(&DB(0,0,i), &ipiv[i*block_size]);
3170 real_t bc_norm_squared = 0.0;
3176 bc_norm_squared += r_entry*r_entry;
3181 if (comm != MPI_COMM_NULL)
3183 double glob_bc_norm_squared = 0.0;
3184 MPI_Reduce(&bc_norm_squared, &glob_bc_norm_squared, 1,
3187 bc_norm_squared = glob_bc_norm_squared;
3189 MPI_Comm_rank(comm, &rank);
3190 print = (rank == 0);
3193 if ((it == 0 ||
final || bc_norm_squared > 0.0) && print)
3195 mfem::out <<
" ResidualBCMonitor : b.c. residual norm = "
3196 << sqrt(bc_norm_squared) << endl;
3201#ifdef MFEM_USE_SUITESPARSE
3226 umfpack_di_free_numeric(&
Numeric);
3230 umfpack_dl_free_numeric(&
Numeric);
3235 MFEM_VERIFY(
mat,
"not a SparseMatrix");
3244 MFEM_VERIFY(
width ==
height,
"not a square matrix");
3257 umfpack_di_report_status(
Control, status);
3259 " umfpack_di_symbolic() failed!");
3267 umfpack_di_report_status(
Control, status);
3269 " umfpack_di_numeric() failed!");
3271 umfpack_di_free_symbolic(&
Symbolic);
3275 SuiteSparse_long status;
3279 AI =
new SuiteSparse_long[
width + 1];
3280 AJ =
new SuiteSparse_long[Ap[
width]];
3281 for (
int i = 0; i <=
width; i++)
3283 AI[i] = (SuiteSparse_long)(Ap[i]);
3285 for (
int i = 0; i < Ap[
width]; i++)
3287 AJ[i] = (SuiteSparse_long)(Ai[i]);
3295 umfpack_dl_report_status(
Control, status);
3297 " umfpack_dl_symbolic() failed!");
3305 umfpack_dl_report_status(
Control, status);
3307 " umfpack_dl_numeric() failed!");
3309 umfpack_dl_free_symbolic(&
Symbolic);
3316 mfem_error(
"UMFPackSolver::Mult : matrix is not set!"
3317 " Call SetOperator first!");
3329 umfpack_di_report_status(
Control, status);
3330 mfem_error(
"UMFPackSolver::Mult : umfpack_di_solve() failed!");
3335 SuiteSparse_long status =
3342 umfpack_dl_report_status(
Control, status);
3343 mfem_error(
"UMFPackSolver::Mult : umfpack_dl_solve() failed!");
3351 mfem_error(
"UMFPackSolver::MultTranspose : matrix is not set!"
3352 " Call SetOperator first!");
3364 umfpack_di_report_status(
Control, status);
3366 " umfpack_di_solve() failed!");
3371 SuiteSparse_long status =
3378 umfpack_dl_report_status(
Control, status);
3380 " umfpack_dl_solve() failed!");
3393 umfpack_di_free_numeric(&
Numeric);
3397 umfpack_dl_free_numeric(&
Numeric);
3412 "Had Numeric pointer in KLU, but not Symbolic");
3420 MFEM_VERIFY(
mat != NULL,
"not a SparseMatrix");
3428 MFEM_VERIFY(
width ==
height,
"not a square matrix");
3440 MFEM_VERIFY(
mat != NULL,
3441 "KLUSolver::Mult : matrix is not set! Call SetOperator first!");
3454 MFEM_VERIFY(
mat != NULL,
3455 "KLUSolver::Mult : matrix is not set! Call SetOperator first!");
3482 for (
int i = 0; i < block_dof.
NumRows(); ++i)
3487 block_solvers[i].SetOperator(sub_A);
3496 for (
int i = 0; i < block_dof.
NumRows(); ++i)
3501 block_solvers[i].Mult(sub_rhs, sub_sol);
3515 add(-1.0, z, 1.0, x, z);
3532 add(-1.0, z, 1.0, x, z);
3541 :
Solver(0, false), global_size(-1)
3549 :
Solver(0, false), mycomm(mycomm_), global_size(-1), parallel(true) { }
3557 MFEM_VERIFY(
height ==
width,
"Solver must be a square Operator!");
3563 MFEM_VERIFY(solver,
"Solver hasn't been set, call SetSolver() first.");
3567 MFEM_VERIFY(
height ==
width,
"Solver must be a square Operator!");
3573 MFEM_VERIFY(solver,
"Solver hasn't been set, call SetSolver() first.");
3575 "solver was modified externally! call SetSolver() again!");
3576 MFEM_VERIFY(
height ==
b.Size(),
"incompatible input Vector size!");
3577 MFEM_VERIFY(
height == x.
Size(),
"incompatible output Vector size!");
3580 Orthogonalize(
b, b_ortho);
3586 solver->
Mult(b_ortho, x);
3589 Orthogonalize(x, x);
3592void OrthoSolver::Orthogonalize(
const Vector &v,
Vector &v_ortho)
const
3594 if (global_size == -1)
3600 MPI_Allreduce(MPI_IN_PLACE, &global_size, 1, HYPRE_MPI_BIG_INT,
3613 MPI_Allreduce(MPI_IN_PLACE, &global_sum, 1, MPITypeMap<real_t>::mpi_type,
3618 real_t ratio = global_sum /
static_cast<real_t>(global_size);
3623 for (
int i = 0; i < v_ortho.
Size(); ++i)
3625 v_ortho(i) = v(i) - ratio;
3632 bool op_is_symmetric,
3634 :
Solver(op.NumRows()), aux_map_(aux_map, own_aux_map)
3636 aux_system_.
Reset(
RAP(&op, aux_map));
3639 aux_smoother_.
As<
HypreSmoother>()->SetOperatorSymmetry(op_is_symmetric);
3642void AuxSpaceSmoother::Mult(
const Vector &x,
Vector &y,
bool transpose)
const
3647 Vector aux_sol(aux_rhs.Size());
3654 aux_smoother_->
Mult(aux_rhs, aux_sol);
3658 aux_map_->
Mult(aux_sol, y);
3662#ifdef MFEM_USE_LAPACK
3665 :
Solver(0), mat(nullptr), const_tol_(1.0e-14), min_nnz_(0),
3666 max_nnz_(0), verbosity_(0), res_change_termination_tol_(1.0e-4),
3667 zero_tol_(1.0e-14), rhs_delta_(1.0e-11), n_outer_(100000),
3668 n_inner_(100000), nStallCheck_(100), normalize_(true),
3675 MFEM_VERIFY(mat,
"NNLSSolver operator must be of type DenseMatrix");
3688 qr_residual_mode_ = qr_residual_mode;
3691 NNLS_qrres_on_ =
true;
3700 MFEM_VERIFY(rhs_lb.
Size() == m && rhs_ub.
Size() == m,
"");
3706 Vector rhs_halfgap = rhs_ub;
3707 rhs_halfgap -= rhs_lb;
3710 Vector rhs_avg_glob = rhs_avg;
3711 Vector rhs_halfgap_glob = rhs_halfgap;
3712 Vector halfgap_target(m);
3713 halfgap_target = 1.0e3 * const_tol_;
3718 for (
int i=0; i<m; ++i)
3720 const real_t s = halfgap_target(i) / rhs_halfgap_glob(i);
3721 row_scaling_[i] = s;
3723 rhs_lb(i) = (rhs_avg(i) * s) - halfgap_target(i);
3724 rhs_ub(i) = (rhs_avg(i) * s) + halfgap_target(i);
3730 MFEM_VERIFY(mat,
"NNLSSolver operator must be of type DenseMatrix");
3732 mat->
Mult(w, rhs_ub);
3733 rhs_ub *= row_scaling_;
3738 for (
int i=0; i<rhs_ub.
Size(); ++i)
3740 rhs_lb(i) -= rhs_delta_;
3741 rhs_ub(i) += rhs_delta_;
3745 Solve(rhs_lb, rhs_ub, sol);
3750 for (
int i=0; i<sol.
Size(); ++i)
3758 mfem::out <<
"Number of nonzeros in NNLSSolver solution: " << nnz
3759 <<
", out of " << sol.
Size() << endl;
3763 mat->
Mult(sol, res);
3764 res *= row_scaling_;
3770 const real_t relNorm = res.
Norml2() / std::max(normGsol, normRHS);
3771 mfem::out <<
"Relative residual norm for NNLSSolver solution of Gs = Gw: "
3782 MFEM_VERIFY(rhs_lb.
Size() == m && rhs_lb.
Size() == m && soln.
Size() == n,
"");
3783 MFEM_VERIFY(n >= m,
"NNLSSolver system cannot be over-determined.");
3795 Vector rhs_halfgap(rhs_ub);
3796 rhs_halfgap -= rhs_lb;
3799 Vector rhs_avg_glob(rhs_avg);
3800 Vector rhs_halfgap_glob(rhs_halfgap);
3809 std::vector<unsigned int> nz_ind(m);
3816 int min_nnz_cap = std::min(
static_cast<int>(min_nnz_), std::min(m,n));
3818 std::vector<real_t> l2_res_hist;
3819 std::vector<unsigned int> stalled_indices;
3820 int stalledFlag = 0;
3821 int num_stalled = 0;
3822 int nz_ind_zero = 0;
3825 Vector soln_nz_glob_up(m);
3828 Vector mat_0_data(m * n);
3829 Vector mat_qr_data(m * n);
3830 Vector submat_data(m * n);
3838 std::vector<real_t> work;
3839 int n_outer_iter = 0;
3840 int n_total_inner_iter = 0;
3847 res_glob = rhs_avg_glob;
3848 Vector qt_rhs_glob = rhs_avg_glob;
3849 Vector qqt_rhs_glob = qt_rhs_glob;
3850 Vector sub_qt = rhs_avg_glob;
3856 Vector rhs_scaled(rhs_halfgap_glob);
3858 rhs_scaled *= row_scaling_;
3861 mu_tol = 1.0e-15 * tmp.
Max();
3867 for (
int oiter = 0; oiter < n_outer_; ++oiter)
3871 rmax = fabs(res_glob(0)) - rhs_halfgap_glob(0);
3872 for (
int i=1; i<m; ++i)
3874 rmax = std::max(rmax, fabs(res_glob(i)) - rhs_halfgap_glob(i));
3877 l2_res_hist.push_back(res_glob.
Norml2());
3881 mfem::out <<
"NNLS " << oiter <<
" " << n_total_inner_iter <<
" " << m
3882 <<
" " << n <<
" " << n_glob <<
" " << rmax <<
" "
3883 << l2_res_hist[oiter] << endl;
3885 if (rmax <= const_tol_ && n_glob >= min_nnz_cap)
3889 mfem::out <<
"NNLS target tolerance met" << endl;
3895 if (n_glob >= max_nnz_)
3899 mfem::out <<
"NNLS target nnz met" << endl;
3909 mfem::out <<
"NNLS system is square... exiting" << endl;
3916 if (oiter > nStallCheck_)
3920 for (
int i=0; i<nStallCheck_/2; ++i)
3922 mean0 += l2_res_hist[oiter - i];
3923 mean1 += l2_res_hist[oiter - (nStallCheck_) - i];
3926 real_t mean_res_change = (mean1 / mean0) - 1.0;
3927 if (std::abs(mean_res_change) < res_change_termination_tol_)
3931 mfem::out <<
"NNLSSolver stall detected... exiting" << endl;
3939 res_glob *= row_scaling_;
3942 for (
int i = 0; i < n_nz_ind; ++i)
3944 mu(nz_ind[i]) = 0.0;
3946 for (
unsigned int i = 0; i < stalled_indices.size(); ++i)
3948 mu(stalled_indices[i]) = 0.0;
3955 num_stalled = stalled_indices.size();
3956 if (num_stalled > 0)
3960 mfem::out <<
"NNLS Lagrange multiplier is below the minimum "
3961 <<
"threshold: mumax = " << mumax <<
", mutol = "
3962 << mu_tol <<
"\n" <<
" Resetting stalled indices "
3963 <<
"vector of size " << num_stalled <<
"\n";
3965 stalled_indices.resize(0);
3969 for (
int i = 0; i < n_nz_ind; ++i)
3971 mu(nz_ind[i]) = 0.0;
3981 for (
int i=1; i<n; ++i)
3992 nz_ind[n_nz_ind] = imax;
3997 mfem::out <<
"Found next index: " << imax <<
" " << mumax << endl;
4000 for (
int i=0; i<m; ++i)
4002 mat_0_data(i + (n_glob*m)) = (*mat)(i,imax) * row_scaling_[i];
4003 mat_qr_data(i + (n_glob*m)) = mat_0_data(i + (n_glob*m));
4006 i_qr_start = n_glob;
4011 mfem::out <<
"Updated matrix with new index" << endl;
4014 for (
int iiter = 0; iiter < n_inner_; ++iiter)
4016 ++n_total_inner_iter;
4019 const bool incremental_update =
true;
4020 n_update = n_glob - i_qr_start;
4021 m_update = m - i_qr_start;
4022 if (incremental_update)
4028 MFEM_LAPACK_PREFIX(ormqr_)(&lside, &
trans, &m, &n_update,
4029 &i_qr_start, mat_qr_data.
GetData(), &m,
4031 mat_qr_data.
GetData() + (i_qr_start * m),
4032 &m, work.data(), &lwork, &info);
4033 MFEM_VERIFY(info == 0,
"");
4034 lwork =
static_cast<int>(work[0]);
4036 MFEM_LAPACK_PREFIX(ormqr_)(&lside, &
trans, &m, &n_update,
4037 &i_qr_start, mat_qr_data.
GetData(), &m,
4039 mat_qr_data.
GetData() + (i_qr_start * m),
4040 &m, work.data(), &lwork, &info);
4041 MFEM_VERIFY(info == 0,
"");
4048 for (
int i=0; i<m_update; ++i)
4049 for (
int j=0; j<n_update; ++j)
4051 submat_data[i + (j * m_update)] =
4052 mat_qr_data[i + i_qr_start + ((j + i_qr_start) * m)];
4056 for (
int j=0; j<n_update; ++j)
4058 sub_tau[j] = tau[i_qr_start + j];
4061 MFEM_LAPACK_PREFIX(geqrf_)(&m_update, &n_update, submat_data.
GetData(),
4062 &m_update, sub_tau.
GetData(), work.data(),
4064 MFEM_VERIFY(info == 0,
"");
4065 lwork =
static_cast<int>(work[0]);
4066 if (lwork == 0) { lwork = 1; }
4068 MFEM_LAPACK_PREFIX(geqrf_)(&m_update, &n_update, submat_data.
GetData(),
4069 &m_update, sub_tau.
GetData(), work.data(),
4071 MFEM_VERIFY(info == 0,
"");
4074 for (
int i=0; i<m_update; ++i)
4075 for (
int j=0; j<n_update; ++j)
4077 mat_qr_data[i + i_qr_start + ((j + i_qr_start)* m)] =
4078 submat_data[i + (j * m_update)];
4081 for (
int j=0; j<n_update; ++j)
4083 tau[i_qr_start + j] = sub_tau[j];
4089 for (
int i=0; i<m; ++i)
4090 for (
int j=0; j<n_glob; ++j)
4092 mat_qr_data(i + (j*m)) = mat_0_data(i + (j*m));
4099 MFEM_LAPACK_PREFIX(geqrf_)(&m, &n_glob, mat_qr_data.
GetData(), &m,
4100 tau.
GetData(), work.data(), &lwork, &info);
4101 MFEM_VERIFY(info == 0,
"");
4102 lwork =
static_cast<int>(work[0]);
4104 MFEM_LAPACK_PREFIX(geqrf_)(&m, &n_glob, mat_qr_data.
GetData(), &m,
4105 tau.
GetData(), work.data(), &lwork, &info);
4106 MFEM_VERIFY(info == 0,
"");
4111 mfem::out <<
"Updated QR " << iiter << endl;
4115 if (incremental_update && iiter == 0)
4125 for (
int i=0; i<m_update; ++i)
4127 submat_data[i] = mat_qr_data[i + i_qr_start + (i_qr_start * m)];
4128 sub_qt[i] = qt_rhs_glob[i + i_qr_start];
4131 sub_tau[0] = tau[i_qr_start];
4133 MFEM_LAPACK_PREFIX(ormqr_)(&lside, &
trans, &m_update, &ione, &ione,
4134 submat_data.
GetData(), &m_update,
4136 &m_update, work.data(), &lwork, &info);
4137 MFEM_VERIFY(info == 0,
"");
4138 lwork =
static_cast<int>(work[0]);
4140 MFEM_LAPACK_PREFIX(ormqr_)(&lside, &
trans, &m_update, &ione, &ione,
4141 submat_data.
GetData(), &m_update,
4143 &m_update, work.data(), &lwork, &info);
4144 MFEM_VERIFY(info == 0,
"");
4146 for (
int i=0; i<m_update; ++i)
4148 qt_rhs_glob[i + i_qr_start] = sub_qt[i];
4154 qt_rhs_glob = rhs_avg_glob;
4157 MFEM_LAPACK_PREFIX(ormqr_)(&lside, &
trans, &m, &ione, &n_glob,
4160 work.data(), &lwork, &info);
4161 MFEM_VERIFY(info == 0,
"");
4162 lwork =
static_cast<int>(work[0]);
4164 MFEM_LAPACK_PREFIX(ormqr_)(&lside, &
trans, &m, &ione, &n_glob,
4167 work.data(), &lwork, &info);
4168 MFEM_VERIFY(info == 0,
"");
4173 mfem::out <<
"Updated rhs " << iiter << endl;
4180 MFEM_LAPACK_PREFIX(trsm_)(&lside, &upper, ¬rans, &nounit,
4181 &n_glob, &ione, &fone,
4187 mfem::out <<
"Solved triangular system " << iiter << endl;
4192 real_t smin = n_glob > 0 ? vec1(0) : 0.0;
4193 for (
int i=0; i<n_glob; ++i)
4195 soln_nz_glob_up(i) = vec1(i);
4196 smin = std::min(smin, soln_nz_glob_up(i));
4199 if (smin > zero_tol_)
4202 for (
int i=0; i<n_glob; ++i)
4204 soln_nz_glob(i) = soln_nz_glob_up(i);
4215 mfem::out <<
"Start pruning " << iiter << endl;
4216 for (
int i = 0; i < n_glob; ++i)
4218 if (soln_nz_glob_up(i) <= zero_tol_)
4220 mfem::out << i <<
" " << n_glob <<
" " << soln_nz_glob_up(i) << endl;
4225 if (soln_nz_glob_up(n_glob - 1) <= zero_tol_)
4232 mfem::out <<
"Detected stall due to adding and removing same "
4233 <<
"column. Switching to QR residual calculation "
4234 <<
"method." << endl;
4238 mfem::out <<
"Detected stall due to adding and removing same"
4239 <<
" column. Exiting now." << endl;
4246 NNLS_qrres_on_ =
true;
4253 for (
int i = 0; i < n_glob; ++i)
4255 if (soln_nz_glob_up(i) <= zero_tol_)
4257 alpha = std::min(
alpha, soln_nz_glob(i)/(soln_nz_glob(i) - soln_nz_glob_up(i)));
4262 for (
int i = 0; i < n_glob; ++i)
4264 soln_nz_glob(i) +=
alpha*(soln_nz_glob_up(i) - soln_nz_glob(i));
4265 if (i == 0 || soln_nz_glob(i) < smin)
4267 smin = soln_nz_glob(i);
4271 while (smin > zero_tol_)
4279 smin = soln_nz_glob(0);
4280 for (
int i = 1; i < n_glob; ++i)
4282 if (soln_nz_glob(i) < smin)
4284 smin = soln_nz_glob(i);
4289 alpha = soln_nz_glob(index_min)/(soln_nz_glob(index_min)
4290 - soln_nz_glob_up(index_min));
4293 for (
int i = 0; i < n_glob; ++i)
4295 soln_nz_glob(i) +=
alpha*(soln_nz_glob_up(i) - soln_nz_glob(i));
4300 i_qr_start = n_glob+1;
4306 smin = n_glob > 0 ? soln_nz_glob(0) : 0.0;
4307 for (
int i=1; i<n_glob; ++i)
4309 smin = std::min(smin, soln_nz_glob(i));
4312 if (smin < zero_tol_)
4321 if (zero_ibool == 0)
4330 for (
int i = 0; i < n_glob; ++i)
4332 if (soln_nz_glob(i) < zero_tol_)
4338 MFEM_VERIFY(ind_zero != -1,
"");
4341 for (
int i = 0; i < ind_zero; ++i)
4348 for (
int i=0; i<m; ++i)
4349 for (
int j=ind_zero; j<n_glob-1; ++j)
4351 mat_qr_data(i + (j*m)) = mat_0_data(i + ((j+1)*m));
4356 for (
int i=0; i<m; ++i)
4357 for (
int j=ind_zero; j<n_glob-1; ++j)
4359 mat_0_data(i + (j*m)) = mat_qr_data(i + (j*m));
4364 for (
int i = nz_ind_zero; i < n_nz_ind-1; ++i)
4366 nz_ind[i] = nz_ind[i+1];
4371 for (
int i = ind_zero; i < n_glob-1; ++i)
4373 soln_nz_glob(i) = soln_nz_glob(i+1);
4376 i_qr_start = std::min(i_qr_start, ind_zero);
4382 mfem::out <<
"Finished pruning " << iiter << endl;
4387 if (stalledFlag == 1)
4391 num_stalled = stalled_indices.size();
4392 stalled_indices.resize(num_stalled + 1);
4393 stalled_indices[num_stalled] = imax;
4396 mfem::out <<
"Adding index " << imax <<
" to stalled index list "
4397 <<
"of size " << num_stalled << endl;
4402 if (!NNLS_qrres_on_)
4404 res_glob = rhs_avg_glob;
4406 MFEM_LAPACK_PREFIX(gemv_)(¬rans, &m, &n_glob, &fmone,
4408 soln_nz_glob.
GetData(), &ione, &fone,
4418 for (
int i=0; i<n_glob; ++i)
4420 qqt_rhs_glob(i) = qt_rhs_glob(i);
4423 MFEM_LAPACK_PREFIX(ormqr_)(&lside, ¬rans, &m, &ione, &n_glob,
4426 work.data(), &lwork, &info);
4428 MFEM_VERIFY(info == 0,
"");
4429 lwork =
static_cast<int>(work[0]);
4431 MFEM_LAPACK_PREFIX(ormqr_)(&lside, ¬rans, &m, &ione, &n_glob,
4434 work.data(), &lwork, &info);
4435 MFEM_VERIFY(info == 0,
"");
4436 res_glob = rhs_avg_glob;
4437 res_glob -= qqt_rhs_glob;
4442 mfem::out <<
"Computed residual" << endl;
4449 MFEM_VERIFY(n_glob == n_nz_ind,
"");
4451 for (
int i = 0; i < n_glob; ++i)
4453 soln(nz_ind[i]) = soln_nz_glob(i);
4458 mfem::out <<
"NNLS solver: m = " << m <<
", n = " << n
4459 <<
", outer_iter = " << n_outer_iter <<
", inner_iter = "
4460 << n_total_inner_iter;
4468 mfem::out << endl <<
"Warning, NNLS convergence stalled: "
4469 << (exit_flag == 2) << endl;
4470 mfem::out <<
"resErr = " << rmax <<
" vs tol = " << const_tol_
4471 <<
"; mumax = " << mumax <<
" vs tol = " << mu_tol << endl;
const T * HostRead() const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), false).
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
int Size() const
Return the logical size of the array.
bool UseDevice() const
Return the device flag of the Memory object used by the Array.
void MakeRef(T *data_, int size_, bool own_data=false)
Make this Array a reference to a pointer.
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
AuxSpaceSmoother(const HypreParMatrix &op, HypreParMatrix *aux_map, bool op_is_symmetric=true, bool own_aux_map=false)
void Mult(const Vector &b, Vector &x) const override
Iterative solution of the linear system using the BiCGSTAB method.
void SetOperator(const Operator &op) override
Set/update the solver for the given operator.
Reordering
The reordering method used by the BlockILU factorization.
BlockILU(int block_size_, Reordering reordering_=Reordering::MINIMUM_DISCARDED_FILL, int k_fill_=0)
void Mult(const Vector &b, Vector &x) const
Solve the system LUx = b, where L and U are the block ILU factors.
void SetOperator(const Operator &op)
Conjugate gradient method.
void Mult(const Vector &b, Vector &x) const override
Iterative solution of the linear system using the Conjugate Gradient method.
void SetOperator(const Operator &op) override
Set/update the solver for the given operator.
Data type dense matrix using column-major storage.
void AddMult_a(real_t a, const Vector &x, Vector &y) const
y += a * A.x
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
void Add(const real_t c, const DenseMatrix &A)
Adds the matrix A multiplied by the number c to the matrix.
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
void Mult(const Vector &x, Vector &y) const override
Direct solution of the block diagonal linear system.
DirectSubBlockSolver(const SparseMatrix &A, const SparseMatrix &block_dof)
void Mult(const Vector &b, Vector &x) const override
Iterative solution of the linear system using the FGMRES method.
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
void SetKDim(int dim)
Set the number of iteration to perform between restarts, default is 50.
void Mult(const Vector &b, Vector &x) const override
Iterative solution of the linear system using the GMRES method.
Wrapper for hypre's ParCSR matrix class.
void GetDiag(Vector &diag) const
Get the local diagonal of the matrix.
Parallel smoothers in hypre.
virtual void MonitorSolution(int it, real_t norm, const Vector &x, bool final)
Monitor the solution vector x.
const class IterativeSolver * iter_solver
The last IterativeSolver to which this controller was attached.
virtual void MonitorResidual(int it, real_t norm, const Vector &r, bool final)
Monitor the solution vector r.
Abstract base class for iterative solver.
real_t abs_tol
Absolute tolerance.
real_t rel_tol
Relative tolerance.
PrintLevel print_options
Output behavior for the iterative solver.
real_t GetFinalNorm() const
Returns the final residual norm after termination of the solver during the last call to Mult().
void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
virtual real_t Dot(const Vector &x, const Vector &y) const
Return the standard (l2, i.e., Euclidean) inner product of x and y.
void SetRelTol(real_t rtol)
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
int print_level
(DEPRECATED) Legacy print level definition, which is left for compatibility with custom iterative sol...
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
int max_iter
Limit for the number of iterations the solver is allowed to do.
void SetMaxIter(int max_it)
bool Monitor(int it, real_t norm, const Vector &r, const Vector &x, bool final=false) const
Monitor both the residual r and the solution x.
MPI_Comm GetComm() const
Return the associated MPI communicator, or MPI_COMM_NULL if no communicator is set.
static int GuessLegacyPrintLevel(PrintLevel)
Use some heuristics to guess a legacy print level corresponding to the given PrintLevel.
IterativeSolverMonitor * monitor
bool GetConverged() const
Returns true if the last call to Mult() converged successfully.
PrintLevel FromLegacyPrintLevel(int)
Convert a legacy print level integer to a PrintLevel object.
void SetAbsTol(real_t atol)
real_t Norm(const Vector &x) const
Return the inner product norm of x, using the inner product defined by Dot()
void MultTranspose(const Vector &b, Vector &x) const override
Direct solution of the transposed linear system using KLU.
void Mult(const Vector &b, Vector &x) const override
Direct solution of the linear system using KLU.
void SetOperator(const Operator &op) override
Set/update the solver for the given operator.
Array< Vector * > skArray
void Mult(const Vector &b, Vector &x) const override
Solve the nonlinear system with right-hand side b.
Array< Vector * > ykArray
void Solve(int m, int n, real_t *X) const override
void Mult(const Vector &b, Vector &x) const override
Iterative solution of the linear system using the MINRES method.
void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
void SetPreconditioner(Solver &pr) override
This should be called before SetOperator.
void Solve(const Vector &rhs_lb, const Vector &rhs_ub, Vector &soln) const
Solve the NNLS problem. Specifically, we find a vector soln, such that rhs_lb < mat*soln < rhs_ub is ...
QRresidualMode
Enumerated types of QRresidual mode. Options are 'off': the residual is calculated normally,...
void Mult(const Vector &w, Vector &sol) const override
Compute the non-negative least squares solution to the underdetermined system.
void SetQRResidualMode(const QRresidualMode qr_residual_mode)
Set the residual calculation mode for the NNLS solver. See QRresidualMode enum above for details.
void SetOperator(const Operator &op) override
The operator must be a DenseMatrix.
void NormalizeConstraints(Vector &rhs_lb, Vector &rhs_ub) const
Normalize the constraints such that the tolerances for each constraint (i.e. (UB - LB)/2) are equal....
void SetAdaptiveLinRtol(const int type=2, const real_t rtol0=0.5, const real_t rtol_max=0.9, const real_t alpha=0.5 *(1.0+sqrt(5.0)), const real_t gamma=1.0)
Enable adaptive linear solver relative tolerance algorithm.
void AdaptiveLinRtolPreSolve(const Vector &x, const int it, const real_t fnorm) const
Method for the adaptive linear solver rtol invoked before the linear solve.
virtual real_t ComputeScalingFactor(const Vector &x, const Vector &b) const
This method can be overloaded in derived classes to implement line search algorithms.
void Mult(const Vector &b, Vector &x) const override
Solve the nonlinear system with right-hand side b.
void AdaptiveLinRtolPostSolve(const Vector &x, const Vector &b, const int it, const real_t fnorm) const
Method for the adaptive linear solver rtol invoked after the linear solve.
virtual void ProcessNewState(const Vector &x) const
This method can be overloaded in derived classes to perform computations that need knowledge of the n...
void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Chebyshev accelerated smoothing with given vector, no matrix necessary.
void Mult(const Vector &x, Vector &y) const
Approach the solution of the linear system by applying Chebyshev smoothing.
OperatorChebyshevSmoother(const Operator &oper_, const Vector &d, const Array< int > &ess_tdof_list, int order, real_t max_eig_estimate)
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get().
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Jacobi smoothing for a given bilinear form (no matrix necessary).
void Setup(const Vector &diag)
OperatorJacobiSmoother(const real_t damping=1.0)
Default constructor: the diagonal will be computed by subsequent calls to SetOperator() using the Ope...
void SetOperator(const Operator &op)
Recompute the diagonal using the method AssembleDiagonal of the given new Operator,...
void Mult(const Vector &x, Vector &y) const
Approach the solution of the linear system by applying Jacobi smoothing.
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
int width
Dimension of the input / number of columns in the matrix.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
int height
Dimension of the output / number of rows in the matrix.
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
int NumCols() const
Get the number of columns (size of input) of the Operator. Synonym with Width().
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
virtual void AssembleDiagonal(Vector &diag) const
Computes the diagonal entries into diag. Typically, this operation only makes sense for linear Operat...
int NumRows() const
Get the number of rows (size of output) of the Operator. Synonym with Height().
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
virtual Operator & GetGradient(const Vector &x) const
Evaluate the gradient operator at the point x. The default behavior in class Operator is to generate ...
const Operator * C
Not owned, some can remain unused (NULL).
void SetSolutionBounds(const Vector &xl, const Vector &xh)
void SetEqualityConstraint(const Vector &c)
int GetNumConstraints() const
OptimizationProblem(int insize, const Operator *C_, const Operator *D_)
In parallel, insize is the number of the local true dofs.
void SetInequalityConstraint(const Vector &dl, const Vector &dh)
const OptimizationProblem * problem
void SetOperator(const Operator &op) override
Set the Operator that is the OrthoSolver is to invert (approximately).
void SetSolver(Solver &s)
Set the solver used by the OrthoSolver.
void Mult(const Vector &b, Vector &x) const override
Perform the action of the OrthoSolver: P * solver * P where P is the projection to the subspace of ve...
PowerMethod helper class to estimate the largest eigenvalue of an operator using the iterative power ...
real_t EstimateLargestEigenvalue(Operator &opr, Vector &v0, int numSteps=10, real_t tolerance=1e-8, int seed=12345)
Returns an estimate of the largest eigenvalue of the operator opr using the iterative power method.
General product operator: x -> (A*B)(x) = A(B(x)).
void Mult(const Vector &x, Vector &y) const override
Solution of the linear system using a product of subsolvers.
void MultTranspose(const Vector &x, Vector &y) const override
Solution of the transposed linear system using a product of subsolvers.
void MonitorResidual(int it, real_t norm, const Vector &r, bool final) override
Monitor the solution vector r.
const Array< int > * ess_dofs_list
Not owned.
void Mult(const Vector &xt, Vector &x) const override
void SetBounds(const Vector &lo_, const Vector &hi_)
void SetLinearConstraint(const Vector &w_, real_t a_)
void print_iteration(int it, real_t r, real_t l) const
void SetOptimizationProblem(const OptimizationProblem &prob) override
real_t solve(real_t l, const Vector &xt, Vector &x, int &nclip) const
Solve QP at fixed lambda.
Stationary linear iteration: x <- x + B (b - A x)
void Mult(const Vector &b, Vector &x) const override
Iterative solution of the linear system using Stationary Linear Iteration.
void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
virtual void SetOperator(const Operator &op)=0
Set/update the solver for the given operator.
const int * HostReadJ() const
void GetDiag(Vector &d) const
Returns the Diagonal of A.
void MultTranspose(const Vector &x, Vector &y) const override
Multiply a vector with the transposed matrix. y = At * x.
void GetSubMatrix(const Array< int > &rows, const Array< int > &cols, DenseMatrix &subm) const
bool Finalized() const
Returns whether or not CSR format has been finalized.
const real_t * HostReadData() const
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
int RowSize(const int i) const
Returns the number of elements in row i.
real_t * HostReadWriteData()
void Mult(const Vector &x, Vector &y) const override
Matrix vector multiplication.
real_t * GetData()
Return the element data, i.e. the array A.
void SortColumnIndices()
Sort the column indices corresponding to each row.
int * GetJ()
Return the array J.
int * GetI()
Return the array I.
const int * HostReadI() const
real_t Control[UMFPACK_CONTROL]
void SetOperator(const Operator &op) override
Factorize the given Operator op which must be a SparseMatrix.
void MultTranspose(const Vector &b, Vector &x) const override
Direct solution of the transposed linear system using UMFPACK.
real_t Info[UMFPACK_INFO]
void Mult(const Vector &b, Vector &x) const override
Direct solution of the linear system using UMFPACK.
virtual const real_t * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
virtual const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
virtual real_t * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
void SetDataAndSize(real_t *d, int s)
Set the Vector data and size.
void AddElementVector(const Array< int > &dofs, const Vector &elemvect)
Add elements of the elemvect Vector to the entries listed in dofs. Negative dof values cause the -dof...
real_t Norml2() const
Returns the l2 norm of the vector.
Vector & Set(const real_t a, const Vector &x)
(*this) = a * x
real_t Max() const
Returns the maximal element of the vector.
int Size() const
Returns the size of the vector.
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
real_t Sum() const
Return the sum of the vector entries.
void SetSize(int s)
Resize the vector to size s.
virtual real_t * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
virtual real_t * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
virtual real_t * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Vector & Add(const real_t a, const Vector &Va)
(*this) += a * Va
void trans(const Vector &u, Vector &x)
int BiCGSTAB(const Operator &A, Vector &x, const Vector &b, Solver &M, int &max_iter, real_t &tol, real_t atol, int printit)
BiCGSTAB method. (tolerances are squared)
void Swap(Array< T > &, Array< T > &)
void mfem_error(const char *msg)
void AddMult_a(real_t alpha, const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += alpha * B * C.
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
void add(const Vector &v1, const Vector &v2, Vector &v)
bool IsFinite(const real_t &val)
real_t InnerProduct(HypreParVector *x, HypreParVector *y)
int GMRES(const Operator &A, Vector &x, const Vector &b, Solver &M, int &max_iter, int m, real_t &tol, real_t atol, int printit)
GMRES method. (tolerances are squared)
void ApplyPlaneRotation(real_t &dx, real_t &dy, real_t &cs, real_t &sn)
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
MemoryType GetMemoryType(MemoryClass mc)
Return a suitable MemoryType for a given MemoryClass.
int aGMRES(const Operator &A, Vector &x, const Vector &b, const Operator &M, int &max_iter, int m_max, int m_min, int m_step, real_t cf, real_t &tol, real_t &atol, int printit)
void SLI(const Operator &A, const Vector &b, Vector &x, int print_iter, int max_num_iter, real_t RTOLERANCE, real_t ATOLERANCE)
Stationary linear iteration. (tolerances are squared)
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, real_t RTOLERANCE, real_t ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
void CG(const Operator &A, const Vector &b, Vector &x, int print_iter, int max_num_iter, real_t RTOLERANCE, real_t ATOLERANCE)
Conjugate gradient method. (tolerances are squared)
void GeneratePlaneRotation(real_t &dx, real_t &dy, real_t &cs, real_t &sn)
void MINRES(const Operator &A, const Vector &b, Vector &x, int print_it, int max_it, real_t rtol, real_t atol)
MINRES method without preconditioner. (tolerances are squared)
void subtract(const Vector &x, const Vector &y, Vector &z)
MemoryType
Memory types supported by MFEM.
void MinimumDiscardedFillOrdering(SparseMatrix &C, Array< int > &p)
void forall(int N, lambda &&body)
void Update(Vector &x, int k, DenseMatrix &h, Vector &s, Array< Vector * > &v)
void Add(const DenseMatrix &A, const DenseMatrix &B, real_t alpha, DenseMatrix &C)
C = A + alpha*B.
real_t p(const Vector &x, real_t t)
Settings for the output behavior of the IterativeSolver.
PrintLevel & Iterations()
bool errors
If a fatal problem has been detected the failure will be reported to mfem::err.
bool iterations
Detailed information about each iteration will be reported to mfem::out.
PrintLevel & FirstAndLast()
bool warnings
If a non-fatal problem has been detected some context-specific information will be reported to mfem::...
bool first_and_last
Information about the first and last iteration will be printed to mfem::out.
bool summary
A summary of the solver process will be reported after the last iteration to mfem::out.
Helper struct to convert a C++ type to an MPI type.