Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 
21 namespace 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
Teuchos::SerialDenseMatrix::values
ScalarType * values() const
Teuchos::DefaultBLASImpl::NRM2
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const
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
Teuchos_LAPACK.hpp
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::Value
This struct is used for storing eigenvalues and Ritz values, as a pair of real values.
Definition: AnasaziTypes.hpp:42
AnasaziConfigDefs.hpp
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::LAPACK::LAPY2
ScalarType LAPY2(const ScalarType &x, const ScalarType &y) const
Teuchos::SerialDenseMatrix::numCols
OrdinalType numCols() const
Teuchos::DefaultBLASImpl::SCAL
void SCAL(const OrdinalType &n, const ScalarType &alpha, ScalarType *x, const OrdinalType &incx) const
AnasaziTypes.hpp
Types and exceptions used within Anasazi solvers and interfaces.
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
Teuchos::SerialDenseMatrix::stride
OrdinalType stride() const
Teuchos::SerialDenseMatrix::numRows
OrdinalType numRows() const
Anasazi::HelperTraits
Class which defines basic traits for working with different scalar types.
Definition: AnasaziHelperTraits.hpp:31
Generated on Fri Dec 20 2024 09:18:33 for Anasazi by   doxygen 1.8.5