00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #ifndef ANASAZI_HELPER_TRAITS_HPP
00030 #define ANASAZI_HELPER_TRAITS_HPP
00031
00036 #include "AnasaziConfigDefs.hpp"
00037 #include "AnasaziTypes.hpp"
00038 #include "Teuchos_LAPACK.hpp"
00039
00040 namespace Anasazi {
00041
00049 template <class ScalarType>
00050 class HelperTraits
00051 {
00052 public:
00053
00055
00058 static void sortRitzValues(
00059 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& rRV,
00060 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
00061 std::vector<Value<ScalarType> >* RV, std::vector<int>* RO, std::vector<int>* RI );
00062
00064
00067 static void scaleRitzVectors(
00068 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
00069 Teuchos::SerialDenseMatrix<int, ScalarType>* S );
00070
00072
00074 static void computeRitzResiduals(
00075 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
00076 const Teuchos::SerialDenseMatrix<int, ScalarType>& S,
00077 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>* RR);
00078
00079 };
00080
00081
00082 template<class ScalarType>
00083 void HelperTraits<ScalarType>::sortRitzValues(
00084 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& rRV,
00085 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
00086 std::vector<Value<ScalarType> >* RV, std::vector<int>* RO, std::vector<int>* RI )
00087 {
00088 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00089 MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
00090
00091 int curDim = (int)rRV.size();
00092 int i = 0;
00093
00094
00095 RI->clear();
00096
00097
00098 while( i < curDim ) {
00099 if ( iRV[i] != MT_ZERO ) {
00100
00101
00102 (*RV)[i].set(rRV[i], iRV[i]);
00103 (*RV)[i+1].set(rRV[i+1], iRV[i+1]);
00104
00105
00106 if ( (*RV)[i].imagpart < MT_ZERO ) {
00107
00108 Anasazi::Value<ScalarType> tmp_ritz( (*RV)[i] );
00109 (*RV)[i] = (*RV)[i+1];
00110 (*RV)[i+1] = tmp_ritz;
00111
00112 int tmp_order = (*RO)[i];
00113 (*RO)[i] = (*RO)[i+1];
00114 (*RO)[i+1] = tmp_order;
00115
00116 }
00117 RI->push_back(1); RI->push_back(-1);
00118 i = i+2;
00119 } else {
00120
00121
00122 (*RV)[i].set(rRV[i], MT_ZERO);
00123 RI->push_back(0);
00124 i++;
00125 }
00126 }
00127 }
00128
00129
00130 template<class ScalarType>
00131 void HelperTraits<ScalarType>::scaleRitzVectors(
00132 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
00133 Teuchos::SerialDenseMatrix<int, ScalarType>* S )
00134 {
00135 ScalarType ST_ONE = Teuchos::ScalarTraits<ScalarType>::one();
00136
00137 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00138 MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
00139
00140 Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
00141 Teuchos::BLAS<int,ScalarType> blas;
00142
00143 int i = 0, curDim = S->numRows();
00144 ScalarType temp;
00145 ScalarType* s_ptr = S->values();
00146 while( i < curDim ) {
00147 if ( iRV[i] != MT_ZERO ) {
00148 temp = lapack_mag.LAPY2( blas.NRM2( curDim, s_ptr+i*curDim, 1 ),
00149 blas.NRM2( curDim, s_ptr+(i+1)*curDim, 1 ) );
00150 blas.SCAL( curDim, ST_ONE/temp, s_ptr+i*curDim, 1 );
00151 blas.SCAL( curDim, ST_ONE/temp, s_ptr+(i+1)*curDim, 1 );
00152 i = i+2;
00153 } else {
00154 temp = blas.NRM2( curDim, s_ptr+i*curDim, 1 );
00155 blas.SCAL( curDim, ST_ONE/temp, s_ptr+i*curDim, 1 );
00156 i++;
00157 }
00158 }
00159 }
00160
00161 template<class ScalarType>
00162 void HelperTraits<ScalarType>::computeRitzResiduals(
00163 const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& iRV,
00164 const Teuchos::SerialDenseMatrix<int, ScalarType>& S,
00165 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>* RR )
00166 {
00167 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00168 MagnitudeType MT_ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
00169
00170 Teuchos::LAPACK<int,MagnitudeType> lapack_mag;
00171 Teuchos::BLAS<int,ScalarType> blas;
00172
00173 int i = 0;
00174 int s_stride = S.stride();
00175 int s_rows = S.numRows();
00176 int s_cols = S.numCols();
00177 ScalarType* s_ptr = S.values();
00178
00179 while( i < s_cols ) {
00180 if ( iRV[i] != MT_ZERO ) {
00181 (*RR)[i] = lapack_mag.LAPY2( blas.NRM2(s_rows, s_ptr + i*s_stride, 1),
00182 blas.NRM2(s_rows, s_ptr + (i+1)*s_stride, 1) );
00183 (*RR)[i+1] = (*RR)[i];
00184 i = i+2;
00185 } else {
00186 (*RR)[i] = blas.NRM2(s_rows, s_ptr + i*s_stride, 1);
00187 i++;
00188 }
00189 }
00190 }
00191
00192 #ifdef HAVE_TEUCHOS_COMPLEX
00193
00194
00202 template <class T>
00203 class HelperTraits<ANSZI_CPLX_CLASS<T> >
00204 {
00205 public:
00206 static void sortRitzValues(
00207 const std::vector<T>& rRV,
00208 const std::vector<T>& iRV,
00209 std::vector<Value<ANSZI_CPLX_CLASS<T> > >* RV,
00210 std::vector<int>* RO, std::vector<int>* RI );
00211
00212 static void scaleRitzVectors(
00213 const std::vector<T>& iRV,
00214 Teuchos::SerialDenseMatrix<int, ANSZI_CPLX_CLASS<T> >* S );
00215
00216 static void computeRitzResiduals(
00217 const std::vector<T>& iRV,
00218 const Teuchos::SerialDenseMatrix<int, ANSZI_CPLX_CLASS<T> >& S,
00219 std::vector<T>* RR );
00220 };
00221
00222 template<class T>
00223 void HelperTraits<ANSZI_CPLX_CLASS<T> >::sortRitzValues(
00224 const std::vector<T>& rRV,
00225 const std::vector<T>& iRV,
00226 std::vector<Value<ANSZI_CPLX_CLASS<T> > >* RV,
00227 std::vector<int>* RO, std::vector<int>* RI )
00228 {
00229 (void)RO;
00230 int curDim = (int)rRV.size();
00231 int i = 0;
00232
00233
00234 RI->clear();
00235
00236
00237 while( i < curDim ) {
00238 (*RV)[i].set(rRV[i], iRV[i]);
00239 RI->push_back(0);
00240 i++;
00241 }
00242 }
00243
00244 template<class T>
00245 void HelperTraits<ANSZI_CPLX_CLASS<T> >::scaleRitzVectors(
00246 const std::vector<T>& iRV,
00247 Teuchos::SerialDenseMatrix<int, ANSZI_CPLX_CLASS<T> >* S )
00248 {
00249 (void)iRV;
00250 typedef ANSZI_CPLX_CLASS<T> ST;
00251 ST ST_ONE = Teuchos::ScalarTraits<ST>::one();
00252
00253 Teuchos::BLAS<int,ST> blas;
00254
00255 int i = 0, curDim = S->numRows();
00256 ST temp;
00257 ST* s_ptr = S->values();
00258 while( i < curDim ) {
00259 temp = blas.NRM2( curDim, s_ptr+i*curDim, 1 );
00260 blas.SCAL( curDim, ST_ONE/temp, s_ptr+i*curDim, 1 );
00261 i++;
00262 }
00263 }
00264
00265 template<class T>
00266 void HelperTraits<ANSZI_CPLX_CLASS<T> >::computeRitzResiduals(
00267 const std::vector<T>& iRV,
00268 const Teuchos::SerialDenseMatrix<int, ANSZI_CPLX_CLASS<T> >& S,
00269 std::vector<T>* RR )
00270 {
00271 (void)iRV;
00272 Teuchos::BLAS<int,ANSZI_CPLX_CLASS<T> > blas;
00273
00274 int s_stride = S.stride();
00275 int s_rows = S.numRows();
00276 int s_cols = S.numCols();
00277 ANSZI_CPLX_CLASS<T>* s_ptr = S.values();
00278
00279 for (int i=0; i<s_cols; ++i ) {
00280 (*RR)[i] = blas.NRM2(s_rows, s_ptr + i*s_stride, 1);
00281 }
00282 }
00283 #endif
00284
00285 }
00286
00287
00288 #endif // ANASAZI_HELPER_TRAITS_HPP