27 for(
int i=1; i<=(2*nsteps); ++i){
28 mat(pr[i], pr[i-1], temp);
33 for(
int i=(2*nsteps+2); i<(4*nsteps+1); ++i){
34 mat(pr[i], pr[i-1], temp);
39 static void print(
const double d[],
int n){
40 for(
int i=0; i<n; ++i){
41 std::cout << d[i] <<
" ";
43 std::cout << std::endl;
47 for(
int i=0; i<n; ++i){
48 std::cout <<
"(" << real(d[i]) <<
"," << imag(d[i]) <<
") ";
50 std::cout << std::endl;
55 static void zero(T d[],
int N){
56 for(
int i=0; i<N; ++i) d[i] = static_cast<T>(0);
60 static void computeGramMatrix(
Complex** G, std::vector<cudaColorSpinorField>& v){
62 const int dim = v.size();
64 for(
int i=0; i<
dim; ++i){
65 for(
int j=0; j<
dim; ++j){
72 static void computeGramVector(
Complex* g, cudaColorSpinorField& r0, std::vector<cudaColorSpinorField>& pr){
74 const int dim = pr.size();
76 for(
int i=0; i<
dim; ++i){
109 static T zip(T a[], T b[],
int dim){
111 for(
int i=0; i<
dim; ++i){
121 for(
int i=0; i<
dim; ++i){
122 for(
int j=0; j<
dim; ++j){
123 result +=
conj(u[i])*v[j]*M[i][j];
139 printfQuda(
"Warning: inverting on zero-field source\n");
177 for(
int i=0; i<(4*s+2); ++i){
187 for(
int i=0; i<(2*s+1); ++i){
205 computeMatrixPowers(PR, p, r, s);
206 computeGramVector(g, r0, PR);
207 computeGramMatrix(G, PR);
210 for(
int i=0; i<(2*s+1); ++i){
213 a[i][i] =
static_cast<Complex>(1);
214 c[i][i + (2*s+1)] =
static_cast<Complex>(1);
223 alpha = zip(g,c[0],4*s+2)/zip(g,a[1],4*s+2);
225 Complex omega_num = computeUdaggerMV(c[0],G,c[1],(4*s+2))
226 - alpha*computeUdaggerMV(c[0],G,a[2],(4*s+2))
227 -
conj(alpha)*computeUdaggerMV(a[1],G,c[1],(4*s+2))
228 +
conj(alpha)*alpha*computeUdaggerMV(a[1],G,a[2],(4*s+2));
230 Complex omega_denom = computeUdaggerMV(c[1],G,c[1],(4*s+2))
231 - alpha*computeUdaggerMV(c[1],G,a[2],(4*s+2))
232 -
conj(alpha)*computeUdaggerMV(a[2],G,c[1],(4*s+2))
233 +
conj(alpha)*alpha*computeUdaggerMV(a[2],G,a[2],(4*s+2));
235 omega = omega_num/omega_denom;
237 for(
int i=0; i<(4*s+2); ++i){
238 e[i] += alpha*a[0][i] +
omega*c[0][i] - alpha*
omega*a[1][i];
242 for(
int k=0; k<=(2*(s - j - 1)); ++k){
243 for(
int i=0; i<(4*s+2); ++i){
244 c_new[k][i] = c[k][i] - alpha*a[k+1][i] -
omega*c[k+1][i] + alpha*
omega*a[k+2][i];
249 beta = (zip(g,c_new[0],(4*s+2))/zip(g,c[0],(4*s+2)))*(alpha/
omega);
251 for(
int k=0; k<=(2*(s - j - 1)); ++k){
252 for(
int i=0; i<(4*s+2); ++i){
253 a_new[k][i] = c_new[k][i] + beta*a[k][i] - beta*
omega*a[k+1][i];
256 for(
int i=0; i<(4*s+2); ++i){
257 a[k][i] = a_new[k][i];
258 c[k][i] = c_new[k][i];
262 for(
int i=0; i<(4*s+2); ++i){
271 for(
int i=0; i<(4*s+2); ++i){
290 for(
int i=0; i<(4*s+2); ++i){
299 for(
int i=0; i<(2*s+1); ++i){
MPBiCGstab(const DiracMatrix &mat, SolverParam ¶m, TimeProfile &profile)
void operator()(ColorSpinorField &out, ColorSpinorField &in)
bool convergence(double r2, double hq2, double r2_tol, double hq_tol)
void PrintSummary(const char *name, int k, double r2, double b2, double r2_tol, double hq_tol)
Prints out the summary of the solver convergence (requires a verbosity of QUDA_SUMMARIZE)....
static double stopping(double tol, double b2, QudaResidualType residual_type)
Set the solver L2 stopping condition.
void PrintStats(const char *name, int k, double r2, double b2, double hq2)
Prints out the running statistics of the solver (requires a verbosity of QUDA_VERBOSE)
void mat(void *out, void **link, void *in, int dagger_bit, int mu, QudaPrecision sPrecision, QudaPrecision gPrecision)
double xmyNorm(ColorSpinorField &x, ColorSpinorField &y)
void zero(ColorSpinorField &a)
double norm2(const ColorSpinorField &a)
void caxpy(const Complex &a, ColorSpinorField &x, ColorSpinorField &y)
Complex cDotProduct(ColorSpinorField &, ColorSpinorField &)
void stop()
Stop profiling.
__host__ __device__ ValueType conj(ValueType x)
__device__ __host__ void zero(double &a)
std::complex< double > Complex
__host__ __device__ ValueType sqrt(ValueType x)
void print(const double d[], int n)
QudaResidualType residual_type