13 #include "../general/globals.hpp"
56 if (dot_prod_type == 0)
72 if (dot_prod_type == 0)
79 MPI_Comm_rank(comm, &rank);
151 double r0, nom, nom0, nomold = 1, cf;
167 nom0 = nom =
Dot(
z,
r);
171 nom0 = nom =
Dot(
r,
r);
175 mfem::out <<
" Iteration : " << setw(3) << 0 <<
" (B r, r) = "
214 cf = sqrt(nom/nomold);
216 mfem::out <<
" Iteration : " << setw(3) << i <<
" (B r, r) = "
217 << nom <<
"\tConv. rate: " << cf <<
'\n';
223 mfem::out <<
"Number of SLI iterations: " << i <<
'\n'
224 <<
"Conv. rate: " << cf <<
'\n';
226 mfem::out <<
"(B r_0, r_0) = " << nom0 <<
'\n'
227 <<
"(B r_N, r_N) = " << nom <<
'\n'
228 <<
"Number of SLI iterations: " << i <<
'\n';
242 mfem::err <<
"SLI: No convergence!" <<
'\n';
243 mfem::out <<
"(B r_0, r_0) = " << nom0 <<
'\n'
244 <<
"(B r_N, r_N) = " << nom <<
'\n'
245 <<
"Number of SLI iterations: " <<
final_iter <<
'\n';
249 mfem::out <<
"Average reduction factor = "
256 int print_iter,
int max_num_iter,
257 double RTOLERANCE,
double ATOLERANCE)
269 int print_iter,
int max_num_iter,
270 double RTOLERANCE,
double ATOLERANCE)
293 double r0, den, nom, nom0, betanom,
alpha, beta;
315 nom0 = nom =
Dot(
d,
r);
316 MFEM_ASSERT(
IsFinite(nom),
"nom = " << nom);
320 mfem::out <<
" Iteration : " << setw(3) << 0 <<
" (B r, r) = "
335 MFEM_ASSERT(
IsFinite(den),
"den = " << den);
340 mfem::out <<
"PCG: The operator is not positive definite. (Ad, d) = "
370 MFEM_ASSERT(
IsFinite(betanom),
"betanom = " << betanom);
374 mfem::out <<
" Iteration : " << setw(3) << i <<
" (B r, r) = "
382 mfem::out <<
"Number of PCG iterations: " << i <<
'\n';
386 mfem::out <<
" Iteration : " << setw(3) << i <<
" (B r, r) = "
410 MFEM_ASSERT(
IsFinite(den),
"den = " << den);
415 mfem::out <<
"PCG: The operator is not positive definite. (Ad, d) = "
432 mfem::out <<
" Iteration : " << setw(3) << 0 <<
" (B r, r) = "
438 mfem::out <<
"PCG: No convergence!" <<
'\n';
442 mfem::out <<
"Average reduction factor = "
443 << pow (betanom/nom0, 0.5/
final_iter) <<
'\n';
449 int print_iter,
int max_num_iter,
450 double RTOLERANCE,
double ATOLERANCE)
462 int print_iter,
int max_num_iter,
463 double RTOLERANCE,
double ATOLERANCE)
477 double &cs,
double &sn)
484 else if (fabs(dy) > fabs(dx))
486 double temp = dx / dy;
487 sn = 1.0 / sqrt( 1.0 + temp*temp );
492 double temp = dy / dx;
493 cs = 1.0 / sqrt( 1.0 + temp*temp );
500 double temp = cs * dx + sn * dy;
501 dy = -sn * dx + cs * dy;
511 for (
int i = k; i >= 0; i--)
514 for (
int j = i - 1; j >= 0; j--)
516 y(j) -= h(j,i) * y(i);
520 for (
int j = 0; j <= k; j++)
573 double beta =
Norm(r);
574 MFEM_ASSERT(
IsFinite(beta),
"beta = " << beta);
589 <<
" Iteration : " << setw(3) << 0
590 <<
" ||B r|| = " << beta << (
print_level == 3 ?
" ...\n" :
"\n");
597 if (v[0] == NULL) { v[0] =
new Vector(n); }
598 v[0]->Set(1.0/beta, r);
599 s = 0.0; s(0) = beta;
601 for (i = 0; i <
m && j <=
max_iter; i++, j++)
613 for (k = 0; k <= i; k++)
615 H(k,i) =
Dot(w, *v[k]);
616 w.
Add(-H(k,i), *v[k]);
620 MFEM_ASSERT(
IsFinite(H(i+1,i)),
"Norm(w) = " << H(i+1,i));
621 if (v[i+1] == NULL) { v[i+1] =
new Vector(n); }
622 v[i+1]->Set(1.0/H(i+1,i), w);
624 for (k = 0; k < i; k++)
633 resid = fabs(s(i+1));
634 MFEM_ASSERT(
IsFinite(resid),
"resid = " << resid);
647 mfem::out <<
" Pass : " << setw(2) << (j-1)/
m+1
648 <<
" Iteration : " << setw(3) << j
649 <<
" ||B r|| = " << resid <<
'\n';
671 MFEM_ASSERT(
IsFinite(beta),
"beta = " << beta);
700 for (i = 0; i < v.
Size(); i++)
725 double beta =
Norm(r);
726 MFEM_ASSERT(
IsFinite(beta),
"beta = " << beta);
740 <<
" Iteration : " << setw(3) << 0
741 <<
" || r || = " << beta << endl;
745 for (i= 0; i<=
m; i++)
754 if (v[0] == NULL) { v[0] =
new Vector(b.
Size()); }
756 v[0] ->
Add (1.0/beta, r);
757 s = 0.0; s(0) = beta;
759 for (i = 0; i <
m && j <=
max_iter; i++, j++)
762 if (z[i] == NULL) { z[i] =
new Vector(b.
Size()); }
775 for (k = 0; k <= i; k++)
777 H(k,i) =
Dot( r, *v[k]);
778 r.Add(-H(k,i), (*v[k]));
782 if (v[i+1] == NULL) { v[i+1] =
new Vector(b.
Size()); }
784 v[i+1] ->
Add (1.0/H(i+1,i), r);
786 for (k = 0; k < i; k++)
795 double resid = fabs(s(i+1));
796 MFEM_ASSERT(
IsFinite(resid),
"resid = " << resid);
798 mfem::out <<
" Pass : " << setw(2) << (j-1)/
m+1
799 <<
" Iteration : " << setw(3) << j
800 <<
" || r || = " << resid << endl;
808 for (i= 0; i<=
m; i++)
810 if (v[i]) {
delete v[i]; }
811 if (z[i]) {
delete z[i]; }
827 MFEM_ASSERT(
IsFinite(beta),
"beta = " << beta);
833 for (i= 0; i<=
m; i++)
835 if (v[i]) {
delete v[i]; }
836 if (z[i]) {
delete z[i]; }
842 for (i = 0; i <=
m; i++)
844 if (v[i]) {
delete v[i]; }
845 if (z[i]) {
delete z[i]; }
854 int &max_iter,
int m,
double &tol,
double atol,
int printit)
871 int print_iter,
int max_num_iter,
int m,
double rtol,
double atol)
873 GMRES(A, x, b, B, max_num_iter, m, rtol, atol, print_iter);
895 double resid, tol_goal;
896 double rho_1, rho_2=1.0,
alpha=1.0, beta, omega=1.0;
911 MFEM_ASSERT(
IsFinite(resid),
"resid = " << resid);
913 mfem::out <<
" Iteration : " << setw(3) << 0
914 <<
" ||r|| = " << resid <<
'\n';
918 if (resid <= tol_goal)
932 mfem::out <<
" Iteration : " << setw(3) << i
933 <<
" ||r|| = " << resid <<
'\n';
945 beta = (rho_1/rho_2) * (
alpha/omega);
961 MFEM_ASSERT(
IsFinite(resid),
"resid = " << resid);
962 if (resid < tol_goal)
966 mfem::out <<
" Iteration : " << setw(3) << i
967 <<
" ||s|| = " << resid <<
'\n';
974 mfem::out <<
" Iteration : " << setw(3) << i
975 <<
" ||s|| = " << resid;
992 MFEM_ASSERT(
IsFinite(resid),
"resid = " << resid);
995 mfem::out <<
" ||r|| = " << resid <<
'\n';
997 if (resid < tol_goal)
1019 int &max_iter,
double &tol,
double atol,
int printit)
1028 bicgstab.
Mult(b, x);
1035 int print_iter,
int max_num_iter,
double rtol,
double atol)
1037 BiCGSTAB(A, x, b, B, max_num_iter, rtol, atol, print_iter);
1063 double beta, eta, gamma0, gamma1, sigma0, sigma1;
1064 double alpha, delta, rho1, rho2, rho3, norm_goal;
1084 eta = beta = sqrt(
Dot(*z,
v1));
1085 MFEM_ASSERT(
IsFinite(eta),
"eta = " << eta);
1086 gamma0 = gamma1 = 1.;
1087 sigma0 = sigma1 = 0.;
1091 if (eta <= norm_goal)
1099 mfem::out <<
"MINRES: iteration " << setw(3) << 0 <<
": ||r||_B = "
1112 MFEM_ASSERT(
IsFinite(alpha),
"alpha = " << alpha);
1119 delta = gamma1*alpha - gamma0*sigma1*beta;
1121 rho2 = sigma1*alpha + gamma0*gamma1*beta;
1131 MFEM_ASSERT(
IsFinite(beta),
"beta = " << beta);
1132 rho1 = hypot(delta, beta);
1136 w0.
Set(1./rho1, *z);
1140 add(1./rho1, *z, -rho2/rho1,
w1,
w0);
1144 add(-rho3/rho1,
w0, -rho2/rho1,
w1,
w0);
1145 w0.
Add(1./rho1, *z);
1149 gamma1 = delta/rho1;
1151 x.
Add(gamma1*eta,
w0);
1157 MFEM_ASSERT(
IsFinite(eta),
"eta = " << eta);
1159 if (fabs(eta) <= norm_goal)
1166 mfem::out <<
"MINRES: iteration " << setw(3) << it <<
": ||r||_B = "
1167 << fabs(eta) <<
'\n';
1202 eta = sqrt(
Dot(*z,
v1));
1203 mfem::out <<
"MINRES: iteration " << setw(3) << it <<
": ||r||_B = "
1204 << eta <<
" (re-computed)" <<
'\n';
1209 mfem::out <<
"MINRES: No convergence!\n";
1214 int max_it,
double rtol,
double atol)
1226 int print_it,
int max_it,
double rtol,
double atol)
1244 MFEM_ASSERT(
height ==
width,
"square Operator is required.");
1252 MFEM_ASSERT(
oper != NULL,
"the Operator is not set (use SetOperator).");
1253 MFEM_ASSERT(
prec != NULL,
"the Solver is not set (use SetSolver).");
1256 double norm0, norm, norm_goal;
1270 norm0 = norm =
Norm(
r);
1276 for (it = 0;
true; it++)
1278 MFEM_ASSERT(
IsFinite(norm),
"norm = " << norm);
1281 mfem::out <<
"Newton iteration " << setw(2) << it
1282 <<
" : ||r|| = " << norm;
1285 mfem::out <<
", ||r||/||r_0|| = " << norm/norm0;
1290 if (norm <= norm_goal)
1312 add(x, -c_scale,
c, x);
1329 int m_max,
int m_min,
int m_step,
double cf,
1330 double &tol,
double &atol,
int printit)
1337 Vector s(m+1), cs(m+1), sn(m+1);
1344 double normb = w.Norml2();
1354 double beta = r.
Norml2();
1356 resid = beta / normb;
1358 if (resid * resid <= tol)
1360 tol = resid * resid;
1367 <<
" Iteration : " << setw(3) << 0
1368 <<
" (r, r) = " << beta*beta <<
'\n';
1370 tol *= (normb*normb);
1371 tol = (atol > tol) ? atol : tol;
1375 for (i= 0; i<=m; i++)
1382 while (j <= max_iter)
1385 v[0] ->
Add (1.0/beta, r);
1386 s = 0.0; s(0) = beta;
1390 for (i = 0; i < m && j <= max_iter; i++)
1395 for (k = 0; k <= i; k++)
1397 H(k,i) = w * (*v[k]);
1398 w.
Add(-H(k,i), (*v[k]));
1401 H(i+1,i) = w.Norml2();
1403 v[i+1] ->
Add (1.0/H(i+1,i), w);
1405 for (k = 0; k < i; k++)
1414 resid = fabs(s(i+1));
1417 <<
" Iteration : " << setw(3) << i+1
1418 <<
" (r, r) = " << resid*resid <<
'\n';
1420 if ( resid*resid < tol)
1423 tol = resid * resid;
1425 for (i= 0; i<=m; i++)
1444 if ( resid*resid < tol)
1446 tol = resid * resid;
1448 for (i= 0; i<=m; i++)
1457 if (m - m_step >= m_min)
1470 tol = resid * resid;
1471 for (i= 0; i<=m; i++)
1493 mfem_error(
"SLBQPOptimizer::SetPreconditioner() : "
1494 "not meaningful for this solver");
1499 mfem_error(
"SLBQPOptimizer::SetOperator() : "
1500 "not meaningful for this solver");
1506 mfem::out <<
"SLBQP iteration " << it <<
": residual = " << r
1507 <<
", lambda = " << l <<
'\n';
1529 const double smin = 0.1;
1536 mfem::out <<
"SLBQP bracketing phase" <<
'\n';
1540 r =
solve(l,xt,x,nclip);
1554 llow = l; rlow = r; l = l + dl;
1557 r =
solve(l,xt,x,nclip);
1560 while ((r < 0) && (nclip <
max_iter))
1564 if (s < smin) { s = smin; }
1569 r =
solve(l,xt,x,nclip);
1577 lupp = l; rupp = r; l = l - dl;
1580 r =
solve(l,xt,x,nclip);
1583 while ((r > 0) && (nclip <
max_iter))
1587 if (s < smin) { s = smin; }
1592 r =
solve(l,xt,x,nclip);
1605 mfem::out <<
"SLBQP secant phase" <<
'\n';
1608 s = 1.0 - rlow/rupp; dl = dl/s; l = lupp - dl;
1611 r =
solve(l,xt,x,nclip);
1614 while ( (fabs(r) > tol) && (nclip <
max_iter) )
1620 lupp = l; rupp = r; s = 1.0 - rlow/rupp;
1621 dl = (lupp - llow)/s; l = lupp - dl;
1626 if (s < smin) { s = smin; }
1628 lnew = 0.75*llow + 0.25*l;
1629 if (lnew < l-dl) { lnew = l-dl; }
1630 lupp = l; rupp = r; l = lnew;
1631 s = (lupp - llow)/(lupp - l);
1639 llow = l; rlow = r; s = 1.0 - rlow/rupp;
1640 dl = (lupp - llow)/s; l = lupp - dl;
1645 if (s < smin) { s = smin; }
1647 lnew = 0.75*lupp + 0.25*l;
1648 if (lnew < l+dl) { lnew = l+dl; }
1649 llow = l; rlow = r; l = lnew;
1650 s = (lupp - llow)/(lupp - l);
1655 r =
solve(l,xt,x,nclip);
1664 mfem::err <<
"SLBQP not converged!" <<
'\n';
1674 mfem::out <<
"SLBQP iterations = " << nclip <<
'\n';
1675 mfem::out <<
"SLBQP lambda = " << l <<
'\n';
1676 mfem::out <<
"SLBQP residual = " << r <<
'\n';
1680 #ifdef MFEM_USE_SUITESPARSE
1707 umfpack_di_free_numeric(&
Numeric);
1711 umfpack_dl_free_numeric(&
Numeric);
1716 MFEM_VERIFY(
mat,
"not a SparseMatrix");
1725 MFEM_VERIFY(
width ==
height,
"not a square matrix");
1733 int status = umfpack_di_symbolic(
width,
width, Ap, Ai, Ax, &Symbolic,
1738 umfpack_di_report_status(
Control, status);
1740 " umfpack_di_symbolic() failed!");
1743 status = umfpack_di_numeric(Ap, Ai, Ax, Symbolic, &
Numeric,
1748 umfpack_di_report_status(
Control, status);
1750 " umfpack_di_numeric() failed!");
1752 umfpack_di_free_symbolic(&Symbolic);
1756 SuiteSparse_long status;
1760 AI =
new SuiteSparse_long[
width + 1];
1761 AJ =
new SuiteSparse_long[Ap[
width]];
1762 for (
int i = 0; i <=
width; i++)
1764 AI[i] = (SuiteSparse_long)(Ap[i]);
1766 for (
int i = 0; i < Ap[
width]; i++)
1768 AJ[i] = (SuiteSparse_long)(Ai[i]);
1771 status = umfpack_dl_symbolic(
width,
width,
AI,
AJ, Ax, &Symbolic,
1776 umfpack_dl_report_status(
Control, status);
1778 " umfpack_dl_symbolic() failed!");
1781 status = umfpack_dl_numeric(
AI,
AJ, Ax, Symbolic, &
Numeric,
1786 umfpack_dl_report_status(
Control, status);
1788 " umfpack_dl_numeric() failed!");
1790 umfpack_dl_free_symbolic(&Symbolic);
1797 mfem_error(
"UMFPackSolver::Mult : matrix is not set!"
1798 " Call SetOperator first!");
1805 umfpack_di_report_info(Control,
Info);
1808 umfpack_di_report_status(Control, status);
1809 mfem_error(
"UMFPackSolver::Mult : umfpack_di_solve() failed!");
1814 SuiteSparse_long status =
1817 umfpack_dl_report_info(Control,
Info);
1820 umfpack_dl_report_status(Control, status);
1821 mfem_error(
"UMFPackSolver::Mult : umfpack_dl_solve() failed!");
1829 mfem_error(
"UMFPackSolver::MultTranspose : matrix is not set!"
1830 " Call SetOperator first!");
1837 umfpack_di_report_info(Control,
Info);
1840 umfpack_di_report_status(Control, status);
1842 " umfpack_di_solve() failed!");
1847 SuiteSparse_long status =
1850 umfpack_dl_report_info(Control,
Info);
1853 umfpack_dl_report_status(Control, status);
1855 " umfpack_dl_solve() failed!");
1868 umfpack_di_free_numeric(&
Numeric);
1872 umfpack_dl_free_numeric(&
Numeric);
1887 "Had Numeric pointer in KLU, but not Symbolic");
1895 MFEM_VERIFY(
mat != NULL,
"not a SparseMatrix");
1903 MFEM_VERIFY(
width ==
height,
"not a square matrix");
1915 MFEM_VERIFY(
mat != NULL,
1916 "KLUSolver::Mult : matrix is not set! Call SetOperator first!");
1929 MFEM_VERIFY(
mat != NULL,
1930 "KLUSolver::Mult : matrix is not set! Call SetOperator first!");
1949 #endif // MFEM_USE_SUITESPARSE
void ApplyPlaneRotation(double &dx, double &dy, double &cs, double &sn)
int Size() const
Logical size of the array.
Conjugate gradient method.
virtual void MultTranspose(const Vector &b, Vector &x) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
double Info[UMFPACK_INFO]
int GetNumIterations() const
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
void SetSize(int s)
Resize the vector to size s.
virtual Operator & GetGradient(const Vector &x) const
Evaluate the gradient operator at the point x. The default behavior in class Operator is to generate ...
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
double Norml2() const
Returns the l2 norm of the vector.
int * GetJ()
Return the array J.
Data type dense matrix using column-major storage.
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
int Size() const
Returns the size of the vector.
int * GetI()
Return the array I.
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
virtual void MultTranspose(const Vector &b, Vector &x) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
double * GetData() const
Return a pointer to the beginning of the Vector data.
void SortColumnIndices()
Sort the column indices corresponding to each row.
void add(const Vector &v1, const Vector &v2, Vector &v)
double Norm(const Vector &x) const
double GetFinalNorm() const
double * GetData()
Return the element data, i.e. the array A.
int BiCGSTAB(const Operator &A, Vector &x, const Vector &b, Solver &M, int &max_iter, double &tol, double atol, int printit)
BiCGSTAB method. (tolerances are squared)
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
bool IsFinite(const double &val)
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
void SetLinearConstraint(const Vector &_w, double _a)
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
void SetPrintLevel(int print_lvl)
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
double solve(double l, const Vector &xt, Vector &x, int &nclip) const
Solve QP at fixed lambda.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
void MINRES(const Operator &A, const Vector &b, Vector &x, int print_it, int max_it, double rtol, double atol)
MINRES method without preconditioner. (tolerances are squared)
void Add(const double c, const DenseMatrix &A)
Adds the matrix A multiplied by the number c to the matrix.
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
void SetMaxIter(int max_it)
void SetKDim(int dim)
Set the number of iteration to perform between restarts, default is 50.
void CG(const Operator &A, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Conjugate gradient method. (tolerances are squared)
void print_iteration(int it, double r, double l) const
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
double Dot(const Vector &x, const Vector &y) const
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Stationary linear iteration: x <- x + B (b - A x)
void Swap(Array< T > &, Array< T > &)
void SetAbsTol(double atol)
void SetRelTol(double rtol)
double Control[UMFPACK_CONTROL]
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
void subtract(const Vector &x, const Vector &y, Vector &z)
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
void Update(Vector &x, int k, DenseMatrix &h, Vector &s, Array< Vector * > &v)
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Vector & Set(const double a, const Vector &x)
(*this) = a * x
int height
Dimension of the output / number of rows in the matrix.
void SetBounds(const Vector &_lo, const Vector &_hi)
double InnerProduct(HypreParVector *x, HypreParVector *y)
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const
This method can be overloaded in derived classes to implement line search algorithms.
void GeneratePlaneRotation(double &dx, double &dy, double &cs, double &sn)
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
virtual void Mult(const Vector &xt, Vector &x) const
Operator application: y=A(x).
virtual void SetOperator(const Operator &op)=0
Set/update the solver for the given operator.
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
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, double cf, double &tol, double &atol, int printit)
virtual void SetPreconditioner(Solver &pr)
These are not currently meaningful for this solver and will error out.
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
int width
Dimension of the input / number of columns in the matrix.
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
void SLI(const Operator &A, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Stationary linear iteration. (tolerances are squared)
int GMRES(const Operator &A, Vector &x, const Vector &b, Solver &M, int &max_iter, int m, double &tol, double atol, int printit)
GMRES method. (tolerances are squared)
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.