MOOCHO (Single Doxygen Collection) Version of the Day
NLPInterfacePack_NLPDirectTester.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) 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 // 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 Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include <assert.h>
00043 #include <math.h>
00044 
00045 #include <ostream>
00046 #include <iomanip>
00047 #include <sstream>
00048 #include <limits>
00049 
00050 #include "NLPInterfacePack_NLPDirectTester.hpp"
00051 #include "NLPInterfacePack_NLPDirect.hpp"
00052 #include "AbstractLinAlgPack_VectorMutable.hpp"
00053 #include "AbstractLinAlgPack_VectorOut.hpp"
00054 #include "AbstractLinAlgPack_VectorSpace.hpp"
00055 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00056 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00057 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00058 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00059 #include "TestingHelperPack_update_success.hpp"
00060 #include "Teuchos_Assert.hpp"
00061 
00062 namespace {
00063 template< class T >
00064 inline
00065 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
00066 } // end namespace
00067 
00068 namespace NLPInterfacePack {
00069 
00070 NLPDirectTester::NLPDirectTester(
00071   const calc_fd_prod_ptr_t  &calc_fd_prod
00072   ,ETestingMethod           Gf_testing_method
00073   ,ETestingMethod           Gc_testing_method
00074   ,value_type               Gf_warning_tol
00075   ,value_type               Gf_error_tol
00076   ,value_type               Gc_warning_tol
00077   ,value_type               Gc_error_tol
00078   ,size_type                num_fd_directions
00079   ,bool                     dump_all
00080   )
00081   :calc_fd_prod_(calc_fd_prod)
00082   ,Gf_testing_method_(Gf_testing_method)
00083   ,Gc_testing_method_(Gc_testing_method)
00084   ,Gf_warning_tol_(Gf_warning_tol)
00085   ,Gf_error_tol_(Gf_error_tol)
00086   ,Gc_warning_tol_(Gc_warning_tol)
00087   ,Gc_error_tol_(Gc_error_tol)
00088   ,num_fd_directions_(num_fd_directions)
00089   ,dump_all_(dump_all)
00090 {
00091   if(calc_fd_prod_ == Teuchos::null)
00092     calc_fd_prod_ = Teuchos::rcp(new CalcFiniteDiffProd());
00093 }
00094 
00095 bool NLPDirectTester::finite_diff_check(
00096   NLPDirect         *nlp
00097   ,const Vector     &xo
00098   ,const Vector     *xl
00099   ,const Vector     *xu
00100   ,const Vector     *c
00101   ,const Vector     *Gf
00102   ,const Vector     *py
00103   ,const Vector     *rGf
00104   ,const MatrixOp   *GcU
00105   ,const MatrixOp   *D
00106   ,const MatrixOp   *Uz
00107   ,bool             print_all_warnings
00108   ,std::ostream     *out
00109   ) const
00110 {
00111 
00112   using std::setw;
00113   using std::endl;
00114   using std::right;
00115 
00116   using AbstractLinAlgPack::sum;
00117   using AbstractLinAlgPack::dot;
00118   using AbstractLinAlgPack::Vp_StV;
00119   using AbstractLinAlgPack::random_vector;
00120   using AbstractLinAlgPack::assert_print_nan_inf;
00121   using LinAlgOpPack::V_StV;
00122   using LinAlgOpPack::V_StMtV;
00123   using LinAlgOpPack::Vp_MtV;
00124   using LinAlgOpPack::M_StM;
00125   using LinAlgOpPack::M_StMtM;
00126 
00127   typedef VectorSpace::vec_mut_ptr_t  vec_mut_ptr_t;
00128 
00129 //  using AbstractLinAlgPack::TestingPack::CompareDenseVectors;
00130 //  using AbstractLinAlgPack::TestingPack::CompareDenseSparseMatrices;
00131 
00132   using TestingHelperPack::update_success;
00133 
00134   bool success = true, preformed_fd;
00135   if(out) {
00136     *out << std::boolalpha
00137        << std::endl
00138        << "*********************************************************\n"
00139        << "*** NLPDirectTester::finite_diff_check(...) ***\n"
00140        << "*********************************************************\n";
00141   }
00142 
00143   const Range1D
00144     var_dep      = nlp->var_dep(),
00145     var_indep    = nlp->var_indep(),
00146     con_decomp   = nlp->con_decomp(),
00147     con_undecomp = nlp->con_undecomp();
00148   NLP::vec_space_ptr_t
00149     space_x = nlp->space_x(),
00150     space_c = nlp->space_c();
00151 
00152   // //////////////////////////////////////////////
00153   // Validate the input
00154 
00155   TEUCHOS_TEST_FOR_EXCEPTION(
00156     py && !c, std::invalid_argument
00157     ,"NLPDirectTester::finite_diff_check(...) : "
00158     "Error, if py!=NULL then c!=NULL must also be true!" );
00159 
00160   const CalcFiniteDiffProd
00161     &fd_deriv_prod = this->calc_fd_prod();
00162 
00163   const value_type
00164     rand_y_l = -1.0, rand_y_u = 1.0,
00165     small_num = ::sqrt(std::numeric_limits<value_type>::epsilon());
00166 
00167   try {
00168 
00169   // ///////////////////////////////////////////////
00170   // (1) Check Gf
00171 
00172   if(Gf) {
00173     switch( Gf_testing_method() ) {
00174       case FD_COMPUTE_ALL: {
00175         // Compute FDGf outright
00176         TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: update above!
00177         break;
00178       }
00179       case FD_DIRECTIONAL: {
00180         // Compute FDGF'*y using random y's
00181         if(out)
00182           *out
00183             << "\nComparing products Gf'*y with finite difference values FDGf'*y for "
00184             << "random y's ...\n";
00185         vec_mut_ptr_t y = space_x->create_member();
00186         value_type max_warning_viol = 0.0;
00187         int num_warning_viol = 0;
00188         const int num_fd_directions_used = ( num_fd_directions() > 0 ? num_fd_directions() : 1 );
00189         for( int direc_i = 1; direc_i <= num_fd_directions_used; ++direc_i ) {
00190           if( num_fd_directions() > 0 ) {
00191             random_vector( rand_y_l, rand_y_u, y.get() );
00192             if(out)
00193               *out
00194                 << "\n****"
00195                 << "\n**** Random directional vector " << direc_i << " ( ||y||_1 / n = "
00196                 << (y->norm_1() / y->dim()) << " )"
00197                 << "\n***\n";
00198           }
00199           else {
00200             *y = 1.0;
00201             if(out)
00202               *out
00203                 << "\n****"
00204                 << "\n**** Ones vector y ( ||y||_1 / n = "<<(y->norm_1()/y->dim())<<" )"
00205                 << "\n***\n";
00206           }
00207           value_type
00208             Gf_y = dot( *Gf, *y ),
00209             FDGf_y;
00210           preformed_fd = fd_deriv_prod.calc_deriv_product(
00211             xo,xl,xu
00212             ,*y,NULL,NULL,true,nlp,&FDGf_y,NULL,out,dump_all(),dump_all()
00213             );
00214           if( !preformed_fd )
00215             goto FD_NOT_PREFORMED;
00216           assert_print_nan_inf(FDGf_y, "FDGf'*y",true,out);
00217           const value_type
00218             calc_err = ::fabs( ( Gf_y - FDGf_y )/( ::fabs(Gf_y) + ::fabs(FDGf_y) + small_num ) );
00219           if( calc_err >= Gf_warning_tol() ) {
00220             max_warning_viol = my_max( max_warning_viol, calc_err );
00221             ++num_warning_viol;
00222           }
00223           if(out)
00224             *out
00225               << "\nrel_err(Gf'*y,FDGf'*y) = "
00226               << "rel_err(" << Gf_y << "," << FDGf_y << ") = "
00227               << calc_err << endl;
00228           if( calc_err >= Gf_error_tol() ) {
00229             if(out) {
00230               *out
00231                 << "Error, above relative error exceeded Gf_error_tol = " << Gf_error_tol() << endl;
00232               if(dump_all()) {
00233                 *out << "\ny =\n" << *y;
00234               }
00235             }
00236           }
00237         }
00238         if(out && num_warning_viol)
00239           *out
00240             << "\nThere were " << num_warning_viol << " warning tolerance "
00241             << "violations out of num_fd_directions = " << num_fd_directions()
00242             << " computations of FDGf'*y\n"
00243             << "and the maximum violation was " << max_warning_viol
00244             << " > Gf_waring_tol = " << Gf_warning_tol() << endl;
00245         break;
00246       }
00247       default:
00248         TEUCHOS_TEST_FOR_EXCEPT(true); // Invalid value
00249     }
00250   }
00251 
00252   // /////////////////////////////////////////////
00253   // (2) Check py = -inv(C)*c
00254   //
00255   // We want to check; 
00256   // 
00257   //  FDC * (inv(C)*c) \approx c
00258   //       \_________/
00259   //         -py
00260   //
00261   // We can compute this as:
00262   //           
00263   // FDC * py = [ FDC, FDN ] * [ -py ; 0 ]
00264   //            \__________/
00265   //                FDA'
00266   // 
00267   // t1 =  [ -py ; 0 ]
00268   // 
00269   // t2 = FDA'*t1
00270   // 
00271   // Compare t2 \approx c
00272   // 
00273   if(py) {
00274     if(out)
00275       *out
00276         << "\nComparing c with finite difference product FDA'*[ -py; 0 ] = -FDC*py ...\n";
00277     // t1 =  [ -py ; 0 ]
00278     VectorSpace::vec_mut_ptr_t
00279       t1 = space_x->create_member();
00280     V_StV( t1->sub_view(var_dep).get(), -1.0, *py );
00281     *t1->sub_view(var_indep) = 0.0;
00282     // t2 = FDA'*t1
00283     VectorSpace::vec_mut_ptr_t
00284       t2 = nlp->space_c()->create_member();
00285     preformed_fd = fd_deriv_prod.calc_deriv_product(
00286       xo,xl,xu
00287       ,*t1,NULL,NULL,true,nlp,NULL,t2.get(),out,dump_all(),dump_all()
00288       );
00289     if( !preformed_fd )
00290       goto FD_NOT_PREFORMED;
00291     const value_type
00292       sum_c  = sum(*c),
00293       sum_t2 = sum(*t2);
00294     assert_print_nan_inf(sum_t2, "sum(-FDC*py)",true,out);
00295     const value_type
00296       calc_err = ::fabs( ( sum_c - sum_t2 )/( ::fabs(sum_c) + ::fabs(sum_t2) + small_num ) );
00297     if(out)
00298       *out
00299         << "\nrel_err(sum(c),sum(-FDC*py)) = "
00300         << "rel_err(" << sum_c << "," << sum_t2 << ") = "
00301         << calc_err << endl;
00302     if( calc_err >= Gc_error_tol() ) {
00303       if(out)
00304         *out
00305           << "Error, above relative error exceeded Gc_error_tol = " << Gc_error_tol() << endl;
00306       if(print_all_warnings)
00307         *out << "\nt1 = [ -py; 0 ] =\n" << *t1
00308            << "\nt2 = FDA'*t1 = -FDC*py =\n"   << *t2;
00309       update_success( false, &success );
00310     }
00311     if( calc_err >= Gc_warning_tol() ) {
00312       if(out)
00313         *out
00314           << "\nWarning, above relative error exceeded Gc_warning_tol = " << Gc_warning_tol() << endl;
00315     }
00316   }
00317 
00318   // /////////////////////////////////////////////
00319   // (3) Check D = -inv(C)*N
00320 
00321   if(D) {
00322     switch( Gc_testing_method() ) {
00323       case FD_COMPUTE_ALL: {
00324         //
00325         // Compute FDN outright and check
00326         // -FDC * D \aprox FDN
00327         // 
00328         // We want to compute:
00329         // 
00330         // FDC * -D = [ FDC, FDN ] * [ -D; 0 ]
00331         //            \__________/
00332         //                FDA'
00333         // 
00334         // To compute the above we perform:
00335         // 
00336         // T = FDA' * [ -D; 0 ] (one column at a time)
00337         // 
00338         // Compare T \approx FDN
00339         //
00340 /*
00341         // FDN
00342         DMatrix FDN(m,n-m);
00343         fd_deriv_computer.calc_deriv( xo, xl, xu
00344           , Range1D(m+1,n), nlp, NULL
00345           , &FDN() ,BLAS_Cpp::trans, out );
00346 
00347         // T = FDA' * [ -D; 0 ] (one column at a time)
00348         DMatrix T(m,n-m);
00349         DVector t(n);
00350         t(m+1,n) = 0.0;
00351         for( int s = 1; s <= n-m; ++s ) {
00352           // t = [ -D(:,s); 0 ]
00353           V_StV( &t(1,m), -1.0, D->col(s) );
00354           // T(:,s) =  FDA' * t
00355           fd_deriv_prod.calc_deriv_product(
00356             xo,xl,xu,t(),NULL,NULL,nlp,NULL,&T.col(s),out);
00357         }        
00358 
00359         // Compare T \approx FDN
00360         if(out)
00361           *out
00362             << "\nChecking the computed D = -inv(C)*N\n"
00363             << "where D(i,j) = (-FDC*D)(i,j), dM(i,j) = FDN(i,j) ...\n";
00364         result = comp_M.comp(
00365           T(), FDN, BLAS_Cpp::no_trans
00366           , CompareDenseSparseMatrices::FULL_MATRIX
00367           , CompareDenseSparseMatrices::REL_ERR_BY_COL
00368           , Gc_warning_tol(), Gc_error_tol()
00369           , print_all_warnings, out );
00370         update_success( result, &success );
00371         if(!result) return false;
00372 */
00373         TEUCHOS_TEST_FOR_EXCEPT(true); // Todo: Implement above!
00374         break;
00375       }
00376       case FD_DIRECTIONAL: {
00377         //
00378         // Compute -FDC * D * v \aprox FDN * v
00379         // for random v's
00380         //
00381         // We will compute this as:
00382         // 
00383         // t1 = [ 0; y ] <: R^(n)
00384         // 
00385         // t2 = FDA' * t1  (  FDN * y ) <: R^(m)
00386         //
00387         // t1 = [ -D * y ; 0 ]  <: R^(n)
00388         // 
00389         // t3 = FDA' * t1  ( -FDC * D * y ) <: R^(m)
00390         // 
00391         // Compare t2 \approx t3
00392         //
00393         if(out)
00394           *out
00395             << "\nComparing finite difference products -FDC*D*y with FDN*y for "
00396               "random vectors y ...\n";
00397         VectorSpace::vec_mut_ptr_t
00398           y  = space_x->sub_space(var_indep)->create_member(),
00399           t1 = space_x->create_member(),
00400           t2 = space_c->create_member(),
00401           t3 = space_c->create_member();
00402         value_type max_warning_viol = 0.0;
00403         int num_warning_viol = 0;
00404         const int num_fd_directions_used = ( num_fd_directions() > 0 ? num_fd_directions() : 1 );
00405         for( int direc_i = 1; direc_i <= num_fd_directions_used; ++direc_i ) {
00406           if( num_fd_directions() > 0 ) {
00407             random_vector( rand_y_l, rand_y_u, y.get() );
00408             if(out)
00409               *out
00410                 << "\n****"
00411                 << "\n**** Random directional vector " << direc_i << " ( ||y||_1 / n = "
00412                 << (y->norm_1() / y->dim()) << " )"
00413                 << "\n***\n";
00414           }
00415           else {
00416             *y = 1.0;
00417             if(out)
00418               *out
00419                 << "\n****"
00420                 << "\n**** Ones vector y ( ||y||_1 / n = "<<(y->norm_1()/y->dim())<<" )"
00421                 << "\n***\n";
00422           }
00423           // t1 = [ 0; y ] <: R^(n)
00424           *t1->sub_view(var_dep)   = 0.0;
00425           *t1->sub_view(var_indep) = *y;
00426           // t2 = FDA' * t1  (  FDN * y ) <: R^(m)
00427           preformed_fd = fd_deriv_prod.calc_deriv_product(
00428             xo,xl,xu
00429             ,*t1,NULL,NULL,true,nlp,NULL,t2.get(),out,dump_all(),dump_all()
00430             );
00431           if( !preformed_fd )
00432             goto FD_NOT_PREFORMED;
00433           // t1 = [ -D * y ; 0 ]  <: R^(n)
00434           V_StMtV( t1->sub_view(var_dep).get(), -1.0, *D, BLAS_Cpp::no_trans, *y );
00435           *t1->sub_view(var_indep) = 0.0;
00436           // t3 = FDA' * t1  ( -FDC * D * y ) <: R^(m)
00437           preformed_fd = fd_deriv_prod.calc_deriv_product(
00438             xo,xl,xu
00439             ,*t1,NULL,NULL,true,nlp,NULL,t3.get(),out,dump_all(),dump_all()
00440             );
00441           // Compare t2 \approx t3
00442           const value_type
00443             sum_t2 = sum(*t2),
00444             sum_t3 = sum(*t3);
00445           const value_type
00446             calc_err = ::fabs( ( sum_t2 - sum_t3 )/( ::fabs(sum_t2) + ::fabs(sum_t3) + small_num ) );
00447           if(out)
00448             *out
00449               << "\nrel_err(sum(-FDC*D*y),sum(FDN*y)) = "
00450               << "rel_err(" << sum_t3 << "," << sum_t2 << ") = "
00451               << calc_err << endl;
00452           if( calc_err >= Gc_warning_tol() ) {
00453             max_warning_viol = my_max( max_warning_viol, calc_err );
00454             ++num_warning_viol;
00455           }
00456           if( calc_err >= Gc_error_tol() ) {
00457             if(out)
00458               *out
00459                 << "Error, above relative error exceeded Gc_error_tol = " << Gc_error_tol() << endl
00460                 << "Stoping the tests!\n";
00461             if(print_all_warnings)
00462               *out << "\ny =\n" << *y
00463                    << "\nt1 = [ -D*y; 0 ] =\n" << *t1
00464                    << "\nt2 =  FDA' * [ 0; y ] = FDN * y =\n" << *t2
00465                    << "\nt3 =  FDA' * t1 = -FDC * D * y =\n" << *t3;
00466             update_success( false, &success );
00467           }
00468         }
00469         if(out && num_warning_viol)
00470           *out
00471             << "\nThere were " << num_warning_viol << " warning tolerance "
00472             << "violations out of num_fd_directions = " << num_fd_directions()
00473             << " computations of sum(FDC*D*y) and sum(FDN*y)\n"
00474             << "and the maximum relative iolation was " << max_warning_viol
00475             << " > Gc_waring_tol = " << Gc_warning_tol() << endl;
00476         break;
00477       }
00478       default:
00479         TEUCHOS_TEST_FOR_EXCEPT(true);
00480     }
00481   }
00482 
00483   // ///////////////////////////////////////////////
00484   // (4) Check rGf = Gf(var_indep) + D'*Gf(var_dep)
00485 
00486   if(rGf) {
00487     if( Gf && D ) {
00488       if(out)
00489         *out
00490           << "\nComparing rGf_tmp = Gf(var_indep) - D'*Gf(var_dep) with rGf ...\n";
00491       VectorSpace::vec_mut_ptr_t
00492         rGf_tmp = space_x->sub_space(var_indep)->create_member();
00493       *rGf_tmp = *Gf->sub_view(var_indep);
00494       Vp_MtV( rGf_tmp.get(), *D, BLAS_Cpp::trans, *Gf->sub_view(var_dep) );
00495       const value_type
00496         sum_rGf_tmp  = sum(*rGf_tmp),
00497         sum_rGf      = sum(*rGf);
00498       const value_type
00499         calc_err = ::fabs( ( sum_rGf_tmp - sum_rGf )/( ::fabs(sum_rGf_tmp) + ::fabs(sum_rGf) + small_num ) );
00500       if(out)
00501         *out
00502           << "\nrel_err(sum(rGf_tmp),sum(rGf)) = "
00503           << "rel_err(" << sum_rGf_tmp << "," << sum_rGf << ") = "
00504           << calc_err << endl;
00505       if( calc_err >= Gc_error_tol() ) {
00506         if(out)
00507           *out
00508             << "Error, above relative error exceeded Gc_error_tol = " << Gc_error_tol() << endl;
00509         if(print_all_warnings)
00510           *out << "\nrGf_tmp =\n" << *rGf_tmp
00511              << "\nrGf =\n"   << *rGf;
00512         update_success( false, &success );
00513       }
00514       if( calc_err >= Gc_warning_tol() ) {
00515         if(out)
00516           *out
00517             << "\nWarning, above relative error exceeded Gc_warning_tol = "
00518             << Gc_warning_tol() << endl;
00519       }
00520     }
00521     else if( D ) {
00522       if(out)
00523         *out
00524           << "\nComparing rGf'*y with the finite difference product"
00525           << " fd_prod(f,[D*y;y]) for random vectors y ...\n";
00526       VectorSpace::vec_mut_ptr_t
00527         y  = space_x->sub_space(var_indep)->create_member(),
00528         t  = space_x->create_member();
00529       value_type max_warning_viol = 0.0;
00530       int num_warning_viol = 0;
00531       const int num_fd_directions_used = ( num_fd_directions() > 0 ? num_fd_directions() : 1 );
00532       for( int direc_i = 1; direc_i <= num_fd_directions_used; ++direc_i ) {
00533         if( num_fd_directions() > 0 ) {
00534           random_vector( rand_y_l, rand_y_u, y.get() );
00535           if(out)
00536             *out
00537               << "\n****"
00538               << "\n**** Random directional vector " << direc_i << " ( ||y||_1 / n = "
00539               << (y->norm_1() / y->dim()) << " )"
00540               << "\n***\n";
00541         }
00542         else {
00543           *y = 1.0;
00544           if(out)
00545             *out
00546               << "\n****"
00547               << "\n**** Ones vector y ( ||y||_1 / n = "<<(y->norm_1()/y->dim())<<" )"
00548               << "\n***\n";
00549         }
00550         // t = [ D*y; y ]
00551         LinAlgOpPack::V_MtV(&*t->sub_view(var_dep),*D,BLAS_Cpp::no_trans,*y);
00552         *t->sub_view(var_indep) = *y;
00553         value_type fd_rGf_y = 0.0;
00554         // fd_Gf_y
00555         preformed_fd = fd_deriv_prod.calc_deriv_product(
00556           xo,xl,xu
00557           ,*t,NULL,NULL,true,nlp,&fd_rGf_y,NULL,out,dump_all(),dump_all()
00558           );
00559         if( !preformed_fd )
00560           goto FD_NOT_PREFORMED;
00561         if(out) *out << "fd_prod(f,[D*y;y]) = " << fd_rGf_y << "\n";
00562         // rGf_y = rGf'*y
00563         const value_type rGf_y = dot(*rGf,*y);
00564         if(out) *out << "rGf'*y = " << rGf_y << "\n";
00565         // Compare fd_rGf_y and rGf*y
00566         const value_type
00567           calc_err = ::fabs( ( rGf_y - fd_rGf_y )/( ::fabs(rGf_y) + ::fabs(fd_rGf_y) + small_num ) );
00568         if( calc_err >= Gc_warning_tol() ) {
00569           max_warning_viol = my_max( max_warning_viol, calc_err );
00570           ++num_warning_viol;
00571         }
00572         if(out)
00573           *out
00574             << "\nrel_err(rGf'*y,fd_prod(f,[D*y;y])) = "
00575             << "rel_err(" << fd_rGf_y << "," << rGf_y << ") = "
00576             << calc_err << endl;
00577         if( calc_err >= Gf_error_tol() ) {
00578           if(out)
00579             *out << "Error, above relative error exceeded Gc_error_tol = " << Gc_error_tol() << endl;
00580           if(print_all_warnings)
00581             *out << "\ny =\n" << *y
00582                  << "\nt = [ D*y; y ] =\n" << *t;
00583           update_success( false, &success );
00584         }
00585       }
00586     }
00587     else {
00588       TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Test rGf without D? (This is not going to be easy!)
00589     }
00590   }
00591   
00592   // ///////////////////////////////////////////////////
00593   // (5) Check GcU, and/or Uz (for undecomposed equalities)
00594 
00595   if(GcU || Uz) {
00596     TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement!
00597   }
00598   
00599 FD_NOT_PREFORMED:
00600 
00601   if(!preformed_fd) {
00602     if(out)
00603       *out
00604         << "\nError, the finite difference computation was not preformed due to cramped bounds\n"
00605         << "Finite difference test failed!\n" << endl;
00606     return false;
00607   }
00608 
00609   } // end try
00610   catch( const AbstractLinAlgPack::NaNInfException& except ) {
00611     if(out)
00612       *out
00613         << "Error, found a NaN or Inf.  Stoping tests\n";
00614     success = false;
00615   }
00616   
00617   if( out ) {
00618     if( success )
00619       *out
00620         << "\nCongradulations, all the finite difference errors where within the\n"
00621         "specified error tolerances!\n";
00622     else
00623       *out
00624         << "\nOh no, at least one of the above finite difference tests failed!\n";
00625   }
00626 
00627   return success;
00628 
00629 }
00630 
00631 }  // end namespace NLPInterfacePack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines