Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziHelperTraits.hpp
1// @HEADER
2// ***********************************************************************
3//
4// Anasazi: Block Eigensolvers Package
5// Copyright 2004 Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef ANASAZI_HELPER_TRAITS_HPP
43#define ANASAZI_HELPER_TRAITS_HPP
44
49#include "AnasaziConfigDefs.hpp"
50#include "AnasaziTypes.hpp"
51#include "Teuchos_LAPACK.hpp"
52
53namespace Anasazi {
54
62 template <class ScalarType>
64 {
65 public:
66
68
71 static void sortRitzValues(
72 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& rRV,
73 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
74 std::vector<Value<ScalarType> >* RV, std::vector<int>* RO, std::vector<int>* RI );
75
77
80 static void scaleRitzVectors(
81 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
82 Teuchos::SerialDenseMatrix<int, ScalarType>* S );
83
85
87 static void computeRitzResiduals(
88 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
89 const Teuchos::SerialDenseMatrix<int, ScalarType>& S,
90 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>* RR);
91
92 };
93
94
95 template<class ScalarType>
97 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& rRV,
98 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
99 std::vector<Value<ScalarType> >* RV, std::vector<int>* RO, std::vector<int>* RI )
100 {
101 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
102 MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
103
104 int curDim = (int)rRV.size();
105 int i = 0;
106
107 // Clear the current index.
108 RI->clear();
109
110 // Place the Ritz values from rRV and iRV into the RV container.
111 while( i < curDim ) {
112 if ( iRV[i] != MT_ZERO ) {
113 //
114 // We will have this situation for real-valued, non-Hermitian matrices.
115 (*RV)[i].set(rRV[i], iRV[i]);
116 (*RV)[i+1].set(rRV[i+1], iRV[i+1]);
117
118 // Make sure that complex conjugate pairs have their positive imaginary part first.
119 if ( (*RV)[i].imagpart < MT_ZERO ) {
120 // The negative imaginary part is first, so swap the order of the ritzValues and ritzOrders.
121 Anasazi::Value<ScalarType> tmp_ritz( (*RV)[i] );
122 (*RV)[i] = (*RV)[i+1];
123 (*RV)[i+1] = tmp_ritz;
124
125 int tmp_order = (*RO)[i];
126 (*RO)[i] = (*RO)[i+1];
127 (*RO)[i+1] = tmp_order;
128
129 }
130 RI->push_back(1); RI->push_back(-1);
131 i = i+2;
132 } else {
133 //
134 // The Ritz value is not complex.
135 (*RV)[i].set(rRV[i], MT_ZERO);
136 RI->push_back(0);
137 i++;
138 }
139 }
140 }
141
142
143 template<class ScalarType>
145 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
146 Teuchos::SerialDenseMatrix<int, ScalarType>* S )
147 {
148 ScalarType ST_ONE = Teuchos::ScalarTraits<ScalarType>::one();
149
150 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
151 MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
152
153 Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
154 Teuchos::BLAS<int,ScalarType> blas;
155
156 int i = 0, curDim = S->numRows();
157 ScalarType temp;
158 ScalarType* s_ptr = S->values();
159 while( i < curDim ) {
160 if ( iRV[i] != MT_ZERO ) {
161 temp = lapack_mag.LAPY2( blas.NRM2( curDim, s_ptr+i*curDim, 1 ),
162 blas.NRM2( curDim, s_ptr+(i+1)*curDim, 1 ) );
163 blas.SCAL( curDim, ST_ONE/temp, s_ptr+i*curDim, 1 );
164 blas.SCAL( curDim, ST_ONE/temp, s_ptr+(i+1)*curDim, 1 );
165 i = i+2;
166 } else {
167 temp = blas.NRM2( curDim, s_ptr+i*curDim, 1 );
168 blas.SCAL( curDim, ST_ONE/temp, s_ptr+i*curDim, 1 );
169 i++;
170 }
171 }
172 }
173
174 template<class ScalarType>
176 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
177 const Teuchos::SerialDenseMatrix<int, ScalarType>& S,
178 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>* RR )
179 {
180 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
181 MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
182
183 Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
184 Teuchos::BLAS<int,ScalarType> blas;
185
186 int i = 0;
187 int s_stride = S.stride();
188 int s_rows = S.numRows();
189 int s_cols = S.numCols();
190 ScalarType* s_ptr = S.values();
191
192 while( i < s_cols ) {
193 if ( iRV[i] != MT_ZERO ) {
194 (*RR)[i] = lapack_mag.LAPY2( blas.NRM2(s_rows, s_ptr + i*s_stride, 1),
195 blas.NRM2(s_rows, s_ptr + (i+1)*s_stride, 1) );
196 (*RR)[i+1] = (*RR)[i];
197 i = i+2;
198 } else {
199 (*RR)[i] = blas.NRM2(s_rows, s_ptr + i*s_stride, 1);
200 i++;
201 }
202 }
203 }
204
205#ifdef HAVE_TEUCHOS_COMPLEX
206 // Partial template specializations for the complex scalar type.
207
215 template <class T>
216 class HelperTraits<std::complex<T> >
217 {
218 public:
219 static void sortRitzValues(
220 const std::vector<T>& rRV,
221 const std::vector<T>& iRV,
222 std::vector<Value<std::complex<T> > >* RV,
223 std::vector<int>* RO, std::vector<int>* RI );
224
225 static void scaleRitzVectors(
226 const std::vector<T>& iRV,
227 Teuchos::SerialDenseMatrix<int, std::complex<T> >* S );
228
229 static void computeRitzResiduals(
230 const std::vector<T>& iRV,
231 const Teuchos::SerialDenseMatrix<int, std::complex<T> >& S,
232 std::vector<T>* RR );
233 };
234
235 template<class T>
236 void HelperTraits<std::complex<T> >::sortRitzValues(
237 const std::vector<T>& rRV,
238 const std::vector<T>& iRV,
239 std::vector<Value<std::complex<T> > >* RV,
240 std::vector<int>* RO, std::vector<int>* RI )
241 {
242 (void)RO;
243 int curDim = (int)rRV.size();
244 int i = 0;
245
246 // Clear the current index.
247 RI->clear();
248
249 // Place the Ritz values from rRV and iRV into the RV container.
250 while( i < curDim ) {
251 (*RV)[i].set(rRV[i], iRV[i]);
252 RI->push_back(0);
253 i++;
254 }
255 }
256
257 template<class T>
258 void HelperTraits<std::complex<T> >::scaleRitzVectors(
259 const std::vector<T>& iRV,
260 Teuchos::SerialDenseMatrix<int, std::complex<T> >* S )
261 {
262 (void)iRV;
263 typedef std::complex<T> ST;
264 ST ST_ONE = Teuchos::ScalarTraits<ST>::one();
265
266 Teuchos::BLAS<int,ST> blas;
267
268 int i = 0, curDim = S->numRows();
269 ST temp;
270 ST* s_ptr = S->values();
271 while( i < curDim ) {
272 temp = blas.NRM2( curDim, s_ptr+i*curDim, 1 );
273 blas.SCAL( curDim, ST_ONE/temp, s_ptr+i*curDim, 1 );
274 i++;
275 }
276 }
277
278 template<class T>
279 void HelperTraits<std::complex<T> >::computeRitzResiduals(
280 const std::vector<T>& iRV,
281 const Teuchos::SerialDenseMatrix<int, std::complex<T> >& S,
282 std::vector<T>* RR )
283 {
284 (void)iRV;
285 Teuchos::BLAS<int,std::complex<T> > blas;
286
287 int s_stride = S.stride();
288 int s_rows = S.numRows();
289 int s_cols = S.numCols();
290 std::complex<T>* s_ptr = S.values();
291
292 for (int i=0; i<s_cols; ++i ) {
293 (*RR)[i] = blas.NRM2(s_rows, s_ptr + i*s_stride, 1);
294 }
295 }
296#endif
297
298} // end Anasazi namespace
299
300
301#endif // ANASAZI_HELPER_TRAITS_HPP
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Types and exceptions used within Anasazi solvers and interfaces.
Class which defines basic traits for working with different scalar types.
static void computeRitzResiduals(const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &iRV, const Teuchos::SerialDenseMatrix< int, ScalarType > &S, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > *RR)
Helper function for correctly computing the Ritz residuals of the projected eigenproblem.
static void scaleRitzVectors(const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &iRV, Teuchos::SerialDenseMatrix< int, ScalarType > *S)
Helper function for correctly scaling the eigenvectors of the projected eigenproblem.
static void sortRitzValues(const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &rRV, const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &iRV, std::vector< Value< ScalarType > > *RV, std::vector< int > *RO, std::vector< int > *RI)
Helper function for correctly storing the Ritz values when the eigenproblem is non-Hermitian.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
This struct is used for storing eigenvalues and Ritz values, as a pair of real values.