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