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

Contents of /branches/diaplayground/cusplibrary/cusp/krylov/detail/cg.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: 3199 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 cg(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::cg(A, x, b, monitor);
42 }
43
44 template <class LinearOperator,
45 class Vector,
46 class Monitor>
47 void cg(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::cg(A, x, b, monitor, M);
58 }
59
60 template <class LinearOperator,
61 class Vector,
62 class Monitor,
63 class Preconditioner>
64 void cg(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 cusp::array1d<ValueType,MemorySpace> z(N);
82 cusp::array1d<ValueType,MemorySpace> r(N);
83 cusp::array1d<ValueType,MemorySpace> p(N);
84
85 // y <- Ax
86 cusp::multiply(A, x, y);
87
88 // r <- b - A*x
89 blas::axpby(b, y, r, ValueType(1), ValueType(-1));
90
91 // z <- M*r
92 cusp::multiply(M, r, z);
93
94 // p <- z
95 blas::copy(z, p);
96
97 // rz = <r^H, z>
98 ValueType rz = blas::dotc(r, z);
99
100 while (!monitor.finished(r))
101 {
102 // y <- Ap
103 cusp::multiply(A, p, y);
104
105 // alpha <- <r,z>/<y,p>
106 ValueType alpha = rz / blas::dotc(y, p);
107
108 // x <- x + alpha * p
109 blas::axpy(p, x, alpha);
110
111 // r <- r - alpha * y
112 blas::axpy(y, r, -alpha);
113
114 // z <- M*r
115 cusp::multiply(M, r, z);
116
117 ValueType rz_old = rz;
118
119 // rz = <r^H, z>
120 rz = blas::dotc(r, z);
121
122 // beta <- <r_{i+1},r_{i+1}>/<r,r>
123 ValueType beta = rz / rz_old;
124
125 // p <- r + beta*p
126 blas::axpby(z, p, p, ValueType(1), beta);
127
128 ++monitor;
129 }
130 }
131
132 } // end namespace krylov
133 } // end namespace cusp
134

  ViewVC Help
Powered by ViewVC 1.1.26