Tpetra Matrix/Vector Services Version of the Day
MatrixMarket_util.hpp
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //               Tpetra: Linear Algebra Services Package 
00005 //                 Copyright 2011 Sandia Corporation
00006 // 
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00038 // 
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 #ifndef __MatrixMarket_util_hpp
00043 #define __MatrixMarket_util_hpp
00044 
00045 #include <string>
00046 
00047 namespace Tpetra {
00048   namespace MatrixMarket {
00049 
00050     namespace details {
00051 
00073       template<class Scalar>
00074       class SetScientific {
00075       public:
00076   typedef Scalar scalar_type;
00077 
00082   SetScientific (std::ostream& out) : 
00083     out_ (out), originalFlags_ (out.flags()) 
00084   {
00085     typedef Teuchos::ScalarTraits<scalar_type> STS;
00086     typedef typename STS::magnitudeType magnitude_type;
00087     typedef Teuchos::ScalarTraits<magnitude_type> STM;
00088 
00089     // Print floating-point values in scientific notation.
00090     out << std::scientific;
00091 
00092     // We're writing decimal digits, so compute the number of
00093     // digits we need to get reasonable accuracy when reading
00094     // values back in.
00095     //
00096     // There is actually an algorithm, due to Guy Steele (yes,
00097     // _that_ Guy Steele) et al., for idempotent printing of
00098     // finite-length floating-point values.  We should actually
00099     // implement that algorithm, but I don't have time for that
00100     // now.  Currently, I just print no more than (one decimal
00101     // digit more than (the number of decimal digits justified
00102     // by the precision of magnitude_type)).
00103 
00104     // FIXME (mfh 06 Apr 2011) This will only work if
00105     // std::log10() is defined for magnitude_type inputs.
00106     // Teuchos::ScalarTraits does not currently have an log10()
00107     // class method.
00108     const magnitude_type numDecDigits = STM::t() * std::log10 (STM::base());
00109 
00110     // Round and add one. Hopefully this doesn't overflow...
00111     const magnitude_type one = STM::one();
00112     const magnitude_type two = one + one;
00113     const int prec = 1 + static_cast<int> ((two*numDecDigits + one) / two);
00114     
00115     // Set the number of (decimal) digits after the decimal
00116     // point to print.
00117     out.precision (prec);
00118   } 
00119 
00125   ~SetScientific () {
00126     out_.flags (originalFlags_);
00127   }
00128 
00129       private:
00131   std::ostream& out_;
00132 
00134   std::ios_base::fmtflags originalFlags_;
00135       };
00136 
00137       // We define a class because template functions can't (in the
00138       // current C++ standard) have default template parameters.
00139       template<class Scalar, bool isComplex=Teuchos::ScalarTraits<Scalar>::isComplex>
00140       class ScalarAssigner {
00141       public:
00142   static void
00143   assign (Scalar& val,
00144     const typename Teuchos::ScalarTraits<Scalar>::magnitudeType& real,
00145     const typename Teuchos::ScalarTraits<Scalar>::magnitudeType& imag);
00146       };
00147 
00148       template<class RealType>
00149       class ScalarAssigner<RealType, false> {
00150       public:
00151   static void
00152   assign (RealType& val,
00153     const typename Teuchos::ScalarTraits<RealType>::magnitudeType& real,
00154     const typename Teuchos::ScalarTraits<RealType>::magnitudeType& imag)
00155   {
00156     // imag had better be zero.  We're ignoring it regardless.
00157     (void) imag;
00158     val = real; 
00159   }
00160       };
00161 
00162 #ifdef HAVE_TEUCHOS_COMPLEX
00163       template<class MagType>
00164       class ScalarAssigner<std::complex<MagType>, true> {
00165       public:
00166   static void
00167   assign (std::complex<MagType>& val,
00168     const typename Teuchos::ScalarTraits<std::complex<MagType> >::magnitudeType& real,
00169     const typename Teuchos::ScalarTraits<std::complex<MagType> >::magnitudeType& imag)
00170   {
00171     val = std::complex<MagType> (real, imag);
00172   }
00173       };
00174 #endif // HAVE_TEUCHOS_COMPLEX
00175 
00176       // \fn assignScalar
00177       // \brief val = S(real, imag).
00178       //
00179       // We have to template it because we don't know that S is a
00180       // complex type; if we write S(real,imag), the compiler will
00181       // complain if S is a real type.
00182       template<class Scalar>
00183       void 
00184       assignScalar (Scalar& val, 
00185         const typename Teuchos::ScalarTraits<Scalar>::magnitudeType& real,
00186         const typename Teuchos::ScalarTraits<Scalar>::magnitudeType& imag)
00187       {
00188   ScalarAssigner<Scalar>::assign (val, real, imag);
00189       }
00190 
00191     } // namespace details
00192 
00193 
00194     static bool isSkew (const std::string& symmType) {
00195       return symmType.size() >= 4 && symmType.substr(0,4) == "skew";
00196     }
00197     static bool isConj (const std::string& symmType) {
00198       return std::string::npos != symmType.find ("hermitian");
00199     }
00200     static bool needsSymmetrization (const std::string& symmType) {
00201       return symmType != "general";
00202     }
00203 
00215     template<class AdderType>
00216     class SymmetrizingAdder {
00217     public:
00218       typedef typename AdderType::index_type index_type;
00219       typedef typename AdderType::value_type value_type;
00220     
00221       SymmetrizingAdder (const Teuchos::RCP<AdderType>& adder, 
00222        const std::string& symmType) :
00223   adder_ (adder),
00224   symmetrize_ (needsSymmetrization (symmType)),
00225   conjugate_ (isConj (symmType)),
00226   skew_ (isSkew (symmType))
00227       {}
00228 
00229 
00230     
00231       void operator() (const index_type i, const index_type j, const value_type& Aij) {
00232   AdderType& theAdder = *adder_;
00233 
00234   theAdder (i, j, Aij);
00235   if (symmetrize_ && i != j) {
00236     typedef Teuchos::ScalarTraits<value_type> STS;
00237     const value_type Aji = skew_ ? 
00238       -(conjugate_ ? STS::conjugate(Aij) : Aij) : 
00239       +(conjugate_ ? STS::conjugate(Aij) : Aij);
00240     theAdder (j, i, Aji);
00241   }
00242       }
00243 
00247       Teuchos::RCP<AdderType> getAdder() const {
00248   return adder_;
00249       }
00250 
00251     private:
00252       Teuchos::RCP<AdderType> adder_;
00253       bool symmetrize_, conjugate_, skew_;
00254     };
00255 
00256   } // namespace MatrixMarket
00257 } // namespace Tpetra
00258 
00259 #endif // __MatrixMarket_util_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines