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