Anasazi Version of the Day
Tsqr_Util.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2010) 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef __TSQR_Tsqr_Util_hpp
00030 #define __TSQR_Tsqr_Util_hpp
00031 
00032 #include "Tsqr_ScalarTraits.hpp"
00033 
00034 #include <algorithm>
00035 #include <complex>
00036 #include <ostream>
00037 
00040 
00041 namespace TSQR {
00042 
00051   template< class Scalar, bool isComplex >
00052   class ScalarPrinter {
00053   public:
00056     void operator() (std::ostream& out, const Scalar& elt) const;
00057   };
00058 
00059   // Partial specialization for real Scalar
00060   template< class Scalar >
00061   class ScalarPrinter< Scalar, false > {
00062   public:
00063     void operator() (std::ostream& out, const Scalar& elt) const {
00064       out << elt;
00065     }
00066   };
00067 
00068   // Partial specialization for complex Scalar
00069   template< class Scalar >
00070   class ScalarPrinter< Scalar, true > {
00071   public:
00072     void operator() (std::ostream& out, const Scalar& elt) const {
00073       typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00074 
00075       const magnitude_type ZERO (0);
00076       const magnitude_type& realPart = std::real (elt);
00077       const magnitude_type& imagPart = std::imag (elt);
00078 
00079       out << realPart;
00080       if (imagPart < ZERO)
00081   out << "-" << ScalarTraits< Scalar >::abs(imagPart) << "*i";
00082       else if (imagPart > ZERO)
00083   out << "+" << imagPart << "*i";
00084     }
00085   };
00086 
00087   template< class LocalOrdinal, class Scalar >
00088   void
00089   print_local_matrix (std::ostream& out,
00090           const LocalOrdinal nrows_local,
00091           const LocalOrdinal ncols,
00092           const Scalar A[],
00093           const LocalOrdinal lda)
00094   {
00095     ScalarPrinter< Scalar, ScalarTraits< Scalar >::is_complex > printer;
00096     for (LocalOrdinal i = 0; i < nrows_local; i++)
00097       {
00098   for (LocalOrdinal j = 0; j < ncols; j++)
00099     {
00100       const Scalar& curElt = A[i + j*lda];
00101       printer (out, curElt);
00102       if (j < ncols - 1)
00103         out << ", ";
00104     }
00105   out << ";" << std::endl;
00106       }
00107   }
00108 
00109   template< class Ordinal, class Scalar >
00110   void
00111   copy_matrix (const Ordinal nrows,
00112          const Ordinal ncols,
00113          Scalar* const A,
00114          const Ordinal lda,
00115          const Scalar* const B,
00116          const Ordinal ldb)
00117   {
00118     for (Ordinal j = 0; j < ncols; j++)
00119       {
00120   Scalar* const A_j = &A[j*lda];
00121   const Scalar* const B_j = &B[j*ldb];
00122   std::copy (B_j, B_j + nrows, A_j);
00123       }
00124   }
00125 
00126   template< class Ordinal, class Scalar >
00127   void
00128   fill_matrix (const Ordinal nrows,
00129          const Ordinal ncols,
00130          Scalar* const A,
00131          const Ordinal lda,
00132          const Scalar& default_val)
00133   {
00134     for (Ordinal j = 0; j < ncols; j++)
00135       {
00136   Scalar* const A_j = &A[j*lda];
00137   std::fill (A_j, A_j + nrows, default_val);
00138       }
00139   }
00140 
00141 
00142   template< class Ordinal, class Scalar, class Generator >
00143   void
00144   generate_matrix (const Ordinal nrows,
00145        const Ordinal ncols,
00146        Scalar* const A,
00147        const Ordinal lda,
00148        Generator gen)
00149   {
00150     for (Ordinal j = 0; j < ncols; j++)
00151       {
00152   Scalar* const A_j = &A[j*lda];
00153   std::generate (A_j, A_j + nrows, gen);
00154       }
00155   }
00156 
00157   template< class Ordinal, class Scalar >
00158   void
00159   copy_upper_triangle (const Ordinal nrows,
00160            const Ordinal ncols,
00161            Scalar* const R_out,
00162            const Ordinal ldr_out,
00163            const Scalar* const R_in,
00164            const Ordinal ldr_in)
00165   {
00166     if (nrows >= ncols)
00167       {
00168   for (Ordinal j = 0; j < ncols; j++)
00169     {
00170       Scalar* const A_j = &R_out[j*ldr_out];
00171       const Scalar* const B_j = &R_in[j*ldr_in];
00172       for (Ordinal i = 0; i <= j; i++)
00173         A_j[i] = B_j[i];
00174     }
00175       }
00176     else
00177       {
00178   copy_upper_triangle (nrows, nrows, R_out, ldr_out, R_in, ldr_in);
00179   for (Ordinal j = nrows; j < ncols; j++)
00180     {
00181       Scalar* const A_j = &R_out[j*ldr_out];
00182       const Scalar* const B_j = &R_in[j*ldr_in];
00183       for (Ordinal i = 0; i < nrows; i++)
00184         A_j[i] = B_j[i];
00185     }
00186       }
00187   }
00188 
00189 
00190   template< class Scalar >
00191   class SumSquare {
00192   public:
00193     Scalar operator() (const Scalar& result, const Scalar& x) const { 
00194       return result + x*x; 
00195     }
00196   };
00197 
00198   // Specialization for complex numbers
00199   template< class Scalar >
00200   class SumSquare< std::complex< Scalar > >  {
00201   public:
00202     Scalar operator() (const std::complex<Scalar>& result, 
00203            const std::complex<Scalar>& x) const { 
00204       const Scalar absval = std::norm (x);
00205       return result + absval * absval; 
00206     }
00207   };
00208 
00209   template< class Ordinal, class Scalar >
00210   void
00211   pack_R_factor (const Ordinal nrows, 
00212      const Ordinal ncols, 
00213      const Scalar R_in[], 
00214      const Ordinal ldr_in,
00215      Scalar buffer[])
00216   {
00217     Ordinal count = 0; // current position in output buffer
00218     if (nrows >= ncols)
00219       for (Ordinal j = 0; j < ncols; j++)
00220   for (Ordinal i = 0; i <= j; i++)
00221     buffer[count++] = R_in[i + j*ldr_in];
00222     else
00223       for (Ordinal j = 0; j < nrows; j++)
00224   for (Ordinal i = 0; i <= j; i++)
00225     buffer[count++] = R_in[i + j*ldr_in];
00226   }
00227 
00228   template< class Ordinal, class Scalar >
00229   void
00230   unpack_R_factor (const Ordinal nrows, 
00231        const Ordinal ncols, 
00232        Scalar R_out[], 
00233        const Ordinal ldr_out,
00234        const Scalar buffer[])
00235   {
00236     Ordinal count = 0; // current position in input buffer
00237     if (nrows >= ncols)
00238       for (Ordinal j = 0; j < ncols; j++)
00239   for (Ordinal i = 0; i <= j; i++)
00240     R_out[i + j*ldr_out] = buffer[count++];
00241     else
00242       for (Ordinal j = 0; j < nrows; j++)
00243   for (Ordinal i = 0; i <= j; i++)
00244     R_out[i + j*ldr_out] = buffer[count++];
00245   }
00246 
00247 } // namespace TSQR
00248 
00249 #endif // __TSQR_Tsqr_Util_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends