Kokkos Node API and Local Linear Algebra Kernels Version of the Day
Tsqr_GlobalVerify.hpp
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //          Kokkos: Node API and Parallel Node Kernels
00005 //              Copyright (2008) 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 __TSQR_Tsqr_GlobalVerify_hpp
00043 #define __TSQR_Tsqr_GlobalVerify_hpp
00044 
00045 #include <Tsqr_LocalVerify.hpp>
00046 #include <Tsqr_MessengerBase.hpp>
00047 #include <Tsqr_Blas.hpp>
00048 #include <Tsqr_Util.hpp>
00049 
00050 #include <Teuchos_ScalarTraits.hpp>
00051 
00052 #include <utility> // std::pair
00053 #include <vector>
00054 
00055 
00056 namespace TSQR {
00057 
00070   template<class Scalar, bool isComplex = Teuchos::ScalarTraits< Scalar >::isComplex>
00071   class GlobalSummer {
00072   public:
00073     typedef Scalar scalar_type;
00074     typedef Teuchos::ScalarTraits< Scalar > STS;
00075     typedef typename STS::magnitudeType magnitude_type;
00076 
00077     static magnitude_type
00078     sum (const Scalar& localSum,
00079    MessengerBase<Scalar>* const messenger);
00080   };
00081 
00082   // Complex-arithmetic "forward declaration"
00083   template<class Scalar>
00084   class GlobalSummer<Scalar, true> {
00085   public:
00086     typedef Scalar scalar_type;
00087     typedef Teuchos::ScalarTraits<Scalar> STS;
00088     typedef typename STS::magnitudeType magnitude_type;
00089 
00090     static magnitude_type
00091     sum (const Scalar& localSum,
00092    MessengerBase<Scalar>* const messenger);
00093   };
00094 
00095   // Real-arithmetic "forward declaration"
00096   template<class Scalar>
00097   class GlobalSummer<Scalar, false> {
00098   public:
00099     typedef Scalar scalar_type;
00100     typedef Teuchos::ScalarTraits<Scalar> STS;
00101     typedef typename STS::magnitudeType magnitude_type;
00102 
00103     static magnitude_type
00104     sum (const Scalar& localSum,
00105    MessengerBase<Scalar>* const messenger);
00106   };
00107 
00108   // Complex-arithmetic case
00109   template<class Scalar>
00110   typename GlobalSummer<Scalar, true>::magnitude_type
00111   GlobalSummer<Scalar, true>::sum (const Scalar& localSum,
00112            MessengerBase<Scalar>* const messenger)
00113   {
00114     // In order to use a MessengerBase<Scalar> on magnitude_type
00115     // values, we have to convert local_result to a Scalar, and then
00116     // convert back the result.  We convert by setting the real
00117     // component of the Scalar to the magnitude_type.  This isn't
00118     // guaranteed to work if magnitude_type has a greater dynamic
00119     // range than Scalar.  That's possible, but that's not how we do
00120     // things with ScalarTraits< std::complex< T > >, and that's not
00121     // how LAPACK does it either, so it's fair to assume that
00122     // magnitude_type and the individual components of Scalar have the
00123     // same dynamic range.
00124     const magnitude_type localSumAbs = STS::magnitude (localSum);
00125     const Scalar localSumAsScalar (localSumAbs, magnitude_type(0));
00126     const Scalar globalSumAsScalar = messenger->globalSum (localSumAsScalar);
00127     const magnitude_type globalSum = STS::magnitude (globalSumAsScalar);
00128     return globalSum;
00129   }
00130 
00131   // Real-arithmetic case
00132   template<class Scalar>
00133   typename GlobalSummer<Scalar, false>::magnitude_type
00134   GlobalSummer<Scalar, false>::sum (const Scalar& localSum,
00135             MessengerBase<Scalar>* const messenger)
00136   {
00137     const Scalar localSumAsScalar (localSum);
00138     const Scalar globalSumAsScalar = messenger->globalSum (localSumAsScalar);
00139     const magnitude_type globalSum = STS::magnitude (globalSumAsScalar);
00140     return globalSum;
00141   }
00142 
00143   template<class LocalOrdinal, class Scalar>
00144   typename Teuchos::ScalarTraits<Scalar>::magnitudeType
00145   global_frobenius_norm (const LocalOrdinal nrows_local, 
00146        const LocalOrdinal ncols,
00147        const Scalar A_local[],
00148        const LocalOrdinal lda_local,
00149        MessengerBase<Scalar>* const messenger)
00150   {
00151     typedef Teuchos::ScalarTraits<Scalar> STS;
00152     typedef typename STS::magnitudeType magnitude_type;
00153 
00154     // FIXME (mfh 20 Apr 2010) This is currently implemented using an
00155     // all-reduction.  This may result in different processors getting
00156     // slightly different answers, due to floating-point arithmetic
00157     // roundoff.  We might not want this if we are using this function
00158     // to test a routine.
00159 
00160     magnitude_type localResult (0);
00161     for (LocalOrdinal j = 0; j < ncols; j++)
00162       {
00163   const Scalar* const cur_col = &A_local[j*lda_local];
00164   for (LocalOrdinal i = 0; i < nrows_local; ++i)
00165     {
00166       const magnitude_type abs_xi = STS::magnitude (cur_col[i]);
00167       localResult = localResult + abs_xi * abs_xi;
00168     }
00169       }
00170     // GlobalSummmer() is a hack to let us use a Scalar - type
00171     // MessengerBase with magnitude_type inputs and outputs.
00172     // Otherwise we would need to carry around a
00173     // MessengerBase<magnitude_type> object as well.
00174     const magnitude_type globalResult = 
00175       GlobalSummer<Scalar, STS::isComplex>::sum (localResult, messenger);
00176     return sqrt (globalResult);
00177   }
00178 
00179   template<class LocalOrdinal, class Scalar>
00180   std::vector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>
00181   global_verify (const LocalOrdinal nrows_local, 
00182      const LocalOrdinal ncols, 
00183      const Scalar A_local[],
00184      const LocalOrdinal lda_local,
00185      const Scalar Q_local[],
00186      const LocalOrdinal ldq_local,
00187      const Scalar R[],
00188      const LocalOrdinal ldr,
00189      MessengerBase<Scalar>* const messenger)
00190   {
00191     typedef Teuchos::ScalarTraits<Scalar> STS;
00192     typedef typename STS::magnitudeType magnitude_type;
00193     using std::make_pair;
00194     using std::pair;
00195     using std::vector;
00196 
00197     const magnitude_type ZERO (0);
00198     const magnitude_type ONE (1);
00199     BLAS<LocalOrdinal, Scalar> blas;
00200 
00201     //
00202     // Compute $\| I - Q^T * Q \|_F$
00203     //
00204 
00205     // Compute Q_local^T * Q_local (this node's component of Q^T*Q)
00206     vector<Scalar> Temp (ncols*ncols, STS::nan());
00207     const LocalOrdinal ld_temp = ncols;
00208 
00209     if (STS::isComplex)
00210       blas.GEMM ("C", "N", ncols, ncols, nrows_local, 
00211      ONE, Q_local, ldq_local, Q_local, ldq_local, 
00212      ZERO, &Temp[0], ld_temp);
00213     else
00214       blas.GEMM ("T", "N", ncols, ncols, nrows_local, 
00215      ONE, Q_local, ldq_local, Q_local, ldq_local, 
00216      ZERO, &Temp[0], ld_temp);
00217   
00218     // Reduce over all the processors to get the global Q^T*Q in Temp2.
00219     vector<Scalar> Temp2 (ncols*ncols, STS::nan());
00220     messenger->globalVectorSum (&Temp[0], &Temp2[0], ncols*ncols);
00221 
00222     // Compute I-(Q^T*Q) redundantly on all processors
00223     for (LocalOrdinal j = 0; j < ncols; j++)
00224       Temp2[j + j*ld_temp] = ONE - Temp2[j + j*ld_temp];
00225   
00226     // Compute the Frobenius norm of I - Q^T*Q, redundantly on all processors.
00227     const magnitude_type Orthog_F = 
00228       local_frobenius_norm (ncols, ncols, &Temp2[0], ld_temp);
00229 
00230     // Compute the Frobenius norm of A.  
00231     const magnitude_type A_F = 
00232       global_frobenius_norm (nrows_local, ncols, &A_local[0], lda_local, messenger);
00233 
00234     //
00235     // Compute $\| A - Q*R \|_F$
00236     //
00237 
00238     vector<Scalar> Resid (nrows_local * ncols, STS::nan());
00239     const LocalOrdinal ld_resid = nrows_local;
00240 
00241     // Resid := A (deep copy)
00242     copy_matrix (nrows_local, ncols, &Resid[0], ld_resid, A_local, lda_local);
00243 
00244     // Resid := Resid - Q*R
00245     blas.GEMM ("N", "N", nrows_local, ncols, ncols, 
00246          -ONE, Q_local, ldq_local, R, ldr, 
00247          ONE, &Resid[0], ld_resid);
00248 
00249     const magnitude_type Resid_F = 
00250       global_frobenius_norm (nrows_local, ncols, &Resid[0], ld_resid, messenger);
00251 
00252     vector<magnitude_type> results (3);
00253     results[0] = Resid_F;
00254     results[1] = Orthog_F;
00255     results[2] = A_F;
00256     return results;
00257   }
00258 
00259 } // namespace TSQR
00260 
00261 #endif // __TSQR_Tsqr_GlobalVerify_hpp
00262 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends