5 #ifndef SPIKE_BICGSTAB_2_H
6 #define SPIKE_BICGSTAB_2_H
10 #include <cusp/blas.h>
11 #include <cusp/print.h>
12 #include <cusp/array1d.h>
20 typedef typename cusp::array1d<int, cusp::host_memory>
IntVectorH;
38 template <
typename SolverVector,
typename PrecVector,
typename IntVector>
42 IntVector& compIndices,
43 IntVector& comp_perms,
44 std::vector<IntVector>& comp_reorderings)
46 int numComponents = comp_reorderings.size();
48 for (
int i = 0; i < numComponents; i++) {
49 int loc_n = comp_reorderings[i].size();
51 PrecVector buffer_rhs(loc_n);
52 PrecVector buffer_sol(loc_n);
54 thrust::scatter_if(rhs.begin(), rhs.end(), comp_perms.begin(), compIndices.begin(), buffer_rhs.begin(),
IsEqualTo<int>(i));
55 precond_pointers[i]->solve(buffer_rhs, buffer_sol);
56 thrust::scatter(buffer_sol.begin(), buffer_sol.end(), comp_reorderings[i].begin(), sol.begin());
68 template <
typename SpmvOperator,
typename SolverVector,
typename PrecVector,
int L>
70 const SolverVector& b,
76 std::vector<IntVectorH>& comp_reorderings)
78 typedef typename SolverVector::value_type SolverValueType;
79 typedef typename SolverVector::memory_space MemorySpace;
81 typedef typename cusp::array1d<int, MemorySpace> IntVector;
86 SolverValueType rou0 = SolverValueType(1);
87 SolverValueType alpha = SolverValueType(0);
88 SolverValueType omega = SolverValueType(1);
97 IntVector loc_compIndices = compIndices;
98 IntVector loc_comp_perms = comp_perms;
99 std::vector<IntVector> loc_comp_reorderings;
101 int numComponents = comp_reorderings.size();
103 for (
int i = 0; i < numComponents; i++)
104 loc_comp_reorderings.push_back(comp_reorderings[i]);
106 std::vector<SolverVector> rr(L+1);
107 std::vector<SolverVector> uu(L+1);
109 for(
int k = 0; k <= L; k++) {
114 SolverValueType tao[L+1][L+1];
115 SolverValueType gamma[L+2];
116 SolverValueType gamma_prime[L+2];
117 SolverValueType gamma_primeprime[L+2];
118 SolverValueType sigma[L+2];
122 cusp::blas::axpby(b, r0, r0, SolverValueType(1), SolverValueType(-1));
125 cusp::blas::copy(r0, r);
130 thrust::copy(thrust::make_zip_iterator(thrust::make_tuple(u.begin(), x.begin(), r.begin())),
131 thrust::make_zip_iterator(thrust::make_tuple(u.end(), x.end(), r.end())),
132 thrust::make_zip_iterator(thrust::make_tuple(uu[0].begin(), xx.begin(), rr[0].begin())));
134 while(!monitor.
done(r)) {
136 rou0 = -omega * rou0;
140 for(
int j = 0; j < L; j++) {
141 rou1 = cusp::blas::dotc(rr[j], r0);
147 SolverValueType beta = alpha * rou1 / rou0;
150 for(
int i = 0; i <= j; i++) {
152 cusp::blas::axpby(rr[i], uu[i], uu[i], SolverValueType(1), -beta);
157 precondSolveWrapper(uu[j], Pv, precond_pointers, loc_compIndices, loc_comp_perms, loc_comp_reorderings);
161 SolverValueType gamma = cusp::blas::dotc(uu[j+1], r0);
165 alpha = rou0 / gamma;
167 for(
int i = 0; i <= j; i++) {
169 cusp::blas::axpy(uu[i+1], rr[i], SolverValueType(-alpha));
174 precondSolveWrapper(rr[j], Pv, precond_pointers, loc_compIndices, loc_comp_perms, loc_comp_reorderings);
178 cusp::blas::axpy(uu[0], xx, alpha);
180 if(monitor.
done(rr[0])) {
182 precondSolveWrapper(xx, x, precond_pointers, loc_compIndices, loc_comp_perms, loc_comp_reorderings);
188 for(
int j = 1; j <= L; j++) {
189 for(
int i = 1; i < j; i++) {
190 tao[i][j] = cusp::blas::dotc(rr[j], rr[i]) / sigma[i];
191 cusp::blas::axpy(rr[i], rr[j], -tao[i][j]);
193 sigma[j] = cusp::blas::dotc(rr[j], rr[j]);
196 gamma_prime[j] = cusp::blas::dotc(rr[j], rr[0]) / sigma[j];
199 gamma[L] = gamma_prime[L];
202 for(
int j = L-1; j > 0; j--) {
203 gamma[j] = gamma_prime[j];
204 for(
int i = j+1; i <= L; i++)
205 gamma[j] -= tao[j][i] * gamma[i];
208 for(
int j = 1; j < L; j++) {
209 gamma_primeprime[j] = gamma[j+1];
210 for(
int i = j+1; i < L; i++)
211 gamma_primeprime[j] += tao[j][i] * gamma[i+1];
217 cusp::blas::axpy(rr[0], xx, gamma[1]);
218 cusp::blas::axpy(rr[L], rr[0], -gamma_prime[L]);
219 cusp::blas::axpy(uu[L], uu[0], -gamma[L]);
223 if (monitor.
done(rr[0])) {
225 precondSolveWrapper(xx, x, precond_pointers, loc_compIndices, loc_comp_perms, loc_comp_reorderings);
234 for(
int j = 1; j < L; j++) {
235 cusp::blas::axpy(uu[j], uu[0], -gamma[j]);
236 cusp::blas::axpy(rr[j], xx, gamma_primeprime[j]);
237 cusp::blas::axpy(rr[j], rr[0], -gamma_prime[j]);
239 if (monitor.
done(rr[0])) {
241 precondSolveWrapper(xx, x, precond_pointers, loc_compIndices, loc_comp_perms, loc_comp_reorderings);
249 thrust::copy(thrust::make_zip_iterator(thrust::make_tuple(uu[0].begin(), xx.begin(), rr[0].begin())),
250 thrust::make_zip_iterator(thrust::make_tuple(uu[0].end(), xx.end(), rr[0].end())),
251 thrust::make_zip_iterator(thrust::make_tuple(u.begin(), x.begin(), r.begin())));
259 template <
typename SpmvOperator,
typename SolverVector,
typename PrecVector>
261 const SolverVector& b,
267 std::vector<IntVectorH>& comp_reorderings)
269 bicgstabl<SpmvOperator, SolverVector, PrecVector, 2>(spmv, b, x, monitor, precond_pointers, compIndices, comp_perms, comp_reorderings);
274 template <
typename SpmvOperator,
typename SolverVector,
typename PrecVector>
276 const SolverVector& b,
282 std::vector<IntVectorH>& comp_reorderings)
284 bicgstabl<SpmvOperator, SolverVector, PrecVector, 4>(spmv, b, x, monitor, precond_pointers, compIndices, comp_perms, comp_reorderings);
void bicgstabl(SpmvOperator &spmv, const SolverVector &b, SolverVector &x, Monitor< SolverVector > &monitor, std::vector< Precond< PrecVector > * > &precond_pointers, IntVectorH &compIndices, IntVectorH &comp_perms, std::vector< IntVectorH > &comp_reorderings)
Preconditioned BiCGStab(L) Krylov method.
Definition: bicgstab2.h:69
void increment(float incr)
Definition: monitor.h:38
void precondSolveWrapper(SolverVector &rhs, SolverVector &sol, std::vector< Precond< PrecVector > * > &precond_pointers, IntVector &compIndices, IntVector &comp_perms, std::vector< IntVector > &comp_reorderings)
Definition: bicgstab2.h:39
IsEqualTo(T val=0)
Definition: bicgstab2.h:28
Definition: bicgstab2.h:24
bool done(const SolverVector &r)
Definition: monitor.h:89
__host__ __device__ bool operator()(const T &val)
Definition: bicgstab2.h:31
T m_val
Definition: bicgstab2.h:26
Spike preconditioner.
Definition: precond.h:37
void bicgstab4(SpmvOperator &spmv, const SolverVector &b, SolverVector &x, Monitor< SolverVector > &monitor, std::vector< Precond< PrecVector > * > &precond_pointers, IntVectorH &compIndices, IntVectorH &comp_perms, std::vector< IntVectorH > &comp_reorderings)
Specializations of the generic spike::bicgstabl function for L=4.
Definition: bicgstab2.h:275
cusp::array1d< int, cusp::host_memory > IntVectorH
Definition: bicgstab2.h:20
Definition of the Spike preconditioner class.
void bicgstab2(SpmvOperator &spmv, const SolverVector &b, SolverVector &x, Monitor< SolverVector > &monitor, std::vector< Precond< PrecVector > * > &precond_pointers, IntVectorH &compIndices, IntVectorH &comp_perms, std::vector< IntVectorH > &comp_reorderings)
Specializations of the generic spike::bicgstabl function for L=2.
Definition: bicgstab2.h:260