EpetraExt Development
EpetraExt_ModelEvaluatorScalingTools.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2011) 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 
00043 #include "EpetraExt_ModelEvaluatorScalingTools.h"
00044 #include "Teuchos_implicit_cast.hpp"
00045 #include "Epetra_RowMatrix.h"
00046 
00047 //
00048 // Here in the implementation of scaling we write all scaling routines to
00049 // specifically and individually handle every input and output object and we
00050 // transfer all objects from one [In,Out]Arg container to another to make sure
00051 // that that object has been specifically addressed.  We could just use the
00052 // [In,Out]Args::setArgs(...) function to copy all arguments be default but
00053 // then that would setup subtle errors where a new quality could be silently
00054 // ignored and not be scaled correctly.  Therefore, we feel it is better to
00055 // have to deal with *every* input and output object specifically with respect
00056 // to scaling in order to avoid overlooking scaling.  In this way, if a new
00057 // input or output object is added to [In,Out]Args but the code in this file
00058 // is not updated, then that object will not be passed through and some type
00059 // of error will be generated right away.  We feel this is the best behavior
00060 // and it justifies having every scaling-related function take both an input
00061 // and an output [In,Out]Args object and transferring the objects specifically.
00062 //
00063 
00064 namespace {
00065 
00066 
00067 const std::string fwdScalingVecName = "fwdScalingVec";
00068 
00069 
00070 // Assert that the input scaling vectors have been set up correctly
00071 void assertModelVarScalings(
00072   const EpetraExt::ModelEvaluator::InArgs &varScalings
00073   )
00074 {
00075   typedef EpetraExt::ModelEvaluator EME;
00076   TEST_FOR_EXCEPTION(
00077     (varScalings.supports(EME::IN_ARG_x) && varScalings.supports(EME::IN_ARG_x_dot))
00078     && (varScalings.get_x() != varScalings.get_x_dot()),
00079     std::logic_error,
00080     "Error, if scaling for x is given and x_dot is supported, then\n"
00081     "the scaling for x_dot must also be set and must be the same scaling\n"
00082     "as is used for x!"
00083     );
00084 }
00085   
00086 
00087 
00088 // Scale a single vector using a templated policy object to take care of what
00089 // vector gets used.
00090 template<class InArgsVectorGetterSetter>
00091 void scaleModelVar(
00092   InArgsVectorGetterSetter vecGetterSetter, // Templated policy object!
00093   const EpetraExt::ModelEvaluator::InArgs &origVars,
00094   const EpetraExt::ModelEvaluator::InArgs &varScalings,
00095   EpetraExt::ModelEvaluator::InArgs *scaledVars,
00096   Teuchos::FancyOStream *out,
00097   Teuchos::EVerbosityLevel verbLevel
00098   )
00099 {
00100 
00101   using Teuchos::null;
00102   using Teuchos::rcp;
00103   using Teuchos::RCP;
00104   using Teuchos::Ptr;
00105   using Teuchos::rcp_const_cast;
00106   typedef EpetraExt::ModelEvaluator EME;
00107 
00108 
00109 #ifdef TEUCHOS_DEBUG
00110   TEST_FOR_EXCEPT(!scaledVars);
00111 #endif
00112 
00113   RCP<const Epetra_Vector>
00114     orig_vec = vecGetterSetter.getVector(origVars);
00115   if ( !is_null(orig_vec) ) {
00116     RCP<const Epetra_Vector>
00117       inv_s_vec = vecGetterSetter.getVector(varScalings);
00118     if ( !is_null(inv_s_vec) ) {
00119       RCP<Epetra_Vector>
00120         scaled_vec = rcp_const_cast<Epetra_Vector>(
00121           vecGetterSetter.getVector(*scaledVars) );
00122       if ( is_null(scaled_vec) )
00123         scaled_vec = rcp(new Epetra_Vector(orig_vec->Map()));
00124       // See if there is a "hidden" forward scaling vector to use
00125       Ptr<const RCP<const Epetra_Vector> > fwd_s_vec =
00126         Teuchos::getOptionalEmbeddedObj<Epetra_Vector, RCP<const Epetra_Vector> >(
00127           inv_s_vec);
00128 /*
00129         Teuchos::get_optional_extra_data<const RCP<const Epetra_Vector> >(
00130           inv_s_vec, fwdScalingVecName );
00131 */
00132       if ( !is_null(fwd_s_vec) ) {
00133         // Use the "hidden" forward scaling vector and multiply
00134         scaled_vec->Multiply( 1.0, **fwd_s_vec, *orig_vec, 0.0 );
00135       }
00136       else {
00137         // Just use the inverse scaling vector and divide
00138         EpetraExt::scaleModelVarsGivenInverseScaling(
00139           *orig_vec, *inv_s_vec, &*scaled_vec );
00140       }
00141       vecGetterSetter.setVector( scaled_vec, scaledVars );
00142     }
00143     else {
00144       vecGetterSetter.setVector( orig_vec, scaledVars );
00145     }
00146   }
00147   else {
00148     vecGetterSetter.setVector( null, scaledVars );
00149   }
00150 
00151 }
00152   
00153 
00154 // Scale variable bounds for a single vector using a templated policy object
00155 // to take care of what vector gets used.
00156 template<class InArgsVectorGetterSetter>
00157 void scaleModelBound(
00158   InArgsVectorGetterSetter vecGetterSetter, // Templated policy object!
00159   const EpetraExt::ModelEvaluator::InArgs &origLowerBounds,
00160   const EpetraExt::ModelEvaluator::InArgs &origUpperBounds,
00161   const double infBnd,
00162   const EpetraExt::ModelEvaluator::InArgs &varScalings,
00163   EpetraExt::ModelEvaluator::InArgs *scaledLowerBounds,
00164   EpetraExt::ModelEvaluator::InArgs *scaledUpperBounds,
00165   Teuchos::FancyOStream *out,
00166   Teuchos::EVerbosityLevel verbLevel
00167   )
00168 {
00169 
00170   using Teuchos::null;
00171   using Teuchos::rcp;
00172   using Teuchos::RCP;
00173   using Teuchos::rcp_const_cast;
00174   typedef EpetraExt::ModelEvaluator EME;
00175 
00176 #ifdef TEUCHOS_DEBUG
00177   TEST_FOR_EXCEPT(!scaledLowerBounds);
00178   TEST_FOR_EXCEPT(!scaledUpperBounds);
00179 #endif
00180 
00181   RCP<const Epetra_Vector>
00182     orig_lower_vec = vecGetterSetter.getVector(origLowerBounds);
00183   if ( !is_null(orig_lower_vec) ) {
00184     RCP<const Epetra_Vector>
00185       inv_s_vec = vecGetterSetter.getVector(varScalings);
00186     if ( !is_null(inv_s_vec) ) {
00187       TEST_FOR_EXCEPT("Can't handle scaling bounds yet!");
00188     }
00189     else {
00190       vecGetterSetter.setVector( orig_lower_vec, scaledLowerBounds );
00191     }
00192   }
00193   else {
00194     vecGetterSetter.setVector( null, scaledLowerBounds );
00195   }
00196 
00197   RCP<const Epetra_Vector>
00198     orig_upper_vec = vecGetterSetter.getVector(origUpperBounds);
00199   if ( !is_null(orig_upper_vec) ) {
00200     RCP<const Epetra_Vector>
00201       inv_s_vec = vecGetterSetter.getVector(varScalings);
00202     if ( !is_null(inv_s_vec) ) {
00203       TEST_FOR_EXCEPT("Can't handle scaling bounds yet!");
00204     }
00205     else {
00206       vecGetterSetter.setVector( orig_upper_vec, scaledUpperBounds );
00207     }
00208   }
00209   else {
00210     vecGetterSetter.setVector( null, scaledUpperBounds );
00211   }
00212 
00213 }
00214 
00215 
00216 // Unscale a single vector using a templated policy object to take care of
00217 // what vector gets used.
00218 template<class InArgsVectorGetterSetter>
00219 void unscaleModelVar(
00220   InArgsVectorGetterSetter vecGetterSetter, // Templated policy object!
00221   const EpetraExt::ModelEvaluator::InArgs &scaledVars,
00222   const EpetraExt::ModelEvaluator::InArgs &varScalings,
00223   EpetraExt::ModelEvaluator::InArgs *origVars,
00224   Teuchos::FancyOStream *out,
00225   Teuchos::EVerbosityLevel verbLevel
00226   )
00227 {
00228 
00229   using Teuchos::null;
00230   using Teuchos::rcp;
00231   using Teuchos::RCP;
00232   using Teuchos::rcp_const_cast;
00233   using Teuchos::includesVerbLevel;
00234   typedef EpetraExt::ModelEvaluator EME;
00235 
00236 
00237 #ifdef TEUCHOS_DEBUG
00238   TEST_FOR_EXCEPT(!origVars);
00239 #endif
00240 
00241   RCP<const Epetra_Vector>
00242     scaled_vec = vecGetterSetter.getVector(scaledVars);
00243   if ( !is_null(scaled_vec) ) {
00244     RCP<const Epetra_Vector>
00245       inv_s_vec = vecGetterSetter.getVector(varScalings);
00246     if ( !is_null(inv_s_vec) ) {
00247       RCP<Epetra_Vector>
00248         orig_vec = rcp_const_cast<Epetra_Vector>(
00249           vecGetterSetter.getVector(*origVars) );
00250       if ( is_null(orig_vec) )
00251         orig_vec = rcp(new Epetra_Vector(scaled_vec->Map()));
00252       EpetraExt::unscaleModelVarsGivenInverseScaling(
00253         *scaled_vec, *inv_s_vec,  &*orig_vec );
00254       if (out && includesVerbLevel(verbLevel,Teuchos::VERB_HIGH)) {
00255         *out << "\nUnscaled vector "<<vecGetterSetter.getName()<<":\n";
00256         Teuchos::OSTab tab(*out);
00257         orig_vec->Print(*out);
00258       }
00259       vecGetterSetter.setVector( orig_vec, origVars );
00260     }
00261     else {
00262       vecGetterSetter.setVector( scaled_vec, origVars );
00263     }
00264   }
00265   else {
00266     vecGetterSetter.setVector( null, origVars );
00267   }
00268 
00269 }
00270 
00271 
00272 // Scale a single vector using a templated policy object to take care of what
00273 // vector gets used.
00274 template<class OutArgsVectorGetterSetter>
00275 void scaleModelFunc(
00276   OutArgsVectorGetterSetter vecGetterSetter, // Templated policy object!
00277   const EpetraExt::ModelEvaluator::OutArgs &origFuncs,
00278   const EpetraExt::ModelEvaluator::OutArgs &funcScalings,
00279   EpetraExt::ModelEvaluator::OutArgs *scaledFuncs,
00280   Teuchos::FancyOStream *out,
00281   Teuchos::EVerbosityLevel verbLevel
00282   )
00283 {
00284   TEST_FOR_EXCEPT(0==scaledFuncs);
00285   Teuchos::RCP<Epetra_Vector>
00286     func = vecGetterSetter.getVector(origFuncs);
00287   if (!is_null(func) ) {
00288     Teuchos::RCP<const Epetra_Vector>
00289       funcScaling = vecGetterSetter.getVector(funcScalings);
00290     if (!is_null(funcScaling) ) {
00291       EpetraExt::scaleModelFuncGivenForwardScaling( *funcScaling, &*func );
00292     }
00293   }
00294   vecGetterSetter.setVector( func, scaledFuncs );
00295 }
00296 
00297 
00298 } // namespace
00299 
00300 
00301 void EpetraExt::gatherModelNominalValues(
00302   const ModelEvaluator &model,
00303   ModelEvaluator::InArgs *nominalValues
00304   )
00305 {
00306 
00307   using Teuchos::implicit_cast;
00308   typedef ModelEvaluator EME;
00309 
00310 #ifdef TEUCHOS_DEBUG
00311   TEST_FOR_EXCEPT(!nominalValues);
00312 #endif
00313 
00314   *nominalValues = model.createInArgs();
00315 
00316   if(nominalValues->supports(EME::IN_ARG_x)) {
00317     nominalValues->set_x(model.get_x_init());
00318   }
00319 
00320   if(nominalValues->supports(EME::IN_ARG_x_dot)) {
00321     nominalValues->set_x_dot(model.get_x_dot_init());
00322   }
00323   
00324   for( int l = 0; l < nominalValues->Np(); ++l ) {
00325     nominalValues->set_p( l, model.get_p_init(l) );
00326   }
00327   
00328   if(nominalValues->supports(EME::IN_ARG_t)) {
00329     nominalValues->set_t(model.get_t_init());
00330   }
00331 
00332 }
00333 
00334 
00335 void EpetraExt::gatherModelBounds(
00336   const ModelEvaluator &model,
00337   ModelEvaluator::InArgs *lowerBounds,
00338   ModelEvaluator::InArgs *upperBounds
00339   )
00340 {
00341 
00342   using Teuchos::implicit_cast;
00343   typedef ModelEvaluator EME;
00344 
00345 #ifdef TEUCHOS_DEBUG
00346   TEST_FOR_EXCEPT(!lowerBounds);
00347   TEST_FOR_EXCEPT(!upperBounds);
00348 #endif
00349 
00350   *lowerBounds = model.createInArgs();
00351   *upperBounds = model.createInArgs();
00352 
00353   if(lowerBounds->supports(EME::IN_ARG_x)) {
00354     lowerBounds->set_x(model.get_x_lower_bounds());
00355     upperBounds->set_x(model.get_x_upper_bounds());
00356   }
00357   
00358   for( int l = 0; l < lowerBounds->Np(); ++l ) {
00359     lowerBounds->set_p( l, model.get_p_lower_bounds(l) );
00360     upperBounds->set_p( l, model.get_p_upper_bounds(l) );
00361   }
00362   
00363   if(lowerBounds->supports(EME::IN_ARG_t)) {
00364     lowerBounds->set_t(model.get_t_lower_bound());
00365     upperBounds->set_t(model.get_t_upper_bound());
00366   }
00367 
00368 }
00369 
00370 
00371 void EpetraExt::scaleModelVars(
00372   const ModelEvaluator::InArgs &origVars,
00373   const ModelEvaluator::InArgs &varScalings,
00374   ModelEvaluator::InArgs *scaledVars,
00375   Teuchos::FancyOStream *out,
00376   Teuchos::EVerbosityLevel verbLevel
00377   )
00378 {
00379   typedef ModelEvaluator EME;
00380 
00381 #ifdef TEUCHOS_DEBUG
00382   TEST_FOR_EXCEPT(!scaledVars);
00383   assertModelVarScalings(varScalings);
00384 #endif
00385 
00386   if (origVars.supports(EME::IN_ARG_x)) {
00387     scaleModelVar( InArgsGetterSetter_x(), origVars, varScalings, scaledVars,
00388       out, verbLevel );
00389   }
00390 
00391   if (origVars.supports(EME::IN_ARG_x_dot)) {
00392     scaleModelVar( InArgsGetterSetter_x_dot(), origVars, varScalings, scaledVars,
00393       out, verbLevel );
00394   }
00395 
00396   const int Np = origVars.Np();
00397   for ( int l = 0; l < Np; ++l ) {
00398     scaleModelVar( InArgsGetterSetter_p(l), origVars, varScalings, scaledVars,
00399       out, verbLevel );
00400   }
00401 
00402   if (origVars.supports(EME::IN_ARG_x_poly)) {
00403     TEST_FOR_EXCEPTION(
00404       !is_null(varScalings.get_x()), std::logic_error,
00405       "Error, can't hanlde scaling of x_poly yet!"
00406       );
00407     scaledVars->set_x_poly(origVars.get_x_poly());
00408   }
00409 
00410   if (origVars.supports(EME::IN_ARG_x_dot_poly)) {
00411     TEST_FOR_EXCEPTION(
00412       !is_null(varScalings.get_x()), std::logic_error,
00413       "Error, can't hanlde scaling of x_dot_poly yet!"
00414       );
00415     scaledVars->set_x_dot_poly(origVars.get_x_dot_poly());
00416   }
00417 
00418   if (origVars.supports(EME::IN_ARG_t)) {
00419     TEST_FOR_EXCEPTION(
00420       varScalings.get_t() > 0.0, std::logic_error,
00421       "Error, can't hanlde scaling of t yet!"
00422       );
00423     scaledVars->set_t(origVars.get_t());
00424   }
00425 
00426   if (origVars.supports(EME::IN_ARG_alpha)) {
00427     TEST_FOR_EXCEPTION(
00428       varScalings.get_alpha() > 0.0, std::logic_error,
00429       "Error, can't hanlde scaling of alpha yet!"
00430       );
00431     scaledVars->set_alpha(origVars.get_alpha());
00432   }
00433 
00434   if (origVars.supports(EME::IN_ARG_beta)) {
00435     TEST_FOR_EXCEPTION(
00436       varScalings.get_beta() > 0.0, std::logic_error,
00437       "Error, can't hanlde scaling of beta yet!"
00438       );
00439     scaledVars->set_beta(origVars.get_beta());
00440   }
00441 
00442   // ToDo: Support other input arguments?
00443 
00444 }
00445 
00446 
00447 void EpetraExt::scaleModelBounds(
00448   const ModelEvaluator::InArgs &origLowerBounds,
00449   const ModelEvaluator::InArgs &origUpperBounds,
00450   const double infBnd,
00451   const ModelEvaluator::InArgs &varScalings,
00452   ModelEvaluator::InArgs *scaledLowerBounds,
00453   ModelEvaluator::InArgs *scaledUpperBounds,
00454   Teuchos::FancyOStream *out,
00455   Teuchos::EVerbosityLevel verbLevel
00456   )
00457 {
00458 
00459   typedef ModelEvaluator EME;
00460 
00461 #ifdef TEUCHOS_DEBUG
00462   TEST_FOR_EXCEPT(!scaledLowerBounds);
00463   TEST_FOR_EXCEPT(!scaledUpperBounds);
00464   assertModelVarScalings(varScalings);
00465 #endif
00466 
00467   if (origLowerBounds.supports(EME::IN_ARG_x)) {
00468     scaleModelBound(
00469       InArgsGetterSetter_x(), origLowerBounds, origUpperBounds, infBnd,
00470       varScalings, scaledLowerBounds, scaledUpperBounds,
00471       out, verbLevel );
00472   }
00473 
00474   if (origLowerBounds.supports(EME::IN_ARG_x_dot)) {
00475     scaleModelBound(
00476       InArgsGetterSetter_x_dot(), origLowerBounds, origUpperBounds, infBnd,
00477       varScalings, scaledLowerBounds, scaledUpperBounds,
00478       out, verbLevel );
00479   }
00480 
00481   const int Np = origLowerBounds.Np();
00482   for ( int l = 0; l < Np; ++l ) {
00483     scaleModelBound(
00484       InArgsGetterSetter_p(l), origLowerBounds, origUpperBounds, infBnd,
00485       varScalings, scaledLowerBounds, scaledUpperBounds,
00486       out, verbLevel );
00487   }
00488 
00489   // ToDo: Support other input arguments?
00490 
00491 }
00492 
00493 
00494 void EpetraExt::unscaleModelVars(
00495   const ModelEvaluator::InArgs &scaledVars,
00496   const ModelEvaluator::InArgs &varScalings,
00497   ModelEvaluator::InArgs *origVars,
00498   Teuchos::FancyOStream *out,
00499   Teuchos::EVerbosityLevel verbLevel
00500   )
00501 {
00502 
00503   using Teuchos::RCP;
00504   using Teuchos::includesVerbLevel;
00505   typedef ModelEvaluator EME;
00506 
00507 #ifdef TEUCHOS_DEBUG
00508   TEST_FOR_EXCEPT(!origVars);
00509   assertModelVarScalings(varScalings);
00510 #endif
00511 
00512   // Print scaling vectors
00513 
00514   if (out && includesVerbLevel(verbLevel,Teuchos::VERB_HIGH)) {
00515     RCP<const Epetra_Vector> inv_s_x;
00516     if ( scaledVars.supports(EME::IN_ARG_x) &&
00517       !is_null(inv_s_x=varScalings.get_x()) )
00518     {
00519       *out << "\nState inverse scaling vector inv_s_x:\n";
00520       Teuchos::OSTab tab(*out);
00521       inv_s_x->Print(*out);
00522     }
00523   }
00524 
00525   // Scal the input varaibles
00526   
00527   if (scaledVars.supports(EME::IN_ARG_x_dot)) {
00528     unscaleModelVar( InArgsGetterSetter_x_dot(), scaledVars, varScalings, origVars,
00529       out, verbLevel );
00530   }
00531 
00532   if (scaledVars.supports(EME::IN_ARG_x)) {
00533     unscaleModelVar( InArgsGetterSetter_x(), scaledVars, varScalings, origVars,
00534       out, verbLevel );
00535   }
00536 
00537   const int Np = scaledVars.Np();
00538   for ( int l = 0; l < Np; ++l ) {
00539     unscaleModelVar( InArgsGetterSetter_p(l), scaledVars, varScalings, origVars,
00540       out, verbLevel );
00541   }
00542 
00543   if (scaledVars.supports(EME::IN_ARG_x_dot_poly)) {
00544     TEST_FOR_EXCEPTION(
00545       !is_null(varScalings.get_x()), std::logic_error,
00546       "Error, can't hanlde unscaling of x_dot_poly yet!"
00547       );
00548     origVars->set_x_dot_poly(scaledVars.get_x_dot_poly());
00549   }
00550 
00551   if (scaledVars.supports(EME::IN_ARG_x_poly)) {
00552     TEST_FOR_EXCEPTION(
00553       !is_null(varScalings.get_x()), std::logic_error,
00554       "Error, can't hanlde unscaling of x_poly yet!"
00555       );
00556     origVars->set_x_poly(scaledVars.get_x_poly());
00557   }
00558 
00559   if (scaledVars.supports(EME::IN_ARG_t)) {
00560     TEST_FOR_EXCEPTION(
00561       varScalings.get_t() > 0.0, std::logic_error,
00562       "Error, can't hanlde unscaling of t yet!"
00563       );
00564     origVars->set_t(scaledVars.get_t());
00565   }
00566 
00567   if (scaledVars.supports(EME::IN_ARG_alpha)) {
00568     TEST_FOR_EXCEPTION(
00569       varScalings.get_alpha() > 0.0, std::logic_error,
00570       "Error, can't hanlde unscaling of alpha yet!"
00571       );
00572     origVars->set_alpha(scaledVars.get_alpha());
00573   }
00574 
00575   if (scaledVars.supports(EME::IN_ARG_beta)) {
00576     TEST_FOR_EXCEPTION(
00577       varScalings.get_beta() > 0.0, std::logic_error,
00578       "Error, can't hanlde unscaling of beta yet!"
00579       );
00580     origVars->set_beta(scaledVars.get_beta());
00581   }
00582 
00583 }
00584 
00585 
00586 void EpetraExt::scaleModelFuncs(
00587   const ModelEvaluator::OutArgs &origFuncs,
00588   const ModelEvaluator::InArgs &varScalings,
00589   const ModelEvaluator::OutArgs &funcScalings,
00590   ModelEvaluator::OutArgs *scaledFuncs,
00591   bool *allFuncsWhereScaled,
00592   Teuchos::FancyOStream *out,
00593   Teuchos::EVerbosityLevel verbLevel
00594   )
00595 {
00596 
00597   using Teuchos::RCP;
00598   typedef ModelEvaluator EME;
00599 
00600   TEST_FOR_EXCEPT(0==allFuncsWhereScaled);
00601 
00602   *allFuncsWhereScaled = true;
00603 
00604   const int Np = origFuncs.Np();
00605   const int Ng = origFuncs.Ng();
00606 
00607   // f
00608   if ( origFuncs.supports(EME::OUT_ARG_f) && !is_null(origFuncs.get_f()) ) {
00609     scaleModelFunc( OutArgsGetterSetter_f(), origFuncs, funcScalings, scaledFuncs,
00610       out, verbLevel );
00611   }
00612 
00613   // f_poly
00614   if (
00615     origFuncs.supports(EME::OUT_ARG_f_poly)
00616     && !is_null(origFuncs.get_f_poly())
00617     )
00618   {
00619     TEST_FOR_EXCEPTION(
00620       !is_null(funcScalings.get_f()), std::logic_error,
00621       "Error, we can't handle scaling of f_poly yet!"
00622       );
00623     scaledFuncs->set_f_poly(origFuncs.get_f_poly());
00624   }
00625 
00626   // g(j)
00627   for ( int j = 0; j < Ng; ++j ) {
00628     scaleModelFunc( OutArgsGetterSetter_g(j), origFuncs, funcScalings, scaledFuncs,
00629       out, verbLevel );
00630   }
00631 
00632   // W
00633   RCP<Epetra_Operator> W;
00634   if ( origFuncs.supports(EME::OUT_ARG_W) && !is_null(W=origFuncs.get_W()) ) {
00635     bool didScaling = false;
00636     scaleModelFuncFirstDerivOp(
00637       varScalings.get_x().get(), funcScalings.get_f().get(),
00638       &*W, &didScaling
00639       );
00640     if(didScaling)
00641       scaledFuncs->set_W(W);
00642     else
00643       *allFuncsWhereScaled = false;
00644   }
00645 
00646   // DfDp(l)
00647   for ( int l = 0; l < Np; ++l ) {
00648     EME::Derivative orig_DfDp_l;
00649     if (
00650       !origFuncs.supports(EME::OUT_ARG_DfDp,l).none()
00651       && !(orig_DfDp_l=origFuncs.get_DfDp(l)).isEmpty()
00652       )
00653     {
00654       EME::Derivative scaled_DfDp_l;
00655       bool didScaling = false;
00656       scaleModelFuncFirstDeriv(
00657         orig_DfDp_l, varScalings.get_p(l).get(), funcScalings.get_f().get(),
00658         &scaled_DfDp_l, &didScaling
00659         );
00660       if(didScaling)
00661         scaledFuncs->set_DfDp(l,scaled_DfDp_l);
00662       else
00663         *allFuncsWhereScaled = false;
00664     }
00665 
00666   }
00667 
00668   // DgDx_dot(j), DgDx(j), and DgDp(j,l)
00669   for ( int j = 0; j < Ng; ++j ) {
00670 
00671     EME::Derivative orig_DgDx_dot_j;
00672     if (
00673       !origFuncs.supports(EME::OUT_ARG_DgDx_dot,j).none()
00674       && !(orig_DgDx_dot_j=origFuncs.get_DgDx_dot(j)).isEmpty()
00675       )
00676     {
00677       EME::Derivative scaled_DgDx_dot_j;
00678       bool didScaling = false;
00679       scaleModelFuncFirstDeriv(
00680         orig_DgDx_dot_j, varScalings.get_x().get(), funcScalings.get_g(j).get(),
00681         &scaled_DgDx_dot_j, &didScaling
00682         );
00683       if(didScaling)
00684         scaledFuncs->set_DgDx_dot(j,scaled_DgDx_dot_j);
00685       else
00686         *allFuncsWhereScaled = false;
00687     }
00688 
00689     EME::Derivative orig_DgDx_j;
00690     if (
00691       !origFuncs.supports(EME::OUT_ARG_DgDx,j).none()
00692       && !(orig_DgDx_j=origFuncs.get_DgDx(j)).isEmpty()
00693       )
00694     {
00695       EME::Derivative scaled_DgDx_j;
00696       bool didScaling = false;
00697       scaleModelFuncFirstDeriv(
00698         orig_DgDx_j, varScalings.get_x().get(), funcScalings.get_g(j).get(),
00699         &scaled_DgDx_j, &didScaling
00700         );
00701       if(didScaling)
00702         scaledFuncs->set_DgDx(j,scaled_DgDx_j);
00703       else
00704         *allFuncsWhereScaled = false;
00705     }
00706 
00707     for ( int l = 0; l < Np; ++l ) {
00708       EME::Derivative orig_DgDp_j_l;
00709       if (
00710         !origFuncs.supports(EME::OUT_ARG_DgDp,j,l).none()
00711         && !(orig_DgDp_j_l=origFuncs.get_DgDp(j,l)).isEmpty()
00712         )
00713       {
00714         EME::Derivative scaled_DgDp_j_l;
00715         bool didScaling = false;
00716         scaleModelFuncFirstDeriv(
00717           orig_DgDp_j_l, varScalings.get_p(l).get(), funcScalings.get_g(j).get(),
00718           &scaled_DgDp_j_l, &didScaling
00719           );
00720         if(didScaling)
00721           scaledFuncs->set_DgDp(j,l,scaled_DgDp_j_l);
00722         else
00723           *allFuncsWhereScaled = false;
00724       }
00725     }
00726   }
00727 
00728 }
00729 
00730 
00731 Teuchos::RCP<const Epetra_Vector>
00732 EpetraExt::createInverseModelScalingVector(
00733   Teuchos::RCP<const Epetra_Vector> const& scalingVector
00734   )
00735 {
00736   Teuchos::RCP<Epetra_Vector> invScalingVector =
00737     Teuchos::rcpWithEmbeddedObj(
00738       new Epetra_Vector(scalingVector->Map()),
00739       scalingVector
00740       );
00741   invScalingVector->Reciprocal(*scalingVector);
00742   return invScalingVector;
00743   // Above, we embedd the forward scaling vector.  This is done in order to
00744   // achieve the exact same numerics as before this refactoring and to improve
00745   // runtime speed and accruacy.
00746 }
00747 
00748 
00749 void EpetraExt::scaleModelVarsGivenInverseScaling(
00750   const Epetra_Vector &origVars,
00751   const Epetra_Vector &invVarScaling,
00752   Epetra_Vector *scaledVars
00753   )
00754 {
00755 #ifdef TEUCHOS_DEBUG
00756   TEST_FOR_EXCEPT(!scaledVars);
00757   TEST_FOR_EXCEPT(!origVars.Map().SameAs(invVarScaling.Map()));
00758   TEST_FOR_EXCEPT(!origVars.Map().SameAs(scaledVars->Map()));
00759 #endif
00760   const int localDim = origVars.Map().NumMyElements();
00761   for ( int i = 0; i < localDim; ++i )
00762     (*scaledVars)[i] = origVars[i] / invVarScaling[i];
00763 }
00764 
00765 
00766 void EpetraExt::scaleModelVarBoundsGivenInverseScaling(
00767   const Epetra_Vector &origLowerBounds,
00768   const Epetra_Vector &origUpperBounds,
00769   const double infBnd,
00770   const Epetra_Vector &invVarScaling,
00771   Epetra_Vector *scaledLowerBounds,
00772   Epetra_Vector *scaledUpperBounds
00773   )
00774 {
00775   TEST_FOR_EXCEPT("ToDo: Implement!");
00776 }
00777 
00778 
00779 void EpetraExt::unscaleModelVarsGivenInverseScaling(
00780   const Epetra_Vector &origVars,
00781   const Epetra_Vector &invVarScaling,
00782   Epetra_Vector *scaledVars
00783   )
00784 {
00785   TEST_FOR_EXCEPT(0==scaledVars);
00786   scaledVars->Multiply( 1.0, invVarScaling, origVars, 0.0 );
00787 }
00788 
00789 
00790 void EpetraExt::scaleModelFuncGivenForwardScaling(
00791   const Epetra_Vector &fwdFuncScaling,
00792   Epetra_Vector *funcs
00793   )
00794 {
00795   TEST_FOR_EXCEPT(0==funcs);
00796   funcs->Multiply( 1.0,  fwdFuncScaling, *funcs, 0.0 );
00797   // Note: Above is what Epetra_LinearProblem does to scale the RHS and LHS
00798   // vectors so this type of argument aliasing must be okay in Epetra!
00799 }
00800 
00801 
00802 void EpetraExt::scaleModelFuncFirstDerivOp(
00803   const Epetra_Vector *invVarScaling,
00804   const Epetra_Vector *fwdFuncScaling,
00805   Epetra_Operator *funcDerivOp,
00806   bool *didScaling
00807   )
00808 {
00809   TEST_FOR_EXCEPT(0==funcDerivOp);
00810   TEST_FOR_EXCEPT(0==didScaling);
00811   *didScaling = false; // Assume not scaled to start
00812   Epetra_RowMatrix *funcDerivRowMatrix
00813     = dynamic_cast<Epetra_RowMatrix*>(funcDerivOp);
00814   if (funcDerivRowMatrix) {
00815     if (fwdFuncScaling)
00816       funcDerivRowMatrix->LeftScale(*fwdFuncScaling);
00817     if (invVarScaling)
00818       funcDerivRowMatrix->RightScale(*invVarScaling);
00819     *didScaling = true;
00820     // Note that above I do the func scaling before the var scaling since that
00821     // is the same order it was done for W in Thyra::EpetraModelEvaluator
00822   }
00823 }
00824 
00825 
00826 void EpetraExt::scaleModelFuncFirstDeriv(
00827   const ModelEvaluator::Derivative &origFuncDeriv,
00828   const Epetra_Vector *invVarScaling,
00829   const Epetra_Vector *fwdFuncScaling,
00830   ModelEvaluator::Derivative *scaledFuncDeriv,
00831   bool *didScaling
00832   )
00833 {
00834   using Teuchos::RCP;
00835   typedef ModelEvaluator EME;
00836   TEST_FOR_EXCEPT(0==scaledFuncDeriv);
00837   TEST_FOR_EXCEPT(0==didScaling);
00838   *didScaling = false;
00839   const RCP<Epetra_MultiVector>
00840     funcDerivMv = origFuncDeriv.getMultiVector();
00841   const EME::EDerivativeMultiVectorOrientation
00842     funcDerivMv_orientation = origFuncDeriv.getMultiVectorOrientation();
00843   if(!is_null(funcDerivMv)) {
00844     if ( funcDerivMv_orientation == EME::DERIV_MV_BY_COL )
00845     {
00846       if (fwdFuncScaling) {
00847         funcDerivMv->Multiply(1.0, *fwdFuncScaling, *funcDerivMv, 0.0);
00848       }
00849       if (invVarScaling) {
00850         TEST_FOR_EXCEPT("ToDo: Scale rows!");
00851         //funcDerivMv->Multiply(1.0, *funcDerivMv, *invVarScaling, 0.0);
00852       }
00853     }
00854     else if ( funcDerivMv_orientation == EME::DERIV_TRANS_MV_BY_ROW )
00855     {
00856       if (invVarScaling) {
00857         funcDerivMv->Multiply(1.0, *invVarScaling, *funcDerivMv, 0.0);
00858       }
00859       if (fwdFuncScaling) {
00860         TEST_FOR_EXCEPT("ToDo: Scale rows!");
00861         //funcDerivMv->Multiply(1.0, *funcDerivMv, *fwdFuncScaling, 0.0);
00862       }
00863     }
00864     else {
00865       TEST_FOR_EXCEPT("Should not get here!");
00866     }
00867     *scaledFuncDeriv = EME::Derivative(funcDerivMv,funcDerivMv_orientation);
00868     *didScaling = true;
00869   }
00870   else {
00871     RCP<Epetra_Operator>
00872       funcDerivOp = origFuncDeriv.getLinearOp();
00873     TEST_FOR_EXCEPT(is_null(funcDerivOp));
00874     scaleModelFuncFirstDerivOp( invVarScaling, fwdFuncScaling,
00875       &*funcDerivOp, didScaling );
00876     if (didScaling)
00877       *scaledFuncDeriv = EME::Derivative(funcDerivOp);
00878   }
00879 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines