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

Contents of /branches/diaplayground/cusplibrary/cusp/krylov/detail/bicg.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: 4290 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 bicg(LinearOperator& A,
34 LinearOperator& At,
35 Vector& x,
36 Vector& b)
37 {
38 typedef typename LinearOperator::value_type ValueType;
39
40 cusp::default_monitor<ValueType> monitor(b);
41
42 cusp::krylov::bicg(A, At, x, b, monitor);
43 }
44
45 template <class LinearOperator,
46 class Vector,
47 class Monitor>
48 void bicg(LinearOperator& A,
49 LinearOperator& At,
50 Vector& x,
51 Vector& b,
52 Monitor& monitor)
53 {
54 typedef typename LinearOperator::value_type ValueType;
55 typedef typename LinearOperator::memory_space MemorySpace;
56
57 cusp::identity_operator<ValueType,MemorySpace> M(A.num_rows, A.num_cols);
58
59 cusp::krylov::bicg(A, At, x, b, monitor, M, M);
60 }
61
62 template <class LinearOperator,
63 class Vector,
64 class Monitor,
65 class Preconditioner>
66 void bicg(LinearOperator& A,
67 LinearOperator& At,
68 Vector& x,
69 Vector& b,
70 Monitor& monitor,
71 Preconditioner& M,
72 Preconditioner& Mt)
73 {
74 CUSP_PROFILE_SCOPED();
75
76 typedef typename LinearOperator::value_type ValueType;
77 typedef typename LinearOperator::memory_space MemorySpace;
78
79 assert(A.num_rows == A.num_cols); // sanity check
80
81 const size_t N = A.num_rows;
82
83 // allocate workspace
84 cusp::array1d<ValueType,MemorySpace> y(N);
85
86 cusp::array1d<ValueType,MemorySpace> p(N);
87 cusp::array1d<ValueType,MemorySpace> p_star(N);
88 cusp::array1d<ValueType,MemorySpace> q(N);
89 cusp::array1d<ValueType,MemorySpace> q_star(N);
90 cusp::array1d<ValueType,MemorySpace> r(N);
91 cusp::array1d<ValueType,MemorySpace> r_star(N);
92 cusp::array1d<ValueType,MemorySpace> z(N);
93 cusp::array1d<ValueType,MemorySpace> z_star(N);
94
95 // y <- Ax
96 cusp::multiply(A, x, y);
97
98 // r <- b - A*x
99 blas::axpby(b, y, r, ValueType(1), ValueType(-1));
100 if(monitor.finished(r)){
101 return;
102 }
103
104 // r_star <- r
105 blas::copy(r, r_star);
106
107 // z = M r
108 cusp::multiply(M, r, z);
109 // z_star = Mt r_star
110 cusp::multiply(Mt, r_star, z_star);
111
112 // rho = (z,r_star)
113 ValueType rho = blas::dotc(z, r_star);
114
115 // p <- z
116 blas::copy(z, p);
117
118 // p_star <- r
119 blas::copy(z_star, p_star);
120
121 while (1)
122 {
123 // q = A p
124 cusp::multiply(A, p, q);
125 // q_star = At p_star
126 cusp::multiply(At, p_star, q_star);
127
128 // alpha = (rho) / (p_star, q)
129 ValueType alpha = rho / blas::dotc(p_star, q);
130
131 // x += alpha*p
132 blas::axpby(x, p, x, ValueType(1), ValueType(alpha));
133
134 // r -= alpha*q
135 blas::axpby(r, q, r, ValueType(1), ValueType(-alpha));
136
137 // r_star -= alpha*q_star
138 blas::axpby(r_star, q_star, r_star, ValueType(1), ValueType(-alpha));
139
140 if (monitor.finished(r)){
141 break;
142 }
143
144 // z = M r
145 cusp::multiply(M, r, z);
146 // z_star = Mt r_star
147 cusp::multiply(Mt, r_star, z_star);
148
149 ValueType prev_rho = rho;
150 // rho = (z,r_star)
151 rho = blas::dotc(z, r_star);
152 if(rho == ValueType(0)){
153 // Failure!
154 // TODO: Make the failure more apparent to the user
155 break;
156 }
157
158 ValueType beta = rho/prev_rho;
159
160 // p = beta*p + z
161 blas::axpby(p, z, p, ValueType(beta), ValueType(1));
162
163 // p_star = beta*p_star + z_star
164 blas::axpby(p_star, z_star, p_star, ValueType(beta), ValueType(1));
165 ++monitor;
166 }
167 }
168
169 } // end namespace krylov
170 } // end namespace cusp
171

  ViewVC Help
Powered by ViewVC 1.1.26