1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
|
/*
* Copyright 2014-2015 The University of Queensland
* http://www.uq.edu.au
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#include <cusp/array1d.h>
#include <cusp/blas.h>
#include <cusp/multiply.h>
#include <cusp/monitor.h>
#include <cusp/linear_operator.h>
namespace blas = cusp::blas;
namespace cusp
{
namespace krylov
{
template <class LinearOperator,
typename ValueType,
class Vector>
void cgls(LinearOperator& A,
Vector& x,
Vector& b,
ValueType shift)
{
cusp::default_monitor<ValueType> monitor(b);
cusp::krylov::cgls(A, x, b, shift, monitor);
}
template <class LinearOperator,
typename ValueType,
class Vector,
class Monitor>
void cgls(LinearOperator& A,
Vector& x,
Vector& b,
ValueType shift,
Monitor& monitor)
{
CUSP_PROFILE_SCOPED();
typedef typename LinearOperator::memory_space MemorySpace;
const size_t N = A.num_rows;
// allocate workspace
cusp::array1d<ValueType,MemorySpace> y(N);
cusp::array1d<ValueType,MemorySpace> z(N);
cusp::array1d<ValueType,MemorySpace> r(N);
cusp::array1d<ValueType,MemorySpace> p(N);
// y <- Ax
cusp::multiply(A, x, y);
// r <- b - A*x
blas::axpby(b, y, r, ValueType(1), ValueType(-1));
// z <- A^T*r - shift*x
cusp::transposed_multiply(A, r, z);
blas::axpy(x, z, -shift);
// p <- z
blas::copy(z, p);
// gamma = <r, r>
ValueType gamma = blas::dotc(z, z);
while (!monitor.finished(z))
{
// y <- Ap
cusp::multiply(A, p, y);
// delta <- <y,y> + shift*<p,p>
ValueType delta = blas::dotc(y,y) + shift*blas::dotc(p,p);
// alpha <- <r,r> / <y,y>
ValueType alpha = gamma / delta;
// x <- x + alpha * p
blas::axpy(p, x, alpha);
// r <- r - alpha * y
blas::axpy(y, r, -alpha);
// z <- A^T*r - shift*x
cusp::transposed_multiply(A, r, z);
blas::axpy(x, z, -shift);
ValueType gamma_old = gamma;
// gamma = <r, r>
gamma = blas::dotc(z, z);
// beta <- <r_{i+1},r_{i+1}> / <r,r>
ValueType beta = gamma / gamma_old;
// p <- r + beta*p
blas::axpby(z, p, p, ValueType(1), beta);
++monitor;
}
}
} // end namespace krylov
} // end namespace cusp
|
