/[escript]/branches/diaplayground/cusplibrary/cusp/krylov/arnoldi.h
ViewVC logotype

Contents of /branches/diaplayground/cusplibrary/cusp/krylov/arnoldi.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4955 - (show annotations)
Tue May 20 04:33:15 2014 UTC (6 years, 2 months ago) by caltinay
File MIME type: text/plain
File size: 3298 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 #pragma once
19
20 #include <cusp/detail/config.h>
21
22 #include <cusp/multiply.h>
23 #include <cusp/array1d.h>
24 #include <cusp/detail/random.h>
25
26 namespace cusp
27 {
28 namespace krylov
29 {
30
31 template <typename Matrix, typename Array2d>
32 void lanczos(const Matrix& A, Array2d& H, size_t k = 10)
33 {
34 typedef typename Matrix::value_type ValueType;
35 typedef typename Matrix::memory_space MemorySpace;
36
37 size_t N = A.num_cols;
38 size_t maxiter = std::min(N, k);
39
40 // allocate workspace
41 cusp::array1d<ValueType,MemorySpace> v0(N);
42 cusp::array1d<ValueType,MemorySpace> v1(N);
43 cusp::array1d<ValueType,MemorySpace> w(N);
44
45 // initialize starting vector to random values in [0,1)
46 cusp::copy(cusp::detail::random_reals<ValueType>(N), v1);
47
48 cusp::blas::scal(v1, ValueType(1) / cusp::blas::nrm2(v1));
49
50 Array2d H_(maxiter + 1, maxiter, 0);
51
52 ValueType alpha = 0.0, beta = 0.0;
53
54 size_t j;
55
56 for(j = 0; j < maxiter; j++)
57 {
58 cusp::multiply(A, v1, w);
59
60 if(j >= 1)
61 {
62 H_(j - 1, j) = beta;
63 cusp::blas::axpy(v0, w, -beta);
64 }
65
66 alpha = cusp::blas::dot(w, v1);
67 H_(j,j) = alpha;
68
69 cusp::blas::axpy(v1, w, -alpha);
70
71 beta = cusp::blas::nrm2(w);
72 H_(j + 1, j) = beta;
73
74 if(beta < 1e-10) break;
75
76 cusp::blas::scal(w, ValueType(1) / beta);
77
78 // [v0 v1 w] - > [v1 w v0]
79 v0.swap(v1);
80 v1.swap(w);
81 }
82
83 H.resize(j,j);
84 for(size_t row = 0; row < j; row++)
85 for(size_t col = 0; col < j; col++)
86 H(row,col) = H_(row,col);
87 }
88
89 template <typename Matrix, typename Array2d>
90 void arnoldi(const Matrix& A, Array2d& H, size_t k = 10)
91 {
92 typedef typename Matrix::value_type ValueType;
93 typedef typename Matrix::memory_space MemorySpace;
94
95 size_t N = A.num_rows;
96
97 size_t maxiter = std::min(N, k);
98
99 Array2d H_(maxiter + 1, maxiter, 0);
100
101 // allocate workspace of k + 1 vectors
102 std::vector< cusp::array1d<ValueType,MemorySpace> > V(maxiter + 1);
103 for (size_t i = 0; i < maxiter + 1; i++)
104 V[i].resize(N);
105
106 // initialize starting vector to random values in [0,1)
107 cusp::copy(cusp::detail::random_reals<ValueType>(N), V[0]);
108
109 // normalize v0
110 cusp::blas::scal(V[0], ValueType(1) / cusp::blas::nrm2(V[0]));
111
112 size_t j;
113
114 for(j = 0; j < maxiter; j++)
115 {
116 cusp::multiply(A, V[j], V[j + 1]);
117
118 for(size_t i = 0; i <= j; i++)
119 {
120 H_(i,j) = cusp::blas::dot(V[i], V[j + 1]);
121
122 cusp::blas::axpy(V[i], V[j + 1], -H_(i,j));
123 }
124
125 H_(j+1,j) = cusp::blas::nrm2(V[j + 1]);
126
127 if(H_(j+1,j) < 1e-10) break;
128
129 cusp::blas::scal(V[j + 1], ValueType(1) / H_(j+1,j));
130 }
131
132 H.resize(j,j);
133 for( size_t row = 0; row < j; row++ )
134 for( size_t col = 0; col < j; col++ )
135 H(row,col) = H_(row,col);
136 }
137
138 } // end namespace krylov
139 } // end namespace cusp
140

  ViewVC Help
Powered by ViewVC 1.1.26