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