29 #ifndef ANASAZI_HELPER_TRAITS_HPP 30 #define ANASAZI_HELPER_TRAITS_HPP 38 #include "Teuchos_LAPACK.hpp" 49 template <
class ScalarType>
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 );
68 const std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
69 Teuchos::SerialDenseMatrix<int, ScalarType>* S );
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);
82 template<
class ScalarType>
84 const std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& rRV,
85 const std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
88 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
89 MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
91 int curDim = (int)rRV.size();
99 if ( iRV[i] != MT_ZERO ) {
102 (*RV)[i].set(rRV[i], iRV[i]);
103 (*RV)[i+1].set(rRV[i+1], iRV[i+1]);
106 if ( (*RV)[i].imagpart < MT_ZERO ) {
109 (*RV)[i] = (*RV)[i+1];
110 (*RV)[i+1] = tmp_ritz;
112 int tmp_order = (*RO)[i];
113 (*RO)[i] = (*RO)[i+1];
114 (*RO)[i+1] = tmp_order;
117 RI->push_back(1); RI->push_back(-1);
122 (*RV)[i].set(rRV[i], MT_ZERO);
130 template<
class ScalarType>
132 const std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
133 Teuchos::SerialDenseMatrix<int, ScalarType>* S )
135 ScalarType ST_ONE = Teuchos::ScalarTraits<ScalarType>::one();
137 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
138 MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
140 Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
141 Teuchos::BLAS<int,ScalarType> blas;
143 int i = 0, curDim = S->numRows();
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 );
154 temp = blas.NRM2( curDim, s_ptr+i*curDim, 1 );
155 blas.SCAL( curDim, ST_ONE/temp, s_ptr+i*curDim, 1 );
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 )
167 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
168 MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
170 Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
171 Teuchos::BLAS<int,ScalarType> blas;
174 int s_stride = S.stride();
175 int s_rows = S.numRows();
176 int s_cols = S.numCols();
177 ScalarType* s_ptr = S.values();
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];
186 (*RR)[i] = blas.NRM2(s_rows, s_ptr + i*s_stride, 1);
192 #ifdef HAVE_TEUCHOS_COMPLEX 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 );
213 const std::vector<T>& iRV,
214 Teuchos::SerialDenseMatrix<
int, ANSZI_CPLX_CLASS<T> >* S );
217 const std::vector<T>& iRV,
218 const Teuchos::SerialDenseMatrix<
int, ANSZI_CPLX_CLASS<T> >& S,
219 std::vector<T>* RR );
223 void HelperTraits<ANSZI_CPLX_CLASS<T> >::sortRitzValues(
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 )
230 int curDim = (int)rRV.size();
237 while( i < curDim ) {
238 (*RV)[i].set(rRV[i], iRV[i]);
245 void HelperTraits<ANSZI_CPLX_CLASS<T> >::scaleRitzVectors(
246 const std::vector<T>& iRV,
247 Teuchos::SerialDenseMatrix<
int, ANSZI_CPLX_CLASS<T> >* S )
250 typedef ANSZI_CPLX_CLASS<T> ST;
251 ST ST_ONE = Teuchos::ScalarTraits<ST>::one();
253 Teuchos::BLAS<int,ST> blas;
255 int i = 0, curDim = S->numRows();
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 );
266 void HelperTraits<ANSZI_CPLX_CLASS<T> >::computeRitzResiduals(
267 const std::vector<T>& iRV,
268 const Teuchos::SerialDenseMatrix<
int, ANSZI_CPLX_CLASS<T> >& S,
272 Teuchos::BLAS<int,ANSZI_CPLX_CLASS<T> > blas;
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();
279 for (
int i=0; i<s_cols; ++i ) {
280 (*RR)[i] = blas.NRM2(s_rows, s_ptr + i*s_stride, 1);
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.