42 #ifndef ANASAZI_HELPER_TRAITS_HPP 43 #define ANASAZI_HELPER_TRAITS_HPP 51 #include "Teuchos_LAPACK.hpp" 62 template <
class ScalarType>
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 );
81 const std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
82 Teuchos::SerialDenseMatrix<int, ScalarType>* S );
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);
95 template<
class ScalarType>
97 const std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& rRV,
98 const std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
101 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
102 MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
104 int curDim = (int)rRV.size();
111 while( i < curDim ) {
112 if ( iRV[i] != MT_ZERO ) {
115 (*RV)[i].set(rRV[i], iRV[i]);
116 (*RV)[i+1].set(rRV[i+1], iRV[i+1]);
119 if ( (*RV)[i].imagpart < MT_ZERO ) {
122 (*RV)[i] = (*RV)[i+1];
123 (*RV)[i+1] = tmp_ritz;
125 int tmp_order = (*RO)[i];
126 (*RO)[i] = (*RO)[i+1];
127 (*RO)[i+1] = tmp_order;
130 RI->push_back(1); RI->push_back(-1);
135 (*RV)[i].set(rRV[i], MT_ZERO);
143 template<
class ScalarType>
145 const std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
146 Teuchos::SerialDenseMatrix<int, ScalarType>* S )
148 ScalarType ST_ONE = Teuchos::ScalarTraits<ScalarType>::one();
150 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
151 MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
153 Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
154 Teuchos::BLAS<int,ScalarType> blas;
156 int i = 0, curDim = S->numRows();
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 );
167 temp = blas.NRM2( curDim, s_ptr+i*curDim, 1 );
168 blas.SCAL( curDim, ST_ONE/temp, s_ptr+i*curDim, 1 );
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 )
180 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
181 MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
183 Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
184 Teuchos::BLAS<int,ScalarType> blas;
187 int s_stride = S.stride();
188 int s_rows = S.numRows();
189 int s_cols = S.numCols();
190 ScalarType* s_ptr = S.values();
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];
199 (*RR)[i] = blas.NRM2(s_rows, s_ptr + i*s_stride, 1);
205 #ifdef HAVE_TEUCHOS_COMPLEX 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 );
226 const std::vector<T>& iRV,
227 Teuchos::SerialDenseMatrix<
int, std::complex<T> >* S );
230 const std::vector<T>& iRV,
231 const Teuchos::SerialDenseMatrix<
int, std::complex<T> >& S,
232 std::vector<T>* RR );
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 )
243 int curDim = (int)rRV.size();
250 while( i < curDim ) {
251 (*RV)[i].set(rRV[i], iRV[i]);
258 void HelperTraits<std::complex<T> >::scaleRitzVectors(
259 const std::vector<T>& iRV,
260 Teuchos::SerialDenseMatrix<
int, std::complex<T> >* S )
263 typedef std::complex<T> ST;
264 ST ST_ONE = Teuchos::ScalarTraits<ST>::one();
266 Teuchos::BLAS<int,ST> blas;
268 int i = 0, curDim = S->numRows();
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 );
279 void HelperTraits<std::complex<T> >::computeRitzResiduals(
280 const std::vector<T>& iRV,
281 const Teuchos::SerialDenseMatrix<
int, std::complex<T> >& S,
285 Teuchos::BLAS<int,std::complex<T> > blas;
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();
292 for (
int i=0; i<s_cols; ++i ) {
293 (*RR)[i] = blas.NRM2(s_rows, s_ptr + i*s_stride, 1);
301 #endif // ANASAZI_HELPER_TRAITS_HPP 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.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
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.
This struct is used for storing eigenvalues and Ritz values, as a pair of real values.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
Types and exceptions used within Anasazi solvers and interfaces.
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.
Class which defines basic traits for working with different scalar types.