EpetraExt Package Browser (Single Doxygen Collection) Development
EpetraExt_Scale_LinearProblem.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2011) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 
00042 #include <EpetraExt_Scale_LinearProblem.h>
00043 
00044 #include <Epetra_LinearProblem.h>
00045 #include <Epetra_CrsMatrix.h>
00046 #include <Epetra_Vector.h>
00047 
00048 namespace EpetraExt {
00049 
00050 LinearProblem_Scale::
00051 ~LinearProblem_Scale()
00052 {
00053   int lsize = (int) lScaleVecs_.size();
00054   for( int i = 0; i < lsize; ++i )
00055     delete lScaleVecs_[i];
00056   int rsize = (int) rScaleVecs_.size();
00057   for( int i = 0; i < rsize; ++i )
00058     delete rScaleVecs_[i];
00059 }
00060 
00061 bool
00062 LinearProblem_Scale::
00063 fwd()
00064 {
00065   Epetra_CrsMatrix & Matrix = *(dynamic_cast<Epetra_CrsMatrix*>(origObj_->GetMatrix()));
00066 
00067   const Epetra_BlockMap & RHSMap = origObj_->GetRHS()->Map();
00068   const Epetra_BlockMap & LHSMap = origObj_->GetLHS()->Map();
00069 
00070   if( iters_ > 0 )
00071   {
00072     if( lScale_ != None && !lScaleVecs_.size() )
00073     {
00074       lScaleVecs_.resize(iters_);
00075       for( int i = 0; i < iters_; ++i )
00076         lScaleVecs_[i] = new Epetra_Vector( RHSMap );
00077     }
00078     if( rScale_ != None && !rScaleVecs_.size() )
00079     {
00080       rScaleVecs_.resize(iters_);
00081       for( int i = 0; i < iters_; ++i )
00082         rScaleVecs_[i] = new Epetra_Vector( LHSMap );
00083     }
00084 
00085     for( int i = 0; i < iters_; ++i )
00086     {
00087       if( lScale_ != None )
00088       {
00089         switch( lScale_ )
00090         {
00091           case Max: Matrix.InvRowMaxs( *(lScaleVecs_[i]) );
00092                     break;
00093           case Sum: Matrix.InvRowSums( *(lScaleVecs_[i]) );
00094                     break;
00095           case Diag: Matrix.ExtractDiagonalCopy( *(lScaleVecs_[i]) );
00096          lScaleVecs_[i]->Reciprocal( *(lScaleVecs_[i]) );
00097          break;
00098           default:  break;
00099         };
00100         if( expFac_ != 1.0 )
00101         {
00102           int numVals = RHSMap.NumMyPoints();
00103           for( int j = 0; j < numVals; ++j ) (*(lScaleVecs_[i]))[j] = pow( (*(lScaleVecs_[i]))[j], expFac_ );
00104         }
00105         newObj_->LeftScale( *lScaleVecs_[i] );
00106       }
00107       if( rScale_ != None )
00108       {
00109         switch( rScale_ )
00110         {
00111           case Max: Matrix.InvColMaxs( *(rScaleVecs_[i]) );
00112                     break;
00113           case Sum: Matrix.InvColSums( *(rScaleVecs_[i]) );
00114                     break;
00115           case Diag: Matrix.ExtractDiagonalCopy( *(rScaleVecs_[i]) );
00116          rScaleVecs_[i]->Reciprocal( *(rScaleVecs_[i]) );
00117          break;
00118           default:  break;
00119         };
00120         if( expFac_ != 1.0 )
00121         {
00122           int numVals = LHSMap.NumMyPoints();
00123           for( int j = 0; j < numVals; ++j ) (*(rScaleVecs_[i]))[j] = pow( (*(rScaleVecs_[i]))[j], expFac_ );
00124         }
00125         newObj_->RightScale( *rScaleVecs_[i] );
00126       }
00127     }
00128   }
00129 
00130   scaled_ = true;
00131 
00132   return true;
00133 }
00134 
00135 bool
00136 LinearProblem_Scale::
00137 rvs()
00138 {
00139   if( !scaled_ ) std::cout << "EpetraExt::LinearProblem_Scale::rvs() : Problem Not Scaled!\n";
00140 
00141   if( iters_ > 0 )
00142   {
00143     for( int i = 0; i < iters_; ++i )
00144     {
00145       int loc = iters_-i-1;
00146       if( rScale_ != None )
00147       {
00148         rScaleVecs_[loc]->Reciprocal( *(rScaleVecs_[loc]) );
00149         newObj_->RightScale( *(rScaleVecs_[loc]) );
00150       }
00151       if( lScale_ != None )
00152       {
00153         lScaleVecs_[loc]->Reciprocal( *(lScaleVecs_[loc]) );
00154         newObj_->LeftScale( *(lScaleVecs_[loc]) );
00155       }
00156     }
00157   }
00158   return true;
00159 }
00160 
00161 } //namespace EpetraExt
00162 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines