EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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>
50
54
55#include <set>
56
57using std::vector;
58
59namespace 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_);
96 for( int i = 0; i < NumBlocks_; ++i )
97 {
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
149bool
151fwd()
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
250bool
252rvs()
253{
254 return true;
255}
256
257} //namespace EpetraExt
View
Copy
std::vector< Epetra_SerialDenseMatrix * > Inverses_
NewTypeRef operator()(OriginalTypeRef orig)
Analysis of transform operation on original object and construction of new object.
bool rvs()
Reverse transfer of data from new object created in the operator() method call to the orig object inp...
bool fwd()
Forward transfer of data from orig object input in the operator() method call to the new object creat...
std::vector< Epetra_SerialDenseSVD * > SVDs_
std::vector< Epetra_SerialDenseMatrix * > RHSBlocks_
std::vector< Epetra_SerialDenseMatrix ** > VbrBlocks_
bool DistributedGlobal() const
int ExtractView(double **A, int *MyLDA) const
bool IndicesAreGlobal() const
int NumMyBlockRows() const
int ExtractMyBlockRowView(int BlockRow, int &RowDim, int &NumBlockEntries, int *&BlockIndices, Epetra_SerialDenseMatrix **&Values) const
const Epetra_BlockMap & RowMap() const
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.