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 (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 #include <EpetraExt_Scale_LinearProblem.h>
00030 
00031 #include <Epetra_LinearProblem.h>
00032 #include <Epetra_CrsMatrix.h>
00033 #include <Epetra_Vector.h>
00034 
00035 namespace EpetraExt {
00036 
00037 LinearProblem_Scale::
00038 ~LinearProblem_Scale()
00039 {
00040   int lsize = (int) lScaleVecs_.size();
00041   for( int i = 0; i < lsize; ++i )
00042     delete lScaleVecs_[i];
00043   int rsize = (int) rScaleVecs_.size();
00044   for( int i = 0; i < rsize; ++i )
00045     delete rScaleVecs_[i];
00046 }
00047 
00048 bool
00049 LinearProblem_Scale::
00050 fwd()
00051 {
00052   Epetra_CrsMatrix & Matrix = *(dynamic_cast<Epetra_CrsMatrix*>(origObj_->GetMatrix()));
00053 
00054   const Epetra_BlockMap & RHSMap = origObj_->GetRHS()->Map();
00055   const Epetra_BlockMap & LHSMap = origObj_->GetLHS()->Map();
00056 
00057   if( iters_ > 0 )
00058   {
00059     if( lScale_ != None && !lScaleVecs_.size() )
00060     {
00061       lScaleVecs_.resize(iters_);
00062       for( int i = 0; i < iters_; ++i )
00063         lScaleVecs_[i] = new Epetra_Vector( RHSMap );
00064     }
00065     if( rScale_ != None && !rScaleVecs_.size() )
00066     {
00067       rScaleVecs_.resize(iters_);
00068       for( int i = 0; i < iters_; ++i )
00069         rScaleVecs_[i] = new Epetra_Vector( LHSMap );
00070     }
00071 
00072     for( int i = 0; i < iters_; ++i )
00073     {
00074       if( lScale_ != None )
00075       {
00076         switch( lScale_ )
00077         {
00078           case Max: Matrix.InvRowMaxs( *(lScaleVecs_[i]) );
00079                     break;
00080           case Sum: Matrix.InvRowSums( *(lScaleVecs_[i]) );
00081                     break;
00082           case Diag: Matrix.ExtractDiagonalCopy( *(lScaleVecs_[i]) );
00083          lScaleVecs_[i]->Reciprocal( *(lScaleVecs_[i]) );
00084          break;
00085           default:  break;
00086         };
00087         if( expFac_ != 1.0 )
00088         {
00089           int numVals = RHSMap.NumMyPoints();
00090           for( int j = 0; j < numVals; ++j ) (*(lScaleVecs_[i]))[j] = pow( (*(lScaleVecs_[i]))[j], expFac_ );
00091         }
00092         newObj_->LeftScale( *lScaleVecs_[i] );
00093       }
00094       if( rScale_ != None )
00095       {
00096         switch( rScale_ )
00097         {
00098           case Max: Matrix.InvColMaxs( *(rScaleVecs_[i]) );
00099                     break;
00100           case Sum: Matrix.InvColSums( *(rScaleVecs_[i]) );
00101                     break;
00102           case Diag: Matrix.ExtractDiagonalCopy( *(rScaleVecs_[i]) );
00103          rScaleVecs_[i]->Reciprocal( *(rScaleVecs_[i]) );
00104          break;
00105           default:  break;
00106         };
00107         if( expFac_ != 1.0 )
00108         {
00109           int numVals = LHSMap.NumMyPoints();
00110           for( int j = 0; j < numVals; ++j ) (*(rScaleVecs_[i]))[j] = pow( (*(rScaleVecs_[i]))[j], expFac_ );
00111         }
00112         newObj_->RightScale( *rScaleVecs_[i] );
00113       }
00114     }
00115   }
00116 
00117   scaled_ = true;
00118 
00119   return true;
00120 }
00121 
00122 bool
00123 LinearProblem_Scale::
00124 rvs()
00125 {
00126   if( !scaled_ ) std::cout << "EpetraExt::LinearProblem_Scale::rvs() : Problem Not Scaled!\n";
00127 
00128   if( iters_ > 0 )
00129   {
00130     for( int i = 0; i < iters_; ++i )
00131     {
00132       int loc = iters_-i-1;
00133       if( rScale_ != None )
00134       {
00135         rScaleVecs_[loc]->Reciprocal( *(rScaleVecs_[loc]) );
00136         newObj_->RightScale( *(rScaleVecs_[loc]) );
00137       }
00138       if( lScale_ != None )
00139       {
00140         lScaleVecs_[loc]->Reciprocal( *(lScaleVecs_[loc]) );
00141         newObj_->LeftScale( *(lScaleVecs_[loc]) );
00142       }
00143     }
00144   }
00145   return true;
00146 }
00147 
00148 } //namespace EpetraExt
00149 

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