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