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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4955 - (hide annotations)
Tue May 20 04:33:15 2014 UTC (6 years, 4 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 caltinay 4955 /*
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