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

Contents of /branches/diaplayground/cusplibrary/cusp/krylov/detail/cr.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: 3813 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 cr(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::cr(A, x, b, monitor);
42 }
43
44 template <class LinearOperator,
45 class Vector,
46 class Monitor>
47 void cr(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::cr(A, x, b, monitor, M);
58 }
59
60 template <class LinearOperator,
61 class Vector,
62 class Monitor,
63 class Preconditioner>
64 void cr(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 const size_t recompute_r = 8; // interval to update r
79
80 // allocate workspace
81 cusp::array1d<ValueType,MemorySpace> y(N);
82 cusp::array1d<ValueType,MemorySpace> z(N);
83 cusp::array1d<ValueType,MemorySpace> r(N);
84 cusp::array1d<ValueType,MemorySpace> p(N);
85 cusp::array1d<ValueType,MemorySpace> Az(N);
86 cusp::array1d<ValueType,MemorySpace> Ax(N);
87
88 // y <- A*x
89 cusp::multiply(A, x, Ax);
90
91 // r <- b - A*x
92 blas::axpby(b, Ax, r, ValueType(1), ValueType(-1));
93
94 // z <- M*r
95 cusp::multiply(M, r, z);
96
97 // p <- z
98 blas::copy(z, p);
99
100 // y <- A*p
101 cusp::multiply(A, p, y);
102
103 // Az <- A*z
104 cusp::multiply(A, z, Az);
105
106 // rz = <r^H, z>
107 ValueType rz = blas::dotc(r, Az);
108
109 while (!monitor.finished(r))
110 {
111 // alpha <- <r,z>/<y,p>
112 ValueType alpha = rz / blas::dotc(y, y);
113
114 // x <- x + alpha * p
115 blas::axpy(p, x, alpha);
116
117 size_t iter = monitor.iteration_count();
118 if( (iter % recompute_r) && (iter > 0) )
119 {
120 // r <- r - alpha * y
121 blas::axpy(y, r, -alpha);
122 }
123 else
124 {
125 // y <- A*x
126 cusp::multiply(A, x, Ax);
127
128 // r <- b - A*x
129 blas::axpby(b, Ax, r, ValueType(1), ValueType(-1));
130 }
131
132 // z <- M*r
133 cusp::multiply(M, r, z);
134
135 // Az <- A*z
136 cusp::multiply(A, z, Az);
137
138 ValueType rz_old = rz;
139
140 // rz = <r^H, z>
141 rz = blas::dotc(r, Az);
142
143 // beta <- <r_{i+1},r_{i+1}>/<r,r>
144 ValueType beta = rz / rz_old;
145
146 // p <- r + beta*p
147 blas::axpby(z, p, p, ValueType(1), beta);
148
149 // y <- z + beta*p
150 blas::axpby(Az, y, y, ValueType(1), beta);
151
152 ++monitor;
153 }
154 }
155
156 } // end namespace krylov
157 } // end namespace cusp
158

  ViewVC Help
Powered by ViewVC 1.1.26