Epetra_VbrMatrix.cpp

Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //               Epetra: 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 "Epetra_VbrMatrix.h"
00030 #include "Epetra_BlockMap.h"
00031 #include "Epetra_Map.h"
00032 #include "Epetra_Import.h"
00033 #include "Epetra_Export.h"
00034 #include "Epetra_Vector.h"
00035 #include "Epetra_MultiVector.h"
00036 #include "Epetra_Comm.h"
00037 #include "Epetra_Distributor.h"
00038 #include "Epetra_SerialDenseMatrix.h"
00039 
00040 //==============================================================================
00041 Epetra_VbrMatrix::Epetra_VbrMatrix(Epetra_DataAccess CV, const Epetra_BlockMap& RowMap, int *NumBlockEntriesPerRow) 
00042   : Epetra_DistObject(RowMap, "Epetra::VbrMatrix"),
00043     Epetra_CompObject(),
00044     Epetra_BLAS(),
00045     Graph_(0),
00046     Allocated_(false),
00047     StaticGraph_(false),
00048     constructedWithFilledGraph_(false),
00049     matrixFillCompleteCalled_(false),
00050     NumMyBlockRows_(RowMap.NumMyElements()),
00051     CV_(CV),
00052     NormInf_(0.0),
00053     NormOne_(0.0),
00054     NormFrob_(0.0),
00055     squareFillCompleteCalled_(false)
00056 {
00057   InitializeDefaults();
00058   Graph_ = new Epetra_CrsGraph(CV, RowMap, NumBlockEntriesPerRow);
00059   int err = Allocate();
00060   assert( err == 0 );
00061 }
00062 
00063 //==============================================================================
00064 Epetra_VbrMatrix::Epetra_VbrMatrix(Epetra_DataAccess CV, const Epetra_BlockMap& RowMap, int NumBlockEntriesPerRow) 
00065   : Epetra_DistObject(RowMap, "Epetra::VbrMatrix"),
00066     Epetra_CompObject(),
00067     Epetra_BLAS(),
00068     Graph_(0),
00069     Allocated_(false),
00070     StaticGraph_(false),
00071     constructedWithFilledGraph_(false),
00072     matrixFillCompleteCalled_(false),
00073     NumMyBlockRows_(RowMap.NumMyElements()),
00074     CV_(CV),
00075     squareFillCompleteCalled_(false)
00076 {
00077   InitializeDefaults();
00078   Graph_ = new Epetra_CrsGraph(CV, RowMap, NumBlockEntriesPerRow);
00079   int err = Allocate();
00080   assert( err == 0 );
00081 }
00082 //==============================================================================
00083 Epetra_VbrMatrix::Epetra_VbrMatrix(Epetra_DataAccess CV, const Epetra_BlockMap& RowMap, 
00084            const Epetra_BlockMap& ColMap, int *NumBlockEntriesPerRow) 
00085   : Epetra_DistObject(RowMap, "Epetra::VbrMatrix"),
00086     Epetra_CompObject(),
00087     Epetra_BLAS(),
00088     Graph_(0),
00089     Allocated_(false),
00090     StaticGraph_(false),
00091     constructedWithFilledGraph_(false),
00092     matrixFillCompleteCalled_(false),
00093     NumMyBlockRows_(RowMap.NumMyElements()),
00094     CV_(CV),
00095     squareFillCompleteCalled_(false)
00096 {
00097   InitializeDefaults();
00098   Graph_ = new Epetra_CrsGraph(CV, RowMap, ColMap, NumBlockEntriesPerRow);
00099   int err = Allocate();
00100   assert( err == 0 );
00101 }
00102 
00103 //==============================================================================
00104 Epetra_VbrMatrix::Epetra_VbrMatrix(Epetra_DataAccess CV, const Epetra_BlockMap& RowMap, 
00105            const Epetra_BlockMap& ColMap, int NumBlockEntriesPerRow) 
00106   : Epetra_DistObject(RowMap, "Epetra::VbrMatrix"),
00107     Epetra_CompObject(),
00108     Epetra_BLAS(),
00109     Graph_(0),
00110     Allocated_(false),
00111     StaticGraph_(false),
00112     constructedWithFilledGraph_(false),
00113     matrixFillCompleteCalled_(false),
00114     NumMyBlockRows_(RowMap.NumMyElements()),
00115     CV_(CV),
00116     squareFillCompleteCalled_(false)
00117 {
00118   InitializeDefaults();
00119   Graph_ = new Epetra_CrsGraph(CV, RowMap, ColMap, NumBlockEntriesPerRow);
00120   int err = Allocate();
00121   assert( err == 0 );
00122 }
00123 //==============================================================================
00124 Epetra_VbrMatrix::Epetra_VbrMatrix(Epetra_DataAccess CV, const Epetra_CrsGraph & Graph) 
00125   : Epetra_DistObject(Graph.RowMap(), "Epetra::VbrMatrix"),
00126     Epetra_CompObject(),
00127     Epetra_BLAS(),
00128     Graph_(new Epetra_CrsGraph(Graph)),
00129     Allocated_(false),
00130     StaticGraph_(true),
00131     constructedWithFilledGraph_(false),
00132     matrixFillCompleteCalled_(false),
00133     NumMyBlockRows_(Graph.RowMap().NumMyElements()),
00134     CV_(CV),
00135     squareFillCompleteCalled_(false)
00136 {
00137   constructedWithFilledGraph_ = Graph.Filled();
00138   InitializeDefaults();
00139   int err = Allocate();
00140   assert(err==0);
00141 }
00142 
00143 //==============================================================================
00144 Epetra_VbrMatrix::Epetra_VbrMatrix(const Epetra_VbrMatrix & Source) 
00145   : Epetra_DistObject(Source),
00146     Epetra_CompObject(),
00147     Epetra_BLAS(),
00148     Graph_(new Epetra_CrsGraph(Source.Graph())),
00149     Allocated_(Source.Allocated_),
00150     StaticGraph_(true),
00151     UseTranspose_(Source.UseTranspose_),
00152     constructedWithFilledGraph_(Source.constructedWithFilledGraph_),
00153     matrixFillCompleteCalled_(Source.matrixFillCompleteCalled_),
00154     NumMyBlockRows_(0),
00155     CV_(Copy),
00156     HavePointObjects_(false),
00157     squareFillCompleteCalled_(false)
00158  {
00159 
00160   InitializeDefaults();
00161   operator=(Source);
00162 }
00163 
00164 //==============================================================================
00165 Epetra_VbrMatrix& Epetra_VbrMatrix::operator=(const Epetra_VbrMatrix& src)
00166 {
00167   if (this == &src) {
00168     return(*this);
00169   }
00170 
00171   DeleteMemory();
00172 
00173   Allocated_ = src.Allocated_;
00174   StaticGraph_ = src.StaticGraph_;
00175   UseTranspose_ = src.UseTranspose_;
00176   NumMyBlockRows_ = src.NumMyBlockRows_;
00177   CV_ = src.CV_;
00178 
00179   InitializeDefaults();
00180 
00181   //our Graph_ member was delete in DeleteMemory() so we don't need to
00182   //delete it now before re-creating it.
00183   Graph_ = new Epetra_CrsGraph(src.Graph());
00184 
00185   int err = Allocate();
00186   assert( err == 0 );
00187 
00188   int i, j;
00189 
00190 #if 0  
00191   bool ConstantShape = true;
00192   const int NOTSETYET = -13 ; 
00193   int MyLda = NOTSETYET ; 
00194   int MyColDim = NOTSETYET ; 
00195   int MyRowDim = NOTSETYET ; 
00196   //  int ierr = Graph_->OptimizeStorage(); // Make sure graph has optimized storage
00197   //  if (ierr) EPETRA_CHK_ERR(ierr);
00198   //  if ( ierr != 0 ) ConstantShape = false ;   // Do we really need this?
00199   assert( ConstantShape ) ; 
00200   for (i=0; i<NumMyBlockRows_; i++) {
00201     int NumBlockEntries = NumBlockEntriesPerRow_[i];
00202     for (j=0; j < NumBlockEntries; j++) {
00203       Epetra_SerialDenseMatrix* ThisBlock = src.Entries_[i][j];
00204       if ( MyLda == NOTSETYET ) {
00205   MyLda = ThisBlock->LDA() ; 
00206         MyColDim = ThisBlock->ColDim() ; 
00207         MyRowDim = ThisBlock->RowDim() ; 
00208       } else {
00209   if ( MyLda !=  ThisBlock->LDA() ) ConstantShape = false ; 
00210   if ( MyRowDim !=  ThisBlock->RowDim() ) ConstantShape = false ; 
00211   if ( MyColDim !=  ThisBlock->ColDim() ) ConstantShape = false ; 
00212       }
00213     }
00214   }  
00215   if ( false &&  ConstantShape ) {
00216 
00217     int NumMyNonzeros = src.Graph_->NumMyNonzeros();
00218 
00219     All_Values_ = new double[NumMyNonzeros];
00220     All_Values_Orig_ = All_Values_ ;
00221     for (i=0; i<NumMyBlockRows_; i++) {
00222       double* Values_ThisBlockRow = All_Values_ ; 
00223       int NumBlockEntries = NumBlockEntriesPerRow_[i];
00224       for (j=0; j < NumBlockEntries; j++) {
00225   double* Values_ThisBlockEntry = All_Values_ ; 
00226   Epetra_SerialDenseMatrix* M_SDM = src.Entries_[i][j] ;
00227   for ( int kk = 0 ; kk < MyColDim ; kk++ ) { 
00228     for ( int ll = 0 ; ll < MyRowDim ; ll++ ) { 
00229       *All_Values_ = (*M_SDM)(ll,kk); 
00230       All_Values_++ ; 
00231     }
00232   }
00233   Entries_[i][j] = new Epetra_SerialDenseMatrix(View, Values_ThisBlockEntry,
00234                   MyLda, MyRowDim, MyColDim );
00235       }
00236     }
00237   } else { 
00238 #endif    
00239     for (i=0; i<NumMyBlockRows_; i++) {
00240       int NumBlockEntries = NumBlockEntriesPerRow_[i];
00241       for (j=0; j < NumBlockEntries; j++) {
00242         //if src.Entries_[i][j] != 0, set Entries_[i][j] to be
00243         //a copy of it.
00244         Entries_[i][j] = src.Entries_[i][j] != 0 ?
00245           new Epetra_SerialDenseMatrix(*(src.Entries_[i][j])) : 0;
00246       }
00247 #if 0
00248     }
00249 #endif
00250   }
00251 
00252   if ( src.StorageOptimized() ) this->OptimizeStorage() ;
00253   
00254   return( *this );
00255 }
00256 
00257 //==============================================================================
00258 void Epetra_VbrMatrix::InitializeDefaults() { // Initialize all attributes that have trivial default values
00259 
00260   UseTranspose_ = false;
00261   Entries_ = 0;
00262   All_Values_ = 0;
00263   All_Values_Orig_ = 0;
00264   NormInf_ = -1.0;
00265   NormOne_ = -1.0;
00266   NormFrob_ = -1.0;
00267   ImportVector_ = 0;
00268 
00269   NumBlockEntriesPerRow_  = 0;
00270   NumAllocatedBlockEntriesPerRow_ = 0;
00271   Indices_ = 0;
00272   ElementSizeList_ = 0;
00273   FirstPointInElementList_ = 0;
00274   
00275   // State variables needed for constructing matrix entry-by-entry
00276 
00277   TempRowDims_ = 0;
00278   TempEntries_ = 0;
00279   LenTemps_ = 0;
00280   CurBlockRow_ = 0;
00281   CurNumBlockEntries_ = 0;
00282   CurBlockIndices_ = 0;
00283   CurEntry_ = -1; // Set to -1 to allow a simple sanity check when submitting entries
00284   CurIndicesAreLocal_ = false;
00285   CurSubmitMode_ = Insert;
00286   
00287   // State variables needed for extracting entries
00288   CurExtractBlockRow_ = 0;
00289   CurExtractEntry_ = -1; // Set to -1 to allow a simple sanity check when extracting entries
00290   CurExtractNumBlockEntries_ = 0;
00291   CurExtractIndicesAreLocal_ = false;
00292   CurExtractView_ = false;
00293   CurRowDim_ = 0;
00294 
00295   // State variable for extracting block diagonal entries
00296   CurBlockDiag_ = -1; // Set to -1 to allow a simple sanity check when extracting entries
00297 
00298   // Atributes that support the Epetra_RowMatrix and Epetra_Operator interfaces
00299   RowMatrixRowMap_ = 0;
00300   RowMatrixColMap_ = 0;
00301   RowMatrixImporter_ = 0;
00302   OperatorDomainMap_ = 0;
00303   OperatorRangeMap_ = 0;
00304   HavePointObjects_ = false;
00305   StorageOptimized_ = false ; 
00306   
00307   OperatorX_ = 0;
00308   OperatorY_ = 0;
00309 
00310   return;
00311 }
00312 
00313 //==============================================================================
00314 int Epetra_VbrMatrix::Allocate() {
00315 
00316   int i, j;
00317   
00318   // Set direct access pointers to graph info (needed for speed)
00319   NumBlockEntriesPerRow_ = Graph_->NumIndicesPerRow();
00320   NumAllocatedBlockEntriesPerRow_ = Graph_->NumAllocatedIndicesPerRow();
00321   Indices_ = Graph_->Indices();
00322 
00323   ElementSizeList_ = RowMap().ElementSizeList();
00324 
00325   FirstPointInElementList_ = RowMap().FirstPointInElementList();
00326   
00327 
00328   // Allocate Entries array
00329   Entries_ = new Epetra_SerialDenseMatrix**[NumMyBlockRows_];
00330   // Allocate and initialize entries
00331   for (i=0; i<NumMyBlockRows_; i++) {
00332     int NumAllocatedBlockEntries = NumAllocatedBlockEntriesPerRow_[i];
00333     
00334     if (NumAllocatedBlockEntries > 0) {
00335       Entries_[i] = new Epetra_SerialDenseMatrix*[NumAllocatedBlockEntries];
00336       for (j=0; j < NumAllocatedBlockEntries; j++) {
00337   Entries_[i][j] = 0;
00338       }
00339     }
00340     else {
00341       Entries_[i] = 0;
00342     }
00343   }
00344   SetAllocated(true);
00345   return(0);
00346 }
00347 //==============================================================================
00348 Epetra_VbrMatrix::~Epetra_VbrMatrix()
00349 {
00350   DeleteMemory();
00351 }
00352 
00353 //==============================================================================
00354 void Epetra_VbrMatrix::DeleteMemory()
00355 {
00356   int i;
00357 
00358   // cout << __FILE__ << " " << __LINE__ << endl ; 
00359 
00360   for (i=0; i<NumMyBlockRows_; i++) {
00361     int NumAllocatedBlockEntries = NumAllocatedBlockEntriesPerRow_[i];
00362     
00363   // cout << __FILE__ << " " << __LINE__ << endl ; 
00364     if (NumAllocatedBlockEntries >0) {
00365   
00366       // cout << __FILE__ << " i=" << i << " linje =" << __LINE__ << endl ; 
00367       for (int j=0; j < NumAllocatedBlockEntries; j++) {
00368   if (Entries_[i][j]!=0) {
00369     delete Entries_[i][j];
00370   }
00371       }
00372       delete [] Entries_[i];
00373     }
00374   }
00375 
00376   // cout << __FILE__ << " " << __LINE__ << endl ; 
00377   if (All_Values_Orig_!=0)   delete [] All_Values_Orig_;
00378   All_Values_Orig_ = 0;
00379 
00380   // cout << __FILE__ << " " << __LINE__ << endl ; 
00381   if (Entries_!=0)       delete [] Entries_;
00382   Entries_ = 0;
00383 
00384   // cout << __FILE__ << " " << __LINE__ << endl ; 
00385   if (ImportVector_!=0) delete ImportVector_;
00386   ImportVector_ = 0;
00387 
00388   // cout << __FILE__ << " " << __LINE__ << endl ; 
00389   NumMyBlockRows_ = 0;
00390 
00391   // cout << __FILE__ << " " << __LINE__ << endl ; 
00392   if (LenTemps_>0) {
00393     delete [] TempRowDims_;
00394     delete [] TempEntries_;
00395   }
00396 
00397   // cout << __FILE__ << " " << __LINE__ << endl ; 
00398   // Delete any objects related to supporting the RowMatrix and Operator interfaces
00399   if (HavePointObjects_) {
00400 #if 1
00401     if ( RowMatrixColMap_ != RowMatrixRowMap_ )        delete RowMatrixColMap_;
00402     if ( OperatorDomainMap_ != RowMatrixRowMap_ )        delete OperatorDomainMap_;
00403     if ( OperatorRangeMap_ != RowMatrixRowMap_ )        delete OperatorRangeMap_;
00404 #else
00405     //  this can fail, see bug # 1114
00406     if (!RowMatrixColMap().SameAs(RowMatrixRowMap())) delete RowMatrixColMap_;
00407     if (!OperatorDomainMap().SameAs(RowMatrixRowMap())) delete OperatorDomainMap_;
00408     if (!OperatorRangeMap().SameAs(RowMatrixRowMap())) delete OperatorRangeMap_;
00409 #endif
00410     delete RowMatrixRowMap_;
00411     delete RowMatrixImporter_;
00412     HavePointObjects_ = false;
00413   }
00414   
00415   // cout << __FILE__ << " " << __LINE__ << endl ; 
00416   if (OperatorX_!=0) {
00417     delete OperatorX_;
00418     delete OperatorY_;
00419   }
00420 
00421   // cout << __FILE__ << " " << __LINE__ << endl ; 
00422   InitializeDefaults(); // Reset all basic pointers to zero
00423   Allocated_ = false;
00424 
00425   // cout << __FILE__ << " " << __LINE__ << endl ; 
00426   delete Graph_; // We created the graph, so must delete it.
00427   Graph_ = 0;
00428   // cout << __FILE__ << " " << __LINE__ << endl ; 
00429 
00430 }
00431 
00432 //==============================================================================
00433 int Epetra_VbrMatrix::PutScalar(double ScalarConstant) 
00434 {
00435   if (!Allocated_) return(0);
00436 
00437   for (int i=0; i<NumMyBlockRows_; i++) {
00438     int NumBlockEntries = NumBlockEntriesPerRow_[i];
00439     int RowDim = ElementSizeList_[i];
00440     for (int j=0; j< NumBlockEntries; j++) {
00441       int LDA = Entries_[i][j]->LDA();
00442       int ColDim = Entries_[i][j]->N();
00443       for (int col=0; col < ColDim; col++) {
00444   double * Entries = Entries_[i][j]->A()+col*LDA;
00445   for (int row=0; row < RowDim; row++)
00446     *Entries++ = ScalarConstant;
00447       }
00448     }
00449   }
00450   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
00451   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
00452   NormFrob_ = -1.0;
00453   return(0);
00454 }
00455 
00456 //==============================================================================
00457 int Epetra_VbrMatrix::Scale(double ScalarConstant) 
00458 {
00459   for (int i=0; i<NumMyBlockRows_; i++) {
00460     int NumBlockEntries = NumBlockEntriesPerRow_[i];
00461     int RowDim = ElementSizeList_[i];
00462     for (int j=0; j< NumBlockEntries; j++) {
00463       int LDA = Entries_[i][j]->LDA();
00464       int ColDim = Entries_[i][j]->N();
00465       for (int col=0; col < ColDim; col++) {
00466   double * Entries = Entries_[i][j]->A()+col*LDA;
00467   for (int row=0; row < RowDim; row++)
00468     *Entries++ *= ScalarConstant;
00469       }
00470     }
00471   }
00472   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
00473   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
00474   NormFrob_ = -1.0;
00475   return(0);
00476 }
00477 //==========================================================================
00478 int Epetra_VbrMatrix::BeginInsertGlobalValues(int BlockRow, int NumBlockEntries, int * BlockIndices) {
00479 
00480   if (IndicesAreLocal()) EPETRA_CHK_ERR(-2); // Cannot insert global values into local graph
00481   Graph_->SetIndicesAreGlobal(true);
00482   int LocalBlockRow = LRID(BlockRow); // Find local row number for this global row index
00483   
00484   bool IndicesAreLocal = false;
00485   EPETRA_CHK_ERR(BeginInsertValues(LocalBlockRow, NumBlockEntries, BlockIndices, IndicesAreLocal));
00486   return(0);
00487 }
00488 
00489 //==========================================================================
00490 int Epetra_VbrMatrix::BeginInsertMyValues(int  BlockRow, int NumBlockEntries,
00491             int * BlockIndices)
00492 {
00493   if (IndicesAreGlobal()) EPETRA_CHK_ERR(-2); // Cannot insert local values until MakeIndicesLocal() is called either directly or through FillComplete()
00494   Graph_->SetIndicesAreLocal(true);
00495   bool IndicesAreLocal = true;
00496   EPETRA_CHK_ERR(BeginInsertValues(BlockRow, NumBlockEntries, BlockIndices, IndicesAreLocal));
00497   return(0);
00498 
00499 }
00500 
00501 //==========================================================================
00502 int Epetra_VbrMatrix::BeginInsertValues(int BlockRow, int NumBlockEntries, 
00503           int * BlockIndices, bool IndicesAreLocal)
00504 {
00505   if (StaticGraph()) EPETRA_CHK_ERR(-2); // If the matrix graph is fully constructed, we cannot insert new values
00506 
00507   int ierr = 0;
00508 
00509   if (BlockRow < 0 || BlockRow >= NumMyBlockRows_) {
00510     EPETRA_CHK_ERR(-1); // Not in BlockRow range    
00511   }
00512   if (CV_==View && Entries_[BlockRow]!=0) ierr = 2; // This row has be defined already. Issue warning.    
00513   if (IndicesAreContiguous()) EPETRA_CHK_ERR(-3); // Indices cannot be individually deleted and new
00514 
00515   // Set up pointers, make sure enough temp space for this round of submits
00516 
00517   Epetra_CombineMode SubmitMode = Insert;
00518 
00519   EPETRA_CHK_ERR(ierr);
00520   EPETRA_CHK_ERR(SetupForSubmits(BlockRow, NumBlockEntries, BlockIndices, IndicesAreLocal, SubmitMode));
00521   return(0);
00522 }
00523 //==========================================================================
00524 int Epetra_VbrMatrix::BeginReplaceGlobalValues(int BlockRow, int NumBlockEntries, int *BlockIndices) {
00525 
00526   BlockRow = LRID(BlockRow); // Normalize row range
00527   bool IndicesAreLocal = false;
00528   EPETRA_CHK_ERR(BeginReplaceValues(BlockRow, NumBlockEntries, BlockIndices, IndicesAreLocal));
00529   return(0);
00530 }
00531 
00532 //==========================================================================
00533 int Epetra_VbrMatrix::BeginReplaceMyValues(int BlockRow, int NumBlockEntries, int *BlockIndices) {
00534 
00535   if (!Graph_->IndicesAreLocal()) EPETRA_CHK_ERR(-1);
00536 
00537   bool IndicesAreLocal = true;
00538   EPETRA_CHK_ERR(BeginReplaceValues(BlockRow, NumBlockEntries, BlockIndices, IndicesAreLocal));
00539   return(0);
00540 }
00541 
00542 //==========================================================================
00543 int Epetra_VbrMatrix::BeginReplaceValues(int BlockRow,
00544            int NumBlockEntries, 
00545            int *BlockIndices,
00546            bool IndicesAreLocal)
00547 {
00548   if (BlockRow < 0 || BlockRow >= NumMyBlockRows_) EPETRA_CHK_ERR(-1); // Not in BlockRow range
00549 
00550   Epetra_CombineMode SubmitMode = Zero; // This is a misuse of Zero mode, fix it later
00551   EPETRA_CHK_ERR(SetupForSubmits(BlockRow, NumBlockEntries, BlockIndices, IndicesAreLocal, SubmitMode));
00552   return(0);
00553 }
00554 
00555 //==========================================================================
00556 int Epetra_VbrMatrix::BeginSumIntoGlobalValues(int BlockRow, int NumBlockEntries, int *BlockIndices) {
00557 
00558   BlockRow = LRID(BlockRow); // Normalize row range
00559   bool IndicesAreLocal = false;
00560   EPETRA_CHK_ERR(BeginSumIntoValues(BlockRow, NumBlockEntries, BlockIndices, IndicesAreLocal));
00561   return(0);
00562 }
00563 
00564 //==========================================================================
00565 int Epetra_VbrMatrix::BeginSumIntoMyValues(int BlockRow, int NumBlockEntries, int *BlockIndices) {
00566 
00567   bool IndicesAreLocal = true;
00568   EPETRA_CHK_ERR(BeginSumIntoValues(BlockRow, NumBlockEntries, BlockIndices, IndicesAreLocal));
00569   return(0);
00570 }
00571 
00572 //==========================================================================
00573 int Epetra_VbrMatrix::BeginSumIntoValues(int BlockRow, int NumBlockEntries, 
00574            int *BlockIndices, bool IndicesAreLocal) {
00575 
00576   if (BlockRow < 0 || BlockRow >= NumMyBlockRows_) EPETRA_CHK_ERR(-1); // Not in BlockRow range
00577 
00578   Epetra_CombineMode SubmitMode = Add;
00579   EPETRA_CHK_ERR(SetupForSubmits(BlockRow, NumBlockEntries, BlockIndices, IndicesAreLocal, SubmitMode));
00580   return(0);
00581 }
00582 
00583 //==========================================================================
00584 int Epetra_VbrMatrix::SetupForSubmits(int BlockRow, int NumBlockEntries,
00585               int * BlockIndices, 
00586               bool IndicesAreLocal,
00587               Epetra_CombineMode SubmitMode) {
00588 
00589   if (NumBlockEntries>LenTemps_) {
00590     if (LenTemps_>0) {
00591       delete [] TempRowDims_;
00592       delete [] TempEntries_;
00593     }
00594     TempRowDims_ = new int[NumBlockEntries];
00595     TempEntries_ = new Epetra_SerialDenseMatrix*[NumBlockEntries];
00596     LenTemps_ = NumBlockEntries;
00597   }
00598 
00599   CurBlockRow_ = BlockRow;
00600   CurNumBlockEntries_ = NumBlockEntries;
00601   CurBlockIndices_ = BlockIndices;
00602   CurEntry_ = 0;
00603   CurIndicesAreLocal_ = IndicesAreLocal;
00604   CurSubmitMode_ = SubmitMode;
00605     
00606   return(0);
00607 }
00608 
00609 //==========================================================================
00610 int Epetra_VbrMatrix::SubmitBlockEntry(double *Values, int LDA,
00611                int NumRows, int NumCols)
00612 {
00613   if (CurEntry_==-1) EPETRA_CHK_ERR(-1); // This means that a Begin routine was not called
00614   if (CurEntry_>=CurNumBlockEntries_) EPETRA_CHK_ERR(-4); // Exceeded the number of entries that can be submitted
00615 
00616   // Fill up temp space with entry
00617 
00618   TempRowDims_[CurEntry_] = NumRows;
00619   TempEntries_[CurEntry_] = new Epetra_SerialDenseMatrix(CV_, Values, LDA,
00620                NumRows, NumCols, false);
00621   CurEntry_++;
00622 
00623   return(0);
00624 }
00625 
00626 //==========================================================================
00627 int Epetra_VbrMatrix::SubmitBlockEntry( Epetra_SerialDenseMatrix &Mat)
00628 {
00629   return SubmitBlockEntry( Mat.A(), Mat.LDA(), Mat.M(), Mat.N() );
00630 }
00631 
00632 //==========================================================================
00633 int Epetra_VbrMatrix::EndSubmitEntries() {
00634 
00635   if (CurEntry_!=CurNumBlockEntries_) EPETRA_CHK_ERR(-6); // Did not submit right number of entries
00636 
00637   if (CurSubmitMode_==Insert) {
00638     EPETRA_CHK_ERR(EndInsertValues());
00639   }
00640   else {
00641     EPETRA_CHK_ERR(EndReplaceSumIntoValues());
00642   }
00643   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
00644   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
00645   NormFrob_ = -1.0;
00646  return(0);
00647 }
00648 
00649 //==========================================================================
00650 int Epetra_VbrMatrix::EndReplaceSumIntoValues()
00651 {
00652   int j;
00653   int ierr = 0;
00654   int Loc;
00655 
00656   int RowDim = ElementSizeList_[CurBlockRow_];
00657 
00658   bool SumInto = (CurSubmitMode_==Add);
00659 
00660   for (j=0; j<CurNumBlockEntries_; j++) {
00661     int BlockIndex = CurBlockIndices_[j];
00662 
00663     bool code = false;
00664     if (CurIndicesAreLocal_) {
00665       code =Graph_->FindMyIndexLoc(CurBlockRow_,BlockIndex,j,Loc);
00666     }
00667     else {
00668       code =Graph_->FindGlobalIndexLoc(CurBlockRow_,BlockIndex,j,Loc);
00669     }
00670 
00671     bool need_to_delete_temp_entry = true;
00672 
00673     if (code) {
00674       Epetra_SerialDenseMatrix* src = TempEntries_[j];
00675 
00676       int ColDim = src->N();
00677 
00678       if (Entries_[CurBlockRow_][Loc]==0) {
00679   Entries_[CurBlockRow_][Loc] = src;
00680         need_to_delete_temp_entry = false;
00681 //    new Epetra_SerialDenseMatrix(RowDim, ColDim, false);
00682       }
00683       else {
00684         Epetra_SerialDenseMatrix* target = Entries_[CurBlockRow_][Loc];
00685 
00686         target->CopyMat(src->A_, src->LDA_, RowDim, ColDim,
00687                         target->A_, target->LDA_, SumInto);
00688       }
00689     }
00690     else {
00691       ierr=2; // Block Discarded, Not Found
00692     }
00693 
00694     if (need_to_delete_temp_entry) {
00695       delete TempEntries_[j];
00696     }
00697   }
00698 
00699   EPETRA_CHK_ERR(ierr);
00700   return(0);
00701 }
00702 
00703 //==========================================================================
00704 int Epetra_VbrMatrix::EndInsertValues()
00705 {
00706   int ierr = 0;
00707   int j;
00708 
00709   int NumValidBlockIndices = CurNumBlockEntries_;
00710   int * ValidBlockIndices = new int[CurNumBlockEntries_];
00711   for ( j=0; j < CurNumBlockEntries_; ++j ) ValidBlockIndices[j] = j;
00712     
00713   if( Graph_->HaveColMap() ) { //test and discard indices not in ColMap
00714     NumValidBlockIndices = 0;
00715     const Epetra_BlockMap& map = Graph_->ColMap();
00716 
00717     for ( j = 0; j < CurNumBlockEntries_; ++j ) {
00718       bool myID = CurIndicesAreLocal_ ?
00719   map.MyLID(CurBlockIndices_[j]) : map.MyGID(CurBlockIndices_[j]);
00720 
00721       if( myID ) {
00722   ValidBlockIndices[ NumValidBlockIndices++ ] = j;
00723       }
00724       else ierr=2; // Discarding a Block not found in ColMap
00725     }
00726   }
00727 
00728   int start = NumBlockEntriesPerRow_[CurBlockRow_];
00729   int stop = start + NumValidBlockIndices;
00730   int NumAllocatedEntries = NumAllocatedBlockEntriesPerRow_[CurBlockRow_];
00731 
00732   if (stop > NumAllocatedEntries) {
00733     if (NumAllocatedEntries==0) { // BlockRow was never allocated, so do it
00734       Entries_[CurBlockRow_] = new Epetra_SerialDenseMatrix*[NumValidBlockIndices];
00735     }
00736     else {
00737       ierr = 1; // Out of room.  Must delete and allocate more space...
00738       Epetra_SerialDenseMatrix ** tmp_Entries =
00739   new Epetra_SerialDenseMatrix*[stop];
00740       for (j=0; j< start; j++) {
00741   tmp_Entries[j] = Entries_[CurBlockRow_][j]; // Copy existing entries
00742       }
00743       delete [] Entries_[CurBlockRow_]; // Delete old storage
00744 
00745       Entries_[CurBlockRow_] = tmp_Entries; // Set pointer to new storage
00746     }
00747   }
00748 
00749   for (j=start; j<stop; j++) {
00750     Epetra_SerialDenseMatrix& mat =
00751       *(TempEntries_[ValidBlockIndices[j-start]]);
00752 
00753     Entries_[CurBlockRow_][j] = new Epetra_SerialDenseMatrix(CV_, mat.A(),
00754                    mat.LDA(),
00755                    mat.M(),
00756                    mat.N());
00757   }
00758 
00759   delete [] ValidBlockIndices;
00760 
00761   for (j=0; j<CurNumBlockEntries_; ++j) {
00762     delete TempEntries_[j];
00763   }
00764 
00765   // Update graph
00766   EPETRA_CHK_ERR(Graph_->InsertIndices(CurBlockRow_, CurNumBlockEntries_, CurBlockIndices_));
00767   EPETRA_CHK_ERR(ierr);
00768 
00769   return(0);
00770 }
00771 
00772 //=============================================================================
00773 int Epetra_VbrMatrix::CopyMat(double * A, int LDA, int NumRows, int NumCols, 
00774             double * B, int LDB, bool SumInto) const
00775 {
00776   int i, j;
00777   double * ptr1 = B;
00778   double * ptr2;
00779 
00780   if (LDB<NumRows) EPETRA_CHK_ERR(-1); // Stride of B is not large enough
00781 
00782   if (SumInto) { // Add to existing values
00783     for (j=0; j<NumCols; j++) {
00784       ptr1 = B + j*LDB;
00785       ptr2 = A + j*LDA;
00786       for (i=0; i<NumRows; i++) *ptr1++ += *ptr2++;
00787     }
00788   }
00789   else {  // Replace values
00790     for (j=0; j<NumCols; j++) {
00791       ptr1 = B + j*LDB;
00792       ptr2 = A + j*LDA;
00793       for (i=0; i<NumRows; i++) *ptr1++ = *ptr2++;
00794     }
00795   }
00796   return(0);
00797 }
00798 //==========================================================================
00799 int Epetra_VbrMatrix::FillComplete() {
00800   squareFillCompleteCalled_ = true;
00801   EPETRA_CHK_ERR(FillComplete(RowMap(), RowMap()));
00802   return(0);
00803 }
00804 
00805 //==========================================================================
00806 int Epetra_VbrMatrix::FillComplete(const Epetra_BlockMap& domain_map,
00807            const Epetra_BlockMap& range_map)
00808 {
00809   int returnValue = 0;
00810 
00811   if (Graph_->Filled()) {
00812     if (!constructedWithFilledGraph_ && !matrixFillCompleteCalled_) {
00813       returnValue = 2;
00814     }
00815   }
00816 
00817  if(!StaticGraph()) {
00818     EPETRA_CHK_ERR(Graph_->MakeIndicesLocal(domain_map, range_map));
00819   }
00820 
00821   SortEntries();  // Sort column entries from smallest to largest
00822   MergeRedundantEntries(); // Get rid of any redundant index values
00823 
00824   if(!StaticGraph()) {
00825     EPETRA_CHK_ERR(Graph_->FillComplete(domain_map, range_map));
00826   }
00827 
00828   // NumMyCols_ = Graph_->NumMyCols(); // Redefine based on local number of cols
00829 
00830   matrixFillCompleteCalled_ = true;
00831 
00832   if (squareFillCompleteCalled_) {
00833     if (DomainMap().NumGlobalElements() != RangeMap().NumGlobalElements()) {
00834       returnValue = 3;
00835     }
00836     squareFillCompleteCalled_ = false;
00837     EPETRA_CHK_ERR(returnValue);
00838   }
00839 
00840   return(returnValue);
00841 }
00842 
00843 //==========================================================================
00844 int Epetra_VbrMatrix::TransformToLocal() {
00845   EPETRA_CHK_ERR(FillComplete());
00846   return(0);
00847 }
00848 
00849 //==========================================================================
00850 int Epetra_VbrMatrix::TransformToLocal(const Epetra_BlockMap* DomainMap, const Epetra_BlockMap* RangeMap) {
00851   EPETRA_CHK_ERR(FillComplete(*DomainMap, *RangeMap));
00852   return(0);
00853 }
00854 
00855 //==========================================================================
00856 int Epetra_VbrMatrix::SortEntries() {
00857 
00858   if (!IndicesAreLocal()) EPETRA_CHK_ERR(-1);
00859   if (Sorted()) return(0);
00860 
00861   // For each row, sort column entries from smallest to largest.
00862   // Use shell sort. Stable sort so it is fast if indices are already sorted.
00863 
00864   
00865   for (int i=0; i<NumMyBlockRows_; i++){
00866 
00867     Epetra_SerialDenseMatrix ** Entries = Entries_[i];
00868     int NumEntries = NumBlockEntriesPerRow_[i];
00869     int * Indices = Indices_[i];
00870     int n = NumEntries;
00871     int m = n/2;
00872     
00873     while (m > 0) {
00874       int max = n - m;
00875       for (int j=0; j<max; j++)
00876         {
00877     for (int k=j; k>=0; k-=m)
00878             {
00879         if (Indices[k+m] >= Indices[k])
00880     break;
00881         Epetra_SerialDenseMatrix *dtemp = Entries[k+m];
00882         Entries[k+m] = Entries[k];
00883         Entries[k] = dtemp;
00884 
00885         int itemp = Indices[k+m];
00886         Indices[k+m] = Indices[k];
00887         Indices[k] = itemp;
00888             }
00889         }
00890       m = m/2;
00891     }
00892   }
00893   Graph_->SetSorted(true); // This also sorted the graph
00894   return(0);
00895 }
00896 
00897 //==========================================================================
00898 int Epetra_VbrMatrix::MergeRedundantEntries()
00899 {
00900 
00901   if (NoRedundancies()) return(0);
00902   if (!Sorted()) EPETRA_CHK_ERR(-1);  // Must have sorted entries
00903 
00904   // For each row, remove column indices that are repeated.
00905   // Also, determine if matrix is upper or lower triangular or has no diagonal
00906   // Note:  This function assumes that SortEntries was already called.
00907 
00908   bool SumInto = true;
00909   for (int i=0; i<NumMyBlockRows_; i++){
00910     int NumEntries = NumBlockEntriesPerRow_[i];
00911     if (NumEntries>1) {
00912       Epetra_SerialDenseMatrix ** const Entries = Entries_[i];
00913       int * const Indices = Indices_[i];
00914       int RowDim = ElementSizeList_[i];
00915       int curEntry =0;
00916       Epetra_SerialDenseMatrix* curBlkEntry = Entries[0];
00917       for (int k=1; k<NumEntries; k++) {
00918   if (Indices[k]==Indices[k-1]) {
00919     CopyMat(Entries[k]->A(), Entries[k]->LDA(), RowDim, Entries[k]->N(),
00920       curBlkEntry->A(), curBlkEntry->LDA(), SumInto);
00921   }
00922   else {
00923     CopyMat(curBlkEntry->A(), curBlkEntry->LDA(), RowDim, curBlkEntry->N(),
00924       Entries[curEntry]->A(), Entries[curEntry]->LDA(), false);
00925     curEntry++;
00926     curBlkEntry = Entries[k];
00927   }
00928       }
00929       CopyMat(curBlkEntry->A(), curBlkEntry->LDA(), RowDim, curBlkEntry->N(),
00930         Entries[curEntry]->A(), Entries[curEntry]->LDA(), false);
00931     }
00932   }
00933     
00934   EPETRA_CHK_ERR(Graph_->RemoveRedundantIndices()); // Remove redundant indices and then return
00935   return(0);
00936 
00937 }
00938 
00939 //==========================================================================
00940 int Epetra_VbrMatrix::OptimizeStorage() {
00941 
00942   if (StorageOptimized()) return(0); // Have we been here before?
00943 
00944   // cout << __FILE__ << " " << __LINE__ << " " << Graph_->NumMyNonzeros() << endl ; 
00945 
00946   bool ConstantShape = true;
00947   int i,j;
00948   const int NOTSETYET = -13 ; 
00949   int MyLda = NOTSETYET ; 
00950   int MyColDim = NOTSETYET ; 
00951   int MyRowDim = NOTSETYET ; 
00952   //
00953   //  I don't know why the underlying graph musthave optimized storage, but 
00954   //  the test for optimized storage depends upon the graph hainvg optimized 
00955   //  storage and that is enough for me.  
00956   //
00957   //  Ideally, we would like to have optimized storage for the graph.  But, 
00958   //  for now, this is commented out until bug 1097 is fixed
00959   //
00960   //  int ierr = Graph_->OptimizeStorage(); // Make sure graph has optimized storage
00961   //  if (ierr) EPETRA_CHK_ERR(ierr);
00962   //  if ( ierr != 0 ) ConstantShape = false ; 
00963   if( ConstantShape ) {
00964     for (i=0; i<NumMyBlockRows_; i++) {
00965       int NumBlockEntries = NumBlockEntriesPerRow_[i];
00966       for (j=0; j < NumBlockEntries; j++) {
00967   Epetra_SerialDenseMatrix* ThisBlock = Entries_[i][j];
00968   if ( MyLda == NOTSETYET ) {
00969     MyLda = ThisBlock->LDA() ; 
00970     MyColDim = ThisBlock->ColDim() ; 
00971     MyRowDim = ThisBlock->RowDim() ; 
00972   } else {
00973     if ( MyLda !=  ThisBlock->LDA() ) ConstantShape = false ; 
00974     if ( MyRowDim !=  ThisBlock->RowDim() ) ConstantShape = false ; 
00975     if ( MyColDim !=  ThisBlock->ColDim() ) ConstantShape = false ; 
00976   }
00977       }
00978     }    
00979   }  
00980   // cout << __FILE__ << " " << __LINE__ << " " << Graph_->NumMyNonzeros() << endl ; 
00981 
00982   if ( ConstantShape ) {
00983 
00984     int NumMyNonzeros = Graph_->NumMyNonzeros();
00985     All_Values_ = new double[NumMyNonzeros];
00986     All_Values_Orig_ = All_Values_ ;
00987     for (i=0; i<NumMyBlockRows_; i++) {
00988       int NumBlockEntries = NumBlockEntriesPerRow_[i];
00989       for (j=0; j < NumBlockEntries; j++) {
00990   double* Values_ThisBlockEntry = All_Values_ ; 
00991   Epetra_SerialDenseMatrix* M_SDM = Entries_[i][j] ;
00992   for ( int kk = 0 ; kk < MyColDim ; kk++ ) { 
00993     for ( int ll = 0 ; ll < MyRowDim ; ll++ ) { 
00994       *All_Values_ = (*M_SDM)(ll,kk) ;
00995       All_Values_++ ; 
00996     }
00997   }
00998   //  Entries_[i][j] = new Epetra_SerialDenseMatrix(*(src.Entries_[i][j]));
00999   delete Entries_[i][j]; 
01000   Entries_[i][j] = new Epetra_SerialDenseMatrix(View, Values_ThisBlockEntry,
01001                   MyLda, MyRowDim, MyColDim );
01002       }
01003     }
01004     StorageOptimized_ = true ; 
01005   }
01006 
01007   // cout << __FILE__ << " " << __LINE__ << " " << Graph_->NumMyNonzeros() << endl ; 
01008 
01009 
01010   /* Work on later...
01011      int i, j;
01012 
01013      // The purpose of this routine is to make the block entries in each row contiguous in memory
01014      // so that a single call to GEMV or GEMM call be used to compute an entire block row.
01015 
01016 
01017      bool Contiguous = true; // Assume contiguous is true
01018      for (i=1; i<NumMyBlockRows_; i++){
01019      int NumEntries = NumBlockEntriesPerRow_[i];
01020      int NumAllocatedEntries = NumAllocatedBlockEntriesPerRow_[i];
01021       
01022      // Check if NumEntries is same as NumAllocatedEntries and 
01023      // check if end of beginning of current row starts immediately after end of previous row.
01024      if ((NumEntries!=NumAllocatedEntries) || (Entries_[i]!=Entries_[i-1]+NumEntries)) {
01025      Contiguous = false;
01026      break;
01027      }
01028      }
01029 
01030      // NOTE:  At the end of the above loop set, there is a possibility that NumEntries and NumAllocatedEntries
01031      //        for the last row could be different, but I don't think it matters.
01032 
01033      
01034      if ((CV_==View) && !Contiguous) EPETRA_CHK_ERR(-1);  // This is user data, it's not contiguous and we can't make it so.
01035 
01036      int ierr = Graph_->OptimizeStorage(); // Make sure graph has optimized storage
01037      if (ierr) EPETRA_CHK_ERR(ierr);
01038 
01039      if (Contiguous) return(0); // Everything is done.  Return
01040 
01041      // Compute Number of Nonzero entries (Done in FillComplete, but we may not have been there yet.)
01042      int NumMyNonzeros = Graph_->NumMyNonzeros();
01043 
01044      // Allocate one big integer array for all index values
01045      All_Values_ = new double[NumMyNonzeros];
01046   
01047      // Set Entries_ to point into All_Entries_
01048   
01049      double * tmp = All_Values_;
01050      for (i=0; i<NumMyBlockRows_; i++) {
01051      int NumEntries = NumEntriesPerBlockRow_[i];
01052      for (j=0; j<NumEntries; j++) tmp[j] = Entries_[i][j];
01053      if (Entries_[i] !=0) delete [] Entries_[i];
01054      Entries_[i] = tmp;
01055      tmp += NumEntries;
01056      }
01057   */
01058   // cout << __FILE__ << " " << __LINE__ << " " << Graph_->NumMyNonzeros() << endl ; 
01059   return(0);
01060 }
01061 //==========================================================================
01062 int Epetra_VbrMatrix::ExtractGlobalRowCopy(int GlobalRow,
01063              int Length, 
01064              int & NumEntries,
01065              double *Values,
01066              int * Indices) const
01067 {
01068   (void)GlobalRow;
01069   (void)Length;
01070   (void)NumEntries;
01071   (void)Values;
01072   (void)Indices;
01073   cout << "Must implement..." << endl;
01074   return(0);
01075 }
01076 //==========================================================================
01077 int Epetra_VbrMatrix::ExtractGlobalBlockRowPointers(int BlockRow, int MaxNumBlockEntries, 
01078                 int & RowDim, int & NumBlockEntries, 
01079                 int * BlockIndices,
01080                 Epetra_SerialDenseMatrix** & Entries) const {
01081 
01082   bool IndicesAreLocal = false;
01083   EPETRA_CHK_ERR(ExtractBlockRowPointers(BlockRow, MaxNumBlockEntries, RowDim, NumBlockEntries, BlockIndices, 
01084            Entries, IndicesAreLocal));
01085   return(0);
01086 }
01087 
01088 //==========================================================================
01089 int Epetra_VbrMatrix::ExtractMyBlockRowPointers(int BlockRow, int MaxNumBlockEntries, 
01090             int & RowDim, int & NumBlockEntries, 
01091             int * BlockIndices,
01092             Epetra_SerialDenseMatrix** & Entries) const {
01093 
01094   bool IndicesAreLocal = true;
01095   EPETRA_CHK_ERR(ExtractBlockRowPointers(BlockRow,MaxNumBlockEntries , RowDim, NumBlockEntries, BlockIndices, 
01096            Entries, IndicesAreLocal));
01097   return(0);
01098 }
01099 
01100 //==========================================================================
01101 int Epetra_VbrMatrix::ExtractBlockRowPointers(int BlockRow, int MaxNumBlockEntries, 
01102                 int & RowDim, int & NumBlockEntries, 
01103                 int * BlockIndices,
01104                 Epetra_SerialDenseMatrix** & Entries, 
01105                 bool IndicesAreLocal) const {
01106   int ierr = 0;
01107   if (!IndicesAreLocal) {
01108     ierr = Graph_->ExtractGlobalRowCopy(BlockRow, MaxNumBlockEntries,
01109           NumBlockEntries, BlockIndices);
01110     BlockRow = LRID(BlockRow);
01111   }
01112   else {
01113     ierr = Graph_->ExtractMyRowCopy(BlockRow, MaxNumBlockEntries,
01114             NumBlockEntries, BlockIndices);
01115   }
01116   if (ierr) EPETRA_CHK_ERR(ierr);
01117 
01118   RowDim = ElementSizeList_[BlockRow];
01119 
01120   Entries = Entries_[BlockRow];
01121 
01122 
01123   EPETRA_CHK_ERR(ierr);
01124   return(0);
01125 }
01126 
01127 //==========================================================================
01128 int Epetra_VbrMatrix::BeginExtractGlobalBlockRowCopy(int BlockRow, int MaxNumBlockEntries, 
01129                  int & RowDim, int & NumBlockEntries, 
01130                  int * BlockIndices, int * ColDims) const {
01131 
01132   bool IndicesAreLocal = false;
01133   EPETRA_CHK_ERR(BeginExtractBlockRowCopy(BlockRow, MaxNumBlockEntries, RowDim, NumBlockEntries, BlockIndices, 
01134             ColDims, IndicesAreLocal));
01135   return(0);
01136 }
01137 
01138 //==========================================================================
01139 int Epetra_VbrMatrix::BeginExtractMyBlockRowCopy(int BlockRow, int MaxNumBlockEntries, 
01140              int & RowDim, int & NumBlockEntries, 
01141              int * BlockIndices, int * ColDims) const {
01142 
01143   bool IndicesAreLocal = true;
01144   EPETRA_CHK_ERR(BeginExtractBlockRowCopy(BlockRow,MaxNumBlockEntries , RowDim, NumBlockEntries, BlockIndices, 
01145             ColDims, IndicesAreLocal));
01146   return(0);
01147 }
01148 
01149 //==========================================================================
01150 int Epetra_VbrMatrix::BeginExtractBlockRowCopy(int BlockRow, int MaxNumBlockEntries, 
01151                  int & RowDim, int & NumBlockEntries, 
01152                  int * BlockIndices, int * ColDims, 
01153                  bool IndicesAreLocal) const  {
01154   int ierr = 0;
01155   if (!IndicesAreLocal)
01156     ierr = Graph_->ExtractGlobalRowCopy(BlockRow, MaxNumBlockEntries, NumBlockEntries, BlockIndices);
01157   else
01158     ierr = Graph_->ExtractMyRowCopy(BlockRow, MaxNumBlockEntries, NumBlockEntries, BlockIndices);
01159   if (ierr) EPETRA_CHK_ERR(ierr);
01160 
01161   bool ExtractView = false;
01162   ierr = SetupForExtracts(BlockRow, RowDim, NumBlockEntries, ExtractView, IndicesAreLocal);
01163   if (ierr) EPETRA_CHK_ERR(ierr);
01164 
01165   EPETRA_CHK_ERR(ExtractBlockDimsCopy(NumBlockEntries, ColDims));
01166   return(0);
01167 }
01168 
01169 //==========================================================================
01170 int Epetra_VbrMatrix::SetupForExtracts(int BlockRow, int & RowDim, int NumBlockEntries, bool ExtractView, 
01171                bool IndicesAreLocal) const {
01172 
01173   if (!IndicesAreLocal) BlockRow = LRID(BlockRow); // Normalize row range
01174   CurExtractBlockRow_ = BlockRow;
01175   CurExtractEntry_ = 0;
01176   CurExtractNumBlockEntries_ = NumBlockEntries;
01177   CurExtractIndicesAreLocal_ = IndicesAreLocal;
01178   CurExtractView_ = ExtractView;
01179   CurRowDim_ = ElementSizeList_[CurExtractBlockRow_];
01180   RowDim = CurRowDim_;
01181     
01182   return(0);
01183 }
01184 
01185 //==========================================================================
01186 int Epetra_VbrMatrix::ExtractBlockDimsCopy(int NumBlockEntries, int * ColDims) const {
01187 
01188   for (int i=0; i<NumBlockEntries; i++) {
01189     ColDims[i] = Entries_[CurExtractBlockRow_][i]->N();
01190   }
01191   return(0);
01192 }
01193 
01194 //==========================================================================
01195 int Epetra_VbrMatrix::ExtractEntryCopy(int SizeOfValues,
01196                double * Values,
01197                int LDA,
01198                bool SumInto) const
01199 {
01200   (void)SumInto;
01201   if (CurExtractEntry_==-1) EPETRA_CHK_ERR(-1); // No BeginCopy routine was called
01202   int CurColDim = Entries_[CurExtractBlockRow_][CurExtractEntry_]->N();
01203   if (LDA*CurColDim>SizeOfValues) EPETRA_CHK_ERR(-2);  // Not enough space
01204 
01205   Epetra_SerialDenseMatrix* CurEntries = Entries_[CurExtractBlockRow_][CurExtractEntry_];
01206   int CurRowDim = CurEntries->M();
01207   int CurLDA = CurEntries->LDA();
01208 
01209   CurExtractEntry_++; // Increment Entry Pointer
01210 
01211   double* vals = CurEntries->A();
01212   if (LDA==CurRowDim && CurLDA==CurRowDim) // Columns of the entry are contiguous, can treat as single array
01213     for (int ii=0; ii<CurRowDim*CurColDim; ++ii)
01214       Values[ii] = vals[ii];
01215   else {
01216     double * CurTargetCol = Values;
01217     double * CurSourceCol = vals;
01218 
01219     for (int jj=0; jj<CurColDim; jj++) {
01220       for (int ii=0; ii<CurRowDim; ++ii) 
01221   CurTargetCol[ii] = CurSourceCol[ii];
01222       CurTargetCol += LDA;
01223       CurSourceCol += CurLDA;
01224     }
01225   }
01226 
01227   return(0);
01228 }
01229 
01230 //==========================================================================
01231 int Epetra_VbrMatrix::BeginExtractGlobalBlockRowView(int BlockRow,
01232                  int & RowDim, int & NumBlockEntries, 
01233                  int * & BlockIndices) const
01234 {
01235 
01236   bool IndicesAreLocal = false;
01237   EPETRA_CHK_ERR(BeginExtractBlockRowView(BlockRow, RowDim, NumBlockEntries,
01238             BlockIndices, IndicesAreLocal));
01239   return(0);
01240 }
01241 
01242 //==========================================================================
01243 int Epetra_VbrMatrix::BeginExtractMyBlockRowView(int BlockRow, int & RowDim, int & NumBlockEntries, 
01244              int * & BlockIndices)  const
01245 {
01246 
01247   bool IndicesAreLocal = true;
01248   EPETRA_CHK_ERR(BeginExtractBlockRowView(BlockRow, RowDim, NumBlockEntries, BlockIndices, 
01249             IndicesAreLocal));
01250   return(0);
01251 }
01252 
01253 //==========================================================================
01254 int Epetra_VbrMatrix::BeginExtractBlockRowView(int BlockRow, int & RowDim, int & NumBlockEntries, 
01255                  int * & BlockIndices,
01256                  bool IndicesAreLocal) const
01257 {
01258   int ierr = 0;
01259   if (!IndicesAreLocal)
01260     ierr = Graph_->ExtractGlobalRowView(BlockRow, NumBlockEntries, BlockIndices);
01261   else
01262     ierr = Graph_->ExtractMyRowView(BlockRow,  NumBlockEntries, BlockIndices);
01263   if (ierr) EPETRA_CHK_ERR(ierr);
01264 
01265   bool ExtractView = true;
01266   ierr = SetupForExtracts(BlockRow, RowDim, NumBlockEntries, ExtractView, IndicesAreLocal);
01267   if (ierr) EPETRA_CHK_ERR(ierr);
01268 
01269   return(0);
01270 }
01271 
01272 //==========================================================================
01273 int Epetra_VbrMatrix::ExtractEntryView(Epetra_SerialDenseMatrix* & entry) const
01274 {
01275   entry = Entries_[CurExtractBlockRow_][CurExtractEntry_];
01276 
01277   CurExtractEntry_++; // Increment Entry Pointer
01278   return(0);
01279 }
01280 
01281 //==========================================================================
01282 int Epetra_VbrMatrix::ExtractGlobalBlockRowView(int BlockRow, int & RowDim, int & NumBlockEntries, 
01283             int * & BlockIndices,
01284             Epetra_SerialDenseMatrix** & Entries) const
01285 {
01286 
01287   Entries = Entries_[LRID(BlockRow)]; // Pointer to Array of pointers for this row's block entries
01288   bool IndicesAreLocal = false;
01289   EPETRA_CHK_ERR(BeginExtractBlockRowView(BlockRow, RowDim, NumBlockEntries,
01290             BlockIndices, 
01291             IndicesAreLocal));
01292   return(0);
01293 }
01294 
01295 //==========================================================================
01296 int Epetra_VbrMatrix::ExtractMyBlockRowView(int BlockRow, int & RowDim, int & NumBlockEntries, 
01297               int * & BlockIndices,
01298               Epetra_SerialDenseMatrix** & Entries) const
01299 {
01300 
01301   Entries = Entries_[BlockRow]; // Pointer to Array of pointers for this row's block entries
01302   bool IndicesAreLocal = true;
01303   EPETRA_CHK_ERR(BeginExtractBlockRowView(BlockRow, RowDim, NumBlockEntries, BlockIndices, 
01304             IndicesAreLocal));
01305   return(0);
01306 }
01307 
01308 //==============================================================================
01309 int Epetra_VbrMatrix::ExtractDiagonalCopy(Epetra_Vector & Diagonal) const {
01310   
01311   if (!Filled()) EPETRA_CHK_ERR(-1); // Can't get diagonal unless matrix is filled
01312   if (!RowMap().SameAs(Diagonal.Map())) EPETRA_CHK_ERR(-2); // Maps must be the same
01313   double * diagptr = Diagonal.Values();
01314   for (int i=0; i<NumMyBlockRows_; i++){
01315     int BlockRow = GRID(i);
01316     int RowDim = ElementSizeList_[i];
01317     int NumEntries = NumBlockEntriesPerRow_[i];
01318     int * Indices = Indices_[i];
01319     for (int j=0; j<NumEntries; j++) {
01320       int BlockCol = GCID(Indices[j]);
01321       if (BlockRow==BlockCol) {
01322   CopyMatDiag(Entries_[i][j]->A(), Entries_[i][j]->LDA(), RowDim,
01323         Entries_[i][j]->N(), diagptr+FirstPointInElementList_[i]);
01324   break;
01325       }
01326     }
01327   }
01328   return(0);
01329 }
01330 //==============================================================================
01331 int Epetra_VbrMatrix::ReplaceDiagonalValues(const Epetra_Vector & Diagonal) {
01332   
01333   if (!Filled()) EPETRA_CHK_ERR(-1); // Can't get diagonal unless matrix is filled
01334   if (!RowMap().SameAs(Diagonal.Map())) EPETRA_CHK_ERR(-2); // Maps must be the same
01335   int ierr = 0;
01336   double * diagptr = (double *) Diagonal.Values(); // Dangerous but being lazy
01337   for (int i=0; i<NumMyBlockRows_; i++){
01338     int BlockRow = GRID(i);
01339     int RowDim = ElementSizeList_[i];
01340     int NumEntries = NumBlockEntriesPerRow_[i];
01341     int * Indices = Indices_[i];
01342     bool DiagMissing = true;
01343     for (int j=0; j<NumEntries; j++) {
01344       int BlockCol = GCID(Indices[j]);
01345       if (BlockRow==BlockCol) {
01346   ReplaceMatDiag(Entries_[i][j]->A(), Entries_[i][j]->LDA(), RowDim, Entries_[i][j]->N(), 
01347            diagptr+FirstPointInElementList_[i]);
01348   DiagMissing = false;
01349   break;
01350       }
01351     }
01352     if (DiagMissing) ierr = 1; // flag a warning error
01353   }
01354   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
01355   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
01356   NormFrob_ = -1.0;
01357   EPETRA_CHK_ERR(ierr);
01358   return(0);
01359 }
01360 //==============================================================================
01361 int Epetra_VbrMatrix::BeginExtractBlockDiagonalCopy(int MaxNumBlockDiagonalEntries, 
01362                 int & NumBlockDiagonalEntries, int * RowColDims ) const{
01363   
01364   if (!Filled()) EPETRA_CHK_ERR(-1); // Can't get diagonal unless matrix is filled
01365   CurBlockDiag_ = 0; // Initialize pointer
01366   NumBlockDiagonalEntries = NumMyBlockRows_;
01367   if (NumBlockDiagonalEntries>MaxNumBlockDiagonalEntries) EPETRA_CHK_ERR(-2);
01368   EPETRA_CHK_ERR(RowMap().ElementSizeList(RowColDims));
01369   return(0);
01370 }
01371 //==============================================================================
01372 int Epetra_VbrMatrix::ExtractBlockDiagonalEntryCopy(int SizeOfValues, double * Values, int LDA, bool SumInto ) const {
01373   
01374   if (CurBlockDiag_==-1) EPETRA_CHK_ERR(-1); // BeginExtractBlockDiagonalCopy was not called
01375   int i = CurBlockDiag_;
01376   int BlockRow = i;
01377   int RowDim = ElementSizeList_[i];
01378   int NumEntries = NumBlockEntriesPerRow_[i];
01379   int * Indices = Indices_[i];
01380   for (int j=0; j<NumEntries; j++) {
01381     int Col = Indices[j];
01382     if (BlockRow==Col) {
01383       int ColDim = Entries_[i][j]->N();
01384       if (LDA*ColDim>SizeOfValues) EPETRA_CHK_ERR(-2); // Not enough room in Values
01385       CopyMat(Entries_[i][j]->A(), Entries_[i][j]->LDA(), RowDim, ColDim, Values, LDA, SumInto);
01386       break;
01387     }
01388   }
01389   CurBlockDiag_++; // Increment counter
01390   return(0);
01391 }
01392 //==============================================================================
01393 int Epetra_VbrMatrix::BeginExtractBlockDiagonalView(int & NumBlockDiagonalEntries, 
01394                 int * & RowColDims ) const {
01395   
01396   if (!Filled()) EPETRA_CHK_ERR(-1); // Can't get diagonal unless matrix is filled
01397   CurBlockDiag_ = 0; // Initialize pointer
01398   NumBlockDiagonalEntries = NumMyBlockRows_;
01399   RowColDims = ElementSizeList_;
01400   return(0);
01401 }
01402 //==============================================================================
01403 int Epetra_VbrMatrix::ExtractBlockDiagonalEntryView(double * & Values, int & LDA) const {
01404   
01405   if (CurBlockDiag_==-1) EPETRA_CHK_ERR(-1); // BeginExtractBlockDiagonalCopy was not called
01406   int i = CurBlockDiag_;
01407   int BlockRow = i;
01408   int NumEntries = NumBlockEntriesPerRow_[i];
01409   int * Indices = Indices_[i];
01410   for (int j=0; j<NumEntries; j++) {
01411     int Col = Indices[j];
01412     if (BlockRow==Col) {
01413       Values = Entries_[i][j]->A();
01414       LDA = Entries_[i][j]->LDA();
01415       break;
01416     }
01417   }
01418   CurBlockDiag_++; // Increment counter
01419   return(0);
01420 }
01421 //=============================================================================
01422 int Epetra_VbrMatrix::CopyMatDiag(double * A, int LDA, int NumRows, int NumCols, 
01423           double * Diagonal) const {
01424 
01425   int i;
01426   double * ptr1 = Diagonal;
01427   double * ptr2;
01428   int ndiags = EPETRA_MIN(NumRows,NumCols);
01429 
01430   for (i=0; i<ndiags; i++) {
01431     ptr2 = A + i*LDA+i;
01432     *ptr1++ = *ptr2;
01433   }
01434   return(0);
01435 }
01436 //=============================================================================
01437 int Epetra_VbrMatrix::ReplaceMatDiag(double * A, int LDA, int NumRows, int NumCols, 
01438              double * Diagonal) {
01439 
01440   int i;
01441   double * ptr1 = Diagonal;
01442   double * ptr2;
01443   int ndiags = EPETRA_MIN(NumRows,NumCols);
01444 
01445   for (i=0; i<ndiags; i++) {
01446     ptr2 = A + i*LDA+i;
01447     *ptr2 = *ptr1++;
01448   }
01449   return(0);
01450 }
01451 //=============================================================================
01452 int Epetra_VbrMatrix::MaxNumEntries() const {
01453 
01454   int outval = 0;
01455 
01456   for (int i=0; i<NumMyBlockRows_; i++){
01457     int NumBlockEntries = NumMyBlockEntries(i);
01458     int NumEntries = 0;
01459     for (int j=0; j<NumBlockEntries; j++) NumEntries += Entries_[i][j]->N();
01460     outval = EPETRA_MAX(outval,NumEntries);
01461   }
01462   return(outval);
01463 }
01464 //=============================================================================
01465 int Epetra_VbrMatrix::NumMyRowEntries(int MyRow, int & NumEntries) const {
01466 
01467   int BlockRow, BlockOffset;
01468   int ierr = RowMap().FindLocalElementID(MyRow, BlockRow, BlockOffset);  if (ierr!=0) EPETRA_CHK_ERR(ierr);
01469 
01470   int NumBlockEntries = NumMyBlockEntries(BlockRow);
01471   NumEntries = 0;
01472   for (int i=0; i<NumBlockEntries; i++) NumEntries += Entries_[BlockRow][i]->N();
01473   return(0);  
01474 }
01475 //=============================================================================
01476 int Epetra_VbrMatrix::ExtractMyRowCopy(int MyRow,
01477                int Length,
01478                int & NumEntries,
01479                double *Values,
01480                int * Indices) const
01481 {
01482   if (!Filled()) EPETRA_CHK_ERR(-1); // Can't extract row unless matrix is filled
01483   if (!IndicesAreLocal()) EPETRA_CHK_ERR(-2);
01484 
01485   int ierr = 0;
01486   int BlockRow, BlockOffset;
01487   ierr = RowMap().FindLocalElementID(MyRow, BlockRow, BlockOffset);
01488   if (ierr!=0) EPETRA_CHK_ERR(ierr);
01489 
01490   int RowDim, NumBlockEntries;
01491   int * BlockIndices;
01492   Epetra_SerialDenseMatrix ** ValBlocks;
01493   ierr = ExtractMyBlockRowView(BlockRow, RowDim, NumBlockEntries,
01494              BlockIndices, ValBlocks);
01495   if (ierr!=0) EPETRA_CHK_ERR(ierr);
01496 
01497   int * ColFirstPointInElementList = FirstPointInElementList_;
01498   if (Importer()!=0) {
01499     ColFirstPointInElementList = ColMap().FirstPointInElementList();
01500   }
01501   NumEntries = 0;
01502   for (int i=0; i<NumBlockEntries; i++) {
01503     int ColDim = ValBlocks[i]->N();
01504     NumEntries += ColDim;
01505     if (NumEntries>Length) EPETRA_CHK_ERR(-3); // Not enough space
01506     int LDA = ValBlocks[i]->LDA();
01507     double * A = ValBlocks[i]->A()+BlockOffset; // Point to first element in row
01508     int Index = ColFirstPointInElementList[BlockIndices[i]];
01509     for (int j=0; j < ColDim; j++) {
01510       *Values++ = *A;
01511       A += LDA;
01512       *Indices++ = Index++;
01513     }
01514   }
01515 
01516   return(0);      
01517 }
01518 //=============================================================================
01519 int Epetra_VbrMatrix::Multiply1(bool TransA, const Epetra_Vector& x, Epetra_Vector& y) const
01520 {
01521   //
01522   // This function forms the product y = A * x or y = A' * x
01523   //
01524 
01525   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
01526   
01527   int i, j;
01528   int * NumBlockEntriesPerRow = NumBlockEntriesPerRow_;
01529   int * FirstPointInElement = FirstPointInElementList_;
01530   int * ElementSize = ElementSizeList_;
01531   int ** Indices = Indices_;
01532   Epetra_SerialDenseMatrix*** Entries = Entries_;
01533 
01534   double * xp = (double*)x.Values();
01535   double *yp = (double*)y.Values();
01536 
01537   int * ColElementSizeList = ColMap().ElementSizeList();
01538   int * ColFirstPointInElementList = ColMap().FirstPointInElementList();
01539 
01540   bool x_and_y_same = false;
01541   Epetra_Vector* ytemp = 0;
01542   if (xp == yp) {
01543     x_and_y_same = true;
01544     ytemp = new Epetra_Vector(y.Map());
01545     yp = (double*)ytemp->Values();
01546   }
01547 
01548   Epetra_MultiVector& yref = x_and_y_same ? *ytemp : y;
01549   //
01550   //Now yref is a reference to ytemp if x_and_y_same == true, otherwise it is a reference
01551   //to y.
01552   //
01553 
01554   if (!TransA) {
01555 
01556     // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
01557     if (Importer()!=0) {
01558       if (ImportVector_!=0) {
01559   if (ImportVector_->NumVectors()!=1) { delete ImportVector_; ImportVector_= 0;}
01560       }
01561       if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),1); // Create import vector if needed
01562       EPETRA_CHK_ERR(ImportVector_->Import(x, *Importer(), Insert));
01563       xp = (double*)ImportVector_->Values();
01564       ColElementSizeList = ColMap().ElementSizeList(); // The Import map will always have an existing ElementSizeList
01565       ColFirstPointInElementList = ColMap().FirstPointInElementList(); // Import map will always have an existing ...
01566     }
01567 
01568 
01569     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
01570     if (Exporter()!=0) {
01571       if (ExportVector_!=0) {
01572   if (ExportVector_->NumVectors()!=1) { delete ExportVector_; ExportVector_= 0;}
01573       }
01574       if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),1); // Create Export vector if needed
01575       yp = (double*)ExportVector_->Values();
01576     }
01577     
01578     // Do actual computation
01579     int NumMyRows_ = NumMyRows();
01580     for (i=0; i< NumMyRows_; i++) yp[i] = 0.0;  // Initialize y
01581     
01582     for (i=0; i < NumMyBlockRows_; i++) {
01583       int      NumEntries = *NumBlockEntriesPerRow++;
01584       int *    BlockRowIndices = *Indices++;
01585       Epetra_SerialDenseMatrix** BlockRowValues  = *Entries++;
01586       double * cury = yp + *FirstPointInElement++;
01587       int      RowDim = *ElementSize++;
01588       for (j=0; j < NumEntries; j++) {
01589   //sum += BlockRowValues[j] * xp[BlockRowIndices[j]];
01590   double * A = BlockRowValues[j]->A();
01591   int LDA = BlockRowValues[j]->LDA();
01592   int Index = BlockRowIndices[j];
01593   double * curx = xp + ColFirstPointInElementList[Index];
01594   int ColDim = ColElementSizeList[Index];
01595   GEMV('N', RowDim, ColDim, 1.0, A, LDA, curx, 1.0, cury);      
01596       }
01597     }
01598     if (Exporter()!=0) {
01599       yref.PutScalar(0.0);
01600       EPETRA_CHK_ERR(yref.Export(*ExportVector_, *Exporter(), Add)); // Fill y with Values from export vector
01601     }
01602     // Handle case of rangemap being a local replicated map
01603     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(yref.Reduce());
01604   }
01605   
01606   else { // Transpose operation
01607     
01608     
01609     // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
01610     
01611     if (Exporter()!=0) {
01612       if (ExportVector_!=0) {
01613   if (ExportVector_->NumVectors()!=1) { delete ExportVector_; ExportVector_= 0;}
01614       }
01615       if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),1); // Create Export vector if needed
01616       EPETRA_CHK_ERR(ExportVector_->Import(x, *Exporter(), Insert));
01617       xp = (double*)ExportVector_->Values();
01618     }
01619 
01620     // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
01621     if (Importer()!=0) {
01622       if (ImportVector_!=0) {
01623   if (ImportVector_->NumVectors()!=1) { delete ImportVector_; ImportVector_= 0;}
01624       }
01625       if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),1); // Create import vector if needed
01626       yp = (double*)ImportVector_->Values();
01627       ColElementSizeList = ColMap().ElementSizeList(); // The Import map will always have an existing ElementSizeList
01628       ColFirstPointInElementList = ColMap().FirstPointInElementList(); // Import map will always have an existing ...
01629     }
01630     
01631     // Do actual computation
01632     int NumMyCols_ = NumMyCols();
01633     for (i=0; i < NumMyCols_; i++) yp[i] = 0.0; // Initialize y for transpose multiply
01634     
01635     for (i=0; i < NumMyBlockRows_; i++) {
01636       int      NumEntries = *NumBlockEntriesPerRow++;
01637       int *    BlockRowIndices = *Indices++;
01638       Epetra_SerialDenseMatrix** BlockRowValues  = *Entries++;
01639       double * curx = xp + *FirstPointInElement++;
01640       int      RowDim = *ElementSize++;
01641       for (j=0; j < NumEntries; j++) {
01642   //yp[BlockRowIndices[j]] += BlockRowValues[j] * xp[i];
01643   double * A = BlockRowValues[j]->A();
01644   int LDA = BlockRowValues[j]->LDA();
01645   int Index = BlockRowIndices[j];
01646   double * cury = yp + ColFirstPointInElementList[Index];
01647   int ColDim = ColElementSizeList[Index];
01648   GEMV('T', RowDim, ColDim, 1.0, A, LDA, curx, 1.0, cury);
01649   
01650       }
01651     }
01652     if (Importer()!=0) {
01653       yref.PutScalar(0.0); // Make sure target is zero
01654       EPETRA_CHK_ERR(yref.Export(*ImportVector_, *Importer(), Add)); // Fill y with Values from export vector
01655     }
01656     // Handle case of rangemap being a local replicated map
01657     if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(yref.Reduce());
01658   }
01659 
01660   if (x_and_y_same == true) {
01661     y = *ytemp;
01662     delete ytemp;
01663   }
01664   
01665   UpdateFlops(2*NumGlobalNonzeros());
01666   return(0);
01667 }
01668 
01669 //=============================================================================
01670 int Epetra_VbrMatrix::Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
01671   //
01672   // This function forms the product Y = A * Y or Y = A' * X
01673   //
01674   
01675   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
01676   
01677   int i;
01678   int * NumBlockEntriesPerRow = NumBlockEntriesPerRow_;
01679   int ** Indices = Indices_;
01680   Epetra_SerialDenseMatrix*** Entries = Entries_;
01681   
01682   int * RowElementSizeList = ElementSizeList_;
01683   int * RowFirstPointInElementList = FirstPointInElementList_;
01684   int * ColElementSizeList = ColMap().ElementSizeList();
01685   int * ColFirstPointInElementList = ColMap().FirstPointInElementList();
01686 
01687    
01688   int NumVectors = X.NumVectors();
01689   double **Xp = (double**)X.Pointers();
01690   double **Yp = (double**)Y.Pointers();
01691 
01692   bool x_and_y_same = false;
01693   Epetra_MultiVector* Ytemp = 0;
01694   if (Xp == Yp) {
01695     x_and_y_same = true;
01696     Ytemp = new Epetra_MultiVector(Y.Map(), NumVectors);
01697     Yp = (double**)Ytemp->Pointers();
01698   }
01699 
01700   Epetra_MultiVector& Yref = x_and_y_same ? *Ytemp : Y;
01701   //
01702   //Now Yref is a reference to Ytemp if x_and_y_same == true, otherwise it is a reference
01703   //to Y.
01704   //
01705 
01706   if (!TransA) {
01707     
01708     // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
01709     if (Importer()!=0) {
01710       if (ImportVector_!=0) {
01711   if (ImportVector_->NumVectors()!=NumVectors) { delete ImportVector_; ImportVector_= 0;}
01712       }
01713       // Create import vector if needed
01714       if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),NumVectors);
01715 
01716       EPETRA_CHK_ERR(ImportVector_->Import(X, *Importer(), Insert));
01717       Xp = (double**)ImportVector_->Pointers();
01718       ColElementSizeList = ColMap().ElementSizeList();
01719       ColFirstPointInElementList = ColMap().FirstPointInElementList();
01720     }
01721 
01722     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
01723     if (Exporter()!=0) {
01724       if (ExportVector_!=0) {
01725   if (ExportVector_->NumVectors()!=NumVectors) { delete ExportVector_; ExportVector_= 0;}
01726       }
01727       // Create Export vector if needed
01728       if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),NumVectors);
01729 
01730       ExportVector_->PutScalar(0.0); // Zero y values
01731       Yp = (double**)ExportVector_->Pointers();
01732       RowElementSizeList = ColMap().ElementSizeList();
01733       RowFirstPointInElementList = ColMap().FirstPointInElementList();
01734     }
01735     else {
01736       Yref.PutScalar(0.0); // Zero y values
01737     }
01738 
01739     // Do actual computation
01740     if ( All_Values_ != 0 ) { // Contiguous Storage 
01741       if ( ! TransA && *RowElementSizeList <= 4 ) { 
01742 
01743   int RowDim = *RowElementSizeList ;
01744 
01745   Epetra_SerialDenseMatrix* Asub = **Entries;
01746   double *A = Asub->A_ ;
01747 
01748   if ( ( NumVectors == 1 ) ) { 
01749 
01750     for (i=0; i < NumMyBlockRows_; i++) {
01751       int      NumEntries = *NumBlockEntriesPerRow++;
01752       int *    BlockRowIndices = *Indices++;
01753 
01754       double * xptr = Xp[0];
01755       double y0 = 0.0;
01756       double y1 = 0.0;
01757       double y2 = 0.0;
01758       double y3 = 0.0;
01759       switch(RowDim) {
01760       case 1:
01761         for (int j=0; j < NumEntries; ++j) {
01762     int xoff = ColFirstPointInElementList[*BlockRowIndices++];
01763     
01764     y0 += A[0]*xptr[xoff];
01765     A += 1; 
01766         }
01767         break;
01768         
01769       case 2:
01770         for (int j=0; j < NumEntries; ++j) {
01771     int xoff = ColFirstPointInElementList[*BlockRowIndices++];
01772     y0 += A[0]*xptr[xoff] + A[2]*xptr[xoff+1];
01773     y1 += A[1]*xptr[xoff] + A[3]*xptr[xoff+1];
01774     A += 4 ;
01775         }
01776         break;
01777         
01778       case 3:
01779         for (int j=0; j < NumEntries; ++j) {
01780     int xoff = ColFirstPointInElementList[*BlockRowIndices++];
01781     y0 += A[0]*xptr[xoff+0] + A[3]*xptr[xoff+1] + A[6]*xptr[xoff+2];
01782     y1 += A[1]*xptr[xoff+0] + A[4]*xptr[xoff+1] + A[7]*xptr[xoff+2];
01783     y2 += A[2]*xptr[xoff+0] + A[5]*xptr[xoff+1] + A[8]*xptr[xoff+2];
01784     A += 9 ;
01785         }
01786         break;
01787 
01788       case 4:
01789         for (int j=0; j < NumEntries; ++j) {
01790     int xoff = ColFirstPointInElementList[*BlockRowIndices++];
01791     y0 += A[0]*xptr[xoff+0] + A[4]*xptr[xoff+1] + A[8]*xptr[xoff+2] +A[12]*xptr[xoff+3];
01792     y1 += A[1]*xptr[xoff+0] + A[5]*xptr[xoff+1] + A[9]*xptr[xoff+2] +A[13]*xptr[xoff+3];
01793     y2 += A[2]*xptr[xoff+0] + A[6]*xptr[xoff+1] + A[10]*xptr[xoff+2] +A[14]*xptr[xoff+3];
01794     y3 += A[3]*xptr[xoff+0] + A[7]*xptr[xoff+1] + A[11]*xptr[xoff+2] +A[15]*xptr[xoff+3];
01795     A += 16 ;
01796         }
01797         break;
01798       default:
01799         assert(false);
01800       }
01801       int  yoff = *RowFirstPointInElementList++;
01802       switch(RowDim) {
01803       case 4:
01804         *(Yp[0]+yoff+3) = y3;
01805       case 3:
01806         *(Yp[0]+yoff+2) = y2;
01807       case 2:
01808         *(Yp[0]+yoff+1) = y1;
01809       case 1:
01810         *(Yp[0]+yoff) = y0;
01811       }
01812     }
01813   } else { // NumVectors != 1 
01814     double *OrigA = A ; 
01815     for (i=0; i < NumMyBlockRows_; i++) {
01816       int      NumEntries = *NumBlockEntriesPerRow++;
01817       int  yoff = *RowFirstPointInElementList++;
01818       int *    BRI = *Indices++;
01819     
01820       for (int k=0; k<NumVectors; k++) {
01821         int *    BlockRowIndices = BRI;
01822         A = OrigA ; 
01823         double * xptr = Xp[k];
01824         double y0 = 0.0;
01825         double y1 = 0.0;
01826         double y2 = 0.0;
01827         double y3 = 0.0;
01828         switch(RowDim) {
01829         case 1:
01830     for (int j=0; j < NumEntries; ++j) {
01831       int xoff = ColFirstPointInElementList[*BlockRowIndices++];
01832       
01833       y0 += A[0]*xptr[xoff];
01834       A += 1; 
01835     }
01836     break;
01837     
01838         case 2:
01839     for (int j=0; j < NumEntries; ++j) {
01840       int xoff = ColFirstPointInElementList[*BlockRowIndices++];
01841       y0 += A[0]*xptr[xoff] + A[2]*xptr[xoff+1];
01842       y1 += A[1]*xptr[xoff] + A[3]*xptr[xoff+1];
01843       A += 4 ;
01844     }
01845     break;
01846     
01847         case 3:
01848     for (int j=0; j < NumEntries; ++j) {
01849       int xoff = ColFirstPointInElementList[*BlockRowIndices++];
01850       y0 += A[0]*xptr[xoff+0] + A[3]*xptr[xoff+1] + A[6]*xptr[xoff+2];
01851       y1 += A[1]*xptr[xoff+0] + A[4]*xptr[xoff+1] + A[7]*xptr[xoff+2];
01852       y2 += A[2]*xptr[xoff+0] + A[5]*xptr[xoff+1] + A[8]*xptr[xoff+2];
01853       A += 9 ;
01854     }
01855     break;
01856       case 4:
01857         for (int j=0; j < NumEntries; ++j) {
01858     int xoff = ColFirstPointInElementList[*BlockRowIndices++];
01859     y0 += A[0]*xptr[xoff+0] + A[4]*xptr[xoff+1] + A[8]*xptr[xoff+2] +A[12]*xptr[xoff+3];
01860     y1 += A[1]*xptr[xoff+0] + A[5]*xptr[xoff+1] + A[9]*xptr[xoff+2] +A[13]*xptr[xoff+3];
01861     y2 += A[2]*xptr[xoff+0] + A[6]*xptr[xoff+1] + A[10]*xptr[xoff+2] +A[14]*xptr[xoff+3];
01862     y3 += A[3]*xptr[xoff+0] + A[7]*xptr[xoff+1] + A[11]*xptr[xoff+2] +A[15]*xptr[xoff+3];
01863     A += 16 ;
01864         }
01865         break;
01866         default:
01867     assert(false);
01868         }
01869         switch(RowDim) {
01870         case 4:
01871     *(Yp[k]+yoff+3) = y3;
01872         case 3:
01873     *(Yp[k]+yoff+2) = y2;
01874         case 2:
01875     *(Yp[k]+yoff+1) = y1;
01876         case 1:
01877     *(Yp[k]+yoff) = y0;
01878         }
01879       }
01880       OrigA = A ; 
01881     }
01882   } // end else NumVectors != 1 
01883 
01884       } else {  
01885   
01886   for (i=0; i < NumMyBlockRows_; i++) {
01887     int      NumEntries = *NumBlockEntriesPerRow++;
01888     int *    BlockRowIndices = *Indices++;
01889     Epetra_SerialDenseMatrix** BlockRowValues  = *Entries++;
01890     int  yoff = *RowFirstPointInElementList++;
01891     int RowDim = *RowElementSizeList++;
01892     
01893     FastBlockRowMultiply(TransA, RowDim, NumEntries, BlockRowIndices, yoff, 
01894              ColFirstPointInElementList, ColElementSizeList, 
01895              BlockRowValues, Xp, Yp, NumVectors);
01896   }
01897       }
01898     } else { 
01899       for (i=0; i < NumMyBlockRows_; i++) {
01900   int      NumEntries = *NumBlockEntriesPerRow++;
01901   int *    BlockRowIndices = *Indices++;
01902   Epetra_SerialDenseMatrix** BlockRowValues  = *Entries++;
01903   int  yoff = *RowFirstPointInElementList++;
01904   int RowDim = *RowElementSizeList++;
01905   BlockRowMultiply(TransA, RowDim, NumEntries, BlockRowIndices, yoff, 
01906        ColFirstPointInElementList, ColElementSizeList, 
01907        BlockRowValues, Xp, Yp, NumVectors);
01908       }
01909     }
01910 
01911     if (Exporter()!=0) {
01912       Yref.PutScalar(0.0);
01913       EPETRA_CHK_ERR(Yref.Export(*ExportVector_, *Exporter(), Add)); // Fill Y with Values from export vector
01914     }
01915     // Handle case of rangemap being a local replicated map
01916     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Yref.Reduce());
01917   }
01918   else { // Transpose operation
01919     
01920     
01921     // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
01922     
01923     if (Exporter()!=0) {
01924       if (ExportVector_!=0) {
01925   if (ExportVector_->NumVectors()!=NumVectors) { delete ExportVector_; ExportVector_= 0;}
01926       }
01927       // Create Export vector if needed
01928       if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),NumVectors);
01929 
01930       EPETRA_CHK_ERR(ExportVector_->Import(X, *Exporter(), Insert));
01931       Xp = (double**)ExportVector_->Pointers();
01932       ColElementSizeList = RowMap().ElementSizeList();
01933       ColFirstPointInElementList = RowMap().FirstPointInElementList();
01934     }
01935   
01936     // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
01937     if (Importer()!=0) {
01938       if (ImportVector_!=0) {
01939   if (ImportVector_->NumVectors()!=NumVectors) { delete ImportVector_; ImportVector_= 0;}
01940       }
01941       // Create import vector if needed
01942       if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),NumVectors);
01943 
01944       ImportVector_->PutScalar(0.0); // Zero y values
01945       Yp = (double**)ImportVector_->Pointers();
01946       RowElementSizeList = ColMap().ElementSizeList();
01947       RowFirstPointInElementList = ColMap().FirstPointInElementList();
01948     }
01949     else
01950       Yref.PutScalar(0.0); // Zero y values
01951 
01952     // Do actual computation
01953 
01954     if ( All_Values_ != 0 ) { // Contiguous Storage 
01955       for (i=0; i < NumMyBlockRows_; i++) {
01956   int      NumEntries = *NumBlockEntriesPerRow++;
01957   int *    BlockRowIndices = *Indices++;
01958   Epetra_SerialDenseMatrix** BlockRowValues  = *Entries++;
01959   int  xoff = *RowFirstPointInElementList++;
01960   int RowDim = *RowElementSizeList++;
01961   FastBlockRowMultiply(TransA, RowDim, NumEntries, BlockRowIndices, xoff, 
01962        ColFirstPointInElementList, ColElementSizeList, 
01963        BlockRowValues, Xp, Yp, NumVectors);
01964       }
01965     } else {
01966       for (i=0; i < NumMyBlockRows_; i++) {
01967   int      NumEntries = *NumBlockEntriesPerRow++;
01968   int *    BlockRowIndices = *Indices++;
01969   Epetra_SerialDenseMatrix** BlockRowValues  = *Entries++;
01970   int  xoff = *RowFirstPointInElementList++;
01971   int RowDim = *RowElementSizeList++;
01972   BlockRowMultiply(TransA, RowDim, NumEntries, BlockRowIndices, xoff, 
01973        ColFirstPointInElementList, ColElementSizeList, 
01974        BlockRowValues, Xp, Yp, NumVectors);
01975       }
01976     }
01977 
01978     if (Importer()!=0) {
01979       Yref.PutScalar(0.0); // Make sure target is zero
01980       EPETRA_CHK_ERR(Yref.Export(*ImportVector_, *Importer(), Add)); // Fill Y with Values from export vector
01981     }
01982     // Handle case of rangemap being a local replicated map
01983     if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1)  EPETRA_CHK_ERR(Yref.Reduce());
01984   }
01985 
01986   UpdateFlops(2*NumVectors*NumGlobalNonzeros());
01987 
01988   if (x_and_y_same == true) {
01989     Y = *Ytemp;
01990     delete Ytemp;
01991   }
01992 
01993   return(0);
01994 }
01995 //=============================================================================
01996 void Epetra_VbrMatrix::BlockRowMultiply(bool TransA,
01997           int RowDim,
01998           int NumEntries, 
01999           int * BlockIndices,
02000           int RowOff,
02001           int * FirstPointInElementList,
02002           int * ElementSizeList,
02003           double Alpha,
02004           Epetra_SerialDenseMatrix** As,
02005           double ** X,
02006           double Beta,
02007           double ** Y,
02008           int NumVectors) const
02009 {
02010   //This overloading of BlockRowMultiply is the same as the one below, except
02011   //that this one accepts Alpha and Beta arguments. This BlockRowMultiply is
02012   //called from within the 'solve' methods.
02013   int j, k;
02014   if (!TransA) {
02015     for (j=0; j < NumEntries; j++) {
02016       Epetra_SerialDenseMatrix* Asub = As[j];
02017       double * A = Asub->A();
02018       int LDA = Asub->LDA();
02019       int BlockIndex = BlockIndices[j];
02020       int xoff = FirstPointInElementList[BlockIndex];
02021       int ColDim = ElementSizeList[BlockIndex];
02022 
02023       for (k=0; k<NumVectors; k++) {
02024   double * curx = X[k] + xoff;
02025   double * cury = Y[k] + RowOff;
02026 
02027   GEMV('N', RowDim, ColDim, Alpha, A, LDA, curx, Beta, cury);
02028       }//end for (k
02029     }//end for (j
02030   }
02031   else { //TransA == true
02032     for (j=0; j < NumEntries; j++) {
02033       double * A = As[j]->A();
02034       int LDA = As[j]->LDA();
02035       int BlockIndex = BlockIndices[j];
02036       int yoff = FirstPointInElementList[BlockIndex];
02037       int ColDim = ElementSizeList[BlockIndex];
02038       for (k=0; k<NumVectors; k++) {
02039   double * curx = X[k] + RowOff;
02040   double * cury = Y[k] + yoff;
02041   GEMV('T', RowDim, ColDim, Alpha, A, LDA, curx, Beta, cury);
02042       }
02043     }
02044   }
02045 
02046   return;
02047 }
02048 //=============================================================================
02049 void Epetra_VbrMatrix::FastBlockRowMultiply(bool TransA,
02050   int RowDim,
02051   int NumEntries, 
02052   int * BlockRowIndices,
02053   int yoff,
02054   int * ColFirstPointInElementList,
02055   int * ColElementSizeList,
02056   Epetra_SerialDenseMatrix** BlockRowValues,
02057   double ** Xp,
02058   double ** Yp,
02059   int NumVectors) const
02060 {
02061   int j, k;
02062   if (!TransA) {
02063     for (k=0; k<NumVectors; k++) {
02064       double * y = Yp[k] + yoff;
02065       double * xptr = Xp[k];
02066 
02067       Epetra_SerialDenseMatrix* Asub = BlockRowValues[0];
02068       double *OrigA = Asub->A_;
02069       int LDA = Asub->LDA_;
02070       int OrigBlockIndex = BlockRowIndices[0];
02071       int ColDim = ColElementSizeList[OrigBlockIndex];
02072  
02073       assert( RowDim == ColDim ) ; 
02074       assert( RowDim == LDA ) ; 
02075       
02076       switch(RowDim) {
02077 #if 0
02078       case 1:
02079   double y0 = y[0];
02080   double y1 = y[0];
02081   double *A = OrigA ;
02082   for (j=0; j < NumEntries; ++j) {
02083     
02084     int BlockIndex = BlockRowIndices[j];
02085     int xoff = ColFirstPointInElementList[BlockIndex];
02086     
02087     double *A = OrigA + j * ColDim * LDA  ;
02088     double *x = xptr + xoff;
02089     y[0] += A[0]*x[0];
02090   }//end for (j
02091   break;
02092   
02093       case 2:
02094   double y0 = y[0];
02095   double y1 = y[0];
02096   double *A = OrigA ;
02097   for (j=0; j < NumEntries; ++j) {
02098     
02099     int BlockIndex = BlockRowIndices[j];
02100     int xoff = ColFirstPointInElementList[BlockIndex];
02101     
02102     y0 += A[0]*xptr[xoff] + A[2]*xptr[xoff+1];
02103     y1 += A[1]*xptr[xoff] + A[3]*xptr[xoff+1];
02104     A += ColDim * LDA  ;
02105   }
02106   y[0] = y0;
02107   y[1] = y1;
02108   break;
02109   
02110       case 3:
02111   for (j=0; j < NumEntries; ++j) {
02112       
02113     int BlockIndex = BlockRowIndices[j];
02114     int xoff = ColFirstPointInElementList[BlockIndex];
02115       
02116     double *A = OrigA + j * ColDim * LDA  ;
02117     double *x = xptr + xoff;
02118     y[0] += A[0]*x[0] + A[3]*x[1] + A[6]*x[2];
02119     y[1] += A[1]*x[0] + A[4]*x[1] + A[7]*x[2];
02120     y[2] += A[2]*x[0] + A[5]*x[1] + A[8]*x[2];
02121   }
02122   break;
02123 
02124       case 4:
02125   for (j=0; j < NumEntries; ++j) {
02126       
02127     int BlockIndex = BlockRowIndices[j];
02128     int xoff = ColFirstPointInElementList[BlockIndex];
02129       
02130     double *A = OrigA + j * ColDim * LDA  ;
02131     double *x = xptr + xoff;
02132     y[0] += A[0]*x[0] + A[4]*x[1] + A[8]*x[2] + A[12]*x[3];
02133     y[1] += A[1]*x[0] + A[5]*x[1] + A[9]*x[2] + A[13]*x[3];
02134     y[2] += A[2]*x[0] + A[6]*x[1] + A[10]*x[2] + A[14]*x[3];
02135     y[3] += A[3]*x[0] + A[7]*x[1] + A[11]*x[2] + A[15]*x[3];
02136   }
02137   break;
02138 #endif
02139       case 5:
02140   for (j=0; j < NumEntries; ++j) {
02141       
02142     int BlockIndex = BlockRowIndices[j];
02143     int xoff = ColFirstPointInElementList[BlockIndex];
02144       
02145     double *A = OrigA + j * ColDim * LDA  ;
02146     double *x = xptr + xoff;
02147     y[0] += A[0]*x[0] + A[5]*x[1] + A[10]*x[2] + A[15]*x[3] + A[20]*x[4];
02148     y[1] += A[1]*x[0] + A[6]*x[1] + A[11]*x[2] + A[16]*x[3] + A[21]*x[4];
02149     y[2] += A[2]*x[0] + A[7]*x[1] + A[12]*x[2] + A[17]*x[3] + A[22]*x[4];
02150     y[3] += A[3]*x[0] + A[8]*x[1] + A[13]*x[2] + A[18]*x[3] + A[23]*x[4];
02151     y[4] += A[4]*x[0] + A[9]*x[1] + A[14]*x[2] + A[19]*x[3] + A[24]*x[4];
02152   }
02153   break;
02154 
02155       case 6:
02156   for (j=0; j < NumEntries; ++j) {
02157       
02158     int BlockIndex = BlockRowIndices[j];
02159     int xoff = ColFirstPointInElementList[BlockIndex];
02160       
02161     double *A = OrigA + j * ColDim * LDA  ;
02162     double *x = xptr + xoff;
02163     y[0] += A[0]*x[0] + A[6]*x[1] + A[12]*x[2] + A[18]*x[3] + A[24]*x[4]
02164       + A[30]*x[5];
02165     y[1] += A[1]*x[0] + A[7]*x[1] + A[13]*x[2] + A[19]*x[3] + A[25]*x[4]
02166       + A[31]*x[5];
02167     y[2] += A[2]*x[0] + A[8]*x[1] + A[14]*x[2] + A[20]*x[3] + A[26]*x[4]
02168       + A[32]*x[5];
02169     y[3] += A[3]*x[0] + A[9]*x[1] + A[15]*x[2] + A[21]*x[3] + A[27]*x[4]
02170       + A[33]*x[5];
02171     y[4] += A[4]*x[0] + A[10]*x[1] + A[16]*x[2] + A[22]*x[3] + A[28]*x[4]
02172       + A[34]*x[5];
02173     y[5] += A[5]*x[0] + A[11]*x[1] + A[17]*x[2] + A[23]*x[3] + A[29]*x[4]
02174       + A[35]*x[5];
02175   }
02176   break;
02177 
02178       default:
02179   for (j=0; j < NumEntries; ++j) {
02180       
02181     int BlockIndex = BlockRowIndices[j];
02182     int xoff = ColFirstPointInElementList[BlockIndex];
02183       
02184     double *A = OrigA + j * ColDim * LDA  ;
02185     double *x = xptr + xoff;
02186     GEMV('N', RowDim, ColDim, 1.0, A, LDA, x, 1.0, y);
02187   }
02188       }//end switch
02189 
02190     }//end for (k
02191   }
02192   else { //TransA == true
02193     for (j=0; j < NumEntries; j++) {
02194       double * A = BlockRowValues[j]->A();
02195       int LDA = BlockRowValues[j]->LDA();
02196       int BlockIndex = BlockRowIndices[j];
02197       int Yyoff = ColFirstPointInElementList[BlockIndex];
02198       int ColDim = ColElementSizeList[BlockIndex];
02199       for (k=0; k<NumVectors; k++) {
02200   double * x = Xp[k] + yoff;
02201   double * y = Yp[k] + Yyoff;
02202   GEMV('T', RowDim, ColDim, 1.0, A, LDA, x, 1.0, y);
02203       }
02204     }
02205   }
02206 
02207 }
02208 void Epetra_VbrMatrix::BlockRowMultiply(bool TransA,
02209           int RowDim,
02210           int NumEntries, 
02211           int * BlockIndices,
02212           int RowOff,
02213           int * FirstPointInElementList,
02214           int * ElementSizeList,
02215           Epetra_SerialDenseMatrix** As,
02216           double ** X,
02217           double ** Y,
02218           int NumVectors) const
02219 {
02220   //This overloading of BlockRowMultiply is the same as the one above, except
02221   //that this one doesn't accept the Alpha and Beta arguments (it assumes that
02222   //they are both 1.0) and contains some inlined unrolled code for certain
02223   //cases (certain block-sizes) rather than calling GEMV every time. This
02224   //BlockRowMultiply is called from within the 'Multiply' methods.
02225   //Note: Scott Hutchinson's Aztec method 'dvbr_sparax_basic' was consulted in
02226   //the optimizing of this method.
02227 
02228   int j, k;
02229   if (!TransA) {
02230     for (k=0; k<NumVectors; k++) {
02231       double * y = Y[k] + RowOff;
02232       double * xptr = X[k];
02233 
02234       for (j=0; j < NumEntries; ++j) {
02235   Epetra_SerialDenseMatrix* Asub = As[j];
02236   double * A = Asub->A_;
02237   int LDA = Asub->LDA_;
02238   int BlockIndex = BlockIndices[j];
02239   int xoff = FirstPointInElementList[BlockIndex];
02240   int ColDim = ElementSizeList[BlockIndex];
02241 
02242   double * x = xptr + xoff;
02243 
02244   //Call GEMV if sub-block is non-square or if LDA != RowDim.
02245   if (LDA != RowDim || ColDim != RowDim) {
02246     GEMV('N', RowDim, ColDim, 1.0, A, LDA, x, 1.0, y);
02247     continue;
02248   }
02249 
02250   //It is a big performance win to use inlined, unrolled code for small
02251   //common block sizes rather than calling GEMV.
02252 
02253   switch(RowDim) {
02254   case 1:
02255     y[0] += A[0]*x[0];
02256     break;
02257 
02258   case 2:
02259     y[0] += A[0]*x[0] + A[2]*x[1];
02260     y[1] += A[1]*x[0] + A[3]*x[1];
02261     break;
02262 
02263   case 3:
02264     y[0] += A[0]*x[0] + A[3]*x[1] + A[6]*x[2];
02265     y[1] += A[1]*x[0] + A[4]*x[1] + A[7]*x[2];
02266     y[2] += A[2]*x[0] + A[5]*x[1] + A[8]*x[2];
02267     break;
02268 
02269   case 4:
02270     y[0] += A[0]*x[0] + A[4]*x[1] + A[8]*x[2] + A[12]*x[3];
02271     y[1] += A[1]*x[0] + A[5]*x[1] + A[9]*x[2] + A[13]*x[3];
02272     y[2] += A[2]*x[0] + A[6]*x[1] + A[10]*x[2] + A[14]*x[3];
02273     y[3] += A[3]*x[0] + A[7]*x[1] + A[11]*x[2] + A[15]*x[3];
02274     break;
02275 
02276   case 5:
02277     y[0] += A[0]*x[0] + A[5]*x[1] + A[10]*x[2] + A[15]*x[3] + A[20]*x[4];
02278     y[1] += A[1]*x[0] + A[6]*x[1] + A[11]*x[2] + A[16]*x[3] + A[21]*x[4];
02279     y[2] += A[2]*x[0] + A[7]*x[1] + A[12]*x[2] + A[17]*x[3] + A[22]*x[4];
02280     y[3] += A[3]*x[0] + A[8]*x[1] + A[13]*x[2] + A[18]*x[3] + A[23]*x[4];
02281     y[4] += A[4]*x[0] + A[9]*x[1] + A[14]*x[2] + A[19]*x[3] + A[24]*x[4];
02282     break;
02283 
02284   case 6:
02285     y[0] += A[0]*x[0] + A[6]*x[1] + A[12]*x[2] + A[18]*x[3] + A[24]*x[4]
02286       + A[30]*x[5];
02287     y[1] += A[1]*x[0] + A[7]*x[1] + A[13]*x[2] + A[19]*x[3] + A[25]*x[4]
02288       + A[31]*x[5];
02289     y[2] += A[2]*x[0] + A[8]*x[1] + A[14]*x[2] + A[20]*x[3] + A[26]*x[4]
02290       + A[32]*x[5];
02291     y[3] += A[3]*x[0] + A[9]*x[1] + A[15]*x[2] + A[21]*x[3] + A[27]*x[4]
02292       + A[33]*x[5];
02293     y[4] += A[4]*x[0] + A[10]*x[1] + A[16]*x[2] + A[22]*x[3] + A[28]*x[4]
02294       + A[34]*x[5];
02295     y[5] += A[5]*x[0] + A[11]*x[1] + A[17]*x[2] + A[23]*x[3] + A[29]*x[4]
02296       + A[35]*x[5];
02297     break;
02298 
02299   default:
02300     GEMV('N', RowDim, ColDim, 1.0, A, LDA, x, 1.0, y);
02301   }//end switch
02302       }//end for (k
02303     }//end for (j
02304   }
02305   else { //TransA == true
02306     for (j=0; j < NumEntries; j++) {
02307       double * A = As[j]->A();
02308       int LDA = As[j]->LDA();
02309       int BlockIndex = BlockIndices[j];
02310       int yoff = FirstPointInElementList[BlockIndex];
02311       int ColDim = ElementSizeList[BlockIndex];
02312       for (k=0; k<NumVectors; k++) {
02313   double * x = X[k] + RowOff;
02314   double * y = Y[k] + yoff;
02315   GEMV('T', RowDim, ColDim, 1.0, A, LDA, x, 1.0, y);
02316       }
02317     }
02318   }
02319 
02320   return;
02321 }
02322 //=============================================================================
02323 int Epetra_VbrMatrix::Solve(bool Upper, bool TransA,
02324           bool UnitDiagonal,
02325           const Epetra_MultiVector& X,
02326           Epetra_MultiVector& Y) const
02327 {
02328   (void)UnitDiagonal;
02329   //
02330   // This function find Y such that LY = X or UY = X or the transpose cases.
02331   //
02332 
02333   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
02334 
02335   if ((Upper) && (!UpperTriangular())) EPETRA_CHK_ERR(-2);
02336   if ((!Upper) && (!LowerTriangular())) EPETRA_CHK_ERR(-3);
02337   if (!NoDiagonal()) EPETRA_CHK_ERR(-4); // We must use UnitDiagonal
02338 
02339   int i;
02340   int * NumBlockEntriesPerRow = NumBlockEntriesPerRow_;
02341   int * FirstPointInElement = FirstPointInElementList_;
02342   int * ElementSize = ElementSizeList_;
02343   int ** Indices = Indices_;
02344   Epetra_SerialDenseMatrix*** Entries = Entries_;
02345 
02346   int * ColElementSizeList = ElementSizeList_;
02347   int * ColFirstPointInElementList = FirstPointInElementList_;
02348 
02349   // If upper, point to last row
02350   if (Upper) {
02351     NumBlockEntriesPerRow += NumMyBlockRows_-1;
02352     FirstPointInElement += NumMyBlockRows_-1;
02353     ElementSize += NumMyBlockRows_-1;
02354     Indices += NumMyBlockRows_-1;
02355     Entries += NumMyBlockRows_-1;
02356   }
02357 
02358   double **Yp = (double**)Y.Pointers();
02359 
02360   int NumVectors = X.NumVectors();
02361 
02362   if (X.Pointers() != Y.Pointers()) Y = X; // Copy X into Y if they are not the same multivector
02363 
02364   bool Case1 = (((!TransA) && Upper) ||  (TransA && !Upper)); 
02365   // Case 2 = (((TransA) && Upper) || (!TransA) && Lower);
02366   if (Case1) {
02367     for (i=0; i < NumMyBlockRows_; i++) {
02368       int      NumEntries = *NumBlockEntriesPerRow--;
02369       int *    BlockRowIndices = *Indices--;
02370       Epetra_SerialDenseMatrix** BlockRowValues  = *Entries--;
02371       int  yoff = *FirstPointInElement--;
02372       int RowDim = *ElementSize--;
02373       BlockRowMultiply(TransA, RowDim, NumEntries, BlockRowIndices, yoff, 
02374            ColFirstPointInElementList, ColElementSizeList, 
02375            1.0, BlockRowValues, Yp, -1.0, Yp, NumVectors);
02376     }
02377   }
02378   else {
02379     for (i=0; i < NumMyBlockRows_; i++) {
02380       int      NumEntries = *NumBlockEntriesPerRow++;
02381       int *    BlockRowIndices = *Indices++;
02382       Epetra_SerialDenseMatrix** BlockRowValues  = *Entries++;
02383       int  yoff = *FirstPointInElement++;
02384       int RowDim = *ElementSize++;
02385       BlockRowMultiply(TransA, RowDim, NumEntries, BlockRowIndices, yoff, 
02386            ColFirstPointInElementList, ColElementSizeList, 
02387            1.0, BlockRowValues, Yp, -1.0, Yp, NumVectors);
02388     }
02389   }
02390 
02391   UpdateFlops(2*NumVectors*NumGlobalNonzeros());
02392   return(0);
02393 }
02394 
02395 //=============================================================================
02396 int Epetra_VbrMatrix::InvRowSums(Epetra_Vector& x) const {
02397   //
02398   // Put inverse of the sum of absolute values of the ith row of A in x[i].
02399   //
02400   EPETRA_CHK_ERR(InverseSums(true, x));
02401   return(0);
02402 }
02403 
02404 //=============================================================================
02405 int Epetra_VbrMatrix::InvColSums(Epetra_Vector& x) const {
02406   //
02407   // Put inverse of the sum of absolute values of the jth column of A in x[j].
02408   //
02409   EPETRA_CHK_ERR(InverseSums(false, x));
02410   return(0);
02411 }
02412 //=============================================================================
02413 int Epetra_VbrMatrix::InverseSums(bool DoRows, Epetra_Vector& x) const {
02414   //
02415   // Put inverse of the sum of absolute values of the ith row of A in x[i].
02416   //
02417 
02418   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
02419   bool hasOperatorMap = false;
02420   if (DoRows) {
02421     if ( !Graph().RangeMap().SameAs(x.Map()) ) {
02422       hasOperatorMap = OperatorRangeMap().SameAs(x.Map());
02423       if( !hasOperatorMap )
02424         EPETRA_CHK_ERR(-2);
02425     }
02426   }
02427   else {
02428     if ( !Graph().DomainMap().SameAs(x.Map()) ) {
02429       hasOperatorMap = OperatorDomainMap().SameAs(x.Map());
02430       if( !hasOperatorMap )
02431         EPETRA_CHK_ERR(-2);
02432     }
02433   }
02434   int ierr = 0;
02435   int * NumBlockEntriesPerRow = NumBlockEntriesPerRow_;
02436   int ** Indices = Indices_;
02437   Epetra_SerialDenseMatrix*** Entries = Entries_;
02438   
02439   int * RowElementSizeList = ElementSizeList_;
02440   int * RowFirstPointInElementList = FirstPointInElementList_;
02441   int * ColElementSizeList = ElementSizeList_;
02442   int * ColFirstPointInElementList = FirstPointInElementList_;
02443   if (Importer()!=0) {
02444     ColElementSizeList = ColMap().ElementSizeList();
02445     ColFirstPointInElementList = ColMap().FirstPointInElementList();
02446   }
02447 
02448   x.PutScalar(0.0); // Zero out result vector
02449 
02450   double * xp = (double*)x.Values();
02451 
02452   // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
02453   Epetra_Vector * x_tmp = 0;
02454   if (!DoRows) {
02455     if (Importer()!=0) {
02456       x_tmp = new Epetra_Vector(ColMap()); // Create import vector if needed
02457       xp = (double*)x_tmp->Values();
02458     }
02459   }
02460 
02461   for (int i=0; i < NumMyBlockRows_; i++) {
02462     int      NumEntries = *NumBlockEntriesPerRow++;
02463     int *    BlockRowIndices = *Indices++;
02464     Epetra_SerialDenseMatrix** BlockRowValues  = *Entries++;
02465     int xoff = *RowFirstPointInElementList++;
02466     int RowDim = *RowElementSizeList++;
02467     if (DoRows) {
02468       for (int ii=0; ii < NumEntries; ii++) {
02469   double * xptr = xp+xoff;
02470   double * A = BlockRowValues[ii]->A();
02471   int LDA = BlockRowValues[ii]->LDA();
02472   int BlockIndex = BlockRowIndices[ii];
02473   int ColDim = ColElementSizeList[BlockIndex];
02474   for (int j=0; j<ColDim; j++) {
02475     double * curEntry = A + j*LDA;
02476     for (int k=0; k<RowDim; k++)
02477       xptr[k] += fabs(*curEntry++);
02478   }
02479       }
02480     }
02481     else {
02482       for (int ii=0; ii < NumEntries; ii++) {
02483   double * A = BlockRowValues[ii]->A();
02484   int LDA = BlockRowValues[ii]->LDA();
02485   int BlockIndex = BlockRowIndices[ii];
02486   int off = ColFirstPointInElementList[BlockIndex];
02487   int ColDim = ColElementSizeList[BlockIndex];
02488   double * curx = xp+off;
02489   for (int j=0; j<ColDim; j++) {
02490     double * curEntry = A + j*LDA;
02491     for (int k=0; k<RowDim; k++)
02492       curx[j] += fabs(*curEntry++);
02493   }
02494       }
02495     }
02496   }
02497 
02498   if (!DoRows) {
02499     if (Importer()!=0){
02500       Epetra_Vector  *x_blocked = 0;
02501       if(hasOperatorMap)
02502         x_blocked = new Epetra_Vector( ::View, DoRows ? Graph().RangeMap() : Graph().DomainMap(), &x[0] );
02503       else
02504         x_blocked = &x;
02505       x_blocked->PutScalar(0.0);
02506       EPETRA_CHK_ERR(x_blocked->Export(*x_tmp, *Importer(), Add)); // Fill x with Values from import vector
02507       if(hasOperatorMap)
02508         delete x_blocked;
02509       delete x_tmp;
02510       xp = (double*) x.Values();
02511     }
02512   }
02513   int NumMyRows_ = NumMyRows();
02514   {for (int i=0; i < NumMyRows_; i++) {
02515     double scale = xp[i];
02516     if (scale<Epetra_MinDouble) {
02517       if (scale==0.0) ierr = 1; // Set error to 1 to signal that zero row/col sum found (supercedes ierr = 2)
02518       else if (ierr!=1) ierr = 2;
02519       xp[i] = Epetra_MaxDouble;
02520     }
02521     else
02522       xp[i] = 1.0/scale;
02523   }}
02524   UpdateFlops(NumGlobalNonzeros());
02525 
02526   EPETRA_CHK_ERR(ierr);
02527   return(0);
02528 }
02529 //=============================================================================
02530 int Epetra_VbrMatrix::LeftScale(const Epetra_Vector& x) {
02531   //
02532   // Multiply the ith row of A by x[i].
02533   //
02534   EPETRA_CHK_ERR(Scale(true, x));
02535   return(0);
02536 }
02537 
02538 //=============================================================================
02539 int Epetra_VbrMatrix::RightScale(const Epetra_Vector& x) {
02540   //
02541   // Multiply the jth column of A by x[j].
02542   //
02543   EPETRA_CHK_ERR(Scale (false, x));
02544   return(0);
02545 }
02546 //=============================================================================
02547 int Epetra_VbrMatrix::Scale(bool DoRows, const Epetra_Vector& x) {
02548 
02549   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
02550   bool hasOperatorMap = false;
02551   if (DoRows) {
02552     if ( !Graph().RangeMap().SameAs(x.Map()) ) {
02553       hasOperatorMap = OperatorRangeMap().SameAs(x.Map());
02554       if( !hasOperatorMap )
02555         EPETRA_CHK_ERR(-2);
02556     }
02557   }
02558   else {
02559     if ( !Graph().DomainMap().SameAs(x.Map()) ) {
02560       hasOperatorMap = OperatorDomainMap().SameAs(x.Map());
02561       if( !hasOperatorMap )
02562         EPETRA_CHK_ERR(-2);
02563     }
02564   }
02565   int ierr = 0;
02566   int * NumBlockEntriesPerRow = NumBlockEntriesPerRow_;
02567   int ** Indices = Indices_;
02568   Epetra_SerialDenseMatrix*** Entries = Entries_;
02569   
02570   int * RowElementSizeList = ElementSizeList_;
02571   int * RowFirstPointInElementList = FirstPointInElementList_;
02572   int * ColElementSizeList = ElementSizeList_;
02573   int * ColFirstPointInElementList = FirstPointInElementList_;
02574   if (Importer()!=0) {
02575     ColElementSizeList = ColMap().ElementSizeList();
02576     ColFirstPointInElementList = ColMap().FirstPointInElementList();
02577   }
02578 
02579   double * xp = (double*)x.Values();
02580 
02581   // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
02582   Epetra_Vector * x_tmp = 0;
02583   if (!DoRows) {
02584     if (Importer()!=0) {
02585       x_tmp = new Epetra_Vector(ColMap()); // Create import vector if needed
02586       EPETRA_CHK_ERR(x_tmp->Import(x,*Importer(), Insert)); // x_tmp will have all the values we need
02587       xp = (double*)x_tmp->Values();
02588     }
02589   }
02590 
02591   for (int i=0; i < NumMyBlockRows_; i++) {
02592     int      NumEntries = *NumBlockEntriesPerRow++;
02593     int *    BlockRowIndices = *Indices++;
02594     Epetra_SerialDenseMatrix** BlockRowValues  = *Entries++;
02595     int xoff = *RowFirstPointInElementList++;
02596     int RowDim = *RowElementSizeList++;
02597     if (DoRows) {
02598       for (int ii=0; ii < NumEntries; ii++) {
02599   double * xptr = xp+xoff;
02600   double * A = BlockRowValues[ii]->A();
02601   int LDA = BlockRowValues[ii]->LDA();
02602   int BlockIndex = BlockRowIndices[ii];
02603   int ColDim = ColElementSizeList[BlockIndex];
02604   for (int j=0; j<ColDim; j++) {
02605     double * curEntry = A + j*LDA;
02606     for (int k=0; k<RowDim; k++)
02607       *curEntry++ *= xptr[k];
02608   }
02609       }
02610     }
02611     else {
02612       for (int ii=0; ii < NumEntries; ii++) {
02613   double * A = BlockRowValues[ii]->A();
02614   int LDA = BlockRowValues[ii]->LDA();
02615   int BlockIndex = BlockRowIndices[ii];
02616   int off = ColFirstPointInElementList[BlockIndex];
02617   int ColDim = ColElementSizeList[BlockIndex];
02618   double * curx = xp+off;
02619   for (int j=0; j<ColDim; j++) {
02620     double * curEntry = A + j*LDA;
02621     for (int k=0; k<RowDim; k++)
02622       *curEntry++ *= curx[j];
02623   }
02624       }
02625     }
02626   }
02627 
02628   if (x_tmp!=0) delete x_tmp;
02629   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
02630   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
02631   NormFrob_ = -1.0;
02632   UpdateFlops(NumGlobalNonzeros());
02633 
02634   EPETRA_CHK_ERR(ierr);
02635   return(0);
02636 }
02637 //=============================================================================
02638 double Epetra_VbrMatrix::NormInf() const {
02639 
02640 #if 0
02641   //
02642   //  Commenting this section out disables caching, ie. 
02643   //  causes the norm to be computed each time NormInf is called.
02644   //  See bug #1151 for a full discussion.  
02645   //
02646   double MinNorm ; 
02647   Comm().MinAll( &NormInf_, &MinNorm, 1 ) ; 
02648 
02649   if( MinNorm >= 0.0) 
02650   //  if( NormInf_ >= 0.0) 
02651     return(NormInf_);
02652 #endif
02653 
02654   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
02655 
02656   int MaxRowDim_ = MaxRowDim();
02657   double * tempv = new double[MaxRowDim_];
02658 
02659   int * NumBlockEntriesPerRow = NumBlockEntriesPerRow_;
02660   int * ElementSize = ElementSizeList_;
02661   Epetra_SerialDenseMatrix*** Entries = Entries_;
02662 
02663   double Local_NormInf = 0.0;
02664   for (int i=0; i < NumMyBlockRows_; i++) {
02665     int      NumEntries = *NumBlockEntriesPerRow++ ;
02666     int RowDim = *ElementSize++;
02667     Epetra_SerialDenseMatrix** BlockRowValues  = *Entries++;
02668     BlockRowNormInf(RowDim, NumEntries, 
02669         BlockRowValues, tempv);
02670     for (int j=0; j < RowDim; j++) Local_NormInf = EPETRA_MAX(Local_NormInf, tempv[j]);
02671   }
02672   Comm().MaxAll(&Local_NormInf, &NormInf_, 1);
02673   delete [] tempv;
02674   UpdateFlops(NumGlobalNonzeros());
02675   return(NormInf_);
02676 }
02677 //=============================================================================
02678 void Epetra_VbrMatrix::BlockRowNormInf(int RowDim, int NumEntries,
02679                Epetra_SerialDenseMatrix** As, 
02680                double * Y) const
02681 {
02682   int i, j, k;
02683   for (k=0; k<RowDim; k++) Y[k] = 0.0;
02684 
02685   for (i=0; i < NumEntries; i++) {
02686     double * A = As[i]->A();
02687     int LDA = As[i]->LDA();
02688     int ColDim = As[i]->N();
02689     for (j=0; j<ColDim; j++) {
02690       for (k=0; k<RowDim; k++) Y[k] += fabs(A[k]);
02691       A += LDA;
02692     }
02693   }
02694   return;
02695 }
02696 //=============================================================================
02697 double Epetra_VbrMatrix::NormOne() const {
02698 
02699 #if 0
02700   //
02701   //  Commenting this section out disables caching, ie. 
02702   //  causes the norm to be computed each time NormOne is called.  
02703   //  See bug #1151 for a full discussion.  
02704   //
02705   double MinNorm ; 
02706   Comm().MinAll( &NormOne_, &MinNorm, 1 ) ; 
02707 
02708   if( MinNorm >= 0.0) 
02709     return(NormOne_);
02710 #endif
02711 
02712   int * ColFirstPointInElementList = FirstPointInElementList_;
02713   if (Importer()!=0) {
02714     ColFirstPointInElementList = ColMap().FirstPointInElementList();
02715   }
02716 
02717   Epetra_Vector * x = new Epetra_Vector(RowMap()); // Need temp vector for column sums
02718   
02719   double * xp = (double*)x->Values();
02720   Epetra_MultiVector * x_tmp = 0;
02721 
02722   // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
02723   if (Importer()!=0) {
02724     x_tmp = new Epetra_Vector(ColMap()); // Create temporary import vector if needed
02725     xp = (double*)x_tmp->Values();
02726   }
02727 
02728   int * NumBlockEntriesPerRow = NumBlockEntriesPerRow_;
02729   int * ElementSize = ElementSizeList_;
02730   int ** Indices = Indices_;
02731   Epetra_SerialDenseMatrix*** Entries = Entries_;
02732 
02733   for (int i=0; i < NumMyBlockRows_; i++) {
02734     int NumEntries = *NumBlockEntriesPerRow++;
02735     int RowDim = *ElementSize++;
02736     int *    BlockRowIndices = *Indices++;
02737     Epetra_SerialDenseMatrix** BlockRowValues  = *Entries++;
02738     BlockRowNormOne(RowDim, NumEntries, BlockRowIndices,
02739         BlockRowValues,  ColFirstPointInElementList, xp);
02740   }
02741   if (Importer()!=0) {
02742     x->PutScalar(0.0);
02743     EPETRA_CHK_ERR(x->Export(*x_tmp, *Importer(), Add));
02744   } // Fill x with Values from temp vector
02745   x->MaxValue(&NormOne_); // Find max
02746   if (x_tmp!=0) delete x_tmp;
02747   delete x;
02748   UpdateFlops(NumGlobalNonzeros());
02749   return(NormOne_);
02750 }
02751 //=============================================================================
02752 double Epetra_VbrMatrix::NormFrobenius() const {
02753 
02754 #if 0
02755   //
02756   //  Commenting this section out disables caching, ie. 
02757   //  causes the norm to be computed each time NormOne is called.  
02758   //  See bug #1151 for a full discussion.  
02759   //
02760   double MinNorm ; 
02761   Comm().MinAll( &NormFrob_, &MinNorm, 1 ) ; 
02762 
02763   if( MinNorm >= 0.0) 
02764     return(NormFrob_);
02765 #endif
02766 
02767   int * NumBlockEntriesPerRow = NumBlockEntriesPerRow_;
02768   int * ElementSize = ElementSizeList_;
02769   Epetra_SerialDenseMatrix*** Entries = Entries_;
02770 
02771   double local_sum = 0.0;
02772 
02773   for (int i=0; i < NumMyBlockRows_; i++) {
02774     int NumEntries = *NumBlockEntriesPerRow++;
02775     int RowDim = *ElementSize++;
02776     Epetra_SerialDenseMatrix** BlockRowValues  = *Entries++;
02777 
02778     for(int ii=0; ii<NumEntries; ++ii) {
02779       double* A = BlockRowValues[ii]->A();
02780       int LDA = BlockRowValues[ii]->LDA();
02781       int colDim = BlockRowValues[ii]->N();
02782       for(int j=0; j<colDim; ++j) {
02783         for(int k=0; k<RowDim; ++k) {
02784           local_sum += A[k]*A[k];
02785         }
02786         A += LDA;
02787       }
02788     }
02789 
02790 //The following commented-out code performs the calculation by running
02791 //all the way across each point-entry row before going to the next.
02792 //This is the order in which the CrsMatrix frobenius norm is calculated,
02793 //but for a VbrMatrix it is faster to run the block-entries in
02794 //column-major order (which is how they're stored), as the above code
02795 //does.
02796 //But if the same matrix is stored in both Vbr and Crs form, and you wish
02797 //to compare the frobenius norms, it is necessary to run through the
02798 //coefficients in the same order to avoid round-off differences.
02799 //
02800 //    for(int k=0; k<RowDim; ++k) {
02801 //      for(int ii=0; ii<NumEntries; ++ii) {
02802 //        double* A = BlockRowValues[ii]->A();
02803 //        int LDA = BlockRowValues[ii]->LDA();
02804 //        int colDim = BlockRowValues[ii]->N();
02805 //        for(int j=0; j<colDim; ++j) {
02806 //          double Ak = A[k+j*LDA];
02807 //          local_sum += Ak*Ak;
02808 //        }
02809 //      }
02810 //    }
02811 
02812   }
02813 
02814   double global_sum = 0.0;
02815   Comm().SumAll(&local_sum, &global_sum, 1);
02816 
02817   NormFrob_ = std::sqrt(global_sum);
02818 
02819   UpdateFlops(NumGlobalNonzeros());
02820 
02821   return(NormFrob_);
02822 }
02823 //=============================================================================
02824 void Epetra_VbrMatrix::BlockRowNormOne(int RowDim, int NumEntries, int * BlockRowIndices,
02825                Epetra_SerialDenseMatrix** As, 
02826                int * ColFirstPointInElementList, double * x) const {
02827   int i, j, k;
02828 
02829   for (i=0; i < NumEntries; i++) {
02830     double * A = As[i]->A();
02831     int LDA = As[i]->LDA();
02832     int ColDim = As[i]->N();
02833     double * curx = x + ColFirstPointInElementList[BlockRowIndices[i]];
02834     for (j=0; j<ColDim; j++) {
02835       for (k=0; k<RowDim; k++) curx[j] += fabs(A[k]);
02836       A += LDA;
02837     }
02838   }
02839   return;
02840 }
02841 //=========================================================================
02842 int Epetra_VbrMatrix::CheckSizes(const Epetra_SrcDistObject & Source) {
02843   const Epetra_VbrMatrix & A = dynamic_cast<const Epetra_VbrMatrix &>(Source);
02844   if (!A.Graph().GlobalConstantsComputed()) EPETRA_CHK_ERR(-1); // Must have global constants to proceed
02845   return(0);
02846 }
02847 //=========================================================================
02848 int Epetra_VbrMatrix::CopyAndPermute(const Epetra_SrcDistObject & Source,
02849                                      int NumSameIDs, 
02850                                      int NumPermuteIDs,
02851                                      int * PermuteToLIDs,
02852                                      int *PermuteFromLIDs,
02853                                      const Epetra_OffsetIndex * Indexor)
02854 {
02855   (void)Indexor; 
02856   const Epetra_VbrMatrix & A = dynamic_cast<const Epetra_VbrMatrix &>(Source);
02857   int i, j;
02858   
02859   int BlockRow, NumBlockEntries;
02860   int * BlockIndices;
02861   int RowDim;
02862   Epetra_SerialDenseMatrix ** Entries;
02863   int FromBlockRow, ToBlockRow;
02864   
02865   // Do copy first
02866   if (NumSameIDs>0) {
02867     int MaxNumBlockEntries = A.MaxNumBlockEntries();
02868     BlockIndices = new int[MaxNumBlockEntries];  // Need some temporary space
02869  
02870 
02871     for (i=0; i<NumSameIDs; i++) {
02872       BlockRow = GRID(i);
02873       EPETRA_CHK_ERR( A.ExtractGlobalBlockRowPointers(BlockRow,
02874                                                       MaxNumBlockEntries,
02875                   RowDim, NumBlockEntries, 
02876                   BlockIndices, Entries)); // Set pointers
02877       // Place into target matrix.  Depends on Epetra_DataAccess copy/view and static/dynamic graph.
02878       if (StaticGraph() || IndicesAreLocal()) {
02879   EPETRA_CHK_ERR(BeginReplaceGlobalValues(BlockRow, NumBlockEntries,
02880                                                 BlockIndices));
02881       }
02882       else {
02883   EPETRA_CHK_ERR(BeginInsertGlobalValues(BlockRow, NumBlockEntries, BlockIndices));
02884       }
02885       // Insert block entries one-at-a-time
02886       for (j=0; j<NumBlockEntries; j++) SubmitBlockEntry(Entries[j]->A(),
02887                Entries[j]->LDA(),
02888                RowDim, Entries[j]->N());
02889       EndSubmitEntries(); // Complete this block row
02890     }
02891     delete [] BlockIndices;
02892   }
02893 
02894   // Do local permutation next
02895   if (NumPermuteIDs>0) {
02896     int MaxNumBlockEntries = A.MaxNumBlockEntries();
02897     BlockIndices = new int[MaxNumBlockEntries];  // Need some temporary space
02898       
02899     for (i=0; i<NumPermuteIDs; i++) {
02900       FromBlockRow = A.GRID(PermuteFromLIDs[i]);
02901       ToBlockRow = GRID(PermuteToLIDs[i]);
02902       EPETRA_CHK_ERR(A.ExtractGlobalBlockRowPointers(FromBlockRow, MaxNumBlockEntries, RowDim, NumBlockEntries, 
02903                BlockIndices, Entries)); // Set pointers
02904       // Place into target matrix.  Depends on Epetra_DataAccess copy/view and static/dynamic graph.
02905       if (StaticGraph() || IndicesAreLocal()) {
02906   EPETRA_CHK_ERR(BeginReplaceGlobalValues(ToBlockRow, NumBlockEntries, BlockIndices));
02907       }
02908       else {
02909   EPETRA_CHK_ERR(BeginInsertGlobalValues(ToBlockRow, NumBlockEntries, BlockIndices));
02910       }
02911       // Insert block entries one-at-a-time
02912       for (j=0; j<NumBlockEntries; j++) SubmitBlockEntry(Entries[j]->A(),
02913                Entries[j]->LDA(), RowDim, Entries[j]->N());
02914       EndSubmitEntries(); // Complete this block row
02915     }
02916     delete [] BlockIndices;
02917   }
02918  
02919   return(0);
02920 }
02921 
02922 //=========================================================================
02923 int Epetra_VbrMatrix::PackAndPrepare(const Epetra_SrcDistObject & Source,
02924                                      int NumExportIDs,
02925                                      int * ExportLIDs,
02926                                      int & LenExports,
02927                                      char * & Exports,
02928                                      int & SizeOfPacket,
02929                                      int * Sizes,
02930                                      bool & VarSizes,
02931                                      Epetra_Distributor & Distor)
02932 {
02933   (void)LenExports;
02934   (void)Sizes;
02935   (void)VarSizes;
02936   (void)Distor; 
02937   const Epetra_VbrMatrix & A = dynamic_cast<const Epetra_VbrMatrix &>(Source);
02938 
02939   double * DoubleExports = 0;
02940   int GlobalMaxNumNonzeros = A.GlobalMaxNumNonzeros();
02941   int GlobalMaxNumBlockEntries = A.GlobalMaxNumBlockEntries();
02942   // Will have GlobalMaxNumEntries doubles, GlobalMaxNumEntries +2 ints, pack them interleaved
02943   int DoublePacketSize = GlobalMaxNumNonzeros +  
02944     (((2*GlobalMaxNumBlockEntries+3)+(int)sizeof(int)-1)*(int)sizeof(int))/(int)sizeof(double);
02945   SizeOfPacket = DoublePacketSize * (int)sizeof(double); 
02946 
02947   if (DoublePacketSize*NumExportIDs>LenExports_) {
02948     if (LenExports_>0) delete [] Exports_;
02949     LenExports_ = DoublePacketSize*NumExportIDs;
02950     DoubleExports = new double[LenExports_];
02951     Exports_ = (char *) DoubleExports;
02952   }
02953 
02954   if (NumExportIDs<=0) return(0); // All done if nothing to pack
02955 
02956   int i, j;
02957   
02958   int NumBlockEntries;
02959   int * BlockIndices;
02960   int RowDim, * ColDims;
02961   double * Entries;
02962   int FromBlockRow;
02963   double * valptr, * dintptr;
02964   int * intptr;
02965   
02966   // Each segment of IntExports will be filled by a packed row of information for each row as follows:
02967   // 1st int: GRID of row where GRID is the global row ID for the source matrix
02968   // next int:  RowDim of Block Row
02969   // next int: NumBlockEntries, Number of indices in row.
02970   // next NumBlockEntries: The actual indices for the row.
02971   // Any remaining space (of length GlobalMaxNumNonzeros - NumBlockEntries ints) will be wasted but we need fixed
02972   //   sized segments for current communication routines.
02973 
02974   // Each segment of Exports will be filled with values.
02975 
02976   valptr = (double *) Exports;
02977   dintptr = valptr + GlobalMaxNumNonzeros;
02978   intptr = (int *) dintptr;
02979   bool NoSumInto = false;
02980   for (i=0; i<NumExportIDs; i++) {
02981     FromBlockRow = A.GRID(ExportLIDs[i]);
02982     BlockIndices = intptr + 3;
02983     ColDims = BlockIndices + GlobalMaxNumBlockEntries;
02984     EPETRA_CHK_ERR(A.BeginExtractGlobalBlockRowCopy(FromBlockRow, GlobalMaxNumBlockEntries, RowDim,
02985               NumBlockEntries, BlockIndices, ColDims));
02986     // Now extract each block entry into send buffer
02987     Entries = valptr;
02988     for (j=0; j<NumBlockEntries; j++) {
02989       int SizeOfValues = RowDim*ColDims[j];
02990       A.ExtractEntryCopy(SizeOfValues, Entries, RowDim, NoSumInto);
02991       Entries += SizeOfValues;
02992     }
02993     // Fill first three slots of intptr with info
02994     intptr[0] = FromBlockRow;
02995     intptr[1] = RowDim;
02996     intptr[2] = NumBlockEntries;
02997     valptr += DoublePacketSize; // Point to next segment
02998     dintptr = valptr + GlobalMaxNumNonzeros;
02999     intptr = (int *) dintptr;
03000   }
03001     
03002   return(0);
03003 }
03004 
03005 //=========================================================================
03006 int Epetra_VbrMatrix::UnpackAndCombine(const Epetra_SrcDistObject & Source, 
03007                                        int NumImportIDs,
03008                                        int * ImportLIDs, 
03009                                        int LenImports,
03010                                        char * Imports,
03011                                        int & SizeOfPacket, 
03012                                        Epetra_Distributor & Distor, 
03013                                        Epetra_CombineMode CombineMode,
03014                                        const Epetra_OffsetIndex * Indexor)
03015 {
03016   (void)LenImports;
03017   (void)SizeOfPacket;
03018   (void)Distor;
03019   (void)Indexor;
03020   if (NumImportIDs<=0) return(0);
03021 
03022   if(   CombineMode != Add
03023   && CombineMode != Zero
03024   && CombineMode != Insert )
03025     EPETRA_CHK_ERR(-1); // CombineMode not supported, default to mode Zero
03026 
03027   const Epetra_VbrMatrix & A = dynamic_cast<const Epetra_VbrMatrix &>(Source);
03028   int NumBlockEntries;
03029   int * BlockIndices;
03030   int RowDim, * ColDims;
03031   double * Values;
03032   int ToBlockRow;
03033   int i, j;
03034   
03035   double * valptr, *dintptr;
03036   int * intptr;
03037   int GlobalMaxNumNonzeros = A.GlobalMaxNumNonzeros();
03038   int GlobalMaxNumBlockEntries = A.GlobalMaxNumBlockEntries();
03039   // Will have GlobalMaxNumEntries doubles, GlobalMaxNumEntries +2 ints, pack them interleaved
03040   int DoublePacketSize = GlobalMaxNumNonzeros +  
03041     (((2*GlobalMaxNumBlockEntries+3)+(int)sizeof(int)-1)*(int)sizeof(int))/(int)sizeof(double);
03042   // Unpack it...
03043 
03044 
03045   // Each segment of IntImports will be filled by a packed row of information for each row as follows:
03046   // 1st int: GRID of row where GRID is the global row ID for the source matrix
03047   // next int:  NumBlockEntries, Number of indices in row.
03048   // next int:  RowDim of Block Row
03049   // next NumBlockEntries: The actual indices for the row.
03050   // Any remaining space (of length GlobalMaxNumNonzeros - NumBlockEntries ints) will be 
03051   //  wasted but we need fixed sized segments for current communication routines.
03052 
03053   valptr = (double *) Imports;
03054   dintptr = valptr + GlobalMaxNumNonzeros;
03055   intptr = (int *) dintptr;
03056     
03057   for (i=0; i<NumImportIDs; i++) {
03058     ToBlockRow = GRID(ImportLIDs[i]);
03059     assert((intptr[0])==ToBlockRow); // Sanity check
03060     RowDim = RowMap().ElementSize(ImportLIDs[i]);
03061     assert((intptr[1])==RowDim); // Sanity check
03062     NumBlockEntries = intptr[2];
03063     BlockIndices = intptr + 3; 
03064     ColDims = BlockIndices + GlobalMaxNumBlockEntries;
03065     if (CombineMode==Add) {
03066       if (StaticGraph() || IndicesAreLocal()) {
03067   // Replace any current values
03068   EPETRA_CHK_ERR(BeginSumIntoGlobalValues(ToBlockRow, NumBlockEntries, BlockIndices));
03069       }
03070       else {
03071   // Insert values
03072   EPETRA_CHK_ERR(BeginInsertGlobalValues(ToBlockRow, NumBlockEntries, BlockIndices));
03073       }
03074     }
03075     else if (CombineMode==Insert) {
03076       if (StaticGraph() || IndicesAreLocal()) {
03077   // Replace any current values
03078   EPETRA_CHK_ERR(BeginReplaceGlobalValues(ToBlockRow, NumBlockEntries, BlockIndices));
03079       }
03080       else {
03081   // Insert values
03082   EPETRA_CHK_ERR(BeginInsertGlobalValues(ToBlockRow, NumBlockEntries, BlockIndices));
03083       }
03084     }
03085     // Now extract each block entry into send buffer
03086     Values = valptr;
03087     for (j=0; j<NumBlockEntries; j++) {
03088       int LDA = RowDim;
03089       int ColDim = ColDims[j];
03090       SubmitBlockEntry(Values, LDA, RowDim, ColDim);
03091       Values += (LDA*ColDim);
03092     }
03093     EndSubmitEntries(); // Done with this block row
03094     valptr += DoublePacketSize; // Point to next segment
03095     dintptr = valptr + GlobalMaxNumNonzeros;
03096     intptr = (int *) dintptr;
03097   }
03098   
03099   return(0);
03100 }
03101 //=========================================================================
03102 int Epetra_VbrMatrix::GeneratePointObjects() const {
03103 
03104   if (HavePointObjects_)  return(0); // Already done
03105 
03106   // Generate a point-wise compatible row map
03107   EPETRA_CHK_ERR(BlockMap2PointMap(RowMap(), RowMatrixRowMap_));
03108 
03109   // For each of the next maps, first check if the corresponding block map
03110   // is the same as the block row map for this matrix.  If so, then we simply
03111   // copy the pointer.  Otherwise we create a new point map.
03112 
03113   // This check can save storage and time since it will often be the case that the
03114   // domain, range and row block maps will be the same.  Also, in the serial case,
03115   // the column block map will also often be the same as the row block map.
03116 
03117   if (ColMap().SameAs(RowMap())) {
03118     RowMatrixColMap_ = RowMatrixRowMap_;
03119   }
03120   else {
03121     EPETRA_CHK_ERR(BlockMap2PointMap(ColMap(), RowMatrixColMap_));
03122   }
03123 
03124   if (DomainMap().SameAs(RowMap())) {
03125     OperatorDomainMap_ = RowMatrixRowMap_;
03126   }
03127   else {
03128     EPETRA_CHK_ERR(BlockMap2PointMap(DomainMap(), OperatorDomainMap_));
03129   }
03130   if (RangeMap().SameAs(RowMap())) {
03131     OperatorRangeMap_ = RowMatrixRowMap_;
03132   }
03133   else {
03134     EPETRA_CHK_ERR(BlockMap2PointMap(RangeMap(), OperatorRangeMap_));
03135   }
03136 
03137   // Finally generate Importer that will migrate needed domain elements to the column map
03138   // layout.
03139   RowMatrixImporter_ = new Epetra_Import(*RowMatrixColMap_, *OperatorDomainMap_);
03140 
03141   HavePointObjects_ = true;
03142   return(0);
03143 }
03144 //=========================================================================
03145 int Epetra_VbrMatrix::BlockMap2PointMap(const Epetra_BlockMap & BlockMap,
03146           Epetra_Map * & PointMap) const
03147 {
03148   // Generate an Epetra_Map that has the same number and distribution of points
03149   // as the input Epetra_BlockMap object.  The global IDs for the output PointMap
03150   // are computed by using the MaxElementSize of the BlockMap.  For variable block
03151   // sizes this will create gaps in the GID space, but that is OK for Epetra_Maps.
03152 
03153   int MaxElementSize = BlockMap.MaxElementSize();
03154   int PtNumMyElements = BlockMap.NumMyPoints();
03155   int * PtMyGlobalElements = 0;
03156   if (PtNumMyElements>0) PtMyGlobalElements = new int[PtNumMyElements];
03157 
03158   int NumMyElements = BlockMap.NumMyElements();
03159 
03160   int curID = 0;
03161   for (int i=0; i<NumMyElements; i++) {
03162     int StartID = BlockMap.GID(i)*MaxElementSize;
03163     int ElementSize = BlockMap.ElementSize(i);
03164     for (int j=0; j<ElementSize; j++) PtMyGlobalElements[curID++] = StartID+j;
03165   }
03166   assert(curID==PtNumMyElements); // Sanity test
03167 
03168   PointMap = new Epetra_Map(-1, PtNumMyElements, PtMyGlobalElements, BlockMap.IndexBase(), BlockMap.Comm());
03169 
03170   if (PtNumMyElements>0) delete [] PtMyGlobalElements;
03171 
03172   if (!BlockMap.PointSameAs(*PointMap)) {EPETRA_CHK_ERR(-1);} // Maps not compatible
03173   return(0);
03174 }
03175 //=========================================================================
03176 int Epetra_VbrMatrix::UpdateOperatorXY(const Epetra_MultiVector& X,
03177                const Epetra_MultiVector& Y) const
03178 {
03179   if (OperatorX_!=0)
03180     if (OperatorX_->NumVectors()!=X.NumVectors()) {delete OperatorX_; OperatorX_ = 0; delete OperatorY_; OperatorY_=0;}
03181   if (OperatorX_==0) {
03182     if (!X.Map().PointSameAs(DomainMap())) EPETRA_CHK_ERR(-1); // X not point-wise compatible with the block domain map
03183     if (!Y.Map().PointSameAs(RangeMap())) EPETRA_CHK_ERR(-2); // Y not point-wise compatible with the block col map
03184     OperatorX_ = new Epetra_MultiVector(View, DomainMap(), X.Pointers(), X.NumVectors());
03185     OperatorY_ = new Epetra_MultiVector(View, RangeMap(), Y.Pointers(), Y.NumVectors());
03186   } 
03187   else {
03188     EPETRA_CHK_ERR(OperatorX_->ResetView(X.Pointers()));
03189     EPETRA_CHK_ERR(OperatorY_->ResetView(Y.Pointers()));
03190   }
03191   return(0);
03192 }
03193 //=========================================================================
03194 int Epetra_VbrMatrix::Apply(const Epetra_MultiVector& X,
03195           Epetra_MultiVector& Y) const
03196 {
03197   if (!Epetra_VbrMatrix::UseTranspose()) {
03198     EPETRA_CHK_ERR(UpdateOperatorXY(X,Y)); // Update X and Y vector whose maps are compatible with the Vbr Matrix
03199     EPETRA_CHK_ERR(Epetra_VbrMatrix::Multiply(Epetra_VbrMatrix::UseTranspose(), *OperatorX_, *OperatorY_));
03200   } 
03201   else { // Swap role of OperatorX_ and OperatorY_ to remain compatible with domain and range spaces.
03202     EPETRA_CHK_ERR(UpdateOperatorXY(Y,X)); // Update X and Y vector whose maps are compatible with the Vbr Matrix
03203     EPETRA_CHK_ERR(Epetra_VbrMatrix::Multiply(Epetra_VbrMatrix::UseTranspose(), *OperatorY_, *OperatorX_));
03204   }
03205   return(0);
03206 }
03207 //=========================================================================
03208 int Epetra_VbrMatrix::ApplyInverse(const Epetra_MultiVector& X,
03209            Epetra_MultiVector& Y) const
03210 {
03211   if (!Epetra_VbrMatrix::UseTranspose()) {
03212     EPETRA_CHK_ERR(UpdateOperatorXY(X,Y)); // Update X and Y vector whose maps are compatible with the Vbr Matrix
03213     EPETRA_CHK_ERR(Solve(UpperTriangular(), Epetra_VbrMatrix::UseTranspose(), NoDiagonal(), *OperatorX_, *OperatorY_));
03214   }
03215   else { // Swap role of OperatorX_ and OperatorY_ to remain compatible with domain and range spaces.
03216     EPETRA_CHK_ERR(UpdateOperatorXY(Y,X)); // Update X and Y vector whose maps are compatible with the Vbr Matrix
03217     EPETRA_CHK_ERR(Solve(UpperTriangular(), Epetra_VbrMatrix::UseTranspose(), NoDiagonal(), *OperatorY_, *OperatorX_));
03218   }
03219   return(0);
03220 }
03221 //=========================================================================
03222 void Epetra_VbrMatrix::Print(ostream& os) const {
03223   int MyPID = RowMap().Comm().MyPID();
03224   int NumProc = RowMap().Comm().NumProc();
03225 
03226   for (int iproc=0; iproc < NumProc; iproc++) {
03227     if (MyPID==iproc) {
03228       if (MyPID==0) {
03229   os <<  "\nNumber of Global Block Rows  = "; os << NumGlobalBlockRows(); os << endl;
03230   os <<    "Number of Global Block Cols  = "; os << NumGlobalBlockCols(); os << endl;
03231   os <<    "Number of Global Block Diags = "; os << NumGlobalBlockDiagonals(); os << endl;
03232   os <<    "Number of Global Blk Entries = "; os << NumGlobalBlockEntries(); os << endl;
03233   os <<    "Global Max Num Block Entries = "; os << GlobalMaxNumBlockEntries(); os << endl;
03234   os <<  "\nNumber of Global Rows        = "; os << NumGlobalRows(); os << endl;
03235   os <<    "Number of Global Cols        = "; os << NumGlobalCols(); os << endl;
03236   os <<    "Number of Global Diagonals   = "; os << NumGlobalDiagonals(); os << endl;
03237   os <<    "Number of Global Nonzeros    = "; os << NumGlobalNonzeros(); os << endl;
03238   os <<    "Global Maximum Num Entries   = "; os << GlobalMaxNumNonzeros(); os << endl;
03239   if (LowerTriangular()) os <<    " ** Matrix is Lower Triangular **"; os << endl;
03240   if (UpperTriangular()) os <<    " ** Matrix is Upper Triangular **"; os << endl;
03241   if (NoDiagonal())      os <<    " ** Matrix has no diagonal     **"; os << endl; os << endl;
03242       }
03243 
03244       os <<  "\nNumber of My Block Rows  = "; os << NumMyBlockRows(); os << endl;
03245       os <<    "Number of My Block Cols  = "; os << NumMyBlockCols(); os << endl;
03246       os <<    "Number of My Block Diags = "; os << NumMyBlockDiagonals(); os << endl;
03247       os <<    "Number of My Blk Entries = "; os << NumMyBlockEntries(); os << endl;
03248       os <<    "My Max Num Block Entries = "; os << MaxNumBlockEntries(); os << endl;
03249       os <<  "\nNumber of My Rows        = "; os << NumMyRows(); os << endl;
03250       os <<    "Number of My Cols        = "; os << NumMyCols(); os << endl;
03251       os <<    "Number of My Diagonals   = "; os << NumMyDiagonals(); os << endl;
03252       os <<    "Number of My Nonzeros    = "; os << NumMyNonzeros(); os << endl;
03253       os <<    "My Maximum Num Entries   = "; os << MaxNumBlockEntries(); os << endl; os << endl;
03254 
03255       os << flush;
03256       
03257     }
03258     // Do a few global ops to give I/O a chance to complete
03259     Comm().Barrier();
03260     Comm().Barrier();
03261     Comm().Barrier();
03262   }
03263 
03264   {for (int iproc=0; iproc < NumProc; iproc++) {
03265     if (MyPID==iproc) {
03266       int NumBlockRows1 = NumMyBlockRows();
03267       int MaxNumBlockEntries1 = MaxNumBlockEntries();
03268       int * BlockIndices1  = new int[MaxNumBlockEntries1];
03269       Epetra_SerialDenseMatrix ** Entries1;
03270       int RowDim1, NumBlockEntries1;
03271       int i, j;
03272 
03273       if (MyPID==0) {
03274   os.width(8);
03275   os <<  "   Processor ";
03276   os.width(10);
03277   os <<  "   Block Row Index ";
03278   os.width(10);
03279   os <<  "   Block Col Index \n";
03280   os.width(20);
03281   os <<  "   Values     ";
03282   os << endl;
03283       }
03284       for (i=0; i<NumBlockRows1; i++) {
03285   int BlockRow1 = GRID(i); // Get global row number
03286   ExtractGlobalBlockRowPointers(BlockRow1, MaxNumBlockEntries1, RowDim1, 
03287               NumBlockEntries1, BlockIndices1,
03288               Entries1);
03289   
03290   for (j = 0; j < NumBlockEntries1 ; j++) {   
03291     os.width(8);
03292     os <<  MyPID ; os << "    ";  
03293     os.width(10);
03294     os <<  BlockRow1 ; os << "    ";  
03295     os.width(10);
03296     os <<  BlockIndices1[j]; os << "    " << endl;
03297     os.width(20);
03298 
03299           if (Entries1[j] == 0) {
03300             os << "Block Entry == NULL"<<endl;
03301             continue;
03302           }
03303 
03304     Epetra_SerialDenseMatrix entry(View, Entries1[j]->A(), Entries1[j]->LDA(),
03305            RowDim1, Entries1[j]->N());
03306     os << entry; os << "    ";
03307     os << endl;
03308   }
03309       }
03310 
03311       delete [] BlockIndices1;
03312       
03313       os << flush;
03314       
03315     }
03316     // Do a few global ops to give I/O a chance to complete
03317     Comm().Barrier();
03318     Comm().Barrier();
03319     Comm().Barrier();
03320   }}
03321 
03322   return;
03323 }

Generated on Thu Sep 18 12:37:58 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1