EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
EpetraExt_ModelEvaluatorScalingTools.cpp
Go to the documentation of this file.
1//@HEADER
2// ***********************************************************************
3//
4// EpetraExt: Epetra Extended - Linear Algebra Services Package
5// Copyright (2011) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41
42
44#include "Teuchos_implicit_cast.hpp"
45#include "Epetra_RowMatrix.h"
46
47//
48// Here in the implementation of scaling we write all scaling routines to
49// specifically and individually handle every input and output object and we
50// transfer all objects from one [In,Out]Arg container to another to make sure
51// that that object has been specifically addressed. We could just use the
52// [In,Out]Args::setArgs(...) function to copy all arguments be default but
53// then that would setup subtle errors where a new quality could be silently
54// ignored and not be scaled correctly. Therefore, we feel it is better to
55// have to deal with *every* input and output object specifically with respect
56// to scaling in order to avoid overlooking scaling. In this way, if a new
57// input or output object is added to [In,Out]Args but the code in this file
58// is not updated, then that object will not be passed through and some type
59// of error will be generated right away. We feel this is the best behavior
60// and it justifies having every scaling-related function take both an input
61// and an output [In,Out]Args object and transferring the objects specifically.
62//
63
64namespace {
65
66
67const std::string fwdScalingVecName = "fwdScalingVec";
68
69#ifdef TEUCHOS_DEBUG
70
71// Assert that the input scaling vectors have been set up correctly.
72//
73// This function ONLY exists if TEUCHOS_DEBUG is defined.
74// This avoids "unused function" warnings with Clang 3.5.
75void assertModelVarScalings(
76 const EpetraExt::ModelEvaluator::InArgs &varScalings
77 )
78{
79 typedef EpetraExt::ModelEvaluator EME;
80 TEUCHOS_TEST_FOR_EXCEPTION(
81 (varScalings.supports(EME::IN_ARG_x) && varScalings.supports(EME::IN_ARG_x_dot))
82 && (varScalings.get_x() != varScalings.get_x_dot()),
83 std::logic_error,
84 "Error, if scaling for x is given and x_dot is supported, then\n"
85 "the scaling for x_dot must also be set and must be the same scaling\n"
86 "as is used for x!"
87 );
88 TEUCHOS_TEST_FOR_EXCEPTION(
89 (varScalings.supports(EME::IN_ARG_x) && varScalings.supports(EME::IN_ARG_x_dotdot))
90 && (varScalings.get_x() != varScalings.get_x_dotdot()),
91 std::logic_error,
92 "Error, if scaling for x is given and x_dotdot is supported, then\n"
93 "the scaling for x_dotdot must also be set and must be the same scaling\n"
94 "as is used for x!"
95 );
96}
97
98#endif // TEUCHOS_DEBUG
99
100
101// Scale a single vector using a templated policy object to take care of what
102// vector gets used.
103template<class InArgsVectorGetterSetter>
104void scaleModelVar(
105 InArgsVectorGetterSetter vecGetterSetter, // Templated policy object!
106 const EpetraExt::ModelEvaluator::InArgs &origVars,
107 const EpetraExt::ModelEvaluator::InArgs &varScalings,
109 Teuchos::FancyOStream * /* out */,
110 Teuchos::EVerbosityLevel /* verbLevel */
111 )
112{
113
114 using Teuchos::null;
115 using Teuchos::rcp;
116 using Teuchos::RCP;
117 using Teuchos::Ptr;
118 using Teuchos::rcp_const_cast;
119
120
121#ifdef TEUCHOS_DEBUG
122 TEUCHOS_TEST_FOR_EXCEPT(!scaledVars);
123#endif
124
125 RCP<const Epetra_Vector>
126 orig_vec = vecGetterSetter.getVector(origVars);
127 if ( !is_null(orig_vec) ) {
128 RCP<const Epetra_Vector>
129 inv_s_vec = vecGetterSetter.getVector(varScalings);
130 if ( !is_null(inv_s_vec) ) {
131 RCP<Epetra_Vector>
132 scaled_vec = rcp_const_cast<Epetra_Vector>(
133 vecGetterSetter.getVector(*scaledVars) );
134 if ( is_null(scaled_vec) )
135 scaled_vec = rcp(new Epetra_Vector(orig_vec->Map()));
136 // See if there is a "hidden" forward scaling vector to use
137 Ptr<const RCP<const Epetra_Vector> > fwd_s_vec =
138 Teuchos::getOptionalEmbeddedObj<Epetra_Vector, RCP<const Epetra_Vector> >(
139 inv_s_vec);
140/*
141 Teuchos::get_optional_extra_data<const RCP<const Epetra_Vector> >(
142 inv_s_vec, fwdScalingVecName );
143*/
144 if ( !is_null(fwd_s_vec) ) {
145 // Use the "hidden" forward scaling vector and multiply
146 scaled_vec->Multiply( 1.0, **fwd_s_vec, *orig_vec, 0.0 );
147 }
148 else {
149 // Just use the inverse scaling vector and divide
151 *orig_vec, *inv_s_vec, &*scaled_vec );
152 }
153 vecGetterSetter.setVector( scaled_vec, scaledVars );
154 }
155 else {
156 vecGetterSetter.setVector( orig_vec, scaledVars );
157 }
158 }
159 else {
160 vecGetterSetter.setVector( null, scaledVars );
161 }
162
163}
164
165
166// Scale variable bounds for a single vector using a templated policy object
167// to take care of what vector gets used.
168template<class InArgsVectorGetterSetter>
169void scaleModelBound(
170 InArgsVectorGetterSetter vecGetterSetter, // Templated policy object!
171 const EpetraExt::ModelEvaluator::InArgs &origLowerBounds,
172 const EpetraExt::ModelEvaluator::InArgs &origUpperBounds,
173 const double /* infBnd */,
174 const EpetraExt::ModelEvaluator::InArgs &varScalings,
175 EpetraExt::ModelEvaluator::InArgs *scaledLowerBounds,
176 EpetraExt::ModelEvaluator::InArgs *scaledUpperBounds,
177 Teuchos::FancyOStream * /* out */,
178 Teuchos::EVerbosityLevel /* verbLevel */
179 )
180{
181
182 using Teuchos::null;
183 using Teuchos::rcp;
184 using Teuchos::RCP;
185 using Teuchos::rcp_const_cast;
186
187#ifdef TEUCHOS_DEBUG
188 TEUCHOS_TEST_FOR_EXCEPT(!scaledLowerBounds);
189 TEUCHOS_TEST_FOR_EXCEPT(!scaledUpperBounds);
190#endif
191
192 RCP<const Epetra_Vector>
193 orig_lower_vec = vecGetterSetter.getVector(origLowerBounds);
194 if ( !is_null(orig_lower_vec) ) {
195 RCP<const Epetra_Vector>
196 inv_s_vec = vecGetterSetter.getVector(varScalings);
197 if ( !is_null(inv_s_vec) ) {
198 TEUCHOS_TEST_FOR_EXCEPT("Can't handle scaling bounds yet!");
199 }
200 else {
201 vecGetterSetter.setVector( orig_lower_vec, scaledLowerBounds );
202 }
203 }
204 else {
205 vecGetterSetter.setVector( null, scaledLowerBounds );
206 }
207
208 RCP<const Epetra_Vector>
209 orig_upper_vec = vecGetterSetter.getVector(origUpperBounds);
210 if ( !is_null(orig_upper_vec) ) {
211 RCP<const Epetra_Vector>
212 inv_s_vec = vecGetterSetter.getVector(varScalings);
213 if ( !is_null(inv_s_vec) ) {
214 TEUCHOS_TEST_FOR_EXCEPT("Can't handle scaling bounds yet!");
215 }
216 else {
217 vecGetterSetter.setVector( orig_upper_vec, scaledUpperBounds );
218 }
219 }
220 else {
221 vecGetterSetter.setVector( null, scaledUpperBounds );
222 }
223
224}
225
226
227// Unscale a single vector using a templated policy object to take care of
228// what vector gets used.
229template<class InArgsVectorGetterSetter>
230void unscaleModelVar(
231 InArgsVectorGetterSetter vecGetterSetter, // Templated policy object!
232 const EpetraExt::ModelEvaluator::InArgs &scaledVars,
233 const EpetraExt::ModelEvaluator::InArgs &varScalings,
235 Teuchos::FancyOStream *out,
236 Teuchos::EVerbosityLevel verbLevel
237 )
238{
239
240 using Teuchos::null;
241 using Teuchos::rcp;
242 using Teuchos::RCP;
243 using Teuchos::rcp_const_cast;
244 using Teuchos::includesVerbLevel;
245
246
247#ifdef TEUCHOS_DEBUG
248 TEUCHOS_TEST_FOR_EXCEPT(!origVars);
249#endif
250
251 RCP<const Epetra_Vector>
252 scaled_vec = vecGetterSetter.getVector(scaledVars);
253 if ( !is_null(scaled_vec) ) {
254 RCP<const Epetra_Vector>
255 inv_s_vec = vecGetterSetter.getVector(varScalings);
256 if ( !is_null(inv_s_vec) ) {
257 RCP<Epetra_Vector>
258 orig_vec = rcp_const_cast<Epetra_Vector>(
259 vecGetterSetter.getVector(*origVars) );
260 if ( is_null(orig_vec) )
261 orig_vec = rcp(new Epetra_Vector(scaled_vec->Map()));
263 *scaled_vec, *inv_s_vec, &*orig_vec );
264 if (out && includesVerbLevel(verbLevel,Teuchos::VERB_HIGH)) {
265 *out << "\nUnscaled vector "<<vecGetterSetter.getName()<<":\n";
266 Teuchos::OSTab tab(*out);
267 orig_vec->Print(*out);
268 }
269 vecGetterSetter.setVector( orig_vec, origVars );
270 }
271 else {
272 vecGetterSetter.setVector( scaled_vec, origVars );
273 }
274 }
275 else {
276 vecGetterSetter.setVector( null, origVars );
277 }
278
279}
280
281
282// Scale a single vector using a templated policy object to take care of what
283// vector gets used.
284template<class OutArgsVectorGetterSetter>
285void scaleModelFunc(
286 OutArgsVectorGetterSetter vecGetterSetter, // Templated policy object!
287 const EpetraExt::ModelEvaluator::OutArgs &origFuncs,
288 const EpetraExt::ModelEvaluator::OutArgs &funcScalings,
290 Teuchos::FancyOStream * /* out */,
291 Teuchos::EVerbosityLevel /* verbLevel */
292 )
293{
294 TEUCHOS_TEST_FOR_EXCEPT(0==scaledFuncs);
295 Teuchos::RCP<Epetra_Vector>
296 func = vecGetterSetter.getVector(origFuncs);
297 if (!is_null(func) ) {
298 Teuchos::RCP<const Epetra_Vector>
299 funcScaling = vecGetterSetter.getVector(funcScalings);
300 if (!is_null(funcScaling) ) {
301 EpetraExt::scaleModelFuncGivenForwardScaling( *funcScaling, &*func );
302 }
303 }
304 vecGetterSetter.setVector( func, scaledFuncs );
305}
306
307
308} // namespace
309
310
312 const ModelEvaluator &model,
313 ModelEvaluator::InArgs *nominalValues
314 )
315{
316
317 using Teuchos::implicit_cast;
318 typedef ModelEvaluator EME;
319
320#ifdef TEUCHOS_DEBUG
321 TEUCHOS_TEST_FOR_EXCEPT(!nominalValues);
322#endif
323
324 *nominalValues = model.createInArgs();
325
326 if(nominalValues->supports(EME::IN_ARG_x)) {
327 nominalValues->set_x(model.get_x_init());
328 }
329
330 if(nominalValues->supports(EME::IN_ARG_x_dot)) {
331 nominalValues->set_x_dot(model.get_x_dot_init());
332 }
333
334 if(nominalValues->supports(EME::IN_ARG_x_dotdot)) {
335 nominalValues->set_x_dotdot(model.get_x_dotdot_init());
336 }
337
338 for( int l = 0; l < nominalValues->Np(); ++l ) {
339 nominalValues->set_p( l, model.get_p_init(l) );
340 }
341
342 if(nominalValues->supports(EME::IN_ARG_t)) {
343 nominalValues->set_t(model.get_t_init());
344 }
345
346}
347
348
350 const ModelEvaluator &model,
351 ModelEvaluator::InArgs *lowerBounds,
352 ModelEvaluator::InArgs *upperBounds
353 )
354{
355
356 using Teuchos::implicit_cast;
357 typedef ModelEvaluator EME;
358
359#ifdef TEUCHOS_DEBUG
360 TEUCHOS_TEST_FOR_EXCEPT(!lowerBounds);
361 TEUCHOS_TEST_FOR_EXCEPT(!upperBounds);
362#endif
363
364 *lowerBounds = model.createInArgs();
365 *upperBounds = model.createInArgs();
366
367 if(lowerBounds->supports(EME::IN_ARG_x)) {
368 lowerBounds->set_x(model.get_x_lower_bounds());
369 upperBounds->set_x(model.get_x_upper_bounds());
370 }
371
372 for( int l = 0; l < lowerBounds->Np(); ++l ) {
373 lowerBounds->set_p( l, model.get_p_lower_bounds(l) );
374 upperBounds->set_p( l, model.get_p_upper_bounds(l) );
375 }
376
377 if(lowerBounds->supports(EME::IN_ARG_t)) {
378 lowerBounds->set_t(model.get_t_lower_bound());
379 upperBounds->set_t(model.get_t_upper_bound());
380 }
381
382}
383
384
386 const ModelEvaluator::InArgs &origVars,
387 const ModelEvaluator::InArgs &varScalings,
388 ModelEvaluator::InArgs *scaledVars,
389 Teuchos::FancyOStream *out,
390 Teuchos::EVerbosityLevel verbLevel
391 )
392{
393 typedef ModelEvaluator EME;
394
395#ifdef TEUCHOS_DEBUG
396 TEUCHOS_TEST_FOR_EXCEPT(!scaledVars);
397 assertModelVarScalings(varScalings);
398#endif
399
400 if (origVars.supports(EME::IN_ARG_x)) {
401 scaleModelVar( InArgsGetterSetter_x(), origVars, varScalings, scaledVars,
402 out, verbLevel );
403 }
404
405 if (origVars.supports(EME::IN_ARG_x_dot)) {
406 scaleModelVar( InArgsGetterSetter_x_dot(), origVars, varScalings, scaledVars,
407 out, verbLevel );
408 }
409
410 if (origVars.supports(EME::IN_ARG_x_dotdot)) {
411 scaleModelVar( InArgsGetterSetter_x_dotdot(), origVars, varScalings, scaledVars,
412 out, verbLevel );
413 }
414
415 const int Np = origVars.Np();
416 for ( int l = 0; l < Np; ++l ) {
417 scaleModelVar( InArgsGetterSetter_p(l), origVars, varScalings, scaledVars,
418 out, verbLevel );
419 }
420
421 if (origVars.supports(EME::IN_ARG_x_poly)) {
422 TEUCHOS_TEST_FOR_EXCEPTION(
423 !is_null(varScalings.get_x()), std::logic_error,
424 "Error, can't hanlde scaling of x_poly yet!"
425 );
426 scaledVars->set_x_poly(origVars.get_x_poly());
427 }
428
429 if (origVars.supports(EME::IN_ARG_x_dot_poly)) {
430 TEUCHOS_TEST_FOR_EXCEPTION(
431 !is_null(varScalings.get_x()), std::logic_error,
432 "Error, can't hanlde scaling of x_dot_poly yet!"
433 );
434 scaledVars->set_x_dot_poly(origVars.get_x_dot_poly());
435 }
436
437 if (origVars.supports(EME::IN_ARG_x_dotdot_poly)) {
438 TEUCHOS_TEST_FOR_EXCEPTION(
439 !is_null(varScalings.get_x()), std::logic_error,
440 "Error, can't hanlde scaling of x_dotdot_poly yet!"
441 );
442 scaledVars->set_x_dotdot_poly(origVars.get_x_dotdot_poly());
443 }
444
445 if (origVars.supports(EME::IN_ARG_t)) {
446 TEUCHOS_TEST_FOR_EXCEPTION(
447 varScalings.get_t() > 0.0, std::logic_error,
448 "Error, can't hanlde scaling of t yet!"
449 );
450 scaledVars->set_t(origVars.get_t());
451 }
452
453 if (origVars.supports(EME::IN_ARG_alpha)) {
454 TEUCHOS_TEST_FOR_EXCEPTION(
455 varScalings.get_alpha() > 0.0, std::logic_error,
456 "Error, can't hanlde scaling of alpha yet!"
457 );
458 scaledVars->set_alpha(origVars.get_alpha());
459 }
460
461 if (origVars.supports(EME::IN_ARG_beta)) {
462 TEUCHOS_TEST_FOR_EXCEPTION(
463 varScalings.get_beta() > 0.0, std::logic_error,
464 "Error, can't hanlde scaling of beta yet!"
465 );
466 scaledVars->set_beta(origVars.get_beta());
467 }
468
469 // ToDo: Support other input arguments?
470
471}
472
473
475 const ModelEvaluator::InArgs &origLowerBounds,
476 const ModelEvaluator::InArgs &origUpperBounds,
477 const double infBnd,
478 const ModelEvaluator::InArgs &varScalings,
479 ModelEvaluator::InArgs *scaledLowerBounds,
480 ModelEvaluator::InArgs *scaledUpperBounds,
481 Teuchos::FancyOStream *out,
482 Teuchos::EVerbosityLevel verbLevel
483 )
484{
485
486 typedef ModelEvaluator EME;
487
488#ifdef TEUCHOS_DEBUG
489 TEUCHOS_TEST_FOR_EXCEPT(!scaledLowerBounds);
490 TEUCHOS_TEST_FOR_EXCEPT(!scaledUpperBounds);
491 assertModelVarScalings(varScalings);
492#endif
493
494 if (origLowerBounds.supports(EME::IN_ARG_x)) {
495 scaleModelBound(
496 InArgsGetterSetter_x(), origLowerBounds, origUpperBounds, infBnd,
497 varScalings, scaledLowerBounds, scaledUpperBounds,
498 out, verbLevel );
499 }
500
501 if (origLowerBounds.supports(EME::IN_ARG_x_dot)) {
502 scaleModelBound(
503 InArgsGetterSetter_x_dot(), origLowerBounds, origUpperBounds, infBnd,
504 varScalings, scaledLowerBounds, scaledUpperBounds,
505 out, verbLevel );
506 }
507
508 if (origLowerBounds.supports(EME::IN_ARG_x_dotdot)) {
509 scaleModelBound(
510 InArgsGetterSetter_x_dotdot(), origLowerBounds, origUpperBounds, infBnd,
511 varScalings, scaledLowerBounds, scaledUpperBounds,
512 out, verbLevel );
513 }
514
515 const int Np = origLowerBounds.Np();
516 for ( int l = 0; l < Np; ++l ) {
517 scaleModelBound(
518 InArgsGetterSetter_p(l), origLowerBounds, origUpperBounds, infBnd,
519 varScalings, scaledLowerBounds, scaledUpperBounds,
520 out, verbLevel );
521 }
522
523 // ToDo: Support other input arguments?
524
525}
526
527
529 const ModelEvaluator::InArgs &scaledVars,
530 const ModelEvaluator::InArgs &varScalings,
531 ModelEvaluator::InArgs *origVars,
532 Teuchos::FancyOStream *out,
533 Teuchos::EVerbosityLevel verbLevel
534 )
535{
536
537 using Teuchos::RCP;
538 using Teuchos::includesVerbLevel;
539 typedef ModelEvaluator EME;
540
541#ifdef TEUCHOS_DEBUG
542 TEUCHOS_TEST_FOR_EXCEPT(!origVars);
543 assertModelVarScalings(varScalings);
544#endif
545
546 // Print scaling vectors
547
548 if (out && includesVerbLevel(verbLevel,Teuchos::VERB_HIGH)) {
549 RCP<const Epetra_Vector> inv_s_x;
550 if ( scaledVars.supports(EME::IN_ARG_x) &&
551 !is_null(inv_s_x=varScalings.get_x()) )
552 {
553 *out << "\nState inverse scaling vector inv_s_x:\n";
554 Teuchos::OSTab tab(*out);
555 inv_s_x->Print(*out);
556 }
557 }
558
559 // Scal the input varaibles
560
561 if (scaledVars.supports(EME::IN_ARG_x_dot)) {
562 unscaleModelVar( InArgsGetterSetter_x_dot(), scaledVars, varScalings, origVars,
563 out, verbLevel );
564 }
565
566 if (scaledVars.supports(EME::IN_ARG_x_dotdot)) {
567 unscaleModelVar( InArgsGetterSetter_x_dotdot(), scaledVars, varScalings, origVars,
568 out, verbLevel );
569 }
570
571 if (scaledVars.supports(EME::IN_ARG_x)) {
572 unscaleModelVar( InArgsGetterSetter_x(), scaledVars, varScalings, origVars,
573 out, verbLevel );
574 }
575
576 const int Np = scaledVars.Np();
577 for ( int l = 0; l < Np; ++l ) {
578 unscaleModelVar( InArgsGetterSetter_p(l), scaledVars, varScalings, origVars,
579 out, verbLevel );
580 }
581
582 if (scaledVars.supports(EME::IN_ARG_x_dot_poly)) {
583 TEUCHOS_TEST_FOR_EXCEPTION(
584 !is_null(varScalings.get_x()), std::logic_error,
585 "Error, can't hanlde unscaling of x_dot_poly yet!"
586 );
587 origVars->set_x_dot_poly(scaledVars.get_x_dot_poly());
588 }
589
590 if (scaledVars.supports(EME::IN_ARG_x_dotdot_poly)) {
591 TEUCHOS_TEST_FOR_EXCEPTION(
592 !is_null(varScalings.get_x()), std::logic_error,
593 "Error, can't hanlde unscaling of x_dotdot_poly yet!"
594 );
595 origVars->set_x_dotdot_poly(scaledVars.get_x_dotdot_poly());
596 }
597
598 if (scaledVars.supports(EME::IN_ARG_x_poly)) {
599 TEUCHOS_TEST_FOR_EXCEPTION(
600 !is_null(varScalings.get_x()), std::logic_error,
601 "Error, can't hanlde unscaling of x_poly yet!"
602 );
603 origVars->set_x_poly(scaledVars.get_x_poly());
604 }
605
606 if (scaledVars.supports(EME::IN_ARG_t)) {
607 TEUCHOS_TEST_FOR_EXCEPTION(
608 varScalings.get_t() > 0.0, std::logic_error,
609 "Error, can't hanlde unscaling of t yet!"
610 );
611 origVars->set_t(scaledVars.get_t());
612 }
613
614 if (scaledVars.supports(EME::IN_ARG_alpha)) {
615 TEUCHOS_TEST_FOR_EXCEPTION(
616 varScalings.get_alpha() > 0.0, std::logic_error,
617 "Error, can't hanlde unscaling of alpha yet!"
618 );
619 origVars->set_alpha(scaledVars.get_alpha());
620 }
621
622 if (scaledVars.supports(EME::IN_ARG_beta)) {
623 TEUCHOS_TEST_FOR_EXCEPTION(
624 varScalings.get_beta() > 0.0, std::logic_error,
625 "Error, can't hanlde unscaling of beta yet!"
626 );
627 origVars->set_beta(scaledVars.get_beta());
628 }
629
630}
631
632
634 const ModelEvaluator::OutArgs &origFuncs,
635 const ModelEvaluator::InArgs &varScalings,
636 const ModelEvaluator::OutArgs &funcScalings,
637 ModelEvaluator::OutArgs *scaledFuncs,
638 bool *allFuncsWhereScaled,
639 Teuchos::FancyOStream *out,
640 Teuchos::EVerbosityLevel verbLevel
641 )
642{
643
644 using Teuchos::RCP;
645 typedef ModelEvaluator EME;
646
647 TEUCHOS_TEST_FOR_EXCEPT(0==allFuncsWhereScaled);
648
649 *allFuncsWhereScaled = true;
650
651 const int Np = origFuncs.Np();
652 const int Ng = origFuncs.Ng();
653
654 // f
655 if ( origFuncs.supports(EME::OUT_ARG_f) && !is_null(origFuncs.get_f()) ) {
656 scaleModelFunc( OutArgsGetterSetter_f(), origFuncs, funcScalings, scaledFuncs,
657 out, verbLevel );
658 }
659
660 // f_poly
661 if (
662 origFuncs.supports(EME::OUT_ARG_f_poly)
663 && !is_null(origFuncs.get_f_poly())
664 )
665 {
666 TEUCHOS_TEST_FOR_EXCEPTION(
667 !is_null(funcScalings.get_f()), std::logic_error,
668 "Error, we can't handle scaling of f_poly yet!"
669 );
670 scaledFuncs->set_f_poly(origFuncs.get_f_poly());
671 }
672
673 // g(j)
674 for ( int j = 0; j < Ng; ++j ) {
675 scaleModelFunc( OutArgsGetterSetter_g(j), origFuncs, funcScalings, scaledFuncs,
676 out, verbLevel );
677 }
678
679 // W
680 RCP<Epetra_Operator> W;
681 if ( origFuncs.supports(EME::OUT_ARG_W) && !is_null(W=origFuncs.get_W()) ) {
682 bool didScaling = false;
684 varScalings.get_x().get(), funcScalings.get_f().get(),
685 &*W, &didScaling
686 );
687 if(didScaling)
688 scaledFuncs->set_W(W);
689 else
690 *allFuncsWhereScaled = false;
691 }
692
693 // DfDp(l)
694 for ( int l = 0; l < Np; ++l ) {
695 EME::Derivative orig_DfDp_l;
696 if (
697 !origFuncs.supports(EME::OUT_ARG_DfDp,l).none()
698 && !(orig_DfDp_l=origFuncs.get_DfDp(l)).isEmpty()
699 )
700 {
701 EME::Derivative scaled_DfDp_l;
702 bool didScaling = false;
704 orig_DfDp_l, varScalings.get_p(l).get(), funcScalings.get_f().get(),
705 &scaled_DfDp_l, &didScaling
706 );
707 if(didScaling)
708 scaledFuncs->set_DfDp(l,scaled_DfDp_l);
709 else
710 *allFuncsWhereScaled = false;
711 }
712
713 }
714
715 // DgDx_dot(j), DgDx(j), and DgDp(j,l)
716 for ( int j = 0; j < Ng; ++j ) {
717
718 EME::Derivative orig_DgDx_dot_j;
719 if (
720 !origFuncs.supports(EME::OUT_ARG_DgDx_dot,j).none()
721 && !(orig_DgDx_dot_j=origFuncs.get_DgDx_dot(j)).isEmpty()
722 )
723 {
724 EME::Derivative scaled_DgDx_dot_j;
725 bool didScaling = false;
727 orig_DgDx_dot_j, varScalings.get_x().get(), funcScalings.get_g(j).get(),
728 &scaled_DgDx_dot_j, &didScaling
729 );
730 if(didScaling)
731 scaledFuncs->set_DgDx_dot(j,scaled_DgDx_dot_j);
732 else
733 *allFuncsWhereScaled = false;
734 }
735
736 EME::Derivative orig_DgDx_dotdot_j;
737 if (
738 !origFuncs.supports(EME::OUT_ARG_DgDx_dotdot,j).none()
739 && !(orig_DgDx_dotdot_j=origFuncs.get_DgDx_dotdot(j)).isEmpty()
740 )
741 {
742 EME::Derivative scaled_DgDx_dotdot_j;
743 bool didScaling = false;
745 orig_DgDx_dotdot_j, varScalings.get_x().get(), funcScalings.get_g(j).get(),
746 &scaled_DgDx_dotdot_j, &didScaling
747 );
748 if(didScaling)
749 scaledFuncs->set_DgDx_dotdot(j,scaled_DgDx_dotdot_j);
750 else
751 *allFuncsWhereScaled = false;
752 }
753
754 EME::Derivative orig_DgDx_j;
755 if (
756 !origFuncs.supports(EME::OUT_ARG_DgDx,j).none()
757 && !(orig_DgDx_j=origFuncs.get_DgDx(j)).isEmpty()
758 )
759 {
760 EME::Derivative scaled_DgDx_j;
761 bool didScaling = false;
763 orig_DgDx_j, varScalings.get_x().get(), funcScalings.get_g(j).get(),
764 &scaled_DgDx_j, &didScaling
765 );
766 if(didScaling)
767 scaledFuncs->set_DgDx(j,scaled_DgDx_j);
768 else
769 *allFuncsWhereScaled = false;
770 }
771
772 for ( int l = 0; l < Np; ++l ) {
773 EME::Derivative orig_DgDp_j_l;
774 if (
775 !origFuncs.supports(EME::OUT_ARG_DgDp,j,l).none()
776 && !(orig_DgDp_j_l=origFuncs.get_DgDp(j,l)).isEmpty()
777 )
778 {
779 EME::Derivative scaled_DgDp_j_l;
780 bool didScaling = false;
782 orig_DgDp_j_l, varScalings.get_p(l).get(), funcScalings.get_g(j).get(),
783 &scaled_DgDp_j_l, &didScaling
784 );
785 if(didScaling)
786 scaledFuncs->set_DgDp(j,l,scaled_DgDp_j_l);
787 else
788 *allFuncsWhereScaled = false;
789 }
790 }
791 }
792
793}
794
795
796Teuchos::RCP<const Epetra_Vector>
798 Teuchos::RCP<const Epetra_Vector> const& scalingVector
799 )
800{
801 Teuchos::RCP<Epetra_Vector> invScalingVector =
802 Teuchos::rcpWithEmbeddedObj(
803 new Epetra_Vector(scalingVector->Map()),
804 scalingVector
805 );
806 invScalingVector->Reciprocal(*scalingVector);
807 return invScalingVector;
808 // Above, we embedd the forward scaling vector. This is done in order to
809 // achieve the exact same numerics as before this refactoring and to improve
810 // runtime speed and accruacy.
811}
812
813
815 const Epetra_Vector &origVars,
816 const Epetra_Vector &invVarScaling,
817 Epetra_Vector *scaledVars
818 )
819{
820#ifdef TEUCHOS_DEBUG
821 TEUCHOS_TEST_FOR_EXCEPT(!scaledVars);
822 TEUCHOS_TEST_FOR_EXCEPT(!origVars.Map().SameAs(invVarScaling.Map()));
823 TEUCHOS_TEST_FOR_EXCEPT(!origVars.Map().SameAs(scaledVars->Map()));
824#endif
825 const int localDim = origVars.Map().NumMyElements();
826 for ( int i = 0; i < localDim; ++i )
827 (*scaledVars)[i] = origVars[i] / invVarScaling[i];
828}
829
830
832 const Epetra_Vector &/* origLowerBounds */,
833 const Epetra_Vector &/* origUpperBounds */,
834 const double /* infBnd */,
835 const Epetra_Vector &/* invVarScaling */,
836 Epetra_Vector * /* scaledLowerBounds */,
837 Epetra_Vector * /* scaledUpperBounds */
838 )
839{
840 TEUCHOS_TEST_FOR_EXCEPT("ToDo: Implement!");
841}
842
843
845 const Epetra_Vector &origVars,
846 const Epetra_Vector &invVarScaling,
847 Epetra_Vector *scaledVars
848 )
849{
850 TEUCHOS_TEST_FOR_EXCEPT(0==scaledVars);
851 scaledVars->Multiply( 1.0, invVarScaling, origVars, 0.0 );
852}
853
854
856 const Epetra_Vector &fwdFuncScaling,
857 Epetra_Vector *funcs
858 )
859{
860 TEUCHOS_TEST_FOR_EXCEPT(0==funcs);
861 funcs->Multiply( 1.0, fwdFuncScaling, *funcs, 0.0 );
862 // Note: Above is what Epetra_LinearProblem does to scale the RHS and LHS
863 // vectors so this type of argument aliasing must be okay in Epetra!
864}
865
866
868 const Epetra_Vector *invVarScaling,
869 const Epetra_Vector *fwdFuncScaling,
870 Epetra_Operator *funcDerivOp,
871 bool *didScaling
872 )
873{
874 TEUCHOS_TEST_FOR_EXCEPT(0==funcDerivOp);
875 TEUCHOS_TEST_FOR_EXCEPT(0==didScaling);
876 *didScaling = false; // Assume not scaled to start
877 Epetra_RowMatrix *funcDerivRowMatrix
878 = dynamic_cast<Epetra_RowMatrix*>(funcDerivOp);
879 if (funcDerivRowMatrix) {
880 if (fwdFuncScaling)
881 funcDerivRowMatrix->LeftScale(*fwdFuncScaling);
882 if (invVarScaling)
883 funcDerivRowMatrix->RightScale(*invVarScaling);
884 *didScaling = true;
885 // Note that above I do the func scaling before the var scaling since that
886 // is the same order it was done for W in Thyra::EpetraModelEvaluator
887 }
888}
889
890
892 const ModelEvaluator::Derivative &origFuncDeriv,
893 const Epetra_Vector *invVarScaling,
894 const Epetra_Vector *fwdFuncScaling,
895 ModelEvaluator::Derivative *scaledFuncDeriv,
896 bool *didScaling
897 )
898{
899 using Teuchos::RCP;
900 typedef ModelEvaluator EME;
901 TEUCHOS_TEST_FOR_EXCEPT(0==scaledFuncDeriv);
902 TEUCHOS_TEST_FOR_EXCEPT(0==didScaling);
903 *didScaling = false;
904 const RCP<Epetra_MultiVector>
905 funcDerivMv = origFuncDeriv.getMultiVector();
906 const EME::EDerivativeMultiVectorOrientation
907 funcDerivMv_orientation = origFuncDeriv.getMultiVectorOrientation();
908 if(!is_null(funcDerivMv)) {
909 if ( funcDerivMv_orientation == EME::DERIV_MV_BY_COL )
910 {
911 if (fwdFuncScaling) {
912 funcDerivMv->Multiply(1.0, *fwdFuncScaling, *funcDerivMv, 0.0);
913 }
914 if (invVarScaling) {
915 TEUCHOS_TEST_FOR_EXCEPT("ToDo: Scale rows!");
916 //funcDerivMv->Multiply(1.0, *funcDerivMv, *invVarScaling, 0.0);
917 }
918 }
919 else if ( funcDerivMv_orientation == EME::DERIV_TRANS_MV_BY_ROW )
920 {
921 if (invVarScaling) {
922 funcDerivMv->Multiply(1.0, *invVarScaling, *funcDerivMv, 0.0);
923 }
924 if (fwdFuncScaling) {
925 TEUCHOS_TEST_FOR_EXCEPT("ToDo: Scale rows!");
926 //funcDerivMv->Multiply(1.0, *funcDerivMv, *fwdFuncScaling, 0.0);
927 }
928 }
929 else {
930 TEUCHOS_TEST_FOR_EXCEPT("Should not get here!");
931 }
932 *scaledFuncDeriv = EME::Derivative(funcDerivMv,funcDerivMv_orientation);
933 *didScaling = true;
934 }
935 else {
936 RCP<Epetra_Operator>
937 funcDerivOp = origFuncDeriv.getLinearOp();
938 TEUCHOS_TEST_FOR_EXCEPT(is_null(funcDerivOp));
939 scaleModelFuncFirstDerivOp( invVarScaling, fwdFuncScaling,
940 &*funcDerivOp, didScaling );
941 if (didScaling)
942 *scaledFuncDeriv = EME::Derivative(funcDerivOp);
943 }
944}
Class that gets and sets p(l) in an InArgs object.
Class that gets and sets x_dot in an InArgs object.
Class that gets and sets x_dotdot in an InArgs object.
Class that gets and sets x in an InArgs object.
Simple aggregate class that stores a derivative object as a general linear operator or as a multi-vec...
Teuchos::RCP< Epetra_MultiVector > getMultiVector() const
Teuchos::RCP< Epetra_Operator > getLinearOp() const
EDerivativeMultiVectorOrientation getMultiVectorOrientation() const
bool supports(EInArgsMembers arg) const
void set_x_dot(const Teuchos::RCP< const Epetra_Vector > &x_dot)
void set_x_dotdot_poly(const Teuchos::RCP< const Teuchos::Polynomial< Epetra_Vector > > &x_dotdot_poly)
Teuchos::RCP< const Teuchos::Polynomial< Epetra_Vector > > get_x_dotdot_poly() const
Teuchos::RCP< const Teuchos::Polynomial< Epetra_Vector > > get_x_poly() const
Get solution vector Taylor polynomial.
Teuchos::RCP< const Epetra_Vector > get_x_dotdot() const
Teuchos::RCP< const Teuchos::Polynomial< Epetra_Vector > > get_x_dot_poly() const
Get time derivative vector Taylor polynomial.
Teuchos::RCP< const Epetra_Vector > get_p(int l) const
void set_x_dotdot(const Teuchos::RCP< const Epetra_Vector > &x_dotdot)
void set_x_dot_poly(const Teuchos::RCP< const Teuchos::Polynomial< Epetra_Vector > > &x_dot_poly)
Set time derivative vector Taylor polynomial.
void set_x(const Teuchos::RCP< const Epetra_Vector > &x)
Teuchos::RCP< const Epetra_Vector > get_x() const
Set solution vector Taylor polynomial.
void set_x_poly(const Teuchos::RCP< const Teuchos::Polynomial< Epetra_Vector > > &x_poly)
Teuchos::RCP< const Epetra_Vector > get_x_dot() const
void set_p(int l, const Teuchos::RCP< const Epetra_Vector > &p_l)
Teuchos::RCP< Epetra_Operator > get_W() const
bool supports(EOutArgsMembers arg) const
void set_f_poly(const Teuchos::RCP< Teuchos::Polynomial< Epetra_Vector > > &f_poly)
Set residual vector Taylor polynomial.
void set_DgDx_dotdot(int j, const Derivative &DgDx_dotdot_j)
void set_DfDp(int l, const Derivative &DfDp_l)
void set_DgDp(int j, int l, const Derivative &DgDp_j_l)
void set_DgDx(int j, const Derivative &DgDx_j)
Evaluation< Epetra_Vector > get_g(int j) const
Get g(j) where 0 <= j && j < this->Ng().
void set_DgDx_dot(int j, const Derivative &DgDx_dot_j)
Teuchos::RCP< Teuchos::Polynomial< Epetra_Vector > > get_f_poly() const
Get residual vector Taylor polynomial.
Evaluation< Epetra_Vector > get_f() const
Derivative get_DgDp(int j, int l) const
void set_W(const Teuchos::RCP< Epetra_Operator > &W)
Base interface for evaluating a stateless "model".
virtual double get_t_upper_bound() const
virtual Teuchos::RCP< const Epetra_Vector > get_p_upper_bounds(int l) const
virtual double get_t_lower_bound() const
virtual Teuchos::RCP< const Epetra_Vector > get_x_lower_bounds() const
virtual Teuchos::RCP< const Epetra_Vector > get_x_dotdot_init() const
virtual Teuchos::RCP< const Epetra_Vector > get_p_lower_bounds(int l) const
virtual Teuchos::RCP< const Epetra_Vector > get_x_upper_bounds() const
virtual Teuchos::RCP< const Epetra_Vector > get_x_dot_init() const
virtual Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
virtual Teuchos::RCP< const Epetra_Vector > get_x_init() const
virtual InArgs createInArgs() const =0
Class that gets and sets f in an OutArgs object.
Class that gets and sets g(j) in an OutArgs object.
bool SameAs(const Epetra_BlockMap &Map) const
int NumMyElements() const
const Epetra_BlockMap & Map() const
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
virtual int RightScale(const Epetra_Vector &x)=0
virtual int LeftScale(const Epetra_Vector &x)=0
void scaleModelFuncFirstDeriv(const ModelEvaluator::Derivative &origFuncDeriv, const Epetra_Vector *invVarScaling, const Epetra_Vector *fwdFuncScaling, ModelEvaluator::Derivative *scaledFuncDeriv, bool *didScaling)
Scale (in place) an output first-order function derivative object given its function and variable sca...
void scaleModelFuncs(const ModelEvaluator::OutArgs &origFuncs, const ModelEvaluator::InArgs &varScalings, const ModelEvaluator::OutArgs &funcScalings, ModelEvaluator::OutArgs *scaledFuncs, bool *allFuncsWhereScaled, Teuchos::FancyOStream *out=0, Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_LOW)
Scale the output functions and their derivative objects.
void scaleModelVars(const ModelEvaluator::InArgs &origVars, const ModelEvaluator::InArgs &varScalings, ModelEvaluator::InArgs *scaledVars, Teuchos::FancyOStream *out=0, Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_LOW)
Scale the original unscaled variables into the scaled variables.
void scaleModelFuncFirstDerivOp(const Epetra_Vector *invVarScaling, const Epetra_Vector *fwdFuncScaling, Epetra_Operator *funcDerivOp, bool *didScaling)
Scale (in place) an output first-order function derivative object represented as an Epetra_Operator g...
void scaleModelVarsGivenInverseScaling(const Epetra_Vector &origVars, const Epetra_Vector &invVarScaling, Epetra_Vector *scaledVars)
Scale a vector given its inverse scaling vector.
void unscaleModelVars(const ModelEvaluator::InArgs &scaledVars, const ModelEvaluator::InArgs &varScalings, ModelEvaluator::InArgs *origVars, Teuchos::FancyOStream *out=0, Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_LOW)
Unscale the scaled variables.
void gatherModelNominalValues(const ModelEvaluator &model, ModelEvaluator::InArgs *nominalValues)
Gather the nominal values from a model evaluator.
void scaleModelBounds(const ModelEvaluator::InArgs &origLowerBounds, const ModelEvaluator::InArgs &origUpperBounds, const double infBnd, const ModelEvaluator::InArgs &varScalings, ModelEvaluator::InArgs *scaledLowerBounds, ModelEvaluator::InArgs *scaledUpperBounds, Teuchos::FancyOStream *out=0, Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_LOW)
Scale the lower and upper model variable bounds.
Teuchos::RCP< const Epetra_Vector > createInverseModelScalingVector(Teuchos::RCP< const Epetra_Vector > const &scalingVector)
Create an inverse scaling vector.
void scaleModelFuncGivenForwardScaling(const Epetra_Vector &fwdFuncScaling, Epetra_Vector *funcs)
Scale (in place) an output function vector given its forward scaling vector.
void gatherModelBounds(const ModelEvaluator &model, ModelEvaluator::InArgs *lowerBounds, ModelEvaluator::InArgs *upperBounds)
Gather the lower and upper bounds from a model evaluator.
void scaleModelVarBoundsGivenInverseScaling(const Epetra_Vector &origLowerBounds, const Epetra_Vector &origUpperBounds, const double infBnd, const Epetra_Vector &invVarScaling, Epetra_Vector *scaledLowerBounds, Epetra_Vector *scaledUpperBounds)
Scale the model variable bounds.
void unscaleModelVarsGivenInverseScaling(const Epetra_Vector &origVars, const Epetra_Vector &invVarScaling, Epetra_Vector *scaledVars)
Unscale a vector given its inverse scaling vector.