35 #ifndef MADNESS_LINALG_GMRES_H__INCLUDED
36 #define MADNESS_LINALG_GMRES_H__INCLUDED
69 virtual void action(
const T &in,
T &out)
const = 0;
101 template <
typename T,
typename real_type,
typename scalar_type>
123 virtual real_type
norm(
const T &)
const = 0;
129 virtual T &
scale(
T &,
const scalar_type &)
const = 0;
135 virtual T &
gaxpy(
T &x,
const scalar_type &
a,
const T &y,
136 const scalar_type &
b)
const = 0;
139 virtual scalar_type
inner(
const T &,
const T &)
const = 0;
151 static inline double real(
double x) {
return x; }
156 static inline double imag(
double x) {
return 0.0; }
159 template <
typename T,
int NDIM>
161 typename TensorTypeData<T>::float_scalar_type, T> {
188 for(i = 1; i <
NDIM; ++i) {
212 for(
int i = 0; i <
NDIM; ++i)
213 x[i] =
a * x[i] +
b * y[i];
223 ret =
conj(l[0]) * r[0];
224 for(
int i = 1; i <
NDIM; ++i)
225 ret +=
conj(l[i]) * r[i];
229 for(
int i = 1; i <
NDIM; ++i)
238 template <
typename T,
int NDIM>
240 typename TensorTypeData<T>::float_scalar_type, T> {
285 template <
typename T,
int VDIM,
int FDIM>
287 Vector<Function<T, FDIM>, VDIM>,
288 typename TensorTypeData<T>::float_scalar_type, T> {
306 for(
int i = 1; i < VDIM; ++i) {
307 temp =
vec[i].norm2();
317 for(
int i = 0; i < VDIM; ++i)
328 for(
int i = 0; i < VDIM; ++i)
338 for(
int i = 0; i < VDIM; ++i)
339 ret += l[i].
inner(r[i]);
345 for(
int i = 0; i < VDIM; ++i)
386 template <
typename T,
typename real_type,
typename scalar_type>
389 real_type &resid_thresh, real_type &update_thresh,
390 const bool outp =
false) {
401 real_type resid,
norm, updatenorm;
411 space.
gaxpy(r, -1.0,
b, 1.0);
412 betae[0] = resid = space.
norm(r);
413 if(outp && world.
rank() == 0)
414 printf(
"itr rnk update_norm resid\n%.3d N/A N/A %.6e\n",
416 if(resid < resid_thresh) {
418 resid_thresh = resid;
422 space.
scale(r, 1.0 / resid);
428 op.applyOp(
V[iter - 1], r);
431 for(i = 0; i < iter; ++i) {
432 H(i, iter-1) = space.
inner(
V[i], r);
433 space.
gaxpy(r, 1.0,
V[i], -H(i, iter-1));
435 H(iter, iter-1) =
norm = space.
norm(r);
443 1.0e-12, y, s, rank, sumsq);
454 updatenorm = y.
normf();
456 scalar_type temp = y[0] - yold[0];
458 for(i = 1; i < iter-1; ++i) {
459 temp = y[i] - yold[i];
460 updatenorm +=
real(temp)*
real(temp) +
463 updatenorm +=
real(y[iter-1]) *
real(y[iter-1]) +
465 updatenorm = sqrt(updatenorm);
475 space.
gaxpy(x, 1.0,
V[0], betae[0]);
477 space.
gaxpy(r, -1.0,
b, 1.0);
478 resid_thresh = space.
norm(r);
479 update_thresh = updatenorm;
484 if(outp && world.
rank() == 0)
485 printf(
"%.3d N/A %.6e %.6e ** Zero Vector Encount" \
486 "ered **\n", iter, updatenorm, resid_thresh);
490 if(outp && world.
rank() == 0)
491 printf(
"%.3d %.3ld %.6e %.6e ** Zero Vector Encount" \
492 "ered **\n", iter, rank, updatenorm, resid);
496 if(outp && world.
rank() == 0) {
497 printf(
"%.3d %.3ld %.6e %.6e", iter, rank, updatenorm, resid);
499 printf(
" ** Questionable Progress **");
504 }
while(iter <= maxiters && resid > resid_thresh &&
505 updatenorm > update_thresh);
508 for(i = 0; i < iter; ++i) {
509 space.
gaxpy(x, 1.0,
V[i], y[i]);
513 resid_thresh = resid;
514 update_thresh = updatenorm;
A generic vector space which provides common operations needed by linear algebra routines (norm,...
Definition: gmres.h:102
AbstractVectorSpace()
Bury the default constructor to prevent its use.
AbstractVectorSpace(World &world)
Make a vector space.
Definition: gmres.h:118
virtual T & gaxpy(T &x, const scalar_type &a, const T &y, const scalar_type &b) const =0
Standard bilinear gaxpy.
virtual T & scale(T &, const scalar_type &) const =0
Scales the vector (in-place) by a constant.
virtual ~AbstractVectorSpace()
Definition: gmres.h:120
virtual scalar_type inner(const T &, const T &) const =0
The inner product between two vectors.
World & world
The world.
Definition: gmres.h:109
virtual void destroy(T &) const
Any special instructions to be executed when a vector is no longer needed.
Definition: gmres.h:145
virtual real_type norm(const T &) const =0
The norm of a vector.
A vector space using MADNESS Functions.
Definition: gmres.h:240
virtual scalar_type inner(const Function< scalar_type, NDIM > &l, const Function< scalar_type, NDIM > &r) const
The inner product between two vectors.
Definition: gmres.h:273
virtual real_type norm(const Function< scalar_type, NDIM > &vec) const
The norm of a vector.
Definition: gmres.h:251
virtual Function< scalar_type, NDIM > & gaxpy(Function< scalar_type, NDIM > &x, const scalar_type &a, const Function< scalar_type, NDIM > &y, const scalar_type &b) const
Standard bilinear gaxpy.
Definition: gmres.h:264
virtual Function< scalar_type, NDIM > & scale(Function< scalar_type, NDIM > &vec, const scalar_type &c) const
Scales the vector (in-place) by a constant.
Definition: gmres.h:257
virtual ~FunctionSpace()
Definition: gmres.h:249
FunctionSpace(World &world)
Definition: gmres.h:246
virtual void destroy(Function< scalar_type, NDIM > &f) const
Any special instructions to be executed when a vector is no longer needed.
Definition: gmres.h:279
T scalar_type
Definition: gmres.h:244
TensorTypeData< T >::float_scalar_type real_type
Definition: gmres.h:243
A multiresolution adaptive numerical function.
Definition: mra.h:122
Function< T, NDIM > & gaxpy(const T &alpha, const Function< Q, NDIM > &other, const R &beta, bool fence=true)
Inplace, general bi-linear operation in wavelet basis. No communication except for optional fence.
Definition: mra.h:981
A generic operator: takes in one T and produces another T.
Definition: gmres.h:62
virtual ~Operator()
Definition: gmres.h:83
virtual void action(const T &in, T &out) const =0
The action of the operator.
T & applyOp(const T &in, T &out) const
Public access to the operator's action, returns out for convenience.
Definition: gmres.h:78
A slice defines a sub-range or patch of a dimension.
Definition: slice.h:103
Traits class to specify support of numeric types.
Definition: type_data.h:56
A tensor is a multidimension array.
Definition: tensor.h:317
float_scalar_type normf() const
Returns the Frobenius norm of the tensor.
Definition: tensor.h:1726
A vector space using MADNESS Vectors of MADNESS Functions.
Definition: gmres.h:288
virtual Vector< Function< scalar_type, FDIM >, VDIM > & scale(Vector< Function< scalar_type, FDIM >, VDIM > &vec, const scalar_type &c) const
Scales the vector (in-place) by a constant.
Definition: gmres.h:313
virtual scalar_type inner(const Vector< Function< scalar_type, FDIM >, VDIM > &l, const Vector< Function< scalar_type, FDIM >, VDIM > &r) const
The inner product between two vectors.
Definition: gmres.h:333
TensorTypeData< T >::float_scalar_type real_type
Definition: gmres.h:291
VectorOfFunctionsSpace(World &world)
Definition: gmres.h:294
virtual ~VectorOfFunctionsSpace()
Definition: gmres.h:298
virtual Vector< Function< scalar_type, FDIM >, VDIM > & gaxpy(Vector< Function< scalar_type, FDIM >, VDIM > &x, const scalar_type &a, const Vector< Function< scalar_type, FDIM >, VDIM > &y, const scalar_type &b) const
Standard bilinear gaxpy.
Definition: gmres.h:322
virtual void destroy(Vector< Function< scalar_type, FDIM >, VDIM > &f) const
Any special instructions to be executed when a vector is no longer needed.
Definition: gmres.h:343
T scalar_type
Definition: gmres.h:292
virtual real_type norm(const Vector< Function< scalar_type, FDIM >, VDIM > &vec) const
The norm of a vector.
Definition: gmres.h:300
A vector space using MADNESS Vectors.
Definition: gmres.h:161
virtual Vector< scalar_type, NDIM > & scale(Vector< scalar_type, NDIM > &vec, const scalar_type &c) const
Scales the vector (in-place) by a constant.
Definition: gmres.h:200
TensorTypeData< T >::float_scalar_type real_type
Definition: gmres.h:167
virtual Vector< scalar_type, NDIM > & gaxpy(Vector< scalar_type, NDIM > &x, const scalar_type &a, const Vector< scalar_type, NDIM > &y, const scalar_type &b) const
Standard bilinear gaxpy.
Definition: gmres.h:207
virtual real_type norm(const Vector< scalar_type, NDIM > &vec) const
The norm of a vector.
Definition: gmres.h:175
VectorSpace(World &world)
Definition: gmres.h:170
virtual scalar_type inner(const Vector< scalar_type, NDIM > &l, const Vector< scalar_type, NDIM > &r) const
The inner product between two vectors.
Definition: gmres.h:217
static const bool iscplx
Definition: gmres.h:164
T scalar_type
Definition: gmres.h:168
virtual ~VectorSpace()
Definition: gmres.h:173
A simple, fixed dimension vector.
Definition: vector.h:64
A parallel world class.
Definition: world.h:132
ProcessID rank() const
Returns the process rank in this World (same as MPI_Comm_rank()).
Definition: world.h:318
auto T(World &world, response_space &f) -> response_space
Definition: global_functions.cc:34
Tensor< double > op(const Tensor< double > &x)
Definition: kain.cc:508
double norm(const T &t)
Definition: adquad.h:42
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
Function< T, NDIM > conj(const Function< T, NDIM > &f, bool fence=true)
Return the complex conjugate of the input function with the same distribution and optional fence.
Definition: mra.h:2046
Function< T, NDIM > copy(const Function< T, NDIM > &f, const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > &pmap, bool fence=true)
Create a new copy of the function with different distribution and optional fence.
Definition: mra.h:2002
void gelss(const Tensor< T > &a, const Tensor< T > &b, double rcond, Tensor< T > &x, Tensor< typename Tensor< T >::scalar_type > &s, long &rank, Tensor< typename Tensor< T >::scalar_type > &sumsq)
Solve Ax = b for general A using the LAPACK *gelss routines.
Definition: lapack.cc:884
NDIM & f
Definition: mra.h:2416
double imag(double x)
Definition: complexfun.h:56
double real(double x)
Definition: complexfun.h:52
void GMRES(const AbstractVectorSpace< T, real_type, scalar_type > &space, const Operator< T > &op, const T &b, T &x, int &maxiters, real_type &resid_thresh, real_type &update_thresh, const bool outp=false)
A GMRES solver routine for linear systems, .
Definition: gmres.h:387
Vector< T, sizeof...(Ts)+1 > vec(T t, Ts... ts)
Factory function for creating a madness::Vector.
Definition: vector.h:711
static const double b
Definition: nonlinschro.cc:119
static const double a
Definition: nonlinschro.cc:118
Defines simple templates for printing to std::cout "a la Python".
static const double c
Definition: relops.cc:10
static double V(const coordT &r)
Definition: tdse.cc:288
Defines and implements most of Tensor.
Prototypes for a partial interface from Tensor to LAPACK.
static const std::size_t NDIM
Definition: testpdiff.cc:42