EpetraExt_BlockJacobi_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_BlockJacobi_LinearProblem.h>
00030 
00031 #include <Epetra_VbrMatrix.h>
00032 #include <Epetra_CrsGraph.h>
00033 #include <Epetra_Map.h>
00034 #include <Epetra_BlockMap.h>
00035 #include <Epetra_MultiVector.h>
00036 #include <Epetra_LinearProblem.h>
00037 
00038 #include <Epetra_SerialDenseMatrix.h>
00039 #include <Epetra_SerialDenseVector.h>
00040 #include <Epetra_SerialDenseSVD.h>
00041 
00042 #include <set>
00043 
00044 using std::vector;
00045 
00046 namespace EpetraExt {
00047 
00048 LinearProblem_BlockJacobi::
00049 ~LinearProblem_BlockJacobi()
00050 {
00051   for( int i = 0; i < NumBlocks_; ++i )
00052   {
00053     if( SVDs_[i] ) delete SVDs_[i];
00054     else if( Inverses_[i] ) delete Inverses_[i];
00055 
00056     if( RHSBlocks_[i] ) delete RHSBlocks_[i];
00057   }
00058 
00059   if( NewProblem_ ) delete NewProblem_;
00060   if( NewMatrix_ ) delete NewMatrix_;
00061 }
00062 
00063 LinearProblem_BlockJacobi::NewTypeRef
00064 LinearProblem_BlockJacobi::
00065 operator()( OriginalTypeRef orig )
00066 {
00067   origObj_ = &orig;
00068 
00069   Epetra_VbrMatrix * OrigMatrix = dynamic_cast<Epetra_VbrMatrix*>( orig.GetMatrix() );
00070 
00071   if( OrigMatrix->RowMap().DistributedGlobal() )
00072   { std::cout << "FAIL for Global!\n"; abort(); }
00073   if( OrigMatrix->IndicesAreGlobal() )
00074   { std::cout << "FAIL for Global Indices!\n"; abort(); }
00075 
00076   NumBlocks_ = OrigMatrix->NumMyBlockRows();
00077 
00078   //extract serial dense matrices from vbr matrix
00079   VbrBlocks_.resize(NumBlocks_);
00080   VbrBlockCnt_.resize(NumBlocks_);
00081   VbrBlockDim_.resize(NumBlocks_);
00082   VbrBlockIndices_.resize(NumBlocks_);
00083   for( int i = 0; i < NumBlocks_; ++i )
00084   {
00085     OrigMatrix->ExtractMyBlockRowView( i, VbrBlockDim_[i], VbrBlockCnt_[i], VbrBlockIndices_[i], VbrBlocks_[i] );
00086   }
00087 
00088   SVDs_.resize(NumBlocks_);
00089   Inverses_.resize(NumBlocks_);
00090   for( int i = 0; i < NumBlocks_; ++i )
00091   {
00092     if( VbrBlockDim_[i] > 1 )
00093     {
00094       SVDs_[i] = new Epetra_SerialDenseSVD();
00095       SVDs_[i]->SetMatrix( *(VbrBlocks_[i][ VbrBlockCnt_[i]-1 ]) );
00096       SVDs_[i]->Factor();
00097       SVDs_[i]->Invert( rthresh_, athresh_ );
00098       Inverses_[i] = SVDs_[i]->InvertedMatrix();
00099     }
00100     else
00101     {
00102       SVDs_[i] = 0;
00103       double inv = 1. / (*(VbrBlocks_[i][ VbrBlockCnt_[i]-1 ]))(0,0);
00104       Inverses_[i] = new Epetra_SerialDenseMatrix( Copy, &inv, 1, 1, 1 );
00105     }
00106   }
00107 
00108   if( verbose_ > 2 )
00109   {
00110     std::cout << "SVDs and Inverses!\n";
00111     for( int i = 0; i < NumBlocks_; ++i )
00112     {
00113       std::cout << "Block: " << i << " Size: " << VbrBlockDim_[i] << std::endl;
00114       if( SVDs_[i] ) SVDs_[i]->Print(std::cout);
00115       std::cout << *(Inverses_[i]) << std::endl;
00116     }
00117   }
00118 
00119   Epetra_MultiVector * RHS = orig.GetRHS();
00120   double * A;
00121   int LDA;
00122   RHS->ExtractView( &A, &LDA );
00123   double * currLoc = A;
00124   RHSBlocks_.resize(NumBlocks_);
00125   for( int i = 0; i < NumBlocks_; ++i )
00126   {
00127     RHSBlocks_[i] = new Epetra_SerialDenseVector( View, currLoc, VbrBlockDim_[i] );
00128     currLoc += VbrBlockDim_[i];
00129   }
00130 
00131   newObj_ = &orig;
00132 
00133   return *newObj_;
00134 }
00135 
00136 bool
00137 LinearProblem_BlockJacobi::
00138 fwd()
00139 {
00140   if( verbose_ > 2 )
00141   {
00142     std::cout << "-------------------\n";
00143     std::cout << "BlockJacobi\n";
00144     std::cout << "-------------------\n";
00145   }
00146   
00147   double MinSV =  1e15;
00148   double MaxSV =  0.0;
00149 
00150   std::multiset<double> SVs;
00151 
00152   for( int i = 0; i < NumBlocks_; ++i )
00153   {
00154     if( VbrBlockDim_[i] > 1 )
00155     {
00156       SVDs_[i]->Factor();
00157       if( SVDs_[i]->S()[0] > MaxSV ) MaxSV = SVDs_[i]->S()[0];
00158       if( SVDs_[i]->S()[VbrBlockDim_[i]-1] < MinSV ) MinSV = SVDs_[i]->S()[VbrBlockDim_[i]-1];
00159       for( int j = 0; j < VbrBlockDim_[i]; ++j ) SVs.insert( SVDs_[i]->S()[j] );
00160     }
00161     else
00162     {
00163       SVs.insert(1.0);
00164       MaxSV = std::max( MaxSV, 1.0 );
00165     }
00166   }
00167 
00168   std::multiset<double>::iterator iterSI = SVs.begin();
00169   std::multiset<double>::iterator endSI = SVs.end();
00170   int i = 0;
00171   if( verbose_ > 2 )
00172   {
00173     std::cout << std::endl;
00174     std::cout << "Singular Values\n";
00175     for( ; iterSI != endSI; ++iterSI, i++ ) std::cout << i << "\t" << *iterSI << std::endl;
00176     std::cout << std::endl;
00177   }
00178 
00179   Epetra_VbrMatrix * OrigMatrix = dynamic_cast<Epetra_VbrMatrix*>( origObj_->GetMatrix() );
00180 
00181   double abs_thresh = athresh_;
00182   double rel_thresh = rthresh_;
00183   if( thresholding_ == 1 )
00184   {
00185     abs_thresh = MaxSV * rel_thresh;
00186     rel_thresh = 0.0;
00187   }
00188       
00189   for( int i = 0; i < NumBlocks_; ++i )
00190   {
00191     if( VbrBlockDim_[i] > 1 )
00192       SVDs_[i]->Invert( rel_thresh, abs_thresh );
00193     else
00194       (*Inverses_[i])(0,0) = 1./(*(VbrBlocks_[i][ VbrBlockCnt_[i]-1 ]))(0,0);
00195   }
00196 
00197   for( int i = 0; i < NumBlocks_; ++i )
00198   {
00199     for( int j = 0; j < (VbrBlockCnt_[i]-1); ++j )
00200     {
00201       Epetra_SerialDenseMatrix tempMat( *(VbrBlocks_[i][j]) );
00202       VbrBlocks_[i][j]->Multiply( false, false, 1.0, *(Inverses_[i]), tempMat, 0.0 );
00203     }
00204 
00205     Epetra_SerialDenseMatrix tempMat2( *(RHSBlocks_[i]) );
00206     RHSBlocks_[i]->Multiply( false, false, 1.0, *(Inverses_[i]), tempMat2, 0.0 );
00207 
00208     if( verbose_ > 2 )
00209     {
00210       std::cout << "DiagBlock: " << i << std::endl;
00211       std::cout << *(VbrBlocks_[i][VbrBlockCnt_[i]-1]);
00212       std::cout << "RHSBlock: " << i << std::endl;
00213       std::cout << *(RHSBlocks_[i]);
00214     }
00215   }
00216 
00217   if( verbose_ > 2 )
00218   {
00219     std::cout << "Block Jacobi'd Matrix!\n";
00220     if( removeDiag_ ) std::cout << *NewMatrix_ << std::endl;
00221     else              std::cout << *(dynamic_cast<Epetra_VbrMatrix*>(origObj_->GetMatrix())) << std::endl;
00222     std::cout << "Block Jacobi'd RHS!\n";
00223     std::cout << *(origObj_->GetRHS());
00224     std::cout << std::endl;
00225   }
00226 
00227   if( verbose_ > 0 )
00228   {
00229     std::cout << "Min Singular Value: " << MinSV << std::endl;
00230     std::cout << "Max Singular Value: " << MaxSV << std::endl;
00231     std::cout << "--------------------\n";
00232   }
00233 
00234   return true;
00235 }
00236 
00237 bool
00238 LinearProblem_BlockJacobi::
00239 rvs()
00240 {
00241   return true;
00242 }
00243 
00244 } //namespace EpetraExt

Generated on Wed May 12 21:24:45 2010 for EpetraExt by  doxygen 1.4.7