Anasazi Version of the Day
AnasaziHelperTraits.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
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         // Clear the current index.
00095         RI->clear();
00096 
00097         // Place the Ritz values from rRV and iRV into the RV container.
00098         while( i < curDim ) {
00099             if ( iRV[i] != MT_ZERO ) {
00100                 //
00101                 // We will have this situation for real-valued, non-Hermitian matrices.
00102                 (*RV)[i].set(rRV[i], iRV[i]);
00103                 (*RV)[i+1].set(rRV[i+1], iRV[i+1]);
00104 
00105                 // Make sure that complex conjugate pairs have their positive imaginary part first.
00106                 if ( (*RV)[i].imagpart < MT_ZERO ) {
00107                     // The negative imaginary part is first, so swap the order of the ritzValues and ritzOrders.
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                 // The Ritz value is not complex.
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     // Partial template specializations for the complex scalar type.
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         // Clear the current index.
00234         RI->clear();
00235 
00236         // Place the Ritz values from rRV and iRV into the RV container.
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 } // end Anasazi namespace
00286 
00287 
00288 #endif // ANASAZI_HELPER_TRAITS_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends