Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziHelperTraits.hpp
1// @HEADER
2// *****************************************************************************
3// Anasazi: Block Eigensolvers Package
4//
5// Copyright 2004 NTESS and the Anasazi contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef ANASAZI_HELPER_TRAITS_HPP
11#define ANASAZI_HELPER_TRAITS_HPP
12
17#include "AnasaziConfigDefs.hpp"
18#include "AnasaziTypes.hpp"
19#include "Teuchos_LAPACK.hpp"
20
21namespace Anasazi {
22
30 template <class ScalarType>
31 class HelperTraits
32 {
33 public:
34
36
39 static void sortRitzValues(
40 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& rRV,
41 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
42 std::vector<Value<ScalarType> >* RV, std::vector<int>* RO, std::vector<int>* RI );
43
45
48 static void scaleRitzVectors(
49 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
50 Teuchos::SerialDenseMatrix<int, ScalarType>* S );
51
53
55 static void computeRitzResiduals(
56 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
57 const Teuchos::SerialDenseMatrix<int, ScalarType>& S,
58 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>* RR);
59
60 };
61
62
63 template<class ScalarType>
64 void HelperTraits<ScalarType>::sortRitzValues(
65 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& rRV,
66 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
67 std::vector<Value<ScalarType> >* RV, std::vector<int>* RO, std::vector<int>* RI )
68 {
69 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
70 MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
71
72 int curDim = (int)rRV.size();
73 int i = 0;
74
75 // Clear the current index.
76 RI->clear();
77
78 // Place the Ritz values from rRV and iRV into the RV container.
79 while( i < curDim ) {
80 if ( iRV[i] != MT_ZERO ) {
81 //
82 // We will have this situation for real-valued, non-Hermitian matrices.
83 (*RV)[i].set(rRV[i], iRV[i]);
84 (*RV)[i+1].set(rRV[i+1], iRV[i+1]);
85
86 // Make sure that complex conjugate pairs have their positive imaginary part first.
87 if ( (*RV)[i].imagpart < MT_ZERO ) {
88 // The negative imaginary part is first, so swap the order of the ritzValues and ritzOrders.
89 Anasazi::Value<ScalarType> tmp_ritz( (*RV)[i] );
90 (*RV)[i] = (*RV)[i+1];
91 (*RV)[i+1] = tmp_ritz;
92
93 int tmp_order = (*RO)[i];
94 (*RO)[i] = (*RO)[i+1];
95 (*RO)[i+1] = tmp_order;
96
97 }
98 RI->push_back(1); RI->push_back(-1);
99 i = i+2;
100 } else {
101 //
102 // The Ritz value is not complex.
103 (*RV)[i].set(rRV[i], MT_ZERO);
104 RI->push_back(0);
105 i++;
106 }
107 }
108 }
109
110
111 template<class ScalarType>
112 void HelperTraits<ScalarType>::scaleRitzVectors(
113 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
114 Teuchos::SerialDenseMatrix<int, ScalarType>* S )
115 {
116 ScalarType ST_ONE = Teuchos::ScalarTraits<ScalarType>::one();
117
118 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
119 MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
120
121 Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
122 Teuchos::BLAS<int,ScalarType> blas;
123
124 int i = 0, curDim = S->numRows();
125 ScalarType temp;
126 ScalarType* s_ptr = S->values();
127 while( i < curDim ) {
128 if ( iRV[i] != MT_ZERO ) {
129 temp = lapack_mag.LAPY2( blas.NRM2( curDim, s_ptr+i*curDim, 1 ),
130 blas.NRM2( curDim, s_ptr+(i+1)*curDim, 1 ) );
131 blas.SCAL( curDim, ST_ONE/temp, s_ptr+i*curDim, 1 );
132 blas.SCAL( curDim, ST_ONE/temp, s_ptr+(i+1)*curDim, 1 );
133 i = i+2;
134 } else {
135 temp = blas.NRM2( curDim, s_ptr+i*curDim, 1 );
136 blas.SCAL( curDim, ST_ONE/temp, s_ptr+i*curDim, 1 );
137 i++;
138 }
139 }
140 }
141
142 template<class ScalarType>
143 void HelperTraits<ScalarType>::computeRitzResiduals(
144 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
145 const Teuchos::SerialDenseMatrix<int, ScalarType>& S,
146 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>* RR )
147 {
148 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
149 MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
150
151 Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
152 Teuchos::BLAS<int,ScalarType> blas;
153
154 int i = 0;
155 int s_stride = S.stride();
156 int s_rows = S.numRows();
157 int s_cols = S.numCols();
158 ScalarType* s_ptr = S.values();
159
160 while( i < s_cols ) {
161 if ( iRV[i] != MT_ZERO ) {
162 (*RR)[i] = lapack_mag.LAPY2( blas.NRM2(s_rows, s_ptr + i*s_stride, 1),
163 blas.NRM2(s_rows, s_ptr + (i+1)*s_stride, 1) );
164 (*RR)[i+1] = (*RR)[i];
165 i = i+2;
166 } else {
167 (*RR)[i] = blas.NRM2(s_rows, s_ptr + i*s_stride, 1);
168 i++;
169 }
170 }
171 }
172
173#ifdef HAVE_TEUCHOS_COMPLEX
174 // Partial template specializations for the complex scalar type.
175
183 template <class T>
184 class HelperTraits<std::complex<T> >
185 {
186 public:
187 static void sortRitzValues(
188 const std::vector<T>& rRV,
189 const std::vector<T>& iRV,
190 std::vector<Value<std::complex<T> > >* RV,
191 std::vector<int>* RO, std::vector<int>* RI );
192
193 static void scaleRitzVectors(
194 const std::vector<T>& iRV,
195 Teuchos::SerialDenseMatrix<int, std::complex<T> >* S );
196
197 static void computeRitzResiduals(
198 const std::vector<T>& iRV,
199 const Teuchos::SerialDenseMatrix<int, std::complex<T> >& S,
200 std::vector<T>* RR );
201 };
202
203 template<class T>
204 void HelperTraits<std::complex<T> >::sortRitzValues(
205 const std::vector<T>& rRV,
206 const std::vector<T>& iRV,
207 std::vector<Value<std::complex<T> > >* RV,
208 std::vector<int>* RO, std::vector<int>* RI )
209 {
210 (void)RO;
211 int curDim = (int)rRV.size();
212 int i = 0;
213
214 // Clear the current index.
215 RI->clear();
216
217 // Place the Ritz values from rRV and iRV into the RV container.
218 while( i < curDim ) {
219 (*RV)[i].set(rRV[i], iRV[i]);
220 RI->push_back(0);
221 i++;
222 }
223 }
224
225 template<class T>
226 void HelperTraits<std::complex<T> >::scaleRitzVectors(
227 const std::vector<T>& iRV,
228 Teuchos::SerialDenseMatrix<int, std::complex<T> >* S )
229 {
230 (void)iRV;
231 typedef std::complex<T> ST;
232 ST ST_ONE = Teuchos::ScalarTraits<ST>::one();
233
234 Teuchos::BLAS<int,ST> blas;
235
236 int i = 0, curDim = S->numRows();
237 ST temp;
238 ST* s_ptr = S->values();
239 while( i < curDim ) {
240 temp = blas.NRM2( curDim, s_ptr+i*curDim, 1 );
241 blas.SCAL( curDim, ST_ONE/temp, s_ptr+i*curDim, 1 );
242 i++;
243 }
244 }
245
246 template<class T>
247 void HelperTraits<std::complex<T> >::computeRitzResiduals(
248 const std::vector<T>& iRV,
249 const Teuchos::SerialDenseMatrix<int, std::complex<T> >& S,
250 std::vector<T>* RR )
251 {
252 (void)iRV;
253 Teuchos::BLAS<int,std::complex<T> > blas;
254
255 int s_stride = S.stride();
256 int s_rows = S.numRows();
257 int s_cols = S.numCols();
258 std::complex<T>* s_ptr = S.values();
259
260 for (int i=0; i<s_cols; ++i ) {
261 (*RR)[i] = blas.NRM2(s_rows, s_ptr + i*s_stride, 1);
262 }
263 }
264#endif
265
266} // end Anasazi namespace
267
268
269#endif // ANASAZI_HELPER_TRAITS_HPP
AnasaziConfigDefs.hpp
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
AnasaziTypes.hpp
Types and exceptions used within Anasazi solvers and interfaces.
Anasazi::HelperTraits
Class which defines basic traits for working with different scalar types.
Definition AnasaziHelperTraits.hpp:32
Anasazi::HelperTraits::computeRitzResiduals
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.
Definition AnasaziHelperTraits.hpp:143
Anasazi::HelperTraits::scaleRitzVectors
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.
Definition AnasaziHelperTraits.hpp:112
Anasazi::HelperTraits::sortRitzValues
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.
Definition AnasaziHelperTraits.hpp:64
Anasazi::Operator
Anasazi's templated virtual class for constructing an operator that can interface with the OperatorTr...
Definition AnasaziOperator.hpp:35
Anasazi::Operator::Operator
Operator()
Default constructor.
Definition AnasaziOperator.hpp:40
Anasazi
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
Anasazi::Value
This struct is used for storing eigenvalues and Ritz values, as a pair of real values.
Definition AnasaziTypes.hpp:42
Generated by doxygen 1.9.8