EpetraExt Development
EpetraExt_BTF_LinearProblem.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2011) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 
00042 #include <EpetraExt_BTF_LinearProblem.h>
00043 
00044 #include <Epetra_CrsMatrix.h>
00045 #include <Epetra_VbrMatrix.h>
00046 #include <Epetra_CrsGraph.h>
00047 #include <Epetra_Map.h>
00048 #include <Epetra_BlockMap.h>
00049 #include <Epetra_MultiVector.h>
00050 #include <Epetra_LinearProblem.h>
00051 
00052 #include <set>
00053 
00054 using std::vector;
00055 using std::map;
00056 using std::set;
00057 
00058 #define MATTRANS_F77 F77_FUNC(mattrans,MATTRANS)
00059 #define GENBTF_F77   F77_FUNC(genbtf,GENBTF)
00060                                                                                                   
00061 extern "C" {
00062 extern void MATTRANS_F77( int*, int*, int*, int*, int*, int* );
00063 extern void GENBTF_F77( int*, int*, int*, int*, int*, int*, int*, int*, int*,
00064                      int*, int*, int*, int*, int*, int*, int*, int*, int*,
00065                      int*, int*, int* );
00066 }
00067 
00068 namespace EpetraExt {
00069 
00070 LinearProblem_BTF::
00071 ~LinearProblem_BTF()
00072 {
00073   deleteNewObjs_();
00074 }
00075 
00076 void
00077 LinearProblem_BTF::
00078 deleteNewObjs_()
00079 {
00080   if( NewProblem_ ) delete NewProblem_;
00081 
00082   if( NewMatrix_ ) delete NewMatrix_;
00083 
00084   if( NewLHS_ ) delete NewLHS_;
00085   if( NewRHS_ ) delete NewRHS_;
00086 
00087   if( NewMap_ ) delete NewMap_;
00088 
00089   for( int i = 0; i < Blocks_.size(); ++i )
00090     for( int j = 0; j < Blocks_[i].size(); ++j )
00091       delete Blocks_[i][j];
00092 }
00093 
00094 LinearProblem_BTF::NewTypeRef
00095 LinearProblem_BTF::
00096 operator()( OriginalTypeRef orig )
00097 {
00098   changedLP_ = false;
00099 
00100   //check if there is a previous analysis and if it is valid
00101   if( &orig == origObj_ && NewProblem_ )
00102   {
00103     int * indices;
00104     double * values;
00105     int currCnt;
00106     int numRows = OrigMatrix_->NumMyRows();
00107 
00108     for( int i = 0; i < numRows; ++i )
00109       if( ZeroElements_[i].size() )
00110       {
00111         int loc = 0;
00112         OrigMatrix_->ExtractMyRowView( i, currCnt, values, indices );
00113         for( int j = 0; j < currCnt; ++j )
00114           if( ZeroElements_[i].count( indices[j] ) )
00115           {
00116             if( values[j] > threshold_ ) changedLP_ = true;
00117             ++loc;
00118           }
00119       }
00120 
00121 //    changedLP_ = true;
00122 
00123     //if changed, dump all the old stuff and start over
00124     if( changedLP_ )
00125       deleteNewObjs_();
00126     else
00127       return *newObj_;
00128   }
00129     
00130   origObj_ = &orig;
00131 
00132   OrigMatrix_ = dynamic_cast<Epetra_CrsMatrix*>(orig.GetMatrix());
00133   OrigLHS_ = orig.GetLHS();
00134   OrigRHS_ = orig.GetRHS();
00135 
00136   if( OrigMatrix_->RowMap().DistributedGlobal() )
00137   { cout << "FAIL for Global!\n"; abort(); }
00138   if( OrigMatrix_->IndicesAreGlobal() )
00139   { cout << "FAIL for Global Indices!\n"; abort(); }
00140 
00141   int n = OrigMatrix_->NumMyRows();
00142   int nnz = OrigMatrix_->NumMyNonzeros();
00143 
00144 //  cout << "Orig Matrix:\n";
00145 //  cout << *OrigMatrix_ << endl;
00146 
00147   //create std CRS format
00148   //also create graph without zero elements
00149   vector<int> ia(n+1,0);
00150   int maxEntries = OrigMatrix_->MaxNumEntries();
00151   vector<int> ja_tmp(maxEntries);
00152   vector<double> jva_tmp(maxEntries);
00153   vector<int> ja(nnz);
00154   int cnt;
00155 
00156   const Epetra_BlockMap & OldRowMap = OrigMatrix_->RowMap();
00157   const Epetra_BlockMap & OldColMap = OrigMatrix_->ColMap();
00158   Epetra_CrsGraph strippedGraph( Copy, OldRowMap, OldColMap, 0 );
00159   ZeroElements_.resize(n);
00160 
00161   for( int i = 0; i < n; ++i )
00162   {
00163     ZeroElements_[i].clear();
00164     OrigMatrix_->ExtractMyRowCopy( i, maxEntries, cnt, &jva_tmp[0], &ja_tmp[0] );
00165     ia[i+1] = ia[i];
00166     for( int j = 0; j < cnt; ++j )
00167     {
00168       if( fabs(jva_tmp[j]) > threshold_ )
00169         ja[ ia[i+1]++ ] = ja_tmp[j];
00170       else
00171         ZeroElements_[i].insert( ja_tmp[j] );
00172     }
00173 
00174     int new_cnt = ia[i+1] - ia[i];
00175     strippedGraph.InsertMyIndices( i, new_cnt, &ja[ ia[i] ] );
00176   }
00177   nnz = ia[n];
00178   strippedGraph.FillComplete();
00179 
00180   if( verbose_ > 2 )
00181   {
00182     cout << "Stripped Graph\n";
00183     cout << strippedGraph;
00184   }
00185 
00186   vector<int> iat(n+1,0);
00187   vector<int> jat(nnz);
00188   for( int i = 0; i < n; ++i )
00189     for( int j = ia[i]; j < ia[i+1]; ++j )
00190       ++iat[ ja[j]+1 ];
00191   for( int i = 0; i < n; ++i )
00192     iat[i+1] += iat[i];
00193   for( int i = 0; i < n; ++i )
00194     for( int j = ia[i]; j < ia[i+1]; ++j )
00195       jat[ iat[ ja[j] ]++ ] = i;
00196   for( int i = 0; i < n; ++i )
00197     iat[n-i] = iat[n-i-1];
00198   iat[0] = 0;
00199 
00200   //convert to Fortran indexing
00201   for( int i = 0; i < n+1; ++i ) ++ia[i];
00202   for( int i = 0; i < nnz; ++i ) ++ja[i];
00203   for( int i = 0; i < n+1; ++i ) ++iat[i];
00204   for( int i = 0; i < nnz; ++i ) ++jat[i];
00205 
00206   if( verbose_ > 2 )
00207   {
00208     cout << "-----------------------------------------\n";
00209     cout << "CRS Format Graph\n";
00210     cout << "-----------------------------------------\n";
00211     for( int i = 0; i < n; ++i )
00212     {
00213       cout << i+1 << ": " << ia[i+1] << ": ";
00214       for( int j = ia[i]-1; j<ia[i+1]-1; ++j )
00215         cout << " " << ja[j];
00216       cout << endl;
00217     }
00218     cout << "-----------------------------------------\n";
00219   }
00220 
00221   if( verbose_ > 2 )
00222   {
00223     cout << "-----------------------------------------\n";
00224     cout << "CCS Format Graph\n";
00225     cout << "-----------------------------------------\n";
00226     for( int i = 0; i < n; ++i )
00227     {
00228       cout << i+1 << ": " << iat[i+1] << ": ";
00229       for( int j = iat[i]-1; j<iat[i+1]-1; ++j )
00230         cout << " " << jat[j];
00231       cout << endl;
00232     }
00233     cout << "-----------------------------------------\n";
00234   }
00235 
00236   vector<int> w(10*n);
00237 
00238   vector<int> rowperm(n);
00239   vector<int> colperm(n);
00240 
00241   //horizontal block
00242   int nhrows, nhcols, hrzcmp;
00243   //square block
00244   int nsrows, sqcmpn;
00245   //vertial block
00246   int nvrows, nvcols, vrtcmp;
00247 
00248   vector<int> rcmstr(n+1);
00249   vector<int> ccmstr(n+1);
00250 
00251   int msglvl = 0;
00252   int output = 6;
00253 
00254   GENBTF_F77( &n, &n, &iat[0], &jat[0], &ia[0], &ja[0], &w[0],
00255           &rowperm[0], &colperm[0], &nhrows, &nhcols,
00256           &hrzcmp, &nsrows, &sqcmpn, &nvrows, &nvcols, &vrtcmp,
00257           &rcmstr[0], &ccmstr[0], &msglvl, &output );
00258 
00259   //convert back to C indexing
00260   for( int i = 0; i < n; ++i )
00261   {
00262     --rowperm[i];
00263     --colperm[i];
00264   }
00265   for( int i = 0; (i<n+1) && (rcmstr[i]!=n+1); ++i )
00266   {
00267     --rcmstr[i];
00268     --ccmstr[i];
00269   }
00270 
00271   if( verbose_ > 0 )
00272   {
00273     cout << "-----------------------------------------\n";
00274     cout << "BTF Output\n";
00275     cout << "-----------------------------------------\n";
00276 //    cout << "RowPerm and ColPerm\n";
00277 //    for( int i = 0; i<n; ++i )
00278 //      cout << rowperm[i] << "\t" << colperm[i] << endl;
00279     if( hrzcmp )
00280     {
00281       cout << "Num Horizontal: Rows, Cols, Comps\n";
00282       cout << nhrows << "\t" << nhcols << "\t" << hrzcmp << endl;
00283     }
00284     cout << "Num Square: Rows, Comps\n";
00285     cout << nsrows << "\t" << sqcmpn << endl;
00286     if( vrtcmp )
00287     {
00288       cout << "Num Vertical: Rows, Cols, Comps\n";
00289       cout << nvrows << "\t" << nvcols << "\t" << vrtcmp << endl;
00290     }
00291 //    cout << "Row, Col of upper left pt in blocks\n";
00292 //    for( int i = 0; (i<n+1) && (rcmstr[i]!=n+1); ++i )
00293 //      cout << i << " " << rcmstr[i] << " " << ccmstr[i] << endl;
00294 //    cout << "-----------------------------------------\n";
00295   }
00296 
00297   if( hrzcmp || vrtcmp )
00298   {
00299     cout << "FAILED! hrz cmp's:" << hrzcmp << " vrtcmp: " << vrtcmp << endl;
00300     exit(0);
00301   }
00302 
00303   rcmstr[sqcmpn] = n;
00304 
00305   //convert rowperm to OLD->NEW
00306   //reverse ordering of permutation to get upper triangular
00307   vector<int> rowperm_t( n );
00308   vector<int> colperm_t( n );
00309   for( int i = 0; i < n; ++i )
00310   {
00311     rowperm_t[i] = rowperm[i];
00312     colperm_t[i] = colperm[i];
00313   }
00314 
00315   //Generate New Domain and Range Maps
00316   //for now, assume they start out as identical
00317   OldGlobalElements_.resize(n);
00318   OldRowMap.MyGlobalElements( &OldGlobalElements_[0] );
00319 
00320   vector<int> newDomainElements( n );
00321   vector<int> newRangeElements( n );
00322   for( int i = 0; i < n; ++i )
00323   {
00324     newDomainElements[ i ] = OldGlobalElements_[ rowperm_t[i] ];
00325     newRangeElements[ i ] = OldGlobalElements_[ colperm_t[i] ];
00326   }
00327 
00328   //Setup New Block Map Info
00329   Blocks_.resize( sqcmpn );
00330   BlockDim_.resize( sqcmpn );
00331   for( int i = 0; i < sqcmpn; ++i )
00332   {
00333     BlockDim_[i] = rcmstr[i+1]-rcmstr[i];
00334     for( int j = rcmstr[i]; j < rcmstr[i+1]; ++j )
00335     {
00336       BlockRowMap_[ newDomainElements[j] ] = i;
00337       SubBlockRowMap_[ newDomainElements[j] ] = j-rcmstr[i];
00338       BlockColMap_[ newRangeElements[j] ] = i;
00339       SubBlockColMap_[ newRangeElements[j] ] = j-rcmstr[i];
00340     }
00341   }
00342 
00343   if( verbose_ > 2 )
00344   {
00345 /*
00346     cout << "Block Mapping!\n";
00347     cout << "--------------------------\n";
00348     for( int i = 0; i < n; ++i )
00349     {
00350       cout << "Row: " << newDomainElements[i] << " " << BlockRowMap_[newDomainElements[i]] << " " <<
00351               SubBlockRowMap_[newDomainElements[i]] << "\t" << "Col: " << newRangeElements[i] << " " <<
00352               BlockColMap_[newRangeElements[i]] << " " << SubBlockColMap_[newRangeElements[i]] << endl;
00353     }
00354     for( int i = 0; i < sqcmpn; ++i )
00355       cout << "BlockDim: " << i << " " << BlockDim_[i] << endl;
00356     cout << "--------------------------\n";
00357 */
00358     int MinSize = 1000000000;
00359     int MaxSize = 0;
00360     for( int i = 0; i < sqcmpn; ++i )
00361     {
00362       if( MinSize > BlockDim_[i] ) MinSize = BlockDim_[i];
00363       if( MaxSize < BlockDim_[i] ) MaxSize = BlockDim_[i];
00364     }
00365     cout << "Min Block Size: " << MinSize << " " << "Max Block Size: " << MaxSize << endl;
00366   }
00367 
00368   vector<int> myBlockElements( sqcmpn );
00369   for( int i = 0; i < sqcmpn; ++i ) myBlockElements[i] = i;
00370   NewMap_ = new Epetra_BlockMap( sqcmpn, sqcmpn, &myBlockElements[0], &BlockDim_[0], 0, OldRowMap.Comm() );
00371 
00372   if( verbose_ > 2 )
00373   {
00374     cout << "New Block Map!\n";
00375     cout << *NewMap_;
00376   }
00377 
00378   //setup new graph
00379   vector< set<int> > crsBlocks( sqcmpn );
00380   BlockCnt_.resize( sqcmpn );
00381   int maxLength = strippedGraph.MaxNumIndices();
00382   vector<int> sIndices( maxLength );
00383   int currLength;
00384   for( int i = 0; i < n; ++i )
00385   {
00386     strippedGraph.ExtractGlobalRowCopy( OldGlobalElements_[i], maxLength, currLength, &sIndices[0] );
00387     for( int j = 0; j < currLength; ++j )
00388       crsBlocks[ BlockRowMap_[ OldGlobalElements_[i] ] ].insert( BlockColMap_[ sIndices[j] ] );
00389   }
00390 
00391   for( int i = 0; i < sqcmpn; ++i )
00392   {
00393     BlockCnt_[i] = crsBlocks[i].size();
00394     Blocks_[i].resize( BlockCnt_[i] );
00395   }
00396 
00397   NewBlockRows_.resize( sqcmpn );
00398   for( int i = 0; i < sqcmpn; ++i )
00399   {
00400     NewBlockRows_[i] = vector<int>( crsBlocks[i].begin(), crsBlocks[i].end() );
00401     for( int j = 0; j < BlockCnt_[i]; ++j )
00402     {
00403       Blocks_[i][j] = new Epetra_SerialDenseMatrix();
00404       Blocks_[i][j]->Shape( BlockDim_[i], BlockDim_[ NewBlockRows_[i][j] ] );
00405     }
00406   }
00407 
00408   //put data in Blocks_ and new LHS and RHS
00409   NewLHS_ = new Epetra_MultiVector( *NewMap_, 1 );
00410   NewRHS_ = new Epetra_MultiVector( *NewMap_, 1 );
00411 
00412   maxLength = OrigMatrix_->MaxNumEntries();
00413   vector<int> indices( maxLength );
00414   vector<double> values( maxLength );
00415   for( int i = 0; i < n; ++i )
00416   {
00417     int BlockRow = BlockRowMap_[ OldGlobalElements_[i] ];
00418     int SubBlockRow = SubBlockRowMap_[ OldGlobalElements_[i] ];
00419     OrigMatrix_->ExtractGlobalRowCopy( OldGlobalElements_[i], maxLength, currLength, &values[0], &indices[0] );
00420     for( int j = 0; j < currLength; ++j )
00421     {
00422       int BlockCol = BlockColMap_[ indices[j] ];
00423       int SubBlockCol = SubBlockColMap_[ indices[j] ];
00424       for( int k = 0; k < BlockCnt_[BlockRow]; ++k )
00425         if( BlockCol == NewBlockRows_[BlockRow][k] )
00426         {
00427           if( values[j] > threshold_ )
00428           {
00429 //          cout << BlockRow << " " << SubBlockRow << " " << BlockCol << " " << SubBlockCol << endl;
00430 //      cout << k << endl;
00431       (*(Blocks_[BlockRow][k]))(SubBlockRow,SubBlockCol) = values[j];
00432             break;
00433           }
00434           else
00435             ZeroElements_[i].erase( OrigMatrix_->RowMap().LID( indices[j] ) );
00436         }
00437     }
00438 
00439 //    NewLHS_->ReplaceGlobalValue( BlockCol, SubBlockCol, 0, (*OrigLHS_)[0][i] );
00440     NewRHS_->ReplaceGlobalValue( BlockRow, SubBlockRow, 0, (*OrigRHS_)[0][i] );
00441   }
00442 
00443   if( verbose_ > 2 )
00444   {
00445     cout << "Zero Elements: \n";
00446     cout << "--------------\n";
00447     int cnt = 0;
00448     for( int i = 0; i < n; ++i )
00449     {
00450       set<int>::iterator iterSI = ZeroElements_[i].begin();
00451       set<int>::iterator endSI = ZeroElements_[i].end();
00452       for( ; iterSI != endSI; ++iterSI )
00453       {
00454         cout << " " << *iterSI;
00455         ++cnt;
00456       }
00457       cout << endl;
00458     }
00459     cout << "ZE Cnt: " << cnt << endl;
00460     cout << "--------------\n";
00461   }
00462 
00463   //setup new matrix
00464   NewMatrix_ = new Epetra_VbrMatrix( View, *NewMap_, &BlockCnt_[0] );
00465   for( int i = 0; i < sqcmpn; ++i )
00466   {
00467     NewMatrix_->BeginInsertGlobalValues( i, BlockCnt_[i], &(NewBlockRows_[i])[0] );
00468     for( int j = 0; j < BlockCnt_[i]; ++j )
00469       NewMatrix_->SubmitBlockEntry( *(Blocks_[i][j]) );
00470     NewMatrix_->EndSubmitEntries();
00471   }
00472   NewMatrix_->FillComplete();
00473 
00474   if( verbose_ > 2 )
00475   {
00476     cout << "New Block Matrix!\n";
00477     cout << *NewMatrix_;
00478     cout << "New Block LHS!\n";
00479     cout << *NewLHS_;
00480     cout << "New Block RHS!\n";
00481     cout << *NewRHS_;
00482   }
00483 
00484   //create new LP
00485   NewProblem_ = new Epetra_LinearProblem( NewMatrix_, NewLHS_, NewRHS_ );
00486   newObj_ = NewProblem_;
00487 
00488   if( verbose_ ) cout << "-----------------------------------------\n";
00489 
00490   return *newObj_;
00491 }
00492 
00493 bool
00494 LinearProblem_BTF::
00495 fwd()
00496 {
00497   //zero out matrix
00498   int NumBlockRows = BlockDim_.size();
00499   for( int i = 0; i < NumBlockRows; ++i )
00500   {
00501     int NumBlocks = BlockCnt_[i];
00502     for( int j = 0; j < NumBlocks; ++j )
00503     {
00504       int Size = BlockDim_[i] * BlockDim_[ NewBlockRows_[i][j] ];
00505       double * p = Blocks_[i][j]->A();
00506       for( int k = 0; k < Size; ++k ) *p++ = 0.0;
00507     }
00508   }
00509 
00510   int maxLength = OrigMatrix_->MaxNumEntries();
00511   int n = OldGlobalElements_.size();
00512   int currLength;
00513   vector<int> indices( maxLength );
00514   vector<double> values( maxLength );
00515   for( int i = 0; i < n; ++i )
00516   {
00517     int BlockRow = BlockRowMap_[ OldGlobalElements_[i] ];
00518     int SubBlockRow = SubBlockRowMap_[ OldGlobalElements_[i] ];
00519     OrigMatrix_->ExtractGlobalRowCopy( OldGlobalElements_[i], maxLength, currLength, &values[0], &indices[0] );
00520     for( int j = 0; j < currLength; ++j )
00521     {
00522       if( fabs(values[j]) > threshold_ )
00523       {
00524         int BlockCol = BlockColMap_[ indices[j] ];
00525         int SubBlockCol = SubBlockColMap_[ indices[j] ];
00526   for( int k = 0; k < BlockCnt_[BlockRow]; ++k )
00527           if( BlockCol == NewBlockRows_[BlockRow][k] )
00528           {
00529       (*(Blocks_[BlockRow][k]))(SubBlockRow,SubBlockCol) = values[j];
00530             break;
00531           }
00532       }
00533     }
00534 
00535     NewRHS_->ReplaceGlobalValue( BlockRow, SubBlockRow, 0, (*OrigRHS_)[0][i] );
00536   }
00537 
00538 /*
00539   //fill matrix
00540   int sqcmpn = BlockDim_.size();
00541   for( int i = 0; i < sqcmpn; ++i )
00542   {
00543     NewMatrix_->BeginReplaceGlobalValues( i, NewBlockRows_[i].size(), &(NewBlockRows_[i])[0] );
00544     for( int j = 0; j < NewBlockRows_[i].size(); ++j )
00545       NewMatrix_->SubmitBlockEntry( Blocks_[i][j]->A(), Blocks_[i][j]->LDA(), Blocks_[i][j]->M(), Blocks_[i][j]->N() );
00546     NewMatrix_->EndSubmitEntries();
00547   }
00548 */
00549 
00550   return true;
00551 }
00552 
00553 bool
00554 LinearProblem_BTF::
00555 rvs()
00556 {
00557   //copy data from NewLHS_ to OldLHS_;
00558   int rowCnt = OrigLHS_->MyLength();
00559   for( int i = 0; i < rowCnt; ++i )
00560   {
00561     int BlockCol = BlockColMap_[ OldGlobalElements_[i] ];
00562     int SubBlockCol = SubBlockColMap_[ OldGlobalElements_[i] ];
00563     (*OrigLHS_)[0][i] = (*NewLHS_)[0][ NewMap_->FirstPointInElement(BlockCol) + SubBlockCol ];
00564   }
00565 
00566   return true;
00567 }
00568 
00569 } //namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines