Anasazi  Version of the Day
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, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
29 #ifndef ANASAZI_HELPER_TRAITS_HPP
30 #define ANASAZI_HELPER_TRAITS_HPP
31 
36 #include "AnasaziConfigDefs.hpp"
37 #include "AnasaziTypes.hpp"
38 #include "Teuchos_LAPACK.hpp"
39 
40 namespace Anasazi {
41 
49  template <class ScalarType>
50  class HelperTraits
51  {
52  public:
53 
55 
58  static void sortRitzValues(
59  const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& rRV,
60  const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
61  std::vector<Value<ScalarType> >* RV, std::vector<int>* RO, std::vector<int>* RI );
62 
64 
67  static void scaleRitzVectors(
68  const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
69  Teuchos::SerialDenseMatrix<int, ScalarType>* S );
70 
72 
74  static void computeRitzResiduals(
75  const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
76  const Teuchos::SerialDenseMatrix<int, ScalarType>& S,
77  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>* RR);
78 
79  };
80 
81 
82  template<class ScalarType>
84  const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& rRV,
85  const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
86  std::vector<Value<ScalarType> >* RV, std::vector<int>* RO, std::vector<int>* RI )
87  {
88  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
89  MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
90 
91  int curDim = (int)rRV.size();
92  int i = 0;
93 
94  // Clear the current index.
95  RI->clear();
96 
97  // Place the Ritz values from rRV and iRV into the RV container.
98  while( i < curDim ) {
99  if ( iRV[i] != MT_ZERO ) {
100  //
101  // We will have this situation for real-valued, non-Hermitian matrices.
102  (*RV)[i].set(rRV[i], iRV[i]);
103  (*RV)[i+1].set(rRV[i+1], iRV[i+1]);
104 
105  // Make sure that complex conjugate pairs have their positive imaginary part first.
106  if ( (*RV)[i].imagpart < MT_ZERO ) {
107  // The negative imaginary part is first, so swap the order of the ritzValues and ritzOrders.
108  Anasazi::Value<ScalarType> tmp_ritz( (*RV)[i] );
109  (*RV)[i] = (*RV)[i+1];
110  (*RV)[i+1] = tmp_ritz;
111 
112  int tmp_order = (*RO)[i];
113  (*RO)[i] = (*RO)[i+1];
114  (*RO)[i+1] = tmp_order;
115 
116  }
117  RI->push_back(1); RI->push_back(-1);
118  i = i+2;
119  } else {
120  //
121  // The Ritz value is not complex.
122  (*RV)[i].set(rRV[i], MT_ZERO);
123  RI->push_back(0);
124  i++;
125  }
126  }
127  }
128 
129 
130  template<class ScalarType>
132  const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
133  Teuchos::SerialDenseMatrix<int, ScalarType>* S )
134  {
135  ScalarType ST_ONE = Teuchos::ScalarTraits<ScalarType>::one();
136 
137  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
138  MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
139 
140  Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
141  Teuchos::BLAS<int,ScalarType> blas;
142 
143  int i = 0, curDim = S->numRows();
144  ScalarType temp;
145  ScalarType* s_ptr = S->values();
146  while( i < curDim ) {
147  if ( iRV[i] != MT_ZERO ) {
148  temp = lapack_mag.LAPY2( blas.NRM2( curDim, s_ptr+i*curDim, 1 ),
149  blas.NRM2( curDim, s_ptr+(i+1)*curDim, 1 ) );
150  blas.SCAL( curDim, ST_ONE/temp, s_ptr+i*curDim, 1 );
151  blas.SCAL( curDim, ST_ONE/temp, s_ptr+(i+1)*curDim, 1 );
152  i = i+2;
153  } else {
154  temp = blas.NRM2( curDim, s_ptr+i*curDim, 1 );
155  blas.SCAL( curDim, ST_ONE/temp, s_ptr+i*curDim, 1 );
156  i++;
157  }
158  }
159  }
160 
161  template<class ScalarType>
163  const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
164  const Teuchos::SerialDenseMatrix<int, ScalarType>& S,
165  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>* RR )
166  {
167  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
168  MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
169 
170  Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
171  Teuchos::BLAS<int,ScalarType> blas;
172 
173  int i = 0;
174  int s_stride = S.stride();
175  int s_rows = S.numRows();
176  int s_cols = S.numCols();
177  ScalarType* s_ptr = S.values();
178 
179  while( i < s_cols ) {
180  if ( iRV[i] != MT_ZERO ) {
181  (*RR)[i] = lapack_mag.LAPY2( blas.NRM2(s_rows, s_ptr + i*s_stride, 1),
182  blas.NRM2(s_rows, s_ptr + (i+1)*s_stride, 1) );
183  (*RR)[i+1] = (*RR)[i];
184  i = i+2;
185  } else {
186  (*RR)[i] = blas.NRM2(s_rows, s_ptr + i*s_stride, 1);
187  i++;
188  }
189  }
190  }
191 
192 #ifdef HAVE_TEUCHOS_COMPLEX
193  // Partial template specializations for the complex scalar type.
194 
202  template <class T>
203  class HelperTraits<ANSZI_CPLX_CLASS<T> >
204  {
205  public:
206  static void sortRitzValues(
207  const std::vector<T>& rRV,
208  const std::vector<T>& iRV,
209  std::vector<Value<ANSZI_CPLX_CLASS<T> > >* RV,
210  std::vector<int>* RO, std::vector<int>* RI );
211 
212  static void scaleRitzVectors(
213  const std::vector<T>& iRV,
214  Teuchos::SerialDenseMatrix<int, ANSZI_CPLX_CLASS<T> >* S );
215 
216  static void computeRitzResiduals(
217  const std::vector<T>& iRV,
218  const Teuchos::SerialDenseMatrix<int, ANSZI_CPLX_CLASS<T> >& S,
219  std::vector<T>* RR );
220  };
221 
222  template<class T>
224  const std::vector<T>& rRV,
225  const std::vector<T>& iRV,
226  std::vector<Value<ANSZI_CPLX_CLASS<T> > >* RV,
227  std::vector<int>* RO, std::vector<int>* RI )
228  {
229  (void)RO;
230  int curDim = (int)rRV.size();
231  int i = 0;
232 
233  // Clear the current index.
234  RI->clear();
235 
236  // Place the Ritz values from rRV and iRV into the RV container.
237  while( i < curDim ) {
238  (*RV)[i].set(rRV[i], iRV[i]);
239  RI->push_back(0);
240  i++;
241  }
242  }
243 
244  template<class T>
246  const std::vector<T>& iRV,
247  Teuchos::SerialDenseMatrix<int, ANSZI_CPLX_CLASS<T> >* S )
248  {
249  (void)iRV;
250  typedef ANSZI_CPLX_CLASS<T> ST;
251  ST ST_ONE = Teuchos::ScalarTraits<ST>::one();
252 
253  Teuchos::BLAS<int,ST> blas;
254 
255  int i = 0, curDim = S->numRows();
256  ST temp;
257  ST* s_ptr = S->values();
258  while( i < curDim ) {
259  temp = blas.NRM2( curDim, s_ptr+i*curDim, 1 );
260  blas.SCAL( curDim, ST_ONE/temp, s_ptr+i*curDim, 1 );
261  i++;
262  }
263  }
264 
265  template<class T>
267  const std::vector<T>& iRV,
268  const Teuchos::SerialDenseMatrix<int, ANSZI_CPLX_CLASS<T> >& S,
269  std::vector<T>* RR )
270  {
271  (void)iRV;
272  Teuchos::BLAS<int,ANSZI_CPLX_CLASS<T> > blas;
273 
274  int s_stride = S.stride();
275  int s_rows = S.numRows();
276  int s_cols = S.numCols();
277  ANSZI_CPLX_CLASS<T>* s_ptr = S.values();
278 
279  for (int i=0; i<s_cols; ++i ) {
280  (*RR)[i] = blas.NRM2(s_rows, s_ptr + i*s_stride, 1);
281  }
282  }
283 #endif
284 
285 } // end Anasazi namespace
286 
287 
288 #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.