Anasazi Version of the Day
Tsqr_Mgs.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_Mgs_hpp
00030 #define __TSQR_Tsqr_Mgs_hpp
00031 
00032 #include <algorithm>
00033 #include <cassert>
00034 #include <cmath>
00035 #include <utility> // std::pair
00036 
00037 #include <Tsqr_MessengerBase.hpp>
00038 #include <Tsqr_ScalarTraits.hpp>
00039 #include <Tsqr_Util.hpp>
00040 
00041 #include <Teuchos_RCP.hpp>
00042 
00043 // #define MGS_DEBUG 1
00044 #ifdef MGS_DEBUG
00045 #  include <iostream>
00046 using std::cerr;
00047 using std::endl;
00048 #endif // MGS_DEBUG
00049 
00052 
00053 namespace TSQR {
00054 
00055   template< class LocalOrdinal, class Scalar >
00056   class MGS {
00057   public:
00058     typedef Scalar scalar_type;
00059     typedef LocalOrdinal ordinal_type;
00060     typedef typename ScalarTraits<Scalar>::magnitude_type magnitude_type;
00061 
00062     MGS (const Teuchos::RCP< MessengerBase< Scalar > >& messenger) : 
00063       messenger_ (messenger) {}
00064 
00067     bool QR_produces_R_factor_with_nonnegative_diagonal () const {
00068       return true;
00069     }
00070     
00071     void 
00072     mgs (const LocalOrdinal nrows_local, 
00073    const LocalOrdinal ncols, 
00074    Scalar A_local[], 
00075    const LocalOrdinal lda_local,
00076    Scalar R[],
00077    const LocalOrdinal ldr);
00078 
00079   private:
00080     Teuchos::RCP< MessengerBase< Scalar > > messenger_;
00081   };
00082 
00085 
00086   namespace details {
00087 
00088     template< class LocalOrdinal, class Scalar >
00089     class MgsOps {
00090     public:
00091       MgsOps (const Teuchos::RCP< MessengerBase< Scalar > >& messenger) : 
00092   messenger_ (messenger) {}
00093 
00094       void
00095       axpy (const LocalOrdinal nrows_local,
00096       const Scalar alpha,
00097       const Scalar x_local[],
00098       Scalar y_local[]) const
00099       {
00100   for (LocalOrdinal i = 0; i < nrows_local; ++i)
00101     y_local[i] = y_local[i] + alpha * x_local[i];
00102       }
00103 
00104       void
00105       scale (const LocalOrdinal nrows_local, 
00106        Scalar x_local[], 
00107        const Scalar denom) const
00108       {
00109   for (LocalOrdinal i = 0; i < nrows_local; ++i)
00110     x_local[i] = x_local[i] / denom;
00111       }
00112 
00115       Scalar
00116       dot (const LocalOrdinal nrows_local, 
00117      const Scalar x_local[], 
00118      const Scalar y_local[])
00119       {
00120   Scalar local_result (0);
00121 
00122 #ifdef MGS_DEBUG
00123   // for (LocalOrdinal k = 0; k != nrows_local; ++k)
00124   //   cerr << "(x[" << k << "], y[" << k << "]) = (" << x_local[k] << "," << y_local[k] << ")" << " ";
00125   //   cerr << endl;
00126 #endif // MGS_DEBUG
00127 
00128   for (LocalOrdinal i = 0; i < nrows_local; ++i)
00129     local_result += x_local[i] * ScalarTraits< double >::conj (y_local[i]);
00130 
00131 #ifdef MGS_DEBUG
00132     // cerr << "-- Final value on this proc = " << local_result << endl;
00133 #endif // MGS_DEBUG
00134 
00135   // FIXME (mfh 23 Apr 2010) Does MPI_SUM do the right thing for
00136   // complex or otherwise general MPI data types?  Perhaps an MPI_Op
00137   // should belong in the MessengerBase...
00138   return messenger_->globalSum (local_result);
00139       }
00140 
00141       typename ScalarTraits< Scalar >::magnitude_type 
00142       norm2 (const LocalOrdinal nrows_local, 
00143        const Scalar x_local[])
00144       {
00145   Scalar local_result (0);
00146 
00147   // Doing the right thing in the complex case requires taking
00148   // an absolute value.  We want to avoid this additional cost
00149   // in the real case, which is why we check is_complex.
00150   if (ScalarTraits< Scalar >::is_complex) 
00151     {
00152       for (LocalOrdinal i = 0; i < nrows_local; ++i)
00153         {
00154     const Scalar xi = ScalarTraits< Scalar >::abs (x_local[i]);
00155     local_result += xi * xi;
00156         }
00157     }
00158   else
00159     {
00160       for (LocalOrdinal i = 0; i < nrows_local; ++i)
00161         {
00162     const Scalar xi = x_local[i];
00163     local_result += xi * xi;
00164         }
00165     }
00166   const Scalar global_result = messenger_->globalSum (local_result);
00167   // sqrt doesn't make sense if the type of Scalar is complex,
00168   // even if the imaginary part of global_result is zero.
00169   return sqrt (ScalarTraits< Scalar >::abs (global_result));
00170       }
00171 
00172       Scalar
00173       project (const LocalOrdinal nrows_local, 
00174          const Scalar q_local[], 
00175          Scalar v_local[])
00176       {
00177   const Scalar coeff = this->dot (nrows_local, v_local, q_local);
00178   this->axpy (nrows_local, -coeff, q_local, v_local);
00179   return coeff;
00180       }
00181 
00182     private:
00183       Teuchos::RCP< MessengerBase< Scalar > > messenger_;
00184     };
00185   } // namespace details
00186 
00189 
00190   template< class LocalOrdinal, class Scalar >
00191   void
00192   MGS< LocalOrdinal, Scalar >::mgs (const LocalOrdinal nrows_local, 
00193             const LocalOrdinal ncols, 
00194             Scalar A_local[], 
00195             const LocalOrdinal lda_local,
00196             Scalar R[],
00197             const LocalOrdinal ldr)
00198   {
00199     details::MgsOps< LocalOrdinal, Scalar > ops (messenger_);
00200     
00201     for (LocalOrdinal j = 0; j < ncols; ++j)
00202       {
00203   Scalar* const v = &A_local[j*lda_local];
00204   for (LocalOrdinal i = 0; i < j; ++i)
00205     {
00206       const Scalar* const q = &A_local[i*lda_local];
00207       R[i + j*ldr] = ops.project (nrows_local, q, v);
00208 #ifdef MGS_DEBUG
00209       if (my_rank == 0)
00210         cerr << "(i,j) = (" << i << "," << j << "): coeff = " << R[i + j*ldr] << endl;
00211 #endif // MGS_DEBUG
00212     }
00213   const magnitude_type denom = ops.norm2 (nrows_local, v);
00214 #ifdef MGS_DEBUG
00215     if (my_rank == 0)
00216       cerr << "j = " << j << ": denom = " << denom << endl;
00217 #endif // MGS_DEBUG
00218 
00219   // FIXME (mfh 29 Apr 2010)
00220   //
00221   // NOTE IMPLICIT CAST.  This should work for complex numbers.
00222   // If it doesn't work for your Scalar data type, it means that
00223   // you need a different data type for the diagonal elements of
00224   // the R factor, than you need for the other elements.  This
00225   // is unlikely if we're comparing MGS against a Householder QR
00226   // factorization; I don't really understand how the latter
00227   // would work (not that it couldn't be given a sensible
00228   // interpretation) in the case of Scalars that aren't plain
00229   // old real or complex numbers.
00230   R[j + j*ldr] = Scalar (denom);
00231   ops.scale (nrows_local, v, denom);
00232       }
00233   }
00234 
00235 } // namespace TSQR
00236 
00237 #endif // __TSQR_Tsqr_Mgs_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends