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 (2009) 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 
00112   template< class Ordinal, class Scalar >
00113   typename ScalarTraits< Scalar >::magnitude_type
00114   localOrthogonality (const Ordinal nrows,
00115           const Ordinal ncols,
00116           const Scalar Q[],
00117           const Ordinal ldq)
00118   {
00119     typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00120     const Scalar ZERO (0);
00121     const Scalar ONE (1);
00122 
00123     BLAS<Ordinal, Scalar> blas;
00124   
00125     std::vector<Scalar> AbsOrthog (ncols * ncols, std::numeric_limits<Scalar>::quiet_NaN());
00126     const Ordinal AbsOrthog_stride = ncols;
00127 
00128     // Compute AbsOrthog := Q' * Q - I.  First, compute Q' * Q:
00129     if (ScalarTraits< Scalar >::is_complex)
00130       blas.GEMM ("C", "N", ncols, ncols, nrows, ONE, Q, ldq, Q, ldq, 
00131      ZERO, &AbsOrthog[0], AbsOrthog_stride);
00132     else
00133       blas.GEMM ("T", "N", ncols, ncols, nrows, ONE, Q, ldq, Q, ldq, 
00134      ZERO, &AbsOrthog[0], AbsOrthog_stride);
00135 
00136     // Now, compute (Q^T*Q) - I.
00137     for (Ordinal j = 0; j < ncols; j++)
00138       AbsOrthog[j + j*AbsOrthog_stride] = AbsOrthog[j + j*AbsOrthog_stride] - ONE;
00139 
00140     // Now AbsOrthog == Q^T * Q - I.  Compute and return its Frobenius norm.
00141     return local_frobenius_norm (ncols, ncols, &AbsOrthog[0], AbsOrthog_stride);
00142   }
00143 
00144 
00145   
00146   template< class Ordinal, class Scalar >
00147   typename ScalarTraits< Scalar >::magnitude_type
00148   local_relative_orthogonality (const Ordinal nrows,
00149         const Ordinal ncols,
00150         const Scalar Q[],
00151         const Ordinal ldq,
00152         const typename ScalarTraits< Scalar >::magnitude_type A_norm_F)
00153   {
00154     typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00155     const Scalar ZERO (0);
00156     const Scalar ONE (1);
00157 
00158     const bool relative = false; // whether to scale $\|I-Q^T*Q\|_F$ by $\|A\|_F$
00159     BLAS<Ordinal, Scalar> blas;
00160   
00161     std::vector<Scalar> AbsOrthog (ncols * ncols, std::numeric_limits<Scalar>::quiet_NaN());
00162     const Ordinal AbsOrthog_stride = ncols;
00163 
00164     // Compute AbsOrthog := Q' * Q - I.  First, compute Q' * Q:
00165     if (ScalarTraits< Scalar >::is_complex)
00166       blas.GEMM ("C", "N", ncols, ncols, nrows, ONE, Q, ldq, Q, ldq, 
00167      ZERO, &AbsOrthog[0], AbsOrthog_stride);
00168     else
00169       blas.GEMM ("T", "N", ncols, ncols, nrows, ONE, Q, ldq, Q, ldq, 
00170      ZERO, &AbsOrthog[0], AbsOrthog_stride);
00171 
00172     // Now, compute (Q^T*Q) - I.
00173     for (Ordinal j = 0; j < ncols; j++)
00174       AbsOrthog[j + j*AbsOrthog_stride] = AbsOrthog[j + j*AbsOrthog_stride] - ONE;
00175 
00176     // Now AbsOrthog == Q^T * Q - I.  Compute its Frobenius norm.
00177     const magnitude_type AbsOrthog_norm_F = 
00178       local_frobenius_norm (ncols, ncols, &AbsOrthog[0], AbsOrthog_stride);
00179 
00180     // Return the orthogonality measure
00181     return relative ? (AbsOrthog_norm_F / A_norm_F) : AbsOrthog_norm_F;
00182   }
00183 
00184 
00185   template< class Ordinal, class Scalar >
00186   typename ScalarTraits< Scalar >::magnitude_type
00187   localResidual (const Ordinal nrows, 
00188      const Ordinal ncols,
00189      const Scalar A[],
00190      const Ordinal lda,
00191      const Scalar Q[],
00192      const Ordinal ldq,
00193      const Scalar R[],
00194      const Ordinal ldr)
00195   {
00196     typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00197 
00198     std::vector< Scalar > AbsResid (nrows * ncols, std::numeric_limits< Scalar >::quiet_NaN());
00199     const Ordinal AbsResid_stride = nrows;
00200     BLAS< Ordinal, Scalar > blas;
00201     const magnitude_type ONE (1);
00202 
00203     // A_copy := A_copy - Q * R
00204     copy_matrix (nrows, ncols, &AbsResid[0], AbsResid_stride, A, lda);
00205     blas.GEMM ("N", "N", nrows, ncols, ncols, -ONE, Q, ldq, R, ldr, 
00206          ONE, &AbsResid[0], AbsResid_stride);
00207 
00208     return local_frobenius_norm (nrows, ncols, &AbsResid[0], AbsResid_stride);
00209   }
00210 
00211 
00212   template< class Ordinal, class Scalar >
00213   typename ScalarTraits< Scalar >::magnitude_type
00214   local_relative_residual (const Ordinal nrows, 
00215          const Ordinal ncols,
00216          const Scalar A[],
00217          const Ordinal lda,
00218          const Scalar Q[],
00219          const Ordinal ldq,
00220          const Scalar R[],
00221          const Ordinal ldr,
00222          const typename ScalarTraits< Scalar >::magnitude_type A_norm_F)
00223   {
00224     typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00225 
00226     std::vector< Scalar > AbsResid (nrows * ncols, std::numeric_limits< Scalar >::quiet_NaN());
00227     const Ordinal AbsResid_stride = nrows;
00228     BLAS< Ordinal, Scalar > blas;
00229     const magnitude_type ONE (1);
00230 
00231     // if (b_debug)
00232     //   cerr << "relative_residual:" << endl;
00233     // if (matrix_contains_nan (nrows, ncols, A, lda))
00234     //   cerr << "relative_residual: matrix A contains a NaN" << endl;
00235     // if (matrix_contains_nan (nrows, ncols, Q, ldq))
00236     //   cerr << "relative_residual: matrix Q contains a NaN" << endl;
00237     // if (matrix_contains_nan (ncols, ncols, R, ldr))
00238     //   cerr << "relative_residual: matrix R contains a NaN" << endl;
00239 
00240     // A_copy := A_copy - Q * R
00241     copy_matrix (nrows, ncols, &AbsResid[0], AbsResid_stride, A, lda);
00242 
00243     // if (NaN_in_matrix (nrows, ncols, AbsResid, AbsResid_stride))
00244     //   cerr << "relative_residual: matrix AbsResid := A contains a NaN" << endl;
00245 
00246     blas.GEMM ("N", "N", nrows, ncols, ncols, -ONE, Q, ldq, R, ldr, 
00247          ONE, &AbsResid[0], AbsResid_stride);
00248 
00249     // if (NaN_in_matrix (nrows, ncols, AbsResid, AbsResid_stride))
00250     //   cerr << "relative_residual: matrix AbsResid := A - Q*R contains a NaN" << endl;
00251 
00252     const magnitude_type absolute_residual = 
00253       local_frobenius_norm (nrows, ncols, &AbsResid[0], AbsResid_stride);
00254 
00255     // if (b_debug)
00256     //   {
00257     //     cerr << "In relative_residual:" << endl;
00258     //     cerr << "||Q||_2 = " << matrix_2norm(nrows, ncols, Q, ldq) << endl;
00259     //     cerr << "||R||_2 = " << matrix_2norm(ncols, ncols, R, ldr) << endl;
00260     //     cerr << "||A - QR||_2 = " << absolute_residual << endl;
00261     //   }
00262 
00263     return absolute_residual / A_norm_F;
00264   }
00265 
00297   template< class Ordinal, class Scalar >
00298   std::vector< typename ScalarTraits< Scalar >::magnitude_type >
00299   local_verify (const Ordinal nrows, 
00300     const Ordinal ncols, 
00301     const Scalar* const A, 
00302     const Ordinal lda,
00303     const Scalar* const Q, 
00304     const Ordinal ldq, 
00305     const Scalar* const R, 
00306     const Ordinal ldr)
00307   {
00308     typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00309     std::vector< magnitude_type > results(3);
00310     // const bool A_contains_NaN = NaN_in_matrix (nrows, ncols, A, lda);
00311     // const bool Q_contains_NaN = NaN_in_matrix (nrows, ncols, Q, ldq);
00312     // const bool R_contains_NaN = NaN_in_matrix (ncols, ncols, R, ldr);
00313 
00314     results[0] = localResidual (nrows, ncols, A, lda, Q, ldq, R, ldr);
00315     results[1] = localOrthogonality (nrows, ncols, Q, ldq);
00316     results[2] = local_frobenius_norm (nrows, ncols, A, lda);
00317 
00318     return results;
00319   }
00320 
00321 } // namespace TSQR
00322 
00323 #endif // __TSQR_Tsqr_LocalVerify_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends