EpetraExt  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EpetraExt_BlockJacobi_LinearProblem.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 
43 
44 #include <Epetra_VbrMatrix.h>
45 #include <Epetra_CrsGraph.h>
46 #include <Epetra_Map.h>
47 #include <Epetra_BlockMap.h>
48 #include <Epetra_MultiVector.h>
49 #include <Epetra_LinearProblem.h>
50 
51 #include <Epetra_SerialDenseMatrix.h>
52 #include <Epetra_SerialDenseVector.h>
53 #include <Epetra_SerialDenseSVD.h>
54 
55 #include <set>
56 
57 using std::vector;
58 
59 namespace EpetraExt {
60 
63 {
64  for( int i = 0; i < NumBlocks_; ++i )
65  {
66  if( SVDs_[i] ) delete SVDs_[i];
67  else if( Inverses_[i] ) delete Inverses_[i];
68 
69  if( RHSBlocks_[i] ) delete RHSBlocks_[i];
70  }
71 
72  if( NewProblem_ ) delete NewProblem_;
73  if( NewMatrix_ ) delete NewMatrix_;
74 }
75 
79 {
80  origObj_ = &orig;
81 
82  Epetra_VbrMatrix * OrigMatrix = dynamic_cast<Epetra_VbrMatrix*>( orig.GetMatrix() );
83 
84  if( OrigMatrix->RowMap().DistributedGlobal() )
85  { std::cout << "FAIL for Global!\n"; abort(); }
86  if( OrigMatrix->IndicesAreGlobal() )
87  { std::cout << "FAIL for Global Indices!\n"; abort(); }
88 
89  NumBlocks_ = OrigMatrix->NumMyBlockRows();
90 
91  //extract serial dense matrices from vbr matrix
92  VbrBlocks_.resize(NumBlocks_);
93  VbrBlockCnt_.resize(NumBlocks_);
94  VbrBlockDim_.resize(NumBlocks_);
95  VbrBlockIndices_.resize(NumBlocks_);
96  for( int i = 0; i < NumBlocks_; ++i )
97  {
98  OrigMatrix->ExtractMyBlockRowView( i, VbrBlockDim_[i], VbrBlockCnt_[i], VbrBlockIndices_[i], VbrBlocks_[i] );
99  }
100 
101  SVDs_.resize(NumBlocks_);
102  Inverses_.resize(NumBlocks_);
103  for( int i = 0; i < NumBlocks_; ++i )
104  {
105  if( VbrBlockDim_[i] > 1 )
106  {
107  SVDs_[i] = new Epetra_SerialDenseSVD();
108  SVDs_[i]->SetMatrix( *(VbrBlocks_[i][ VbrBlockCnt_[i]-1 ]) );
109  SVDs_[i]->Factor();
110  SVDs_[i]->Invert( rthresh_, athresh_ );
111  Inverses_[i] = SVDs_[i]->InvertedMatrix();
112  }
113  else
114  {
115  SVDs_[i] = 0;
116  double inv = 1. / (*(VbrBlocks_[i][ VbrBlockCnt_[i]-1 ]))(0,0);
117  Inverses_[i] = new Epetra_SerialDenseMatrix( Copy, &inv, 1, 1, 1 );
118  }
119  }
120 
121  if( verbose_ > 2 )
122  {
123  std::cout << "SVDs and Inverses!\n";
124  for( int i = 0; i < NumBlocks_; ++i )
125  {
126  std::cout << "Block: " << i << " Size: " << VbrBlockDim_[i] << std::endl;
127  if( SVDs_[i] ) SVDs_[i]->Print(std::cout);
128  std::cout << *(Inverses_[i]) << std::endl;
129  }
130  }
131 
132  Epetra_MultiVector * RHS = orig.GetRHS();
133  double * A;
134  int LDA;
135  RHS->ExtractView( &A, &LDA );
136  double * currLoc = A;
137  RHSBlocks_.resize(NumBlocks_);
138  for( int i = 0; i < NumBlocks_; ++i )
139  {
140  RHSBlocks_[i] = new Epetra_SerialDenseVector( View, currLoc, VbrBlockDim_[i] );
141  currLoc += VbrBlockDim_[i];
142  }
143 
144  newObj_ = &orig;
145 
146  return *newObj_;
147 }
148 
149 bool
152 {
153  if( verbose_ > 2 )
154  {
155  std::cout << "-------------------\n";
156  std::cout << "BlockJacobi\n";
157  std::cout << "-------------------\n";
158  }
159 
160  double MinSV = 1e15;
161  double MaxSV = 0.0;
162 
163  std::multiset<double> SVs;
164 
165  for( int i = 0; i < NumBlocks_; ++i )
166  {
167  if( VbrBlockDim_[i] > 1 )
168  {
169  SVDs_[i]->Factor();
170  if( SVDs_[i]->S()[0] > MaxSV ) MaxSV = SVDs_[i]->S()[0];
171  if( SVDs_[i]->S()[VbrBlockDim_[i]-1] < MinSV ) MinSV = SVDs_[i]->S()[VbrBlockDim_[i]-1];
172  for( int j = 0; j < VbrBlockDim_[i]; ++j ) SVs.insert( SVDs_[i]->S()[j] );
173  }
174  else
175  {
176  SVs.insert(1.0);
177  MaxSV = std::max( MaxSV, 1.0 );
178  }
179  }
180 
181  std::multiset<double>::iterator iterSI = SVs.begin();
182  std::multiset<double>::iterator endSI = SVs.end();
183  int i = 0;
184  if( verbose_ > 2 )
185  {
186  std::cout << std::endl;
187  std::cout << "Singular Values\n";
188  for( ; iterSI != endSI; ++iterSI, i++ ) std::cout << i << "\t" << *iterSI << std::endl;
189  std::cout << std::endl;
190  }
191 
192  Epetra_VbrMatrix * OrigMatrix = dynamic_cast<Epetra_VbrMatrix*>( origObj_->GetMatrix() );
193 
194  double abs_thresh = athresh_;
195  double rel_thresh = rthresh_;
196  if( thresholding_ == 1 )
197  {
198  abs_thresh = MaxSV * rel_thresh;
199  rel_thresh = 0.0;
200  }
201 
202  for( int i = 0; i < NumBlocks_; ++i )
203  {
204  if( VbrBlockDim_[i] > 1 )
205  SVDs_[i]->Invert( rel_thresh, abs_thresh );
206  else
207  (*Inverses_[i])(0,0) = 1./(*(VbrBlocks_[i][ VbrBlockCnt_[i]-1 ]))(0,0);
208  }
209 
210  for( int i = 0; i < NumBlocks_; ++i )
211  {
212  for( int j = 0; j < (VbrBlockCnt_[i]-1); ++j )
213  {
214  Epetra_SerialDenseMatrix tempMat( *(VbrBlocks_[i][j]) );
215  VbrBlocks_[i][j]->Multiply( false, false, 1.0, *(Inverses_[i]), tempMat, 0.0 );
216  }
217 
218  Epetra_SerialDenseMatrix tempMat2( *(RHSBlocks_[i]) );
219  RHSBlocks_[i]->Multiply( false, false, 1.0, *(Inverses_[i]), tempMat2, 0.0 );
220 
221  if( verbose_ > 2 )
222  {
223  std::cout << "DiagBlock: " << i << std::endl;
224  std::cout << *(VbrBlocks_[i][VbrBlockCnt_[i]-1]);
225  std::cout << "RHSBlock: " << i << std::endl;
226  std::cout << *(RHSBlocks_[i]);
227  }
228  }
229 
230  if( verbose_ > 2 )
231  {
232  std::cout << "Block Jacobi'd Matrix!\n";
233  if( removeDiag_ ) std::cout << *NewMatrix_ << std::endl;
234  else std::cout << *(dynamic_cast<Epetra_VbrMatrix*>(origObj_->GetMatrix())) << std::endl;
235  std::cout << "Block Jacobi'd RHS!\n";
236  std::cout << *(origObj_->GetRHS());
237  std::cout << std::endl;
238  }
239 
240  if( verbose_ > 0 )
241  {
242  std::cout << "Min Singular Value: " << MinSV << std::endl;
243  std::cout << "Max Singular Value: " << MaxSV << std::endl;
244  std::cout << "--------------------\n";
245  }
246 
247  return true;
248 }
249 
250 bool
253 {
254  return true;
255 }
256 
257 } //namespace EpetraExt
bool fwd()
Forward transfer of data from orig object input in the operator() method call to the new object creat...
bool rvs()
Reverse transfer of data from new object created in the operator() method call to the orig object inp...
NewTypeRef operator()(OriginalTypeRef orig)
Analysis of transform operation on original object and construction of new object.