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 using std::cout;
00058 using std::endl;
00059 
00060 #define MATTRANS_F77 F77_FUNC(mattrans,MATTRANS)
00061 #define GENBTF_F77   F77_FUNC(genbtf,GENBTF)
00062                                                                                                   
00063 extern "C" {
00064 extern void MATTRANS_F77( int*, int*, int*, int*, int*, int* );
00065 extern void GENBTF_F77( int*, int*, int*, int*, int*, int*, int*, int*, int*,
00066                      int*, int*, int*, int*, int*, int*, int*, int*, int*,
00067                      int*, int*, int* );
00068 }
00069 
00070 namespace EpetraExt {
00071 
00072 LinearProblem_BTF::
00073 ~LinearProblem_BTF()
00074 {
00075   deleteNewObjs_();
00076 }
00077 
00078 void
00079 LinearProblem_BTF::
00080 deleteNewObjs_()
00081 {
00082   if( NewProblem_ ) delete NewProblem_;
00083 
00084   if( NewMatrix_ ) delete NewMatrix_;
00085 
00086   if( NewLHS_ ) delete NewLHS_;
00087   if( NewRHS_ ) delete NewRHS_;
00088 
00089   if( NewMap_ ) delete NewMap_;
00090 
00091   for( int i = 0; i < Blocks_.size(); ++i )
00092     for( int j = 0; j < Blocks_[i].size(); ++j )
00093       delete Blocks_[i][j];
00094 }
00095 
00096 LinearProblem_BTF::NewTypeRef
00097 LinearProblem_BTF::
00098 operator()( OriginalTypeRef orig )
00099 {
00100   changedLP_ = false;
00101 
00102   //check if there is a previous analysis and if it is valid
00103   if( &orig == origObj_ && NewProblem_ )
00104   {
00105     int * indices;
00106     double * values;
00107     int currCnt;
00108     int numRows = OrigMatrix_->NumMyRows();
00109 
00110     for( int i = 0; i < numRows; ++i )
00111       if( ZeroElements_[i].size() )
00112       {
00113         int loc = 0;
00114         OrigMatrix_->ExtractMyRowView( i, currCnt, values, indices );
00115         for( int j = 0; j < currCnt; ++j )
00116           if( ZeroElements_[i].count( indices[j] ) )
00117           {
00118             if( values[j] > threshold_ ) changedLP_ = true;
00119             ++loc;
00120           }
00121       }
00122 
00123 //    changedLP_ = true;
00124 
00125     //if changed, dump all the old stuff and start over
00126     if( changedLP_ )
00127       deleteNewObjs_();
00128     else
00129       return *newObj_;
00130   }
00131     
00132   origObj_ = &orig;
00133 
00134   OrigMatrix_ = dynamic_cast<Epetra_CrsMatrix*>(orig.GetMatrix());
00135   OrigLHS_ = orig.GetLHS();
00136   OrigRHS_ = orig.GetRHS();
00137 
00138   if( OrigMatrix_->RowMap().DistributedGlobal() )
00139   { cout << "FAIL for Global!\n"; abort(); }
00140   if( OrigMatrix_->IndicesAreGlobal() )
00141   { cout << "FAIL for Global Indices!\n"; abort(); }
00142 
00143   int n = OrigMatrix_->NumMyRows();
00144   int nnz = OrigMatrix_->NumMyNonzeros();
00145 
00146 //  cout << "Orig Matrix:\n";
00147 //  cout << *OrigMatrix_ << endl;
00148 
00149   //create std CRS format
00150   //also create graph without zero elements
00151   vector<int> ia(n+1,0);
00152   int maxEntries = OrigMatrix_->MaxNumEntries();
00153   vector<int> ja_tmp(maxEntries);
00154   vector<double> jva_tmp(maxEntries);
00155   vector<int> ja(nnz);
00156   int cnt;
00157 
00158   const Epetra_BlockMap & OldRowMap = OrigMatrix_->RowMap();
00159   const Epetra_BlockMap & OldColMap = OrigMatrix_->ColMap();
00160   Epetra_CrsGraph strippedGraph( Copy, OldRowMap, OldColMap, 0 );
00161   ZeroElements_.resize(n);
00162 
00163   for( int i = 0; i < n; ++i )
00164   {
00165     ZeroElements_[i].clear();
00166     OrigMatrix_->ExtractMyRowCopy( i, maxEntries, cnt, &jva_tmp[0], &ja_tmp[0] );
00167     ia[i+1] = ia[i];
00168     for( int j = 0; j < cnt; ++j )
00169     {
00170       if( fabs(jva_tmp[j]) > threshold_ )
00171         ja[ ia[i+1]++ ] = ja_tmp[j];
00172       else
00173         ZeroElements_[i].insert( ja_tmp[j] );
00174     }
00175 
00176     int new_cnt = ia[i+1] - ia[i];
00177     strippedGraph.InsertMyIndices( i, new_cnt, &ja[ ia[i] ] );
00178   }
00179   nnz = ia[n];
00180   strippedGraph.FillComplete();
00181 
00182   if( verbose_ > 2 )
00183   {
00184     cout << "Stripped Graph\n";
00185     cout << strippedGraph;
00186   }
00187 
00188   vector<int> iat(n+1,0);
00189   vector<int> jat(nnz);
00190   for( int i = 0; i < n; ++i )
00191     for( int j = ia[i]; j < ia[i+1]; ++j )
00192       ++iat[ ja[j]+1 ];
00193   for( int i = 0; i < n; ++i )
00194     iat[i+1] += iat[i];
00195   for( int i = 0; i < n; ++i )
00196     for( int j = ia[i]; j < ia[i+1]; ++j )
00197       jat[ iat[ ja[j] ]++ ] = i;
00198   for( int i = 0; i < n; ++i )
00199     iat[n-i] = iat[n-i-1];
00200   iat[0] = 0;
00201 
00202   //convert to Fortran indexing
00203   for( int i = 0; i < n+1; ++i ) ++ia[i];
00204   for( int i = 0; i < nnz; ++i ) ++ja[i];
00205   for( int i = 0; i < n+1; ++i ) ++iat[i];
00206   for( int i = 0; i < nnz; ++i ) ++jat[i];
00207 
00208   if( verbose_ > 2 )
00209   {
00210     cout << "-----------------------------------------\n";
00211     cout << "CRS Format Graph\n";
00212     cout << "-----------------------------------------\n";
00213     for( int i = 0; i < n; ++i )
00214     {
00215       cout << i+1 << ": " << ia[i+1] << ": ";
00216       for( int j = ia[i]-1; j<ia[i+1]-1; ++j )
00217         cout << " " << ja[j];
00218       cout << endl;
00219     }
00220     cout << "-----------------------------------------\n";
00221   }
00222 
00223   if( verbose_ > 2 )
00224   {
00225     cout << "-----------------------------------------\n";
00226     cout << "CCS Format Graph\n";
00227     cout << "-----------------------------------------\n";
00228     for( int i = 0; i < n; ++i )
00229     {
00230       cout << i+1 << ": " << iat[i+1] << ": ";
00231       for( int j = iat[i]-1; j<iat[i+1]-1; ++j )
00232         cout << " " << jat[j];
00233       cout << endl;
00234     }
00235     cout << "-----------------------------------------\n";
00236   }
00237 
00238   vector<int> w(10*n);
00239 
00240   vector<int> rowperm(n);
00241   vector<int> colperm(n);
00242 
00243   //horizontal block
00244   int nhrows, nhcols, hrzcmp;
00245   //square block
00246   int nsrows, sqcmpn;
00247   //vertial block
00248   int nvrows, nvcols, vrtcmp;
00249 
00250   vector<int> rcmstr(n+1);
00251   vector<int> ccmstr(n+1);
00252 
00253   int msglvl = 0;
00254   int output = 6;
00255 
00256   GENBTF_F77( &n, &n, &iat[0], &jat[0], &ia[0], &ja[0], &w[0],
00257           &rowperm[0], &colperm[0], &nhrows, &nhcols,
00258           &hrzcmp, &nsrows, &sqcmpn, &nvrows, &nvcols, &vrtcmp,
00259           &rcmstr[0], &ccmstr[0], &msglvl, &output );
00260 
00261   //convert back to C indexing
00262   for( int i = 0; i < n; ++i )
00263   {
00264     --rowperm[i];
00265     --colperm[i];
00266   }
00267   for( int i = 0; (i<n+1) && (rcmstr[i]!=n+1); ++i )
00268   {
00269     --rcmstr[i];
00270     --ccmstr[i];
00271   }
00272 
00273   if( verbose_ > 0 )
00274   {
00275     cout << "-----------------------------------------\n";
00276     cout << "BTF Output\n";
00277     cout << "-----------------------------------------\n";
00278 //    cout << "RowPerm and ColPerm\n";
00279 //    for( int i = 0; i<n; ++i )
00280 //      cout << rowperm[i] << "\t" << colperm[i] << endl;
00281     if( hrzcmp )
00282     {
00283       cout << "Num Horizontal: Rows, Cols, Comps\n";
00284       cout << nhrows << "\t" << nhcols << "\t" << hrzcmp << endl;
00285     }
00286     cout << "Num Square: Rows, Comps\n";
00287     cout << nsrows << "\t" << sqcmpn << endl;
00288     if( vrtcmp )
00289     {
00290       cout << "Num Vertical: Rows, Cols, Comps\n";
00291       cout << nvrows << "\t" << nvcols << "\t" << vrtcmp << endl;
00292     }
00293 //    cout << "Row, Col of upper left pt in blocks\n";
00294 //    for( int i = 0; (i<n+1) && (rcmstr[i]!=n+1); ++i )
00295 //      cout << i << " " << rcmstr[i] << " " << ccmstr[i] << endl;
00296 //    cout << "-----------------------------------------\n";
00297   }
00298 
00299   if( hrzcmp || vrtcmp )
00300   {
00301     cout << "FAILED! hrz cmp's:" << hrzcmp << " vrtcmp: " << vrtcmp << endl;
00302     exit(0);
00303   }
00304 
00305   rcmstr[sqcmpn] = n;
00306 
00307   //convert rowperm to OLD->NEW
00308   //reverse ordering of permutation to get upper triangular
00309   vector<int> rowperm_t( n );
00310   vector<int> colperm_t( n );
00311   for( int i = 0; i < n; ++i )
00312   {
00313     rowperm_t[i] = rowperm[i];
00314     colperm_t[i] = colperm[i];
00315   }
00316 
00317   //Generate New Domain and Range Maps
00318   //for now, assume they start out as identical
00319   OldGlobalElements_.resize(n);
00320   OldRowMap.MyGlobalElements( &OldGlobalElements_[0] );
00321 
00322   vector<int> newDomainElements( n );
00323   vector<int> newRangeElements( n );
00324   for( int i = 0; i < n; ++i )
00325   {
00326     newDomainElements[ i ] = OldGlobalElements_[ rowperm_t[i] ];
00327     newRangeElements[ i ] = OldGlobalElements_[ colperm_t[i] ];
00328   }
00329 
00330   //Setup New Block Map Info
00331   Blocks_.resize( sqcmpn );
00332   BlockDim_.resize( sqcmpn );
00333   for( int i = 0; i < sqcmpn; ++i )
00334   {
00335     BlockDim_[i] = rcmstr[i+1]-rcmstr[i];
00336     for( int j = rcmstr[i]; j < rcmstr[i+1]; ++j )
00337     {
00338       BlockRowMap_[ newDomainElements[j] ] = i;
00339       SubBlockRowMap_[ newDomainElements[j] ] = j-rcmstr[i];
00340       BlockColMap_[ newRangeElements[j] ] = i;
00341       SubBlockColMap_[ newRangeElements[j] ] = j-rcmstr[i];
00342     }
00343   }
00344 
00345   if( verbose_ > 2 )
00346   {
00347 /*
00348     cout << "Block Mapping!\n";
00349     cout << "--------------------------\n";
00350     for( int i = 0; i < n; ++i )
00351     {
00352       cout << "Row: " << newDomainElements[i] << " " << BlockRowMap_[newDomainElements[i]] << " " <<
00353               SubBlockRowMap_[newDomainElements[i]] << "\t" << "Col: " << newRangeElements[i] << " " <<
00354               BlockColMap_[newRangeElements[i]] << " " << SubBlockColMap_[newRangeElements[i]] << endl;
00355     }
00356     for( int i = 0; i < sqcmpn; ++i )
00357       cout << "BlockDim: " << i << " " << BlockDim_[i] << endl;
00358     cout << "--------------------------\n";
00359 */
00360     int MinSize = 1000000000;
00361     int MaxSize = 0;
00362     for( int i = 0; i < sqcmpn; ++i )
00363     {
00364       if( MinSize > BlockDim_[i] ) MinSize = BlockDim_[i];
00365       if( MaxSize < BlockDim_[i] ) MaxSize = BlockDim_[i];
00366     }
00367     cout << "Min Block Size: " << MinSize << " " << "Max Block Size: " << MaxSize << endl;
00368   }
00369 
00370   vector<int> myBlockElements( sqcmpn );
00371   for( int i = 0; i < sqcmpn; ++i ) myBlockElements[i] = i;
00372   NewMap_ = new Epetra_BlockMap( sqcmpn, sqcmpn, &myBlockElements[0], &BlockDim_[0], 0, OldRowMap.Comm() );
00373 
00374   if( verbose_ > 2 )
00375   {
00376     cout << "New Block Map!\n";
00377     cout << *NewMap_;
00378   }
00379 
00380   //setup new graph
00381   vector< set<int> > crsBlocks( sqcmpn );
00382   BlockCnt_.resize( sqcmpn );
00383   int maxLength = strippedGraph.MaxNumIndices();
00384   vector<int> sIndices( maxLength );
00385   int currLength;
00386   for( int i = 0; i < n; ++i )
00387   {
00388     strippedGraph.ExtractGlobalRowCopy( OldGlobalElements_[i], maxLength, currLength, &sIndices[0] );
00389     for( int j = 0; j < currLength; ++j )
00390       crsBlocks[ BlockRowMap_[ OldGlobalElements_[i] ] ].insert( BlockColMap_[ sIndices[j] ] );
00391   }
00392 
00393   for( int i = 0; i < sqcmpn; ++i )
00394   {
00395     BlockCnt_[i] = crsBlocks[i].size();
00396     Blocks_[i].resize( BlockCnt_[i] );
00397   }
00398 
00399   NewBlockRows_.resize( sqcmpn );
00400   for( int i = 0; i < sqcmpn; ++i )
00401   {
00402     NewBlockRows_[i] = vector<int>( crsBlocks[i].begin(), crsBlocks[i].end() );
00403     for( int j = 0; j < BlockCnt_[i]; ++j )
00404     {
00405       Blocks_[i][j] = new Epetra_SerialDenseMatrix();
00406       Blocks_[i][j]->Shape( BlockDim_[i], BlockDim_[ NewBlockRows_[i][j] ] );
00407     }
00408   }
00409 
00410   //put data in Blocks_ and new LHS and RHS
00411   NewLHS_ = new Epetra_MultiVector( *NewMap_, 1 );
00412   NewRHS_ = new Epetra_MultiVector( *NewMap_, 1 );
00413 
00414   maxLength = OrigMatrix_->MaxNumEntries();
00415   vector<int> indices( maxLength );
00416   vector<double> values( maxLength );
00417   for( int i = 0; i < n; ++i )
00418   {
00419     int BlockRow = BlockRowMap_[ OldGlobalElements_[i] ];
00420     int SubBlockRow = SubBlockRowMap_[ OldGlobalElements_[i] ];
00421     OrigMatrix_->ExtractGlobalRowCopy( OldGlobalElements_[i], maxLength, currLength, &values[0], &indices[0] );
00422     for( int j = 0; j < currLength; ++j )
00423     {
00424       int BlockCol = BlockColMap_[ indices[j] ];
00425       int SubBlockCol = SubBlockColMap_[ indices[j] ];
00426       for( int k = 0; k < BlockCnt_[BlockRow]; ++k )
00427         if( BlockCol == NewBlockRows_[BlockRow][k] )
00428         {
00429           if( values[j] > threshold_ )
00430           {
00431 //          cout << BlockRow << " " << SubBlockRow << " " << BlockCol << " " << SubBlockCol << endl;
00432 //      cout << k << endl;
00433       (*(Blocks_[BlockRow][k]))(SubBlockRow,SubBlockCol) = values[j];
00434             break;
00435           }
00436           else
00437             ZeroElements_[i].erase( OrigMatrix_->RowMap().LID( indices[j] ) );
00438         }
00439     }
00440 
00441 //    NewLHS_->ReplaceGlobalValue( BlockCol, SubBlockCol, 0, (*OrigLHS_)[0][i] );
00442     NewRHS_->ReplaceGlobalValue( BlockRow, SubBlockRow, 0, (*OrigRHS_)[0][i] );
00443   }
00444 
00445   if( verbose_ > 2 )
00446   {
00447     cout << "Zero Elements: \n";
00448     cout << "--------------\n";
00449     int cnt = 0;
00450     for( int i = 0; i < n; ++i )
00451     {
00452       set<int>::iterator iterSI = ZeroElements_[i].begin();
00453       set<int>::iterator endSI = ZeroElements_[i].end();
00454       for( ; iterSI != endSI; ++iterSI )
00455       {
00456         cout << " " << *iterSI;
00457         ++cnt;
00458       }
00459       cout << endl;
00460     }
00461     cout << "ZE Cnt: " << cnt << endl;
00462     cout << "--------------\n";
00463   }
00464 
00465   //setup new matrix
00466   NewMatrix_ = new Epetra_VbrMatrix( View, *NewMap_, &BlockCnt_[0] );
00467   for( int i = 0; i < sqcmpn; ++i )
00468   {
00469     NewMatrix_->BeginInsertGlobalValues( i, BlockCnt_[i], &(NewBlockRows_[i])[0] );
00470     for( int j = 0; j < BlockCnt_[i]; ++j )
00471       NewMatrix_->SubmitBlockEntry( *(Blocks_[i][j]) );
00472     NewMatrix_->EndSubmitEntries();
00473   }
00474   NewMatrix_->FillComplete();
00475 
00476   if( verbose_ > 2 )
00477   {
00478     cout << "New Block Matrix!\n";
00479     cout << *NewMatrix_;
00480     cout << "New Block LHS!\n";
00481     cout << *NewLHS_;
00482     cout << "New Block RHS!\n";
00483     cout << *NewRHS_;
00484   }
00485 
00486   //create new LP
00487   NewProblem_ = new Epetra_LinearProblem( NewMatrix_, NewLHS_, NewRHS_ );
00488   newObj_ = NewProblem_;
00489 
00490   if( verbose_ ) cout << "-----------------------------------------\n";
00491 
00492   return *newObj_;
00493 }
00494 
00495 bool
00496 LinearProblem_BTF::
00497 fwd()
00498 {
00499   //zero out matrix
00500   int NumBlockRows = BlockDim_.size();
00501   for( int i = 0; i < NumBlockRows; ++i )
00502   {
00503     int NumBlocks = BlockCnt_[i];
00504     for( int j = 0; j < NumBlocks; ++j )
00505     {
00506       int Size = BlockDim_[i] * BlockDim_[ NewBlockRows_[i][j] ];
00507       double * p = Blocks_[i][j]->A();
00508       for( int k = 0; k < Size; ++k ) *p++ = 0.0;
00509     }
00510   }
00511 
00512   int maxLength = OrigMatrix_->MaxNumEntries();
00513   int n = OldGlobalElements_.size();
00514   int currLength;
00515   vector<int> indices( maxLength );
00516   vector<double> values( maxLength );
00517   for( int i = 0; i < n; ++i )
00518   {
00519     int BlockRow = BlockRowMap_[ OldGlobalElements_[i] ];
00520     int SubBlockRow = SubBlockRowMap_[ OldGlobalElements_[i] ];
00521     OrigMatrix_->ExtractGlobalRowCopy( OldGlobalElements_[i], maxLength, currLength, &values[0], &indices[0] );
00522     for( int j = 0; j < currLength; ++j )
00523     {
00524       if( fabs(values[j]) > threshold_ )
00525       {
00526         int BlockCol = BlockColMap_[ indices[j] ];
00527         int SubBlockCol = SubBlockColMap_[ indices[j] ];
00528   for( int k = 0; k < BlockCnt_[BlockRow]; ++k )
00529           if( BlockCol == NewBlockRows_[BlockRow][k] )
00530           {
00531       (*(Blocks_[BlockRow][k]))(SubBlockRow,SubBlockCol) = values[j];
00532             break;
00533           }
00534       }
00535     }
00536 
00537     NewRHS_->ReplaceGlobalValue( BlockRow, SubBlockRow, 0, (*OrigRHS_)[0][i] );
00538   }
00539 
00540 /*
00541   //fill matrix
00542   int sqcmpn = BlockDim_.size();
00543   for( int i = 0; i < sqcmpn; ++i )
00544   {
00545     NewMatrix_->BeginReplaceGlobalValues( i, NewBlockRows_[i].size(), &(NewBlockRows_[i])[0] );
00546     for( int j = 0; j < NewBlockRows_[i].size(); ++j )
00547       NewMatrix_->SubmitBlockEntry( Blocks_[i][j]->A(), Blocks_[i][j]->LDA(), Blocks_[i][j]->M(), Blocks_[i][j]->N() );
00548     NewMatrix_->EndSubmitEntries();
00549   }
00550 */
00551 
00552   return true;
00553 }
00554 
00555 bool
00556 LinearProblem_BTF::
00557 rvs()
00558 {
00559   //copy data from NewLHS_ to OldLHS_;
00560   int rowCnt = OrigLHS_->MyLength();
00561   for( int i = 0; i < rowCnt; ++i )
00562   {
00563     int BlockCol = BlockColMap_[ OldGlobalElements_[i] ];
00564     int SubBlockCol = SubBlockColMap_[ OldGlobalElements_[i] ];
00565     (*OrigLHS_)[0][i] = (*NewLHS_)[0][ NewMap_->FirstPointInElement(BlockCol) + SubBlockCol ];
00566   }
00567 
00568   return true;
00569 }
00570 
00571 } //namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines