EpetraExt_ModelEvaluatorScalingTools.cpp

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

Generated on Tue Oct 20 12:45:30 2009 for EpetraExt by doxygen 1.4.7