Kokkos Node API and Local Linear Algebra Kernels Version of the Day
Tsqr_LocalVerify.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_LocalVerify_hpp
00043 #define __TSQR_Tsqr_LocalVerify_hpp
00044 
00045 #include <Tsqr_ScalarTraits.hpp>
00046 #include <Tsqr_Blas.hpp>
00047 #include <Tsqr_Util.hpp>
00048 
00049 #include <cmath>
00050 #include <limits>
00051 #include <utility> // std::pair, std::make_pair
00052 #include <vector>
00053 
00056 
00057 namespace TSQR {
00058 
00059   template< class Ordinal, class Scalar >
00060   typename ScalarTraits< Scalar >::magnitude_type
00061   local_frobenius_norm (const Ordinal nrows_local, 
00062       const Ordinal ncols,
00063       const Scalar  A_local[],
00064       const Ordinal lda_local)
00065   {
00066     typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00067     
00068     // FIXME (mfh 22 Apr 2010) This function does no scaling of
00069     // intermediate quantities, so it might overflow unnecessarily.
00070     magnitude_type result (0);
00071     for (Ordinal j = 0; j < ncols; j++)
00072       {
00073   const Scalar* const cur_col = &A_local[j*lda_local];
00074   for (Ordinal i = 0; i < nrows_local; ++i)
00075     {
00076       const magnitude_type abs_xi = ScalarTraits< Scalar >::abs (cur_col[i]);
00077       result = result + abs_xi * abs_xi;
00078     }
00079       }
00080     return sqrt (result);
00081   }
00082 
00083 
00084   template< class Ordinal, class Scalar >
00085   bool
00086   NaN_in_matrix (const Ordinal nrows, 
00087      const Ordinal ncols, 
00088      const Scalar A[], 
00089      const Ordinal lda)
00090   {
00091     // Testing whether a NaN is present in A only makes sense if it is
00092     // possible for NaNs not to signal.  Otherwise the NaNs would have
00093     // signalled and we wouldn't need to be here.  Of course perhaps
00094     // one could change the signal state at runtime, but has_quiet_NaN
00095     // refers to the possibility of quiet NaNs being able to exist at
00096     // all.
00097     if (std::numeric_limits<Scalar>::has_quiet_NaN)
00098       {
00099   for (Ordinal j = 0; j < ncols; j++)
00100     for (Ordinal i = 0; i < nrows; i++)
00101       {
00102         if (std::isnan (A[i + j*lda]))
00103     return true;
00104       }
00105   return false;
00106       }
00107     else
00108       return false;
00109   }
00110 
00111 
00112   template< class Ordinal, class Scalar >
00113   bool
00114   NaN_in_matrix (const Ordinal nrows, 
00115      const Ordinal ncols, 
00116      const std::vector<Scalar>& A, 
00117      const Ordinal lda)
00118   {
00119     const Scalar* const A_ptr = &A[0];
00120     return NaN_in_matrix (nrows, ncols, A_ptr, lda);
00121   }
00122 
00123 
00124 
00125   template< class Ordinal, class Scalar >
00126   typename ScalarTraits< Scalar >::magnitude_type
00127   localOrthogonality (const Ordinal nrows,
00128           const Ordinal ncols,
00129           const Scalar Q[],
00130           const Ordinal ldq)
00131   {
00132     typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00133     const Scalar ZERO (0);
00134     const Scalar ONE (1);
00135 
00136     BLAS<Ordinal, Scalar> blas;
00137   
00138     std::vector<Scalar> AbsOrthog (ncols * ncols, std::numeric_limits<Scalar>::quiet_NaN());
00139     const Ordinal AbsOrthog_stride = ncols;
00140 
00141     // Compute AbsOrthog := Q' * Q - I.  First, compute Q' * Q:
00142     if (ScalarTraits< Scalar >::is_complex)
00143       blas.GEMM ("C", "N", ncols, ncols, nrows, ONE, Q, ldq, Q, ldq, 
00144      ZERO, &AbsOrthog[0], AbsOrthog_stride);
00145     else
00146       blas.GEMM ("T", "N", ncols, ncols, nrows, ONE, Q, ldq, Q, ldq, 
00147      ZERO, &AbsOrthog[0], AbsOrthog_stride);
00148 
00149     // Now, compute (Q^T*Q) - I.
00150     for (Ordinal j = 0; j < ncols; j++)
00151       AbsOrthog[j + j*AbsOrthog_stride] = AbsOrthog[j + j*AbsOrthog_stride] - ONE;
00152 
00153     // Now AbsOrthog == Q^T * Q - I.  Compute and return its Frobenius norm.
00154     return local_frobenius_norm (ncols, ncols, &AbsOrthog[0], AbsOrthog_stride);
00155   }
00156 
00157 
00158   
00159   template< class Ordinal, class Scalar >
00160   typename ScalarTraits< Scalar >::magnitude_type
00161   local_relative_orthogonality (const Ordinal nrows,
00162         const Ordinal ncols,
00163         const Scalar Q[],
00164         const Ordinal ldq,
00165         const typename ScalarTraits< Scalar >::magnitude_type A_norm_F)
00166   {
00167     typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00168     const Scalar ZERO (0);
00169     const Scalar ONE (1);
00170 
00171     const bool relative = false; // whether to scale $\|I-Q^T*Q\|_F$ by $\|A\|_F$
00172     BLAS<Ordinal, Scalar> blas;
00173   
00174     std::vector<Scalar> AbsOrthog (ncols * ncols, std::numeric_limits<Scalar>::quiet_NaN());
00175     const Ordinal AbsOrthog_stride = ncols;
00176 
00177     // Compute AbsOrthog := Q' * Q - I.  First, compute Q' * Q:
00178     if (ScalarTraits< Scalar >::is_complex)
00179       blas.GEMM ("C", "N", ncols, ncols, nrows, ONE, Q, ldq, Q, ldq, 
00180      ZERO, &AbsOrthog[0], AbsOrthog_stride);
00181     else
00182       blas.GEMM ("T", "N", ncols, ncols, nrows, ONE, Q, ldq, Q, ldq, 
00183      ZERO, &AbsOrthog[0], AbsOrthog_stride);
00184 
00185     // Now, compute (Q^T*Q) - I.
00186     for (Ordinal j = 0; j < ncols; j++)
00187       AbsOrthog[j + j*AbsOrthog_stride] = AbsOrthog[j + j*AbsOrthog_stride] - ONE;
00188 
00189     // Now AbsOrthog == Q^T * Q - I.  Compute its Frobenius norm.
00190     const magnitude_type AbsOrthog_norm_F = 
00191       local_frobenius_norm (ncols, ncols, &AbsOrthog[0], AbsOrthog_stride);
00192 
00193     // Return the orthogonality measure
00194     return relative ? (AbsOrthog_norm_F / A_norm_F) : AbsOrthog_norm_F;
00195   }
00196 
00197 
00198   template< class Ordinal, class Scalar >
00199   typename ScalarTraits< Scalar >::magnitude_type
00200   localResidual (const Ordinal nrows, 
00201      const Ordinal ncols,
00202      const Scalar A[],
00203      const Ordinal lda,
00204      const Scalar Q[],
00205      const Ordinal ldq,
00206      const Scalar R[],
00207      const Ordinal ldr)
00208   {
00209     typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00210 
00211     std::vector< Scalar > AbsResid (nrows * ncols, std::numeric_limits< Scalar >::quiet_NaN());
00212     const Ordinal AbsResid_stride = nrows;
00213     BLAS< Ordinal, Scalar > blas;
00214     const magnitude_type ONE (1);
00215 
00216     // A_copy := A_copy - Q * R
00217     copy_matrix (nrows, ncols, &AbsResid[0], AbsResid_stride, A, lda);
00218     blas.GEMM ("N", "N", nrows, ncols, ncols, -ONE, Q, ldq, R, ldr, 
00219          ONE, &AbsResid[0], AbsResid_stride);
00220 
00221     return local_frobenius_norm (nrows, ncols, &AbsResid[0], AbsResid_stride);
00222   }
00223 
00224 
00225   template< class Ordinal, class Scalar >
00226   typename ScalarTraits< Scalar >::magnitude_type
00227   local_relative_residual (const Ordinal nrows, 
00228          const Ordinal ncols,
00229          const Scalar A[],
00230          const Ordinal lda,
00231          const Scalar Q[],
00232          const Ordinal ldq,
00233          const Scalar R[],
00234          const Ordinal ldr,
00235          const typename ScalarTraits< Scalar >::magnitude_type A_norm_F)
00236   {
00237     typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00238 
00239     std::vector< Scalar > AbsResid (nrows * ncols, std::numeric_limits< Scalar >::quiet_NaN());
00240     const Ordinal AbsResid_stride = nrows;
00241     BLAS< Ordinal, Scalar > blas;
00242     const magnitude_type ONE (1);
00243 
00244     // if (b_debug)
00245     //   cerr << "relative_residual:" << endl;
00246     // if (matrix_contains_nan (nrows, ncols, A, lda))
00247     //   cerr << "relative_residual: matrix A contains a NaN" << endl;
00248     // if (matrix_contains_nan (nrows, ncols, Q, ldq))
00249     //   cerr << "relative_residual: matrix Q contains a NaN" << endl;
00250     // if (matrix_contains_nan (ncols, ncols, R, ldr))
00251     //   cerr << "relative_residual: matrix R contains a NaN" << endl;
00252 
00253     // A_copy := A_copy - Q * R
00254     copy_matrix (nrows, ncols, &AbsResid[0], AbsResid_stride, A, lda);
00255 
00256     // if (NaN_in_matrix (nrows, ncols, AbsResid, AbsResid_stride))
00257     //   cerr << "relative_residual: matrix AbsResid := A contains a NaN" << endl;
00258 
00259     blas.GEMM ("N", "N", nrows, ncols, ncols, -ONE, Q, ldq, R, ldr, 
00260          ONE, &AbsResid[0], AbsResid_stride);
00261 
00262     // if (NaN_in_matrix (nrows, ncols, AbsResid, AbsResid_stride))
00263     //   cerr << "relative_residual: matrix AbsResid := A - Q*R contains a NaN" << endl;
00264 
00265     const magnitude_type absolute_residual = 
00266       local_frobenius_norm (nrows, ncols, &AbsResid[0], AbsResid_stride);
00267 
00268     // if (b_debug)
00269     //   {
00270     //     cerr << "In relative_residual:" << endl;
00271     //     cerr << "||Q||_2 = " << matrix_2norm(nrows, ncols, Q, ldq) << endl;
00272     //     cerr << "||R||_2 = " << matrix_2norm(ncols, ncols, R, ldr) << endl;
00273     //     cerr << "||A - QR||_2 = " << absolute_residual << endl;
00274     //   }
00275 
00276     return absolute_residual / A_norm_F;
00277   }
00278 
00310   template< class Ordinal, class Scalar >
00311   std::vector< typename ScalarTraits< Scalar >::magnitude_type >
00312   local_verify (const Ordinal nrows, 
00313     const Ordinal ncols, 
00314     const Scalar* const A, 
00315     const Ordinal lda,
00316     const Scalar* const Q, 
00317     const Ordinal ldq, 
00318     const Scalar* const R, 
00319     const Ordinal ldr)
00320   {
00321     typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00322     std::vector< magnitude_type > results(3);
00323     // const bool A_contains_NaN = NaN_in_matrix (nrows, ncols, A, lda);
00324     // const bool Q_contains_NaN = NaN_in_matrix (nrows, ncols, Q, ldq);
00325     // const bool R_contains_NaN = NaN_in_matrix (ncols, ncols, R, ldr);
00326 
00327     results[0] = localResidual (nrows, ncols, A, lda, Q, ldq, R, ldr);
00328     results[1] = localOrthogonality (nrows, ncols, Q, ldq);
00329     results[2] = local_frobenius_norm (nrows, ncols, A, lda);
00330 
00331     return results;
00332   }
00333 
00334 } // namespace TSQR
00335 
00336 #endif // __TSQR_Tsqr_LocalVerify_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends