Anasazi Version of the Day
Tsqr_LocalVerify.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_LocalVerify_hpp
00030 #define __TSQR_Tsqr_LocalVerify_hpp
00031 
00032 #include <Tsqr_ScalarTraits.hpp>
00033 #include <Tsqr_Blas.hpp>
00034 #include <Tsqr_Util.hpp>
00035 
00036 #include <cmath>
00037 #include <limits>
00038 #include <utility> // std::pair, std::make_pair
00039 #include <vector>
00040 
00043 
00044 namespace TSQR {
00045 
00046   template< class Ordinal, class Scalar >
00047   typename ScalarTraits< Scalar >::magnitude_type
00048   local_frobenius_norm (const Ordinal nrows_local, 
00049       const Ordinal ncols,
00050       const Scalar  A_local[],
00051       const Ordinal lda_local)
00052   {
00053     typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00054     
00055     // FIXME (mfh 22 Apr 2010) This function does no scaling of
00056     // intermediate quantities, so it might overflow unnecessarily.
00057     magnitude_type result (0);
00058     for (Ordinal j = 0; j < ncols; j++)
00059       {
00060   const Scalar* const cur_col = &A_local[j*lda_local];
00061   for (Ordinal i = 0; i < nrows_local; ++i)
00062     {
00063       const magnitude_type abs_xi = ScalarTraits< Scalar >::abs (cur_col[i]);
00064       result = result + abs_xi * abs_xi;
00065     }
00066       }
00067     return sqrt (result);
00068   }
00069 
00070 
00071   template< class Ordinal, class Scalar >
00072   bool
00073   NaN_in_matrix (const Ordinal nrows, 
00074      const Ordinal ncols, 
00075      const Scalar A[], 
00076      const Ordinal lda)
00077   {
00078     // Testing whether a NaN is present in A only makes sense if it is
00079     // possible for NaNs not to signal.  Otherwise the NaNs would have
00080     // signalled and we wouldn't need to be here.  Of course perhaps
00081     // one could change the signal state at runtime, but has_quiet_NaN
00082     // refers to the possibility of quiet NaNs being able to exist at
00083     // all.
00084     if (std::numeric_limits<Scalar>::has_quiet_NaN)
00085       {
00086   for (Ordinal j = 0; j < ncols; j++)
00087     for (Ordinal i = 0; i < nrows; i++)
00088       {
00089         if (std::isnan (A[i + j*lda]))
00090     return true;
00091       }
00092   return false;
00093       }
00094     else
00095       return false;
00096   }
00097 
00098 
00099   template< class Ordinal, class Scalar >
00100   bool
00101   NaN_in_matrix (const Ordinal nrows, 
00102      const Ordinal ncols, 
00103      const std::vector<Scalar>& A, 
00104      const Ordinal lda)
00105   {
00106     const Scalar* const A_ptr = &A[0];
00107     return NaN_in_matrix (nrows, ncols, A_ptr, lda);
00108   }
00109 
00110   
00111   template< class Ordinal, class Scalar >
00112   typename ScalarTraits< Scalar >::magnitude_type
00113   local_relative_orthogonality (const Ordinal nrows,
00114         const Ordinal ncols,
00115         const Scalar Q[],
00116         const Ordinal ldq,
00117         const typename ScalarTraits< Scalar >::magnitude_type A_norm_F)
00118   {
00119     typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00120     const Scalar ZERO (0);
00121     const Scalar ONE (1);
00122 
00123     const bool relative = false; // whether to scale $\|I-Q^T*Q\|_F$ by $\|A\|_F$
00124     BLAS<Ordinal, Scalar> blas;
00125   
00126     std::vector<Scalar> AbsOrthog (ncols * ncols, std::numeric_limits<Scalar>::quiet_NaN());
00127     const Ordinal AbsOrthog_stride = ncols;
00128 
00129     // Compute AbsOrthog := Q' * Q - I.  First, compute Q' * Q:
00130     if (ScalarTraits< Scalar >::is_complex)
00131       blas.GEMM ("C", "N", ncols, ncols, nrows, ONE, Q, ldq, Q, ldq, 
00132      ZERO, &AbsOrthog[0], AbsOrthog_stride);
00133     else
00134       blas.GEMM ("T", "N", ncols, ncols, nrows, ONE, Q, ldq, Q, ldq, 
00135      ZERO, &AbsOrthog[0], AbsOrthog_stride);
00136 
00137     // Now, compute (Q^T*Q) - I.
00138     for (Ordinal j = 0; j < ncols; j++)
00139       AbsOrthog[j + j*AbsOrthog_stride] = AbsOrthog[j + j*AbsOrthog_stride] - ONE;
00140 
00141     // Now AbsOrthog == Q^T * Q - I.  Compute its Frobenius norm.
00142     const magnitude_type AbsOrthog_norm_F = 
00143       local_frobenius_norm (ncols, ncols, &AbsOrthog[0], AbsOrthog_stride);
00144 
00145     // Return the orthogonality measure
00146     return relative ? (AbsOrthog_norm_F / A_norm_F) : AbsOrthog_norm_F;
00147   }
00148 
00149 
00150   template< class Ordinal, class Scalar >
00151   typename ScalarTraits< Scalar >::magnitude_type
00152   local_relative_residual (const Ordinal nrows, 
00153          const Ordinal ncols,
00154          const Scalar A[],
00155          const Ordinal lda,
00156          const Scalar Q[],
00157          const Ordinal ldq,
00158          const Scalar R[],
00159          const Ordinal ldr,
00160          const typename ScalarTraits< Scalar >::magnitude_type A_norm_F)
00161   {
00162     typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00163 
00164     std::vector< Scalar > AbsResid (nrows * ncols, std::numeric_limits< Scalar >::quiet_NaN());
00165     const Ordinal AbsResid_stride = nrows;
00166     BLAS< Ordinal, Scalar > blas;
00167     const magnitude_type ONE (1);
00168 
00169     // if (b_debug)
00170     //   cerr << "relative_residual:" << endl;
00171     // if (matrix_contains_nan (nrows, ncols, A, lda))
00172     //   cerr << "relative_residual: matrix A contains a NaN" << endl;
00173     // if (matrix_contains_nan (nrows, ncols, Q, ldq))
00174     //   cerr << "relative_residual: matrix Q contains a NaN" << endl;
00175     // if (matrix_contains_nan (ncols, ncols, R, ldr))
00176     //   cerr << "relative_residual: matrix R contains a NaN" << endl;
00177 
00178     // A_copy := A_copy - Q * R
00179     copy_matrix (nrows, ncols, &AbsResid[0], AbsResid_stride, A, lda);
00180 
00181     // if (NaN_in_matrix (nrows, ncols, AbsResid, AbsResid_stride))
00182     //   cerr << "relative_residual: matrix AbsResid := A contains a NaN" << endl;
00183 
00184     blas.GEMM ("N", "N", nrows, ncols, ncols, -ONE, Q, ldq, R, ldr, 
00185          ONE, &AbsResid[0], AbsResid_stride);
00186 
00187     // if (NaN_in_matrix (nrows, ncols, AbsResid, AbsResid_stride))
00188     //   cerr << "relative_residual: matrix AbsResid := A - Q*R contains a NaN" << endl;
00189 
00190     const magnitude_type absolute_residual = 
00191       local_frobenius_norm (nrows, ncols, &AbsResid[0], AbsResid_stride);
00192 
00193     // if (b_debug)
00194     //   {
00195     //     cerr << "In relative_residual:" << endl;
00196     //     cerr << "||Q||_2 = " << matrix_2norm(nrows, ncols, Q, ldq) << endl;
00197     //     cerr << "||R||_2 = " << matrix_2norm(ncols, ncols, R, ldr) << endl;
00198     //     cerr << "||A - QR||_2 = " << absolute_residual << endl;
00199     //   }
00200 
00201     return absolute_residual / A_norm_F;
00202   }
00203 
00234   template< class Ordinal, class Scalar >
00235   std::pair< typename ScalarTraits< Scalar >::magnitude_type, typename ScalarTraits< Scalar >::magnitude_type >
00236   local_verify (const Ordinal nrows, 
00237     const Ordinal ncols, 
00238     const Scalar* const A, 
00239     const Ordinal lda,
00240     const Scalar* const Q, 
00241     const Ordinal ldq, 
00242     const Scalar* const R, 
00243     const Ordinal ldr)
00244   {
00245     typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00246     // const bool A_contains_NaN = NaN_in_matrix (nrows, ncols, A, lda);
00247     // const bool Q_contains_NaN = NaN_in_matrix (nrows, ncols, Q, ldq);
00248     // const bool R_contains_NaN = NaN_in_matrix (ncols, ncols, R, ldr);
00249     
00250     const magnitude_type A_norm_F = local_frobenius_norm (nrows, ncols, A, lda);
00251     const magnitude_type rel_orthog = local_relative_orthogonality (nrows, ncols, Q, ldq, A_norm_F);
00252     const magnitude_type rel_resid = local_relative_residual (nrows, ncols, A, lda, Q, ldq, R, ldr, A_norm_F);
00253     return std::make_pair (rel_resid, rel_orthog);
00254   }
00255 
00256 } // namespace TSQR
00257 
00258 #endif // __TSQR_Tsqr_LocalVerify_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends