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(); }
89 NumBlocks_ = OrigMatrix->NumMyBlockRows();
92 VbrBlocks_.resize(NumBlocks_);
93 VbrBlockCnt_.resize(NumBlocks_);
94 VbrBlockDim_.resize(NumBlocks_);
95 VbrBlockIndices_.resize(NumBlocks_);
96 for(
int i = 0; i < NumBlocks_; ++i )
98 OrigMatrix->ExtractMyBlockRowView( i, VbrBlockDim_[i], VbrBlockCnt_[i], VbrBlockIndices_[i], VbrBlocks_[i] );
101 SVDs_.resize(NumBlocks_);
102 Inverses_.resize(NumBlocks_);
103 for(
int i = 0; i < NumBlocks_; ++i )
105 if( VbrBlockDim_[i] > 1 )
108 SVDs_[i]->SetMatrix( *(VbrBlocks_[i][ VbrBlockCnt_[i]-1 ]) );
110 SVDs_[i]->Invert( rthresh_, athresh_ );
111 Inverses_[i] = SVDs_[i]->InvertedMatrix();
116 double inv = 1. / (*(VbrBlocks_[i][ VbrBlockCnt_[i]-1 ]))(0,0);
123 std::cout <<
"SVDs and Inverses!\n";
124 for(
int i = 0; i < NumBlocks_; ++i )
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;
136 double * currLoc = A;
137 RHSBlocks_.resize(NumBlocks_);
138 for(
int i = 0; i < NumBlocks_; ++i )
141 currLoc += VbrBlockDim_[i];
155 std::cout <<
"-------------------\n";
156 std::cout <<
"BlockJacobi\n";
157 std::cout <<
"-------------------\n";
163 std::multiset<double> SVs;
165 for(
int i = 0; i < NumBlocks_; ++i )
167 if( VbrBlockDim_[i] > 1 )
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] );
177 MaxSV = std::max( MaxSV, 1.0 );
181 std::multiset<double>::iterator iterSI = SVs.begin();
182 std::multiset<double>::iterator endSI = SVs.end();
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;
194 double abs_thresh = athresh_;
195 double rel_thresh = rthresh_;
196 if( thresholding_ == 1 )
198 abs_thresh = MaxSV * rel_thresh;
202 for(
int i = 0; i < NumBlocks_; ++i )
204 if( VbrBlockDim_[i] > 1 )
205 SVDs_[i]->Invert( rel_thresh, abs_thresh );
207 (*Inverses_[i])(0,0) = 1./(*(VbrBlocks_[i][ VbrBlockCnt_[i]-1 ]))(0,0);
210 for(
int i = 0; i < NumBlocks_; ++i )
212 for(
int j = 0; j < (VbrBlockCnt_[i]-1); ++j )
215 VbrBlocks_[i][j]->Multiply(
false,
false, 1.0, *(Inverses_[i]), tempMat, 0.0 );
219 RHSBlocks_[i]->Multiply(
false,
false, 1.0, *(Inverses_[i]), tempMat2, 0.0 );
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]);
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";
237 std::cout << std::endl;
242 std::cout <<
"Min Singular Value: " << MinSV << std::endl;
243 std::cout <<
"Max Singular Value: " << MaxSV << std::endl;
244 std::cout <<
"--------------------\n";