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