Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 
53 namespace Anasazi {
54 
62  template <class ScalarType>
63  class HelperTraits
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,
83 
85 
87  static void computeRitzResiduals(
88  const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
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,
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 
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,
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 
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 
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;
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
ScalarType * values() const
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const
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.
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...
ScalarType LAPY2(const ScalarType &x, const ScalarType &y) const
OrdinalType numCols() const
void SCAL(const OrdinalType &n, const ScalarType &alpha, ScalarType *x, const OrdinalType &incx) const
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.
OrdinalType stride() const
OrdinalType numRows() const
Class which defines basic traits for working with different scalar types.