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

Generated on Wed May 12 21:40:37 2010 for EpetraExt by  doxygen 1.4.7