9 orthogonal(orthogonal),
28 const int N = q.size();
34 std::vector<ColorSpinorField*> Q;
35 for (
int i=0; i<N; i++) Q.push_back(q[i]);
52 for (
int i=0; i<N; i++) {
53 phi(i) = A_[i*(N+1)+N];
54 for (
int j=0; j<N; j++) {
55 A(i,j) = A_[i*(N+1)+j];
64 LDLT<matrix> cholesky(A);
65 psi = cholesky.solve(phi);
70 for (
int i=0; i<N; i++) psi_[i] = psi(i);
86 std::vector<ColorSpinorField*> p, std::vector<ColorSpinorField*> q) {
91 const int N = p.size();
94 printfQuda(
"Constructing minimum residual extrapolation with basis size %d\n", N);
116 for (
int i=0; i<N; i++) {
122 std::vector<ColorSpinorField*> Pi;
124 std::vector<ColorSpinorField*> P;
125 for (
int j=i+1; j<N; j++) P.push_back(p[j]);
127 for (
int j=i+1; j<N; j++) alpha[j] = -alpha[j];
132 std::vector<ColorSpinorField*> X;
134 std::vector<ColorSpinorField*> Y;
135 for (
int j=i+1; j<N; j++) { Y.push_back(q[j]); }
143 if (
apply_mat)
for (
int i=0; i<N; i++)
mat(*q[i], *p[i]);
148 std::vector<ColorSpinorField*> X;
154 for (
int i=0; i<N; i++) alpha[i] = -alpha[i];
155 std::vector<ColorSpinorField*> B;
160 printfQuda(
"MinResExt: N = %d, |res| / |src| = %e\n", N, rsd);
170 std::vector<ColorSpinorField*> p(basis.size()), q(basis.size());
171 for (
unsigned int i=0; i<basis.size(); i++) { p[i] = basis[i].first; q[i] = basis[i].second; }
TimeProfile & profile
whether A is hermitian ot not
bool apply_mat
Whether to construct an orthogonal basis or not.
void operator()(ColorSpinorField &x, ColorSpinorField &b, std::vector< std::pair< ColorSpinorField *, ColorSpinorField * > > basis)
bool hermitian
Whether to compute q = Ap or assume it is provided.
MinResExt(const DiracMatrix &mat, bool orthogonal, bool apply_mat, bool hermitian, TimeProfile &profile)
void solve(Complex *psi_, std::vector< ColorSpinorField * > &p, std::vector< ColorSpinorField * > &q, ColorSpinorField &b, bool hermitian)
Solve the equation A p_k psi_k = q_k psi_k = b by minimizing the residual and using Eigen's SVD algor...
bool isRunning(QudaProfileType idx)
void mat(void *out, void **link, void *in, int dagger_bit, int mu, QudaPrecision sPrecision, QudaPrecision gPrecision)
void ax(double a, ColorSpinorField &x)
void zero(ColorSpinorField &a)
double norm2(const ColorSpinorField &a)
void caxpy(const Complex &a, ColorSpinorField &x, ColorSpinorField &y)
void copy(ColorSpinorField &dst, const ColorSpinorField &src)
Complex cDotProduct(ColorSpinorField &, ColorSpinorField &)
std::complex< double > Complex
__host__ __device__ ValueType sqrt(ValueType x)
QudaVerbosity getVerbosity()