/[escript]/branches/diaplayground/cusplibrary/cusp/krylov/detail/bicgstab.inl
ViewVC logotype

Contents of /branches/diaplayground/cusplibrary/cusp/krylov/detail/bicgstab.inl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4955 - (show annotations)
Tue May 20 04:33:15 2014 UTC (6 years, 4 months ago) by caltinay
File size: 4225 byte(s)
added pristine copy of cusplibrary (apache license) to be used by ripley.

1 /*
2 * Copyright 2008-2009 NVIDIA Corporation
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16
17
18 #include <cusp/array1d.h>
19 #include <cusp/blas.h>
20 #include <cusp/multiply.h>
21 #include <cusp/monitor.h>
22 #include <cusp/linear_operator.h>
23
24 namespace blas = cusp::blas;
25
26 namespace cusp
27 {
28 namespace krylov
29 {
30
31 template <class LinearOperator,
32 class Vector>
33 void bicgstab(LinearOperator& A,
34 Vector& x,
35 Vector& b)
36 {
37 typedef typename LinearOperator::value_type ValueType;
38
39 cusp::default_monitor<ValueType> monitor(b);
40
41 cusp::krylov::bicgstab(A, x, b, monitor);
42 }
43
44 template <class LinearOperator,
45 class Vector,
46 class Monitor>
47 void bicgstab(LinearOperator& A,
48 Vector& x,
49 Vector& b,
50 Monitor& monitor)
51 {
52 typedef typename LinearOperator::value_type ValueType;
53 typedef typename LinearOperator::memory_space MemorySpace;
54
55 cusp::identity_operator<ValueType,MemorySpace> M(A.num_rows, A.num_cols);
56
57 cusp::krylov::bicgstab(A, x, b, monitor, M);
58 }
59
60 template <class LinearOperator,
61 class Vector,
62 class Monitor,
63 class Preconditioner>
64 void bicgstab(LinearOperator& A,
65 Vector& x,
66 Vector& b,
67 Monitor& monitor,
68 Preconditioner& M)
69 {
70 CUSP_PROFILE_SCOPED();
71
72 typedef typename LinearOperator::value_type ValueType;
73 typedef typename LinearOperator::memory_space MemorySpace;
74
75 assert(A.num_rows == A.num_cols); // sanity check
76
77 const size_t N = A.num_rows;
78
79 // allocate workspace
80 cusp::array1d<ValueType,MemorySpace> y(N);
81
82 cusp::array1d<ValueType,MemorySpace> p(N);
83 cusp::array1d<ValueType,MemorySpace> r(N);
84 cusp::array1d<ValueType,MemorySpace> r_star(N);
85 cusp::array1d<ValueType,MemorySpace> s(N);
86 cusp::array1d<ValueType,MemorySpace> Mp(N);
87 cusp::array1d<ValueType,MemorySpace> AMp(N);
88 cusp::array1d<ValueType,MemorySpace> Ms(N);
89 cusp::array1d<ValueType,MemorySpace> AMs(N);
90
91 // y <- Ax
92 cusp::multiply(A, x, y);
93
94 // r <- b - A*x
95 blas::axpby(b, y, r, ValueType(1), ValueType(-1));
96
97 // p <- r
98 blas::copy(r, p);
99
100 // r_star <- r
101 blas::copy(r, r_star);
102
103 ValueType r_r_star_old = blas::dotc(r_star, r);
104
105 while (!monitor.finished(r))
106 {
107 // Mp = M*p
108 cusp::multiply(M, p, Mp);
109
110 // AMp = A*Mp
111 cusp::multiply(A, Mp, AMp);
112
113 // alpha = (r_j, r_star) / (A*M*p, r_star)
114 ValueType alpha = r_r_star_old / blas::dotc(r_star, AMp);
115
116 // s_j = r_j - alpha * AMp
117 blas::axpby(r, AMp, s, ValueType(1), ValueType(-alpha));
118
119 if (monitor.finished(s)){
120 // x += alpha*M*p_j
121 blas::axpby(x, Mp, x, ValueType(1), ValueType(alpha));
122 break;
123 }
124
125 // Ms = M*s_j
126 cusp::multiply(M, s, Ms);
127
128 // AMs = A*Ms
129 cusp::multiply(A, Ms, AMs);
130
131 // omega = (AMs, s) / (AMs, AMs)
132 ValueType omega = blas::dotc(AMs, s) / blas::dotc(AMs, AMs);
133
134 // x_{j+1} = x_j + alpha*M*p_j + omega*M*s_j
135 blas::axpbypcz(x, Mp, Ms, x, ValueType(1), alpha, omega);
136
137 // r_{j+1} = s_j - omega*A*M*s
138 blas::axpby(s, AMs, r, ValueType(1), -omega);
139
140 // beta_j = (r_{j+1}, r_star) / (r_j, r_star) * (alpha/omega)
141 ValueType r_r_star_new = blas::dotc(r_star, r);
142 ValueType beta = (r_r_star_new / r_r_star_old) * (alpha / omega);
143 r_r_star_old = r_r_star_new;
144
145 // p_{j+1} = r_{j+1} + beta*(p_j - omega*A*M*p)
146 blas::axpbypcz(r, p, AMp, p, ValueType(1), beta, -beta*omega);
147
148 ++monitor;
149 }
150 }
151
152 } // end namespace krylov
153 } // end namespace cusp
154

  ViewVC Help
Powered by ViewVC 1.1.26