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