Epetra Package Browser (Single Doxygen Collection) Development
Epetra_CrsMatrix.cpp
Go to the documentation of this file.
00001 /*
00002 //@HEADER
00003 // ************************************************************************
00004 //
00005 //               Epetra: Linear Algebra Services Package
00006 //                 Copyright 2011 Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ************************************************************************
00041 //@HEADER
00042 */
00043 
00044 #include "Epetra_ConfigDefs.h"
00045 #include "Epetra_CrsMatrix.h"
00046 #include "Epetra_Map.h"
00047 #include "Epetra_Import.h"
00048 #include "Epetra_Export.h"
00049 #include "Epetra_Vector.h"
00050 #include "Epetra_MultiVector.h"
00051 #include "Epetra_Comm.h"
00052 #include "Epetra_SerialComm.h"
00053 #include "Epetra_Distributor.h"
00054 #include "Epetra_OffsetIndex.h"
00055 #include "Epetra_BLAS_wrappers.h"
00056 
00057 #include "Epetra_IntSerialDenseVector.h"
00058 #include "Epetra_CrsGraphData.h"
00059 #include "Epetra_HashTable.h"
00060 #include "Epetra_Util.h"
00061 #include "Epetra_Import_Util.h"
00062 #include "Epetra_IntVector.h"
00063 
00064 #include <cstdlib>
00065 #include <typeinfo>
00066 
00067 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
00068 # include "Teuchos_VerboseObject.hpp"
00069 bool Epetra_CrsMatrixTraceDumpMultiply = false;
00070 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
00071 
00072 #ifdef HAVE_EPETRA_TEUCHOS
00073 // Define this macro to see some timers for some of these functions
00074 # define EPETRA_CRSMATRIX_TEUCHOS_TIMERS
00075 #endif
00076 
00077 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS
00078 #  include "Teuchos_TimeMonitor.hpp"
00079 #endif
00080 
00081 #if defined(Epetra_ENABLE_MKL_SPARSE) && defined(Epetra_ENABLE_CASK)
00082 #error Error: Epetra_ENABLE_MKL_SPARSE and Epetra_ENABLE_CASK both defined.
00083 #endif
00084 
00085 #ifdef Epetra_ENABLE_MKL_SPARSE
00086 #include "mkl_spblas.h"
00087 #endif
00088 
00089 //==============================================================================
00090 Epetra_CrsMatrix::Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_Map& rowMap, const int* NumEntriesPerRow, bool StaticProfile)
00091   : Epetra_DistObject(rowMap, "Epetra::CrsMatrix"),
00092     Epetra_CompObject(),
00093     Epetra_BLAS(),
00094     Graph_(CV, rowMap, NumEntriesPerRow, StaticProfile),
00095     Allocated_(false),
00096     StaticGraph_(false),
00097     UseTranspose_(false),
00098     constructedWithFilledGraph_(false),
00099     matrixFillCompleteCalled_(false),
00100     StorageOptimized_(false),
00101     Values_(0),
00102     Values_alloc_lengths_(0),
00103     All_Values_(0),
00104     NormInf_(0.0),
00105     NormOne_(0.0),
00106     NormFrob_(0.0),
00107     NumMyRows_(rowMap.NumMyPoints()),
00108     ImportVector_(0),
00109     ExportVector_(0),
00110     CV_(CV),
00111     squareFillCompleteCalled_(false)
00112 {
00113   InitializeDefaults();
00114   Allocate();
00115 }
00116 
00117 //==============================================================================
00118 Epetra_CrsMatrix::Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_Map& rowMap, int NumEntriesPerRow, bool StaticProfile)
00119   : Epetra_DistObject(rowMap, "Epetra::CrsMatrix"),
00120     Epetra_CompObject(),
00121     Epetra_BLAS(),
00122     Graph_(CV, rowMap, NumEntriesPerRow, StaticProfile),
00123     Allocated_(false),
00124     StaticGraph_(false),
00125     UseTranspose_(false),
00126     constructedWithFilledGraph_(false),
00127     matrixFillCompleteCalled_(false),
00128     StorageOptimized_(false),
00129     Values_(0),
00130     Values_alloc_lengths_(0),
00131     All_Values_(0),
00132     NormInf_(0.0),
00133     NormOne_(0.0),
00134     NormFrob_(0.0),
00135     NumMyRows_(rowMap.NumMyPoints()),
00136     ImportVector_(0),
00137     ExportVector_(0),
00138     CV_(CV),
00139     squareFillCompleteCalled_(false)
00140 {
00141   InitializeDefaults();
00142   Allocate();
00143 }
00144 //==============================================================================
00145 Epetra_CrsMatrix::Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_Map& rowMap,
00146            const Epetra_Map& colMap, const int* NumEntriesPerRow, bool StaticProfile)
00147   : Epetra_DistObject(rowMap, "Epetra::CrsMatrix"),
00148     Epetra_CompObject(),
00149     Epetra_BLAS(),
00150     Graph_(CV, rowMap, colMap, NumEntriesPerRow, StaticProfile),
00151     Allocated_(false),
00152     StaticGraph_(false),
00153     UseTranspose_(false),
00154     constructedWithFilledGraph_(false),
00155     matrixFillCompleteCalled_(false),
00156     StorageOptimized_(false),
00157     Values_(0),
00158     Values_alloc_lengths_(0),
00159     All_Values_(0),
00160     NormInf_(0.0),
00161     NormOne_(0.0),
00162     NormFrob_(0.0),
00163     NumMyRows_(rowMap.NumMyPoints()),
00164     ImportVector_(0),
00165     ExportVector_(0),
00166     CV_(CV),
00167     squareFillCompleteCalled_(false)
00168 {
00169   InitializeDefaults();
00170   Allocate();
00171 }
00172 
00173 //==============================================================================
00174 Epetra_CrsMatrix::Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_Map& rowMap,
00175            const Epetra_Map& colMap, int NumEntriesPerRow, bool StaticProfile)
00176   : Epetra_DistObject(rowMap, "Epetra::CrsMatrix"),
00177     Epetra_CompObject(),
00178     Epetra_BLAS(),
00179     Graph_(CV, rowMap, colMap,  NumEntriesPerRow, StaticProfile),
00180     Allocated_(false),
00181     StaticGraph_(false),
00182     UseTranspose_(false),
00183     constructedWithFilledGraph_(false),
00184     matrixFillCompleteCalled_(false),
00185     StorageOptimized_(false),
00186     Values_(0),
00187     Values_alloc_lengths_(0),
00188     All_Values_(0),
00189     NormInf_(0.0),
00190     NormOne_(0.0),
00191     NormFrob_(0.0),
00192     NumMyRows_(rowMap.NumMyPoints()),
00193     ImportVector_(0),
00194     ExportVector_(0),
00195     CV_(CV),
00196     squareFillCompleteCalled_(false)
00197 {
00198   if(!rowMap.GlobalIndicesTypeMatch(colMap))
00199     throw ReportError("Epetra_CrsMatrix::Epetra_CrsMatrix: cannot be called with different indices types for rowMap and colMap", -1);
00200 
00201   InitializeDefaults();
00202   Allocate();
00203 }
00204 //==============================================================================
00205 Epetra_CrsMatrix::Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_CrsGraph& graph)
00206   : Epetra_DistObject(graph.Map(), "Epetra::CrsMatrix"),
00207     Epetra_CompObject(),
00208     Epetra_BLAS(),
00209     Graph_(graph),
00210     Allocated_(false),
00211     StaticGraph_(true),
00212     UseTranspose_(false),
00213     constructedWithFilledGraph_(false),
00214     matrixFillCompleteCalled_(false),
00215     StorageOptimized_(false),
00216     Values_(0),
00217     Values_alloc_lengths_(0),
00218     All_Values_(0),
00219     NormInf_(0.0),
00220     NormOne_(0.0),
00221     NormFrob_(0.0),
00222     NumMyRows_(graph.NumMyRows()),
00223     ImportVector_(0),
00224     ExportVector_(0),
00225     CV_(CV),
00226     squareFillCompleteCalled_(false)
00227 {
00228   constructedWithFilledGraph_ = graph.Filled();
00229   InitializeDefaults();
00230   Allocate();
00231 }
00232 
00233 //==============================================================================
00234 Epetra_CrsMatrix::Epetra_CrsMatrix(const Epetra_CrsMatrix& Matrix)
00235   : Epetra_DistObject(Matrix),
00236     Epetra_CompObject(Matrix),
00237     Epetra_BLAS(),
00238     Graph_(Matrix.Graph()),
00239     Allocated_(false),
00240     StaticGraph_(true),
00241     UseTranspose_(Matrix.UseTranspose_),
00242     constructedWithFilledGraph_(false),
00243     matrixFillCompleteCalled_(false),
00244     StorageOptimized_(false),
00245     Values_(0),
00246     Values_alloc_lengths_(0),
00247     All_Values_(0),
00248     NormInf_(0.0),
00249     NormOne_(0.0),
00250     NormFrob_(0.0),
00251     NumMyRows_(Matrix.NumMyRows()),
00252     ImportVector_(0),
00253     ExportVector_(0),
00254     CV_(Copy),
00255     squareFillCompleteCalled_(false)
00256 {
00257   InitializeDefaults();
00258   operator=(Matrix);
00259 }
00260 
00261 //==============================================================================
00262 Epetra_CrsMatrix& Epetra_CrsMatrix::operator=(const Epetra_CrsMatrix& src)
00263 {
00264   if (this == &src) {
00265     return( *this );
00266   }
00267 
00268   if (!src.Filled()) throw ReportError("Copying an Epetra_CrsMatrix requires source matrix to have Filled()==true", -1);
00269 
00270   Graph_ = src.Graph_; // Copy graph
00271 
00272   DeleteMemory();
00273 
00274   StaticGraph_ = true;
00275   UseTranspose_ = src.UseTranspose_;
00276   constructedWithFilledGraph_ = src.constructedWithFilledGraph_;
00277   matrixFillCompleteCalled_ = src.matrixFillCompleteCalled_;
00278   Values_ = 0;
00279   Values_alloc_lengths_ = 0;
00280   All_Values_ = 0;
00281   NormInf_ = -1.0;
00282   NormOne_ = -1.0;
00283   NormFrob_ = -1.0;
00284   NumMyRows_ = src.NumMyRows_;
00285   ImportVector_ = 0;
00286   ExportVector_ = 0;
00287 
00288   CV_ = Copy;
00289 
00290   StorageOptimized_ = src.StorageOptimized_;
00291   if (src.StorageOptimized()) { // Special copy for case where storage is optimized
00292 
00293     int numMyNonzeros = Graph().NumMyEntries();
00294     if (numMyNonzeros>0) All_Values_ = new double[numMyNonzeros];
00295     double * srcValues = src.All_Values();
00296     for (int i=0; i<numMyNonzeros; ++i) All_Values_[i] = srcValues[i];
00297     Allocated_ = true;
00298 #ifdef Epetra_ENABLE_CASK
00299     if( matrixFillCompleteCalled_  ) {
00300       cask = cask_handler_copy(src.cask);
00301     }
00302 #endif
00303   }
00304   else { // copy for non-optimized storage
00305 
00306     Allocate();
00307     for (int i=0; i<NumMyRows_; i++) {
00308       int NumEntries = src.NumMyEntries(i);
00309       double * const srcValues = src.Values(i);
00310       double * targValues = Values(i);
00311       for (int j=0; j< NumEntries; j++) targValues[j] = srcValues[j];
00312     }
00313   }
00314 
00315   return( *this );
00316 }
00317 
00318 //==============================================================================
00319 void Epetra_CrsMatrix::InitializeDefaults() { // Initialize all attributes that have trivial default values
00320 
00321   UseTranspose_ = false;
00322   Values_ = 0;
00323   Values_alloc_lengths_ = 0;
00324   All_Values_ = 0;
00325   NormInf_ = -1.0;
00326   NormOne_ = -1.0;
00327   NormFrob_ = -1.0;
00328   ImportVector_ = 0;
00329   ExportVector_ = 0;
00330 
00331   return;
00332 }
00333 
00334 //==============================================================================
00335 int Epetra_CrsMatrix::Allocate() {
00336 
00337   int i, j;
00338 
00339   // Allocate Values array
00340   Values_ = NumMyRows_ > 0 ? new double*[NumMyRows_] : NULL;
00341   Values_alloc_lengths_ = NumMyRows_ > 0 ? new int[NumMyRows_] : NULL;
00342   if (NumMyRows_ > 0) {
00343     for(j=0; j<NumMyRows_; ++j) Values_alloc_lengths_[j] = 0;
00344   }
00345 
00346   // Allocate and initialize entries if we are copying data
00347   if (CV_==Copy) {
00348     if (Graph().StaticProfile() || Graph().StorageOptimized()) {
00349       int numMyNonzeros = Graph().NumMyEntries();
00350       if (numMyNonzeros>0) All_Values_ = new double[numMyNonzeros];
00351       if(Graph().StorageOptimized()){
00352         StorageOptimized_ = true;
00353       }
00354     }
00355     double * all_values = All_Values_;
00356     for (i=0; i<NumMyRows_; i++) {
00357       int NumAllocatedEntries = Graph().NumAllocatedMyIndices(i);
00358 
00359       if (NumAllocatedEntries > 0) {
00360         if (Graph().StaticProfile() || Graph().StorageOptimized()) {
00361           Values_[i] = all_values;
00362           all_values += NumAllocatedEntries;
00363         }
00364         else {
00365           Values_[i] = new double[NumAllocatedEntries];
00366           Values_alloc_lengths_[i] = NumAllocatedEntries;
00367         }
00368       }
00369       else
00370         Values_[i] = 0;
00371 
00372       for(j=0; j< NumAllocatedEntries; j++)
00373         Values_[i][j] = 0.0; // Fill values with zero
00374     }
00375   }
00376   else {
00377     for (i=0; i<NumMyRows_; i++) {
00378       Values_[i] = 0;
00379     }
00380   }
00381   SetAllocated(true);
00382 #ifdef Epetra_ENABLE_CASK
00383   cask=NULL;
00384 #endif
00385   return(0);
00386 }
00387 //==============================================================================
00388 Epetra_CrsMatrix::~Epetra_CrsMatrix()
00389 {
00390   DeleteMemory();
00391 }
00392 
00393 //==============================================================================
00394 void Epetra_CrsMatrix::DeleteMemory()
00395 {
00396   int i;
00397 
00398   if (CV_==Copy) {
00399     if (All_Values_!=0)
00400       delete [] All_Values_;
00401     else if (Values_!=0)
00402       for (i=0; i<NumMyRows_; i++)
00403         if (Graph().NumAllocatedMyIndices(i) >0)
00404           delete [] Values_[i];
00405   }
00406 
00407   if (ImportVector_!=0)
00408     delete ImportVector_;
00409   ImportVector_=0;
00410 
00411   if (ExportVector_!=0)
00412     delete ExportVector_;
00413   ExportVector_=0;
00414 
00415   delete [] Values_;
00416   Values_ = 0;
00417 
00418   delete [] Values_alloc_lengths_;
00419   Values_alloc_lengths_ = 0;
00420 
00421   NumMyRows_ = 0;
00422 
00423 #ifdef Epetra_ENABLE_CASK
00424   if( StorageOptimized_  )
00425   {
00426     if( cask != NULL )
00427       cask_handler_destroy(cask);
00428 
00429     cask = NULL;
00430   }
00431 #endif
00432 
00433   Allocated_ = false;
00434 }
00435 
00436 //==============================================================================
00437 int Epetra_CrsMatrix::ReplaceRowMap(const Epetra_BlockMap& newmap)
00438 {
00439   int err = Graph_.ReplaceRowMap(newmap);
00440   if (err == 0) {
00441     //update export vector.
00442 
00443     if (ExportVector_ != 0) {
00444       delete ExportVector_;
00445       ExportVector_= 0;
00446     }
00447 
00448     ExportVector_ = new Epetra_MultiVector(RowMap(),1);
00449   }
00450   return(err);
00451 }
00452 
00453 //==============================================================================
00454 int Epetra_CrsMatrix::ReplaceColMap(const Epetra_BlockMap& newmap)
00455 {
00456   int err = Graph_.ReplaceColMap(newmap);
00457   if (err == 0) {
00458     //update import vector.
00459 
00460     if (ImportVector_ != 0) {
00461       delete ImportVector_;
00462       ImportVector_= 0;
00463     }
00464 
00465     ImportVector_ = new Epetra_MultiVector(ColMap(),1);
00466   }
00467   return(err);
00468 }
00469 
00470 //==============================================================================
00471 int Epetra_CrsMatrix::ReplaceDomainMapAndImporter(const Epetra_Map& NewDomainMap, const Epetra_Import * NewImporter) {
00472   return Graph_.ReplaceDomainMapAndImporter(NewDomainMap,NewImporter);
00473 }
00474 
00475 //==============================================================================
00476 int Epetra_CrsMatrix::RemoveEmptyProcessesInPlace(const Epetra_BlockMap * NewMap) {
00477   // Epetra_DistObject things
00478   if(NewMap) {
00479     Map_  = *NewMap;
00480     Comm_ = &NewMap->Comm();
00481   }
00482 
00483   return Graph_.RemoveEmptyProcessesInPlace(NewMap);
00484 }
00485 
00486 //==============================================================================
00487 int Epetra_CrsMatrix::PutScalar(double ScalarConstant)
00488 {
00489   if (StorageOptimized()) {
00490     int length = NumMyNonzeros();
00491     for (int i=0; i<length; ++i) All_Values_[i] = ScalarConstant;
00492   }
00493   else {
00494     for(int i=0; i<NumMyRows_; i++) {
00495       int NumEntries = Graph().NumMyIndices(i);
00496       double * targValues = Values(i);
00497       for(int j=0; j< NumEntries; j++)
00498   targValues[j] = ScalarConstant;
00499     }
00500   }
00501   return(0);
00502 }
00503 //==============================================================================
00504 int Epetra_CrsMatrix::Scale(double ScalarConstant)
00505 {
00506   if (StorageOptimized()) {
00507     int length = NumMyNonzeros();
00508     for (int i=0; i<length; ++i) All_Values_[i] *= ScalarConstant;
00509   }
00510   else {
00511     for(int i=0; i<NumMyRows_; i++) {
00512       int NumEntries = Graph().NumMyIndices(i);
00513       double * targValues = Values(i);
00514       for(int j=0; j< NumEntries; j++)
00515   targValues[j] *= ScalarConstant;
00516     }
00517   }
00518   return(0);
00519 }
00520 
00521 //==========================================================================
00522 template<typename int_type>
00523 int Epetra_CrsMatrix::TInsertGlobalValues(int_type Row, int NumEntries,
00524            const double* values,
00525            const int_type* Indices)
00526 {
00527   if(IndicesAreLocal())
00528     EPETRA_CHK_ERR(-2); // Cannot insert global values into local graph
00529   if(IndicesAreContiguous())
00530     EPETRA_CHK_ERR(-3); // Indices cannot be individually deleted and newed
00531   Graph_.SetIndicesAreGlobal(true);
00532   int locRow = Graph_.LRID(Row); // Find local row number for this global row index
00533 
00534   EPETRA_CHK_ERR( InsertValues(locRow, NumEntries, values, Indices) );
00535 
00536   return(0);
00537 }
00538 
00539 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00540 int Epetra_CrsMatrix::InsertGlobalValues(int Row, int NumEntries,
00541            const double* values,
00542            const int* Indices)
00543 {
00544   if(RowMap().GlobalIndicesInt())
00545     return TInsertGlobalValues<int>(Row, NumEntries, values, Indices);
00546   else
00547     throw ReportError("Epetra_CrsMatrix::InsertGlobalValues int version called for a matrix that is not int.", -1);
00548 }
00549 #endif
00550 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00551 int Epetra_CrsMatrix::InsertGlobalValues(long long Row, int NumEntries,
00552            const double* values,
00553            const long long* Indices)
00554 {
00555   if(RowMap().GlobalIndicesLongLong())
00556     return TInsertGlobalValues<long long>(Row, NumEntries, values, Indices);
00557   else
00558     throw ReportError("Epetra_CrsMatrix::InsertGlobalValues long long version called for a matrix that is not long long.", -1);
00559 }
00560 #endif
00561 
00562 //==========================================================================
00563 template<typename int_type>
00564 int Epetra_CrsMatrix::TInsertGlobalValues(int_type Row, int NumEntries,
00565            double* values,
00566            int_type* Indices)
00567 {
00568   if(IndicesAreLocal())
00569     EPETRA_CHK_ERR(-2); // Cannot insert global values into local graph
00570   if(IndicesAreContiguous())
00571     EPETRA_CHK_ERR(-3); // Indices cannot be individually deleted and newed
00572   Graph_.SetIndicesAreGlobal(true);
00573   int locRow = Graph_.LRID(Row); // Find local row number for this global row index
00574 
00575   EPETRA_CHK_ERR( InsertValues(locRow, NumEntries, values, Indices) );
00576 
00577   return(0);
00578 }
00579 
00580 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00581 int Epetra_CrsMatrix::InsertGlobalValues(int Row, int NumEntries,
00582            double* values,
00583            int* Indices)
00584 {
00585   if(RowMap().GlobalIndicesInt())
00586     return TInsertGlobalValues<int>(Row, NumEntries, values, Indices);
00587   else
00588     throw ReportError("Epetra_CrsMatrix::InsertGlobalValues int version called for a matrix that is not int.", -1);
00589 }
00590 #endif
00591 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00592 int Epetra_CrsMatrix::InsertGlobalValues(long long Row, int NumEntries,
00593            double* values,
00594            long long* Indices)
00595 {
00596   if(RowMap().GlobalIndicesLongLong()) {
00597     return TInsertGlobalValues<long long>(Row, NumEntries, values, Indices);
00598   }
00599   else
00600     throw ReportError("Epetra_CrsMatrix::InsertGlobalValues long long version called for a matrix that is not long long.", -1);
00601 }
00602 #endif
00603 //==========================================================================
00604 int Epetra_CrsMatrix::InsertMyValues(int Row, int NumEntries,
00605              const double* values,
00606              const int* Indices)
00607 {
00608   if(IndicesAreGlobal())
00609     EPETRA_CHK_ERR(-2); // Cannot insert global values into filled graph
00610   if(IndicesAreContiguous() && CV_==Copy)
00611     EPETRA_CHK_ERR(-3); // Indices cannot be individually deleted and new
00612   Graph_.SetIndicesAreLocal(true);
00613 
00614 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
00615   EPETRA_CHK_ERR( InsertValues(Row, NumEntries, values, Indices) );
00616 #else
00617   throw ReportError("Epetra_CrsMatrix::InsertMyValues: Failure because neither 32 bit nor 64 bit indices insertable.", -1);
00618 #endif
00619 
00620   return(0);
00621 
00622 }
00623 
00624 //==========================================================================
00625 int Epetra_CrsMatrix::InsertMyValues(int Row, int NumEntries,
00626              double* values,
00627              int* Indices)
00628 {
00629   if(IndicesAreGlobal())
00630     EPETRA_CHK_ERR(-2); // Cannot insert global values into filled graph
00631   if(IndicesAreContiguous() && CV_==Copy)
00632     EPETRA_CHK_ERR(-3); // Indices cannot be individually deleted and new
00633   Graph_.SetIndicesAreLocal(true);
00634 
00635 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
00636   EPETRA_CHK_ERR( InsertValues(Row, NumEntries, values, Indices) );
00637 #else
00638     throw ReportError("Epetra_CrsMatrix::InsertMyValues: Failure because neither 32 bit nor 64 bit indices insertable.", -1);
00639 #endif
00640 
00641   return(0);
00642 
00643 }
00644 
00645 //==========================================================================
00646 template<typename int_type>
00647 int Epetra_CrsMatrix::InsertValues(int Row, int NumEntries,
00648            const double* values,
00649            const int_type* Indices)
00650 {
00651   if(CV_ == View){
00652     //cannot allow View mode with const pointers
00653     EPETRA_CHK_ERR(-4);
00654   }
00655   else{
00656     //to avoid code duplication I am using a cheap tactic of removing the constness of
00657     //values and Indices. Since this is only called in copy mode the passed in variables
00658     //will not be modified.
00659     return(InsertValues(Row, NumEntries, const_cast<double*>(values), const_cast<int_type*>(Indices)));
00660   }
00661   return 0;
00662 }
00663 
00664 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00665 int Epetra_CrsMatrix::InsertValues(int Row, int NumEntries,
00666            const double* values,
00667            const int* Indices)
00668 {
00669   if(RowMap().GlobalIndicesInt() || (RowMap().GlobalIndicesLongLong() && IndicesAreLocal()))
00670     return InsertValues<int>(Row, NumEntries, values, Indices);
00671   else
00672     throw ReportError("Epetra_CrsMatrix::InsertValues int version called for a matrix that is not int.", -1);
00673 }
00674 #endif
00675 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00676 int Epetra_CrsMatrix::InsertValues(int Row, int NumEntries,
00677            const double* values,
00678            const long long* Indices)
00679 {
00680   if(RowMap().GlobalIndicesLongLong())
00681     return InsertValues<long long>(Row, NumEntries, values, Indices);
00682   else
00683     throw ReportError("Epetra_CrsMatrix::InsertValues long long version called for a matrix that is not long long.", -1);
00684 }
00685 #endif
00686 //==========================================================================
00687 template<typename int_type>
00688 int Epetra_CrsMatrix::InsertValues(int Row, int NumEntries,
00689            double* values,
00690            int_type* Indices)
00691 {
00692   int j;
00693   int ierr = 0;
00694 
00695   if(Row < 0 || Row >= NumMyRows_)
00696     EPETRA_CHK_ERR(-1); // Not in Row range
00697 
00698   if(CV_ == View) {
00699     //test indices in static graph
00700     if(StaticGraph()) {
00701       int testNumEntries;
00702       int* testIndices;
00703       int testRow = Row;
00704       if(IndicesAreGlobal())
00705         testRow = Graph_.LRID( Row );
00706       EPETRA_CHK_ERR(Graph_.ExtractMyRowView(testRow, testNumEntries, testIndices));
00707 
00708       bool match = true;
00709       if(NumEntries != testNumEntries)
00710         match = false;
00711       for(int i = 0; i < NumEntries; ++i)
00712         match = match && (Indices[i]==testIndices[i]);
00713 
00714       if(!match)
00715         ierr = -3;
00716     }
00717 
00718     if(Values_[Row] != 0)
00719       ierr = 2; // This row has been defined already.  Issue warning.
00720     Values_[Row] = values;
00721   }
00722   else {
00723     if(StaticGraph())
00724       EPETRA_CHK_ERR(-2); // If the matrix graph is fully constructed, we cannot insert new values
00725 
00726     int tmpNumEntries = NumEntries;
00727 
00728     if(Graph_.HaveColMap()) { //must insert only valid indices, values
00729       const double* tmpValues = values;
00730       values = new double[NumEntries];
00731       int loc = 0;
00732       if(IndicesAreLocal()) {
00733         for(int i = 0; i < NumEntries; ++i)
00734           if(Graph_.ColMap().MyLID(static_cast<int>(Indices[i])))
00735             values[loc++] = tmpValues[i];
00736       }
00737       else {
00738         for(int i = 0; i < NumEntries; ++i)
00739           if(Graph_.ColMap().MyGID(Indices[i]))
00740             values[loc++] = tmpValues[i];
00741       }
00742       if(NumEntries != loc)
00743         ierr = 2; //Some columns excluded
00744       NumEntries = loc;
00745     }
00746 
00747     int start = Graph().NumMyIndices(Row);
00748     int stop = start + NumEntries;
00749     int NumAllocatedEntries = Values_alloc_lengths_[Row];
00750     if(stop > NumAllocatedEntries) {
00751       if (Graph().StaticProfile() && stop > Graph().NumAllocatedMyIndices(Row)) {
00752         EPETRA_CHK_ERR(-2); // Cannot expand graph storage if graph created using StaticProfile
00753       }
00754       if(NumAllocatedEntries == 0) {
00755         Values_[Row] = new double[NumEntries]; // Row was never allocated, so do it
00756         Values_alloc_lengths_[Row] = NumEntries;
00757       }
00758       else {
00759         ierr = 1; // Out of room.  Must delete and allocate more space...
00760         double* tmp_Values = new double[stop];
00761         for(j = 0; j < start; j++)
00762           tmp_Values[j] = Values_[Row][j]; // Copy existing entries
00763         delete[] Values_[Row]; // Delete old storage
00764         Values_[Row] = tmp_Values; // Set pointer to new storage
00765         Values_alloc_lengths_[Row] = stop;
00766       }
00767     }
00768 
00769     for(j = start; j < stop; j++)
00770       Values_[Row][j] = values[j-start];
00771 
00772 
00773     NumEntries = tmpNumEntries;
00774     if(Graph_.HaveColMap())
00775       delete[] values;
00776   }
00777 
00778   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
00779   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
00780   NormFrob_ = -1.0;
00781 
00782   if(!StaticGraph()) {
00783     EPETRA_CHK_ERR(Graph_.InsertIndices(Row, NumEntries, Indices));
00784   }
00785 
00786   EPETRA_CHK_ERR(ierr);
00787   return(0);
00788 }
00789 
00790 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00791 int Epetra_CrsMatrix::InsertValues(int Row, int NumEntries,
00792            double* values,
00793            int* Indices)
00794 {
00795   if(RowMap().GlobalIndicesInt() || (RowMap().GlobalIndicesLongLong() && IndicesAreLocal()))
00796   return InsertValues<int>(Row, NumEntries, values, Indices);
00797   else
00798   throw ReportError("Epetra_CrsMatrix::InsertValues int version called for a matrix that is not int.", -1);
00799 }
00800 #endif
00801 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00802 int Epetra_CrsMatrix::InsertValues(int Row, int NumEntries,
00803            double* values,
00804            long long* Indices)
00805 {
00806   if(RowMap().GlobalIndicesLongLong())
00807   return InsertValues<long long>(Row, NumEntries, values, Indices);
00808   else
00809   throw ReportError("Epetra_CrsMatrix::InsertValues long long version called for a matrix that is not long long.", -1);
00810 }
00811 #endif
00812 
00813 //==========================================================================
00814 int Epetra_CrsMatrix::InsertOffsetValues(long long Row, int NumEntries,
00815            double* values,
00816            int* Indices)
00817 {
00818   return ReplaceOffsetValues(Row, NumEntries, values, Indices);
00819 }
00820 
00821 //==========================================================================
00822 template<typename int_type>
00823 int Epetra_CrsMatrix::TReplaceGlobalValues(int_type Row, int NumEntries, const double * srcValues, const int_type *Indices) {
00824 
00825   int j;
00826   int ierr = 0;
00827   int Loc;
00828 
00829   int locRow = Graph_.LRID(Row); // Normalize row range
00830 
00831   if (locRow < 0 || locRow >= NumMyRows_) {
00832     EPETRA_CHK_ERR(-1); // Not in Row range
00833   }
00834   double * targValues = Values(locRow);
00835   for (j=0; j<NumEntries; j++) {
00836     int_type Index = Indices[j];
00837     if (Graph_.FindGlobalIndexLoc(locRow,Index,j,Loc))
00838       targValues[Loc] = srcValues[j];
00839     else
00840       ierr = 2; // Value Excluded
00841   }
00842 
00843   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
00844   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
00845   NormFrob_ = -1.0;
00846 
00847   EPETRA_CHK_ERR(ierr);
00848   return(0);
00849 }
00850 
00851 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00852 int Epetra_CrsMatrix::ReplaceGlobalValues(int Row, int NumEntries, const double * srcValues, const int *Indices) {
00853   if(RowMap().GlobalIndicesInt())
00854   return TReplaceGlobalValues<int>(Row, NumEntries, srcValues, Indices);
00855   else
00856   throw ReportError("Epetra_CrsMatrix::ReplaceGlobalValues int version called for a matrix that is not int.", -1);
00857 }
00858 #endif
00859 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00860 int Epetra_CrsMatrix::ReplaceGlobalValues(long long Row, int NumEntries, const double * srcValues, const long long *Indices) {
00861   if(RowMap().GlobalIndicesLongLong())
00862   return TReplaceGlobalValues<long long>(Row, NumEntries, srcValues, Indices);
00863   else
00864   throw ReportError("Epetra_CrsMatrix::ReplaceGlobalValues long long version called for a matrix that is not long long.", -1);
00865 }
00866 #endif
00867 
00868 //==========================================================================
00869 int Epetra_CrsMatrix::ReplaceMyValues(int Row, int NumEntries, const double * srcValues, const int *Indices) {
00870 
00871   if (!IndicesAreLocal())
00872     EPETRA_CHK_ERR(-4); // Indices must be local.
00873 
00874   int j;
00875   int ierr = 0;
00876   int Loc;
00877 
00878   if (Row < 0 || Row >= NumMyRows_) {
00879     EPETRA_CHK_ERR(-1); // Not in Row range
00880   }
00881 
00882   double* RowValues = Values(Row);
00883   for (j=0; j<NumEntries; j++) {
00884     int Index = Indices[j];
00885     if (Graph_.FindMyIndexLoc(Row,Index,j,Loc))
00886       RowValues[Loc] = srcValues[j];
00887     else
00888       ierr = 2; // Value Excluded
00889   }
00890 
00891   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
00892   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
00893   NormFrob_ = -1.0;
00894 
00895   EPETRA_CHK_ERR(ierr);
00896   return(0);
00897 }
00898 
00899 //==========================================================================
00900 int Epetra_CrsMatrix::ReplaceOffsetValues(long long Row, int NumEntries,
00901             const double * srcValues, const int *Offsets)
00902 {
00903   int j;
00904   int ierr = 0;
00905 
00906   // Normalize row range
00907 #ifdef EPETRA_NO_64BIT_GLOBAL_INDICES
00908   int locRow = LRID((int) Row);
00909 #else
00910   int locRow = LRID(Row);
00911 #endif
00912 
00913   if (locRow < 0 || locRow >= NumMyRows_) {
00914     EPETRA_CHK_ERR(-1); // Not in Row range
00915   }
00916 
00917   double* RowValues = Values(locRow);
00918   for(j=0; j<NumEntries; j++) {
00919     if( Offsets[j] != -1 )
00920       RowValues[Offsets[j]] = srcValues[j];
00921   }
00922 
00923   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
00924   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
00925   NormFrob_ = -1.0;
00926 
00927   EPETRA_CHK_ERR(ierr);
00928   return(0);
00929 }
00930 
00931 //==========================================================================
00932 template<typename int_type>
00933 int Epetra_CrsMatrix::TSumIntoGlobalValues(int_type Row,
00934             int NumEntries,
00935             const double * srcValues,
00936             const int_type *Indices)
00937 {
00938   int j;
00939   int ierr = 0;
00940   int Loc = 0;
00941 
00942 
00943   int locRow = Graph_.LRID(Row); // Normalize row range
00944 
00945   if (locRow < 0 || locRow >= NumMyRows_) {
00946     EPETRA_CHK_ERR(-1); // Not in Row range
00947   }
00948 
00949   if (StaticGraph() && !Graph_.HaveColMap()) {
00950     EPETRA_CHK_ERR(-1);
00951   }
00952 
00953   double * RowValues = Values(locRow);
00954 
00955   if (!StaticGraph()) {
00956     for (j=0; j<NumEntries; j++) {
00957       int_type Index = Indices[j];
00958       if (Graph_.FindGlobalIndexLoc(locRow,Index,j,Loc))
00959 #ifdef EPETRA_HAVE_OMP
00960 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE
00961 #pragma omp atomic
00962 #endif
00963 #endif
00964         RowValues[Loc] += srcValues[j];
00965       else
00966         ierr = 2; // Value Excluded
00967     }
00968   }
00969   else {
00970     const Epetra_BlockMap& colmap = Graph_.ColMap();
00971     int NumColIndices = Graph_.NumMyIndices(locRow);
00972     const int* ColIndices = Graph_.Indices(locRow);
00973 
00974     if (Graph_.Sorted()) {
00975       int insertPoint;
00976       for (j=0; j<NumEntries; j++) {
00977         int Index = colmap.LID(Indices[j]);
00978 
00979         // Check whether the next added element is the subsequent element in
00980         // the graph indices, then we can skip the binary search
00981         if (Loc < NumColIndices && Index == ColIndices[Loc])
00982 #ifdef EPETRA_HAVE_OMP
00983 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE
00984 #pragma omp atomic
00985 #endif
00986 #endif
00987           RowValues[Loc] += srcValues[j];
00988         else {
00989           Loc = Epetra_Util_binary_search(Index, ColIndices, NumColIndices, insertPoint);
00990           if (Loc > -1)
00991 #ifdef EPETRA_HAVE_OMP
00992 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE
00993 #pragma omp atomic
00994 #endif
00995 #endif
00996             RowValues[Loc] += srcValues[j];
00997           else
00998             ierr = 2; // Value Excluded
00999         }
01000         ++Loc;
01001       }
01002     }
01003     else
01004       for (j=0; j<NumEntries; j++) {
01005         int Index = colmap.LID(Indices[j]);
01006         if (Graph_.FindMyIndexLoc(NumColIndices,ColIndices,Index,j,Loc))
01007 #ifdef EPETRA_HAVE_OMP
01008 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE
01009 #pragma omp atomic
01010 #endif
01011 #endif
01012           RowValues[Loc] += srcValues[j];
01013         else
01014           ierr = 2; // Value Excluded
01015       }
01016   }
01017 
01018   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
01019   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
01020   NormFrob_ = -1.0;
01021 
01022   EPETRA_CHK_ERR(ierr);
01023 
01024   return(0);
01025 }
01026 
01027 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01028 int Epetra_CrsMatrix::SumIntoGlobalValues(int Row,
01029             int NumEntries,
01030             const double * srcValues,
01031             const int *Indices)
01032 {
01033   if(RowMap().GlobalIndicesInt())
01034   return TSumIntoGlobalValues<int>(Row, NumEntries, srcValues, Indices);
01035   else
01036   throw ReportError("Epetra_CrsMatrix::SumIntoGlobalValues int version called for a matrix that is not int.", -1);
01037 }
01038 #endif
01039 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01040 int Epetra_CrsMatrix::SumIntoGlobalValues(long long Row,
01041             int NumEntries,
01042             const double * srcValues,
01043             const long long *Indices)
01044 {
01045   if(RowMap().GlobalIndicesLongLong())
01046   return TSumIntoGlobalValues<long long>(Row, NumEntries, srcValues, Indices);
01047   else
01048   throw ReportError("Epetra_CrsMatrix::SumIntoGlobalValues long long version called for a matrix that is not long long.", -1);
01049 }
01050 #endif
01051 
01052 //==========================================================================
01053 int Epetra_CrsMatrix::SumIntoMyValues(int Row, int NumEntries, const double * srcValues, const int *Indices) {
01054 
01055   if (!IndicesAreLocal())
01056     EPETRA_CHK_ERR(-4); // Indices must be local.
01057 
01058   int j;
01059   int ierr = 0;
01060   int Loc = 0;
01061   int insertPoint;
01062 
01063   if (Row < 0 || Row >= NumMyRows_) {
01064     EPETRA_CHK_ERR(-1); // Not in Row range
01065   }
01066 
01067   double* RowValues = Values(Row);
01068   int NumColIndices = Graph_.NumMyIndices(Row);
01069   const int* ColIndices = Graph_.Indices(Row);
01070   if (Graph_.Sorted()) {
01071     for (j=0; j<NumEntries; j++) {
01072       int Index = Indices[j];
01073 
01074       // Check whether the next added element is the subsequent element in
01075       // the graph indices.
01076       if (Loc < NumColIndices && Index == ColIndices[Loc])
01077         RowValues[Loc] += srcValues[j];
01078       else {
01079         Loc = Epetra_Util_binary_search(Index, ColIndices, NumColIndices, insertPoint);
01080         if (Loc > -1)
01081           RowValues[Loc] += srcValues[j];
01082         else
01083           ierr = 2; // Value Excluded
01084       }
01085       ++Loc;
01086     }
01087   }
01088   else {
01089     for (j=0; j<NumEntries; j++) {
01090       int Index = Indices[j];
01091       if (Graph_.FindMyIndexLoc(Row,Index,j,Loc))
01092         RowValues[Loc] += srcValues[j];
01093       else
01094         ierr = 2; // Value Excluded
01095     }
01096   }
01097 
01098   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
01099   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
01100   NormFrob_ = -1.0;
01101 
01102   EPETRA_CHK_ERR(ierr);
01103 
01104   return(0);
01105 }
01106 
01107 //==========================================================================
01108 int Epetra_CrsMatrix::SumIntoOffsetValues(long long Row, int NumEntries, const double * srcValues, const int *Offsets) {
01109 
01110   int j;
01111   int ierr = 0;
01112 
01113   // Normalize row range
01114 #ifdef EPETRA_NO_64BIT_GLOBAL_INDICES
01115   int locRow = LRID((int) Row);
01116 #else
01117   int locRow = LRID(Row);
01118 #endif
01119 
01120   if (locRow < 0 || locRow >= NumMyRows_) {
01121     EPETRA_CHK_ERR(-1); // Not in Row range
01122   }
01123 
01124   double* RowValues = Values(locRow);
01125   for (j=0; j<NumEntries; j++) {
01126     if( Offsets[j] != -1 )
01127       RowValues[Offsets[j]] += srcValues[j];
01128   }
01129 
01130   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
01131   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
01132   NormFrob_ = -1.0;
01133 
01134   EPETRA_CHK_ERR(ierr);
01135 
01136   return(0);
01137 }
01138 
01139 //==========================================================================
01140 int Epetra_CrsMatrix::FillComplete(bool OptimizeDataStorage) {
01141   squareFillCompleteCalled_ = true;
01142   EPETRA_CHK_ERR(FillComplete(RowMap(), RowMap(), OptimizeDataStorage));
01143   return(0);
01144 }
01145 
01146 //==========================================================================
01147 int Epetra_CrsMatrix::FillComplete(const Epetra_Map& domain_map,
01148            const Epetra_Map& range_map, bool OptimizeDataStorage)
01149 {
01150   int returnValue = 0;
01151 
01152   if (Graph_.Filled()) {
01153     if (!constructedWithFilledGraph_ && !matrixFillCompleteCalled_) {
01154       returnValue = 2;
01155     }
01156   }
01157 
01158   if (!StaticGraph()) {
01159     if (Graph_.MakeIndicesLocal(domain_map, range_map) < 0) {
01160       return(-1);
01161     }
01162   }
01163   SortEntries();  // Sort column entries from smallest to largest
01164   MergeRedundantEntries(); // Get rid of any redundant index values
01165   if (!StaticGraph()) {
01166     if (Graph_.FillComplete(domain_map, range_map) < 0) {
01167       return(-2);
01168     }
01169   }
01170 
01171   matrixFillCompleteCalled_ = true;
01172 
01173   if (squareFillCompleteCalled_) {
01174     if (DomainMap().NumGlobalElements64() != RangeMap().NumGlobalElements64()) {
01175       returnValue = 3;
01176     }
01177     squareFillCompleteCalled_ = false;
01178     EPETRA_CHK_ERR(returnValue);
01179   }
01180 
01181   if (OptimizeDataStorage) { EPETRA_CHK_ERR(OptimizeStorage()); }
01182 
01183   return(returnValue);
01184 }
01185 
01186 //==========================================================================
01187 int Epetra_CrsMatrix::TransformToLocal() {
01188   EPETRA_CHK_ERR(FillComplete());
01189   return(0);
01190 }
01191 
01192 //==========================================================================
01193 int Epetra_CrsMatrix::TransformToLocal(const Epetra_Map* domainMap, const Epetra_Map* rangeMap) {
01194   EPETRA_CHK_ERR(FillComplete(*domainMap, *rangeMap));
01195   return(0);
01196 }
01197 
01198 //==========================================================================
01199 int Epetra_CrsMatrix::SortEntries() {
01200 
01201   if(!IndicesAreLocal())
01202     EPETRA_CHK_ERR(-1);
01203   if(Sorted())
01204     return(0);
01205 
01206   // For each row, sort column entries from smallest to largest.
01207   // Use shell sort. Stable sort so it is fast if indices are already sorted.
01208 
01209 
01210   for(int i = 0; i < NumMyRows_; i++){
01211 
01212     double* locValues = Values(i);
01213     int NumEntries = Graph().NumMyIndices(i);
01214     int* locIndices = Graph().Indices(i);
01215 
01216     int n = NumEntries;
01217     int m = n/2;
01218 
01219     while(m > 0) {
01220       int max = n - m;
01221       for(int j = 0; j < max; j++) {
01222   for(int k = j; k >= 0; k-=m) {
01223     if(locIndices[k+m] >= locIndices[k])
01224       break;
01225     double dtemp = locValues[k+m];
01226     locValues[k+m] = locValues[k];
01227     locValues[k] = dtemp;
01228     int itemp = locIndices[k+m];
01229     locIndices[k+m] = locIndices[k];
01230     locIndices[k] = itemp;
01231   }
01232       }
01233       m = m/2;
01234     }
01235   }
01236   Graph_.SetSorted(true); // This also sorted the graph
01237   return(0);
01238 }
01239 
01240 //==========================================================================
01241 int Epetra_CrsMatrix::MergeRedundantEntries() {
01242 
01243   int i;
01244 
01245   if(NoRedundancies())
01246     return(0);
01247   if(!Sorted())
01248     EPETRA_CHK_ERR(-1);  // Must have sorted entries
01249 
01250   // For each row, remove column indices that are repeated.
01251   // Also, determine if matrix is upper or lower triangular or has no diagonal (Done in graph)
01252   // Note:  This function assumes that SortEntries was already called.
01253 
01254   for(i = 0; i<NumMyRows_; i++) {
01255     int NumEntries = Graph().NumMyIndices(i);
01256     if(NumEntries > 1) {
01257       double* const locValues = Values(i);
01258       int* const locIndices = Graph().Indices(i);
01259       int curEntry =0;
01260       double curValue = locValues[0];
01261       for(int k = 1; k < NumEntries; k++) {
01262   if(locIndices[k] == locIndices[k-1])
01263     curValue += locValues[k];
01264   else {
01265     locValues[curEntry++] = curValue;
01266     curValue = locValues[k];
01267   }
01268       }
01269       locValues[curEntry] = curValue;
01270 
01271     }
01272   }
01273 
01274   EPETRA_CHK_ERR(Graph_.RemoveRedundantIndices()); // Remove redundant indices and then return
01275   return(0);
01276 }
01277 
01278 //==========================================================================
01279 int Epetra_CrsMatrix::OptimizeStorage() {
01280 
01281 
01282   if (StorageOptimized())
01283     return(0); // Have we been here before?
01284   if (!Filled()) EPETRA_CHK_ERR(-1); // Cannot optimize storage before calling FillComplete()
01285 
01286 
01287   int ierr = Graph_.OptimizeStorage();
01288   if (ierr!=0) EPETRA_CHK_ERR(ierr);  // In order for OptimizeStorage to make sense for the matrix, it must work on the graph.
01289 
01290   bool Contiguous = true; // Assume contiguous is true
01291   for (int i=1; i<NumMyRows_; i++){
01292     int NumEntries = Graph().NumMyIndices(i-1);
01293 
01294     // check if end of beginning of current row starts immediately after end of previous row.
01295     if (Values_[i]!=Values_[i-1]+NumEntries) {
01296       Contiguous = false;
01297       break;
01298     }
01299   }
01300 
01301   // NOTE:  At the end of the above loop set, there is a possibility that NumEntries and NumAllocatedEntries
01302   //        for the last row could be different, but I don't think it matters.
01303 
01304 
01305   if ((CV_==View) && !Contiguous)
01306     EPETRA_CHK_ERR(-1);  // This is user data, it's not contiguous and we can't make it so.
01307 
01308 
01309   if(!Contiguous) { // Must pack indices if not already contiguous.  Pack values into All_values_
01310 
01311     if (All_Values_==0) {
01312       // Compute Number of Nonzero entries (Done in FillComplete, but we may not have been there yet.)
01313       int numMyNonzeros = Graph_.NumMyNonzeros();
01314 
01315       // Allocate one big array for all values
01316       All_Values_ = new double[numMyNonzeros];
01317       if(All_Values_ == 0) throw ReportError("Error with All_Values_ allocation.", -99);
01318 
01319       const int *const  IndexOffset = Graph().IndexOffset();
01320       const int numMyRows = NumMyRows_;
01321       double ** Values_s = Values_;
01322       double * All_Values_s = All_Values_;
01323 #ifdef EPETRA_HAVE_OMP
01324 #pragma omp parallel for default(none) shared(Values_s,All_Values_s)
01325 #endif
01326       for (int i=0; i<numMyRows; i++) {
01327         int NumEntries = Graph().NumMyIndices(i);
01328         int curOffset = IndexOffset[i];
01329         double * values = Values_s[i];
01330         double * newValues = All_Values_s+curOffset;
01331   for (int j=0; j<NumEntries; j++) newValues[j] = values[j];
01332       }
01333     }
01334     else { // Static Profile, so just pack into existing storage (can't be threaded)
01335       double * tmp = All_Values_;
01336       for (int i=0; i<NumMyRows_; i++) {
01337         int NumEntries = Graph().NumMyIndices(i);
01338         double * values = Values_[i];
01339         if (tmp!=values) // Copy values if not pointing to same location
01340           for (int j=0; j<NumEntries; j++) tmp[j] = values[j];
01341         tmp += NumEntries;
01342       }
01343     }
01344 
01345     // Free Values_ arrays
01346     for (int i=0;i<NumMyRows_; ++i) {
01347       if (Values_alloc_lengths_[i] != 0) delete [] Values_[i];
01348     }
01349 
01350     delete [] Values_alloc_lengths_; Values_alloc_lengths_ = 0;
01351   } // End of !Contiguous section
01352   else {
01353     //if already contiguous, we'll simply set All_Values_ to be
01354     //a copy of Values_[0].
01355     if (All_Values_!=0 && All_Values_!=Values_[0]) delete [] All_Values_;
01356     All_Values_ = NumMyRows_ > 0 ? Values_[0] : NULL;
01357   }
01358 
01359   // Delete unneeded storage
01360   delete [] Values_; Values_=0;
01361 
01362 #ifdef Epetra_ENABLE_CASK
01363   if (cask == NULL  && Graph().StorageOptimized() )  {
01364      int * Indices = Graph().All_Indices();
01365      int * IndexOffset = Graph().IndexOffset();
01366      int NumMyCols_ = NumMyCols();
01367      cask_handler_initialize(&cask);
01368      cask_csr_analysis(NumMyRows_, NumMyCols_, IndexOffset, Indices,
01369                        NumGlobalNonzeros64(),cask);
01370   }
01371 #endif
01372 
01373 
01374   StorageOptimized_ = true;
01375 
01376 
01377   return(0);
01378 }
01379 
01380 //==========================================================================
01381 template<typename int_type>
01382 int Epetra_CrsMatrix::ExtractGlobalRowCopy(int_type Row, int Length, int & NumEntries, double * values,
01383              int_type * Indices) const
01384 {
01385 
01386   int ierr = Graph_.ExtractGlobalRowCopy(Row, Length, NumEntries, Indices);
01387   if (ierr)
01388     EPETRA_CHK_ERR(ierr);
01389 
01390   EPETRA_CHK_ERR(ExtractGlobalRowCopy(Row, Length, NumEntries, values));
01391   return(0);
01392 }
01393 
01394 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01395 int Epetra_CrsMatrix::ExtractGlobalRowCopy(int Row, int Length, int & NumEntries, double * values,
01396              int * Indices) const
01397 {
01398   if(RowMap().GlobalIndicesInt())
01399     return ExtractGlobalRowCopy<int>(Row, Length, NumEntries, values, Indices);
01400   else
01401     throw ReportError("Epetra_CrsMatrix::ExtractGlobalRowCopy int version called for a matrix that is not int.", -1);
01402 }
01403 #endif
01404 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01405 int Epetra_CrsMatrix::ExtractGlobalRowCopy(long long Row, int Length, int & NumEntries, double * values,
01406              long long * Indices) const
01407 {
01408   if(RowMap().GlobalIndicesLongLong())
01409     return ExtractGlobalRowCopy<long long>(Row, Length, NumEntries, values, Indices);
01410   else
01411     throw ReportError("Epetra_CrsMatrix::ExtractGlobalRowCopy long long version called for a matrix that is not long long.", -1);
01412 }
01413 #endif
01414 
01415 //==========================================================================
01416 int Epetra_CrsMatrix::ExtractMyRowCopy(int Row, int Length, int & NumEntries, double * values,
01417                int * Indices) const
01418 {
01419 
01420   int ierr = Graph_.ExtractMyRowCopy(Row, Length, NumEntries, Indices);
01421   if (ierr)
01422     EPETRA_CHK_ERR(ierr);
01423 
01424   EPETRA_CHK_ERR(ExtractMyRowCopy(Row, Length, NumEntries, values));
01425   return(0);
01426 }
01427 //==========================================================================
01428 int Epetra_CrsMatrix::NumMyRowEntries(int Row, int & NumEntries) const
01429 {
01430 
01431   if (!MyLRID(Row))
01432     EPETRA_CHK_ERR(-1); // Not in the range of local rows
01433   NumEntries = NumMyEntries(Row);
01434   return(0);
01435 }
01436 
01437 //==========================================================================
01438 template<typename int_type>
01439 int Epetra_CrsMatrix::ExtractGlobalRowCopy(int_type Row, int Length, int & NumEntries, double * values) const
01440 {
01441   int Row0 = Graph_.RowMap().LID(Row); // Normalize row range
01442 
01443   EPETRA_CHK_ERR(ExtractMyRowCopy(Row0, Length, NumEntries, values));
01444   return(0);
01445 }
01446 
01447 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01448 int Epetra_CrsMatrix::ExtractGlobalRowCopy(int Row, int Length, int & NumEntries, double * values) const
01449 {
01450   if(RowMap().GlobalIndicesInt())
01451     return ExtractGlobalRowCopy<int>(Row, Length, NumEntries, values);
01452   else
01453     throw ReportError("Epetra_CrsMatrix::ExtractGlobalRowCopy int version called for a matrix that is not int.", -1);
01454 }
01455 #endif
01456 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01457 int Epetra_CrsMatrix::ExtractGlobalRowCopy(long long Row, int Length, int & NumEntries, double * values) const
01458 {
01459   if(RowMap().GlobalIndicesLongLong())
01460     return ExtractGlobalRowCopy<long long>(Row, Length, NumEntries, values);
01461   else
01462     throw ReportError("Epetra_CrsMatrix::ExtractGlobalRowCopy long long version called for a matrix that is not long long.", -1);
01463 }
01464 #endif
01465 
01466 //==========================================================================
01467 int Epetra_CrsMatrix::ExtractMyRowCopy(int Row, int Length, int & NumEntries, double * targValues) const
01468 {
01469   int j;
01470 
01471   if (Row < 0 || Row >= NumMyRows_)
01472     EPETRA_CHK_ERR(-1); // Not in Row range
01473 
01474   NumEntries = Graph().NumMyIndices(Row);
01475   if (Length < NumEntries)
01476     EPETRA_CHK_ERR(-2); // Not enough space for copy. Needed size is passed back in NumEntries
01477 
01478   double * srcValues = Values(Row);
01479 
01480   for(j=0; j<NumEntries; j++)
01481     targValues[j] = srcValues[j];
01482 
01483   return(0);
01484 }
01485 
01486 //==============================================================================
01487 int Epetra_CrsMatrix::ExtractDiagonalCopy(Epetra_Vector & Diagonal) const {
01488 
01489   if(!Filled())
01490     EPETRA_CHK_ERR(-1); // Can't get diagonal unless matrix is filled (and in local index space)
01491   if(!RowMap().SameAs(Diagonal.Map()))
01492     EPETRA_CHK_ERR(-2); // Maps must be the same
01493 
01494   for(int i = 0; i < NumMyRows_; i++) {
01495     long long ii = GRID64(i);
01496     int NumEntries = Graph().NumMyIndices(i);
01497     int* Indices = Graph().Indices(i);
01498     double * srcValues = Values(i);
01499 
01500     Diagonal[i] = 0.0;
01501     for(int j = 0; j < NumEntries; j++) {
01502       if(ii == GCID64(Indices[j])) {
01503   Diagonal[i] = srcValues[j];
01504   break;
01505       }
01506     }
01507   }
01508   return(0);
01509 }
01510 //==============================================================================
01511 int Epetra_CrsMatrix::ReplaceDiagonalValues(const Epetra_Vector & Diagonal) {
01512 
01513   if(!Filled())
01514     EPETRA_CHK_ERR(-1); // Can't replace diagonal unless matrix is filled (and in local index space)
01515   if(!RowMap().SameAs(Diagonal.Map()))
01516     EPETRA_CHK_ERR(-2); // Maps must be the same
01517 
01518   int ierr = 0;
01519   for(int i = 0; i < NumMyRows_; i++) {
01520     long long ii = GRID64(i);
01521     int NumEntries = Graph().NumMyIndices(i);
01522     int* Indices = Graph().Indices(i);
01523     double * targValues = Values(i);
01524     bool DiagMissing = true;
01525     for(int j = 0; j < NumEntries; j++) {
01526       if(ii == GCID64(Indices[j])) {
01527   targValues[j] = Diagonal[i];
01528   DiagMissing = false;
01529   break;
01530       }
01531     }
01532     if(DiagMissing)
01533       ierr = 1; // flag a warning error
01534   }
01535 
01536   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
01537   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
01538   NormFrob_ = -1.0;
01539 
01540   EPETRA_CHK_ERR(ierr);
01541 
01542   return(0);
01543 }
01544 
01545 //==========================================================================
01546 template<typename int_type>
01547 int Epetra_CrsMatrix::ExtractGlobalRowView(int_type Row, int & NumEntries, double *& values, int_type *& Indices) const
01548 {
01549   int ierr = Graph_.ExtractGlobalRowView(Row, NumEntries, Indices);
01550   if (ierr)
01551     EPETRA_CHK_ERR(ierr);
01552 
01553   EPETRA_CHK_ERR(ExtractGlobalRowView(Row, NumEntries, values));
01554   return(0);
01555 }
01556 
01557 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01558 int Epetra_CrsMatrix::ExtractGlobalRowView(int Row, int & NumEntries, double *& values, int *& Indices) const
01559 {
01560   if(RowMap().GlobalIndicesInt())
01561     return ExtractGlobalRowView<int>(Row, NumEntries, values, Indices);
01562   else
01563     throw ReportError("Epetra_CrsMatrix::ExtractGlobalRowView int version called for a matrix that is not int.", -1);
01564 }
01565 #endif
01566 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01567 int Epetra_CrsMatrix::ExtractGlobalRowView(long long Row, int & NumEntries, double *& values, long long *& Indices) const
01568 {
01569   if(RowMap().GlobalIndicesLongLong())
01570     return ExtractGlobalRowView<long long>(Row, NumEntries, values, Indices);
01571   else
01572     throw ReportError("Epetra_CrsMatrix::ExtractGlobalRowView long long version called for a matrix that is not long long.", -1);
01573 }
01574 #endif
01575 
01576 //==========================================================================
01577 int Epetra_CrsMatrix::ExtractMyRowView(int Row, int & NumEntries, double *& values, int *& Indices) const
01578 {
01579   int ierr = Graph_.ExtractMyRowView(Row, NumEntries, Indices);
01580   if (ierr)
01581     EPETRA_CHK_ERR(ierr);
01582 
01583   EPETRA_CHK_ERR(ExtractMyRowView(Row, NumEntries, values));
01584   return(0);
01585 }
01586 //==========================================================================
01587 template<typename int_type>
01588 int Epetra_CrsMatrix::ExtractGlobalRowView(int_type Row, int & NumEntries, double *& values) const
01589 {
01590   int Row0 = Graph_.RowMap().LID(Row); // Normalize row range
01591 
01592   EPETRA_CHK_ERR(ExtractMyRowView(Row0, NumEntries, values));
01593   return(0);
01594 }
01595 
01596 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
01597 int Epetra_CrsMatrix::ExtractGlobalRowView(int Row, int & NumEntries, double *& values) const
01598 {
01599   if(RowMap().GlobalIndicesInt())
01600     return ExtractGlobalRowView<int>(Row, NumEntries, values);
01601   else
01602     throw ReportError("Epetra_CrsMatrix::ExtractGlobalRowView int version called for a matrix that is not int.", -1);
01603 }
01604 #endif
01605 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
01606 int Epetra_CrsMatrix::ExtractGlobalRowView(long long Row, int & NumEntries, double *& values) const
01607 {
01608   if(RowMap().GlobalIndicesLongLong())
01609     return ExtractGlobalRowView<long long>(Row, NumEntries, values);
01610   else
01611     throw ReportError("Epetra_CrsMatrix::ExtractGlobalRowView long long version called for a matrix that is not long long.", -1);
01612 }
01613 #endif
01614 
01615 //==========================================================================
01616 int Epetra_CrsMatrix::ExtractMyRowView(int Row, int & NumEntries, double *& targValues) const
01617 {
01618   if (Row < 0 || Row >= NumMyRows_)
01619     EPETRA_CHK_ERR(-1); // Not in Row range
01620 
01621   NumEntries = Graph().NumMyIndices(Row);
01622 
01623   targValues = Values(Row);
01624 
01625   return(0);
01626 }
01627 
01628 //=============================================================================
01629 int Epetra_CrsMatrix::Solve(bool Upper, bool Trans, bool UnitDiagonal,
01630           const Epetra_Vector& x, Epetra_Vector& y) const
01631 {
01632 
01633 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS
01634   TEUCHOS_FUNC_TIME_MONITOR("Epetra_CrsMatrix::Solve(Upper,Trans,UnitDiag,x,y)");
01635 #endif
01636 
01637   //
01638   // This function finds y such that Ly = x or Uy = x or the transpose cases.
01639   //
01640 
01641   // First short-circuit to the pre-5.0 version if no storage optimization was performed
01642   if (!StorageOptimized() && !Graph().StorageOptimized()) {
01643     EPETRA_CHK_ERR(Solve1(Upper, Trans, UnitDiagonal, x, y));
01644     return(0);
01645   }
01646 
01647   if (!Filled()) {
01648     EPETRA_CHK_ERR(-1); // Matrix must be filled.
01649   }
01650 
01651   if ((Upper) && (!UpperTriangular()))
01652     EPETRA_CHK_ERR(-2);
01653   if ((!Upper) && (!LowerTriangular()))
01654     EPETRA_CHK_ERR(-3);
01655   if ((!UnitDiagonal) && (NoDiagonal()))
01656     EPETRA_CHK_ERR(-4); // If matrix has no diagonal, we must use UnitDiagonal
01657   if ((!UnitDiagonal) && (NumMyDiagonals()<NumMyRows_))
01658     EPETRA_CHK_ERR(-5); // Need each row to have a diagonal
01659 
01660   double *xp = (double*)x.Values();
01661   double *yp = (double*)y.Values();
01662 
01663 
01664   GeneralSV(Upper, Trans, UnitDiagonal, xp, yp);
01665 
01666   UpdateFlops(2*NumGlobalNonzeros64());
01667   return(0);
01668 }
01669 
01670 //=============================================================================
01671 int Epetra_CrsMatrix::Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
01672 
01673 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS
01674   TEUCHOS_FUNC_TIME_MONITOR("Epetra_CrsMatrix::Solve(Upper,Trans,UnitDiag,X,Y)");
01675 #endif
01676 
01677   //
01678   // This function find Y such that LY = X or UY = X or the transpose cases.
01679   //
01680 
01681   // First short-circuit to the pre-5.0 version if no storage optimization was performed
01682   if (!StorageOptimized() && !Graph().StorageOptimized()) {
01683     EPETRA_CHK_ERR(Solve1(Upper, Trans, UnitDiagonal, X, Y));
01684     return(0);
01685   }
01686   if(!Filled())
01687     EPETRA_CHK_ERR(-1); // Matrix must be filled.
01688 
01689   if((Upper) && (!UpperTriangular()))
01690     EPETRA_CHK_ERR(-2);
01691   if((!Upper) && (!LowerTriangular()))
01692     EPETRA_CHK_ERR(-3);
01693   if((!UnitDiagonal) && (NoDiagonal()))
01694     EPETRA_CHK_ERR(-4); // If matrix has no diagonal, we must use UnitDiagonal
01695   if((!UnitDiagonal) && (NumMyDiagonals()<NumMyRows_))
01696     EPETRA_CHK_ERR(-5); // Need each row to have a diagonal
01697 
01698   double** Xp = (double**) X.Pointers();
01699   double** Yp = (double**) Y.Pointers();
01700   int LDX = X.ConstantStride() ? X.Stride() : 0;
01701   int LDY = Y.ConstantStride() ? Y.Stride() : 0;
01702   int NumVectors = X.NumVectors();
01703 
01704 
01705     // Do actual computation
01706   if (NumVectors==1)
01707     GeneralSV(Upper, Trans, UnitDiagonal, *Xp, *Yp);
01708   else
01709     GeneralSM(Upper, Trans, UnitDiagonal, Xp, LDX, Yp, LDY, NumVectors);
01710 
01711   UpdateFlops(2 * NumVectors * NumGlobalNonzeros64());
01712   return(0);
01713 }
01714 
01715 //=============================================================================
01716 int Epetra_CrsMatrix::InvRowSums(Epetra_Vector& x) const {
01717   //
01718   // Put inverse of the sum of absolute values of the ith row of A in x[i].
01719   //
01720 
01721   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
01722   int ierr = 0;
01723   int i, j;
01724   x.PutScalar(0.0); // Make sure we sum into a vector of zeros.
01725   double * xp = (double*)x.Values();
01726   if (Graph().RangeMap().SameAs(x.Map()) && Exporter() != 0) {
01727     Epetra_Vector x_tmp(RowMap());
01728     x_tmp.PutScalar(0.0);
01729     double * x_tmp_p = (double*)x_tmp.Values();
01730     for (i=0; i < NumMyRows_; i++) {
01731       int      NumEntries = NumMyEntries(i);
01732       double * RowValues  = Values(i);
01733       for (j=0; j < NumEntries; j++)  x_tmp_p[i] += std::abs(RowValues[j]);
01734     }
01735     EPETRA_CHK_ERR(x.Export(x_tmp, *Exporter(), Add)); //Export partial row sums to x.
01736     int myLength = x.MyLength();
01737     for (i=0; i<myLength; i++) {
01738       if (xp[i]<Epetra_MinDouble) {
01739         if (xp[i]==0.0) ierr = 1; // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)
01740         else if (ierr!=1) ierr = 2;
01741         xp[i] = Epetra_MaxDouble;
01742       }
01743       else
01744         xp[i] = 1.0/xp[i];
01745     }
01746   }
01747   else if (Graph().RowMap().SameAs(x.Map())) {
01748     for (i=0; i < NumMyRows_; i++) {
01749       int      NumEntries = NumMyEntries(i);
01750       double * RowValues  = Values(i);
01751       double scale = 0.0;
01752       for (j=0; j < NumEntries; j++) scale += std::abs(RowValues[j]);
01753       if (scale<Epetra_MinDouble) {
01754         if (scale==0.0) ierr = 1; // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)
01755         else if (ierr!=1) ierr = 2;
01756         xp[i] = Epetra_MaxDouble;
01757       }
01758       else
01759         xp[i] = 1.0/scale;
01760     }
01761   }
01762   else { // x.Map different than both Graph().RowMap() and Graph().RangeMap()
01763     EPETRA_CHK_ERR(-2); // The map of x must be the RowMap or RangeMap of A.
01764   }
01765   UpdateFlops(NumGlobalNonzeros64());
01766   EPETRA_CHK_ERR(ierr);
01767   return(0);
01768 }
01769 
01770 //=============================================================================
01771 int Epetra_CrsMatrix::InvRowMaxs(Epetra_Vector& x) const {
01772   //
01773   // Put inverse of the max of absolute values of the ith row of A in x[i].
01774   //
01775 
01776   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
01777   int ierr = 0;
01778   int i, j;
01779   bool needExport = false;
01780   double * xp = (double*)x.Values();
01781   Epetra_Vector* x_tmp = 0;
01782   if (Graph().RangeMap().SameAs(x.Map())) {
01783     if (Exporter() != 0) {
01784       needExport = true; //Having this information later avoids a .SameAs
01785       x_tmp = new Epetra_Vector(RowMap()); // Create import vector if needed
01786       xp = (double*)x_tmp->Values();
01787     }
01788   }
01789   else if (!Graph().RowMap().SameAs(x.Map())) {
01790     EPETRA_CHK_ERR(-2); // The map of x must be the RowMap or RangeMap of A.
01791   }
01792   for (i=0; i < NumMyRows_; i++) {
01793     int      NumEntries = NumMyEntries(i);
01794     double * RowValues  = Values(i);
01795     double scale = 0.0;
01796     for (j=0; j < NumEntries; j++) scale = EPETRA_MAX(std::abs(RowValues[j]),scale);
01797     if (scale<Epetra_MinDouble) {
01798       if (scale==0.0) ierr = 1; // Set error to 1 to signal that zero rowmax found (supercedes ierr = 2)
01799       else if (ierr!=1) ierr = 2;
01800       xp[i] = Epetra_MaxDouble;
01801     }
01802     else
01803       xp[i] = 1.0/scale;
01804   }
01805   if(needExport) {
01806     x.PutScalar(0.0);
01807     EPETRA_CHK_ERR(x.Export(*x_tmp, *Exporter(), Insert)); // Fill x with values from temp vector
01808     delete x_tmp;
01809   }
01810   UpdateFlops(NumGlobalNonzeros64());
01811   EPETRA_CHK_ERR(ierr);
01812   return(0);
01813 }
01814 
01815 //=============================================================================
01816 int Epetra_CrsMatrix::InvColSums(Epetra_Vector& x) const {
01817   //
01818   // Put inverse of the sum of absolute values of the jth column of A in x[j].
01819   //
01820 
01821   if(!Filled())  EPETRA_CHK_ERR(-1); // Matrix must be filled.
01822   int ierr = 0;
01823   int i, j;
01824   int MapNumMyElements = x.Map().NumMyElements();
01825   x.PutScalar(0.0); // Make sure we sum into a vector of zeros.
01826   double* xp = (double*)x.Values();
01827   if(Graph().DomainMap().SameAs(x.Map()) && Importer() != 0) {
01828     Epetra_Vector x_tmp(ColMap());
01829     x_tmp.PutScalar(0.0);
01830     double * x_tmp_p = (double*)x_tmp.Values();
01831     for(i = 0; i < NumMyRows_; i++) {
01832       int     NumEntries = NumMyEntries(i);
01833       int*    ColIndices = Graph().Indices(i);
01834       double* RowValues  = Values(i);
01835       for(j = 0; j < NumEntries; j++)
01836         x_tmp_p[ColIndices[j]] += std::abs(RowValues[j]);
01837     }
01838     EPETRA_CHK_ERR(x.Export(x_tmp, *Importer(), Add)); // Fill x with partial column sums
01839   }
01840   else if(Graph().ColMap().SameAs(x.Map())) {
01841     for(i = 0; i < NumMyRows_; i++) {
01842       int     NumEntries = NumMyEntries(i);
01843       int*    ColIndices = Graph().Indices(i);
01844       double* RowValues  = Values(i);
01845       for(j = 0; j < NumEntries; j++)
01846         xp[ColIndices[j]] += std::abs(RowValues[j]);
01847     }
01848   }
01849   else { //x.Map different than both Graph().ColMap() and Graph().DomainMap()
01850     EPETRA_CHK_ERR(-2); // x must have the same distribution as the domain of A
01851   }
01852 
01853   // Invert values, don't allow them to get too large
01854   for(i = 0; i < MapNumMyElements; i++) {
01855     double scale = xp[i];
01856     if(scale < Epetra_MinDouble) {
01857       if(scale == 0.0)
01858   ierr = 1; // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)
01859       else if(ierr != 1)
01860   ierr = 2;
01861       xp[i] = Epetra_MaxDouble;
01862     }
01863     else
01864       xp[i] = 1.0 / scale;
01865   }
01866 
01867   UpdateFlops(NumGlobalNonzeros64());
01868   EPETRA_CHK_ERR(ierr);
01869   return(0);
01870 }
01871 
01872 //=============================================================================
01873 int Epetra_CrsMatrix::InvColMaxs(Epetra_Vector& x) const {
01874   //
01875   // Put inverse of the max of absolute values of the jth column of A in x[j].
01876   //
01877 
01878   if(!Filled())  EPETRA_CHK_ERR(-1); // Matrix must be filled.
01879   int ierr = 0;
01880   int i, j;
01881   int MapNumMyElements = x.Map().NumMyElements();
01882   x.PutScalar(0.0); // Make sure we sum into a vector of zeros.
01883   double* xp = (double*)x.Values();
01884   if(Graph().DomainMap().SameAs(x.Map()) && Importer() != 0) {
01885     Epetra_Vector x_tmp(ColMap());
01886     x_tmp.PutScalar(0.0);
01887     double * x_tmp_p = (double*)x_tmp.Values();
01888     for(i = 0; i < NumMyRows_; i++) {
01889       int     NumEntries = NumMyEntries(i);
01890       int*    ColIndices = Graph().Indices(i);
01891       double* RowValues  = Values(i);
01892       for(j = 0; j < NumEntries; j++)
01893         x_tmp_p[ColIndices[j]] = EPETRA_MAX(std::abs(RowValues[j]),x_tmp_p[ColIndices[j]]);
01894     }
01895     EPETRA_CHK_ERR(x.Export(x_tmp, *Importer(), AbsMax)); // Fill x with partial column sums
01896   }
01897   else if(Graph().ColMap().SameAs(x.Map())) {
01898     for(i = 0; i < NumMyRows_; i++) {
01899       int     NumEntries = NumMyEntries(i);
01900       int*    ColIndices = Graph().Indices(i);
01901       double* RowValues  = Values(i);
01902       for(j = 0; j < NumEntries; j++)
01903         xp[ColIndices[j]] = EPETRA_MAX(std::abs(RowValues[j]),xp[ColIndices[j]]);
01904     }
01905   }
01906   else { //x.Map different than both Graph().ColMap() and Graph().DomainMap()
01907     EPETRA_CHK_ERR(-2); // x must have the same distribution as the domain of A
01908   }
01909 
01910   // Invert values, don't allow them to get too large
01911   for(i = 0; i < MapNumMyElements; i++) {
01912     double scale = xp[i];
01913     if(scale < Epetra_MinDouble) {
01914       if(scale == 0.0)
01915   ierr = 1; // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)
01916       else if(ierr != 1)
01917   ierr = 2;
01918       xp[i] = Epetra_MaxDouble;
01919     }
01920     else
01921       xp[i] = 1.0 / scale;
01922   }
01923 
01924   UpdateFlops(NumGlobalNonzeros64());
01925   EPETRA_CHK_ERR(ierr);
01926   return(0);
01927 }
01928 
01929 //=============================================================================
01930 int Epetra_CrsMatrix::LeftScale(const Epetra_Vector& x) {
01931   //
01932   // This function scales the ith row of A by x[i].
01933   //
01934 
01935   if(!Filled())
01936     EPETRA_CHK_ERR(-1); // Matrix must be filled.
01937   double* xp = 0;
01938   if(Graph().RangeMap().SameAs(x.Map()))
01939     // If we have a non-trivial exporter, we must import elements that are
01940     // permuted or are on other processors.  (We will use the exporter to
01941     // perform the import.)
01942     if(Exporter() != 0) {
01943       UpdateExportVector(1);
01944       EPETRA_CHK_ERR(ExportVector_->Import(x,*Exporter(), Insert));
01945       xp = (double*) ExportVector_->Values();
01946     }
01947     else
01948       xp = (double*)x.Values();
01949   else if (Graph().RowMap().SameAs(x.Map()))
01950     xp = (double*)x.Values();
01951   else {
01952     EPETRA_CHK_ERR(-2); // The Map of x must be the RowMap or RangeMap of A.
01953   }
01954   int i, j;
01955 
01956   for(i = 0; i < NumMyRows_; i++) {
01957     int      NumEntries = NumMyEntries(i);
01958     double* RowValues  = Values(i);
01959     double scale = xp[i];
01960     for(j = 0; j < NumEntries; j++)
01961       RowValues[j] *= scale;
01962   }
01963   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
01964   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
01965   NormFrob_ = -1.0;
01966 
01967   UpdateFlops(NumGlobalNonzeros64());
01968 
01969   return(0);
01970 }
01971 
01972 //=============================================================================
01973 int Epetra_CrsMatrix::RightScale(const Epetra_Vector& x) {
01974   //
01975   // This function scales the jth column of A by x[j].
01976   //
01977 
01978   if(!Filled())
01979     EPETRA_CHK_ERR(-1); // Matrix must be filled.
01980   double* xp = 0;
01981   if(Graph().DomainMap().SameAs(x.Map()))
01982     // If we have a non-trivial exporter, we must import elements that are
01983     // permuted or are on other processors.
01984     if(Importer() != 0) {
01985       UpdateImportVector(1);
01986       EPETRA_CHK_ERR(ImportVector_->Import(x, *Importer(), Insert));
01987       xp = (double*) ImportVector_->Values();
01988     }
01989     else
01990       xp = (double*)x.Values();
01991   else if(Graph().ColMap().SameAs(x.Map()))
01992     xp = (double*)x.Values();
01993   else
01994     EPETRA_CHK_ERR(-2); // The Map of x must be the RowMap or RangeMap of A.
01995   int i, j;
01996 
01997   for(i = 0; i < NumMyRows_; i++) {
01998     int     NumEntries = NumMyEntries(i);
01999     int*    ColIndices = Graph().Indices(i);
02000     double* RowValues  = Values(i);
02001     for(j = 0; j < NumEntries; j++)
02002       RowValues[j] *=  xp[ColIndices[j]];
02003   }
02004   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
02005   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
02006   NormFrob_ = -1.0;
02007 
02008   UpdateFlops(NumGlobalNonzeros64());
02009   return(0);
02010 }
02011 
02012 //=============================================================================
02013 double Epetra_CrsMatrix::NormInf() const {
02014 
02015 #if 0
02016   //
02017   //  Commenting this section out disables caching, ie.
02018   //  causes the norm to be computed each time NormInf is called.
02019   //  See bug #1151 for a full discussion.
02020   //
02021   double MinNorm ;
02022   Comm().MinAll( &NormInf_, &MinNorm, 1 ) ;
02023 
02024   if( MinNorm >= 0.0)
02025     return(NormInf_);
02026 #endif
02027 
02028   if(!Filled())
02029     EPETRA_CHK_ERR(-1); // Matrix must be filled.
02030 
02031   Epetra_Vector x(RangeMap()); // Need temp vector for row sums
02032   double* xp = (double*)x.Values();
02033   Epetra_MultiVector* x_tmp = 0;
02034 
02035   // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
02036   if(Exporter() != 0) {
02037     x_tmp = new Epetra_Vector(RowMap()); // Create temporary import vector if needed
02038     xp = (double*)x_tmp->Values();
02039   }
02040   int i, j;
02041 
02042   for(i = 0; i < NumMyRows_; i++) {
02043     xp[i] = 0.0;
02044     int     NumEntries = NumMyEntries(i);
02045     double* RowValues  = Values(i);
02046     for(j = 0; j < NumEntries; j++)
02047       xp[i] += std::abs(RowValues[j]);
02048   }
02049   if(Exporter() != 0) {
02050     x.PutScalar(0.0);
02051     EPETRA_CHK_ERR(x.Export(*x_tmp, *Exporter(), Add)); // Fill x with Values from temp vector
02052   }
02053   x.MaxValue(&NormInf_); // Find max
02054   if(x_tmp != 0)
02055     delete x_tmp;
02056   UpdateFlops(NumGlobalNonzeros64());
02057   return(NormInf_);
02058 }
02059 //=============================================================================
02060 double Epetra_CrsMatrix::NormOne() const {
02061 
02062 #if 0
02063   //
02064   //  Commenting this section out disables caching, ie.
02065   //  causes the norm to be computed each time NormOne is called.
02066   //  See bug #1151 for a full discussion.
02067   //
02068   double MinNorm ;
02069   Comm().MinAll( &NormOne_, &MinNorm, 1 ) ;
02070 
02071   if( MinNorm >= 0.0)
02072     return(NormOne_);
02073 #endif
02074 
02075   if(!Filled())
02076     EPETRA_CHK_ERR(-1); // Matrix must be filled.
02077 
02078   Epetra_Vector x(DomainMap()); // Need temp vector for column sums
02079 
02080   double* xp = (double*)x.Values();
02081   Epetra_MultiVector* x_tmp = 0;
02082   int NumCols = NumMyCols();
02083 
02084 
02085   // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
02086   if(Importer() != 0) {
02087     x_tmp = new Epetra_Vector(ColMap()); // Create temporary import vector if needed
02088     xp = (double*)x_tmp->Values();
02089   }
02090   int i, j;
02091 
02092   for(i = 0; i < NumCols; i++)
02093     xp[i] = 0.0;
02094 
02095   for(i = 0; i < NumMyRows_; i++) {
02096     int     NumEntries = NumMyEntries(i);
02097     int*    ColIndices = Graph().Indices(i);
02098     double* RowValues  = Values(i);
02099     for(j = 0; j < NumEntries; j++)
02100       xp[ColIndices[j]] += std::abs(RowValues[j]);
02101   }
02102   if(Importer() != 0) {
02103     x.PutScalar(0.0);
02104     EPETRA_CHK_ERR(x.Export(*x_tmp, *Importer(), Add)); // Fill x with Values from temp vector
02105   }
02106   x.MaxValue(&NormOne_); // Find max
02107   if(x_tmp != 0)
02108     delete x_tmp;
02109   UpdateFlops(NumGlobalNonzeros64());
02110   return(NormOne_);
02111 }
02112 //=============================================================================
02113 double Epetra_CrsMatrix::NormFrobenius() const {
02114 
02115 #if 0
02116   //
02117   //  Commenting this section out disables caching, ie.
02118   //  causes the norm to be computed each time NormFrobenius is called.
02119   //  See bug #1151 for a full discussion.
02120   //
02121   double MinNorm ;
02122   Comm().MinAll( &NormFrob_, &MinNorm, 1 ) ;
02123 
02124   if( MinNorm >= 0.0)
02125     return(NormFrob_);
02126 #endif
02127 
02128   if(!Filled())
02129     EPETRA_CHK_ERR(-1); // Matrix must be filled.
02130 
02131   double local_sum = 0.0;
02132 
02133   for(int i = 0; i < NumMyRows_; i++) {
02134     int     NumEntries = NumMyEntries(i);
02135     double* RowValues  = Values(i);
02136     for(int j = 0; j < NumEntries; j++) {
02137       local_sum += RowValues[j]*RowValues[j];
02138     }
02139   }
02140 
02141   double global_sum = 0.0;
02142   Comm().SumAll(&local_sum, &global_sum, 1);
02143 
02144   NormFrob_ = std::sqrt(global_sum);
02145 
02146   UpdateFlops(NumGlobalNonzeros64());
02147 
02148   return(NormFrob_);
02149 }
02150 //=========================================================================
02151 int Epetra_CrsMatrix::CheckSizes(const Epetra_SrcDistObject & Source) {
02152   try {
02153     const Epetra_CrsMatrix & A = dynamic_cast<const Epetra_CrsMatrix &>(Source);
02154     if (!A.Graph().GlobalConstantsComputed()) EPETRA_CHK_ERR(-1); // Must have global constants to proceed
02155   }
02156   catch (...) {
02157     return(0); // No error at this point, object could be a RowMatrix
02158   }
02159   return(0);
02160 }
02161 
02162 //=========================================================================
02163 int Epetra_CrsMatrix::CopyAndPermute(const Epetra_SrcDistObject & Source,
02164              int NumSameIDs,
02165              int NumPermuteIDs,
02166                                      int * PermuteToLIDs,
02167              int *PermuteFromLIDs,
02168                                      const Epetra_OffsetIndex * Indexor ) {
02169 
02170   if(!Source.Map().GlobalIndicesTypeMatch(RowMap()))
02171     throw ReportError("Epetra_CrsMatrix::CopyAndPermute: Incoming global index type does not match the one for *this",-1);
02172 
02173   try {
02174     const Epetra_CrsMatrix & A = dynamic_cast<const Epetra_CrsMatrix &>(Source);
02175     EPETRA_CHK_ERR(CopyAndPermuteCrsMatrix(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs,
02176              PermuteFromLIDs,Indexor));
02177   }
02178   catch (...) {
02179     try {
02180       const Epetra_RowMatrix & A = dynamic_cast<const Epetra_RowMatrix &>(Source);
02181       EPETRA_CHK_ERR(CopyAndPermuteRowMatrix(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs,
02182                PermuteFromLIDs,Indexor));
02183     }
02184     catch (...) {
02185       EPETRA_CHK_ERR(-1); // Incompatible SrcDistObject
02186     }
02187   }
02188 
02189   return(0);
02190 }
02191 
02192 //=========================================================================
02193 template<typename int_type>
02194 int Epetra_CrsMatrix::TCopyAndPermuteCrsMatrix(const Epetra_CrsMatrix & A,
02195                                               int NumSameIDs,
02196                 int NumPermuteIDs,
02197                                               int * PermuteToLIDs,
02198                 int *PermuteFromLIDs,
02199                                               const Epetra_OffsetIndex * Indexor) {
02200 
02201   int i, ierr;
02202 
02203   int_type Row;
02204   int NumEntries;
02205   int maxNumEntries = A.MaxNumEntries();
02206   int_type * Indices = 0;
02207   double * values = 0;
02208 
02209   if (maxNumEntries>0 && A.IndicesAreLocal() ) { //Need Temp Space
02210     Indices = new int_type[maxNumEntries];
02211     values = new double[maxNumEntries];
02212   }
02213 
02214   // Do copy first
02215   if (NumSameIDs>0) {
02216     if (A.IndicesAreLocal()) {
02217       if (StaticGraph() || IndicesAreLocal()) {
02218         if(Indexor) {
02219           for (i=0; i<NumSameIDs; i++) {
02220       Row = (int_type) GRID64(i);
02221       EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(Row, maxNumEntries, NumEntries, values, Indices)); // Set pointers
02222             ierr = ReplaceOffsetValues(Row, NumEntries, values, Indexor->SameOffsets()[i]);
02223             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
02224           }
02225         }
02226         else {
02227           for (i=0; i<NumSameIDs; i++) {
02228       Row = (int_type) GRID64(i);
02229       EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(Row, maxNumEntries, NumEntries, values, Indices)); // Set pointers
02230             ierr = ReplaceGlobalValues(Row, NumEntries, values, Indices);
02231             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
02232           }
02233         }
02234       }
02235       else {
02236         if(Indexor) {
02237           for (i=0; i<NumSameIDs; i++) {
02238       Row = (int_type) GRID64(i);
02239       EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(Row, maxNumEntries, NumEntries, values, Indices)); // Set pointers
02240             ierr = InsertOffsetValues(Row, NumEntries, values, Indexor->SameOffsets()[i]);
02241             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
02242           }
02243         }
02244         else {
02245           for (i=0; i<NumSameIDs; i++) {
02246       Row = (int_type) GRID64(i);
02247       EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(Row, maxNumEntries, NumEntries, values, Indices)); // Set pointers
02248             ierr = InsertGlobalValues(Row, NumEntries, values, Indices);
02249             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
02250           }
02251         }
02252       }
02253     }
02254     else { // A.IndicesAreGlobal()
02255       if (StaticGraph() || IndicesAreLocal()) {
02256         if(Indexor) {
02257           for (i=0; i<NumSameIDs; i++) {
02258       Row = (int_type) GRID64(i);
02259       EPETRA_CHK_ERR(A.ExtractGlobalRowView(Row, NumEntries, values, Indices)); // Set pointers
02260             ierr = ReplaceOffsetValues(Row, NumEntries, values, Indexor->SameOffsets()[i]);
02261             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
02262           }
02263         }
02264         else {
02265           for (i=0; i<NumSameIDs; i++) {
02266       Row = (int_type) GRID64(i);
02267       EPETRA_CHK_ERR(A.ExtractGlobalRowView(Row, NumEntries, values, Indices)); // Set pointers
02268             ierr = ReplaceGlobalValues(Row, NumEntries, values, Indices);
02269             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
02270           }
02271         }
02272       }
02273       else {
02274         if(Indexor) {
02275           for (i=0; i<NumSameIDs; i++) {
02276       Row = (int_type) GRID64(i);
02277       EPETRA_CHK_ERR(A.ExtractGlobalRowView(Row, NumEntries, values, Indices)); // Set pointers
02278             ierr = InsertOffsetValues(Row, NumEntries, values, Indexor->SameOffsets()[i]);
02279             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
02280           }
02281         }
02282         else {
02283           for (i=0; i<NumSameIDs; i++) {
02284       Row = (int_type) GRID64(i);
02285       EPETRA_CHK_ERR(A.ExtractGlobalRowView(Row, NumEntries, values, Indices)); // Set pointers
02286             ierr = InsertGlobalValues(Row, NumEntries, values, Indices);
02287             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
02288           }
02289         }
02290       }
02291     }
02292   }
02293 
02294   // Do local permutation next
02295   int_type FromRow, ToRow;
02296   if (NumPermuteIDs>0) {
02297     if (A.IndicesAreLocal()) {
02298       if (StaticGraph() || IndicesAreLocal()) {
02299         if(Indexor) {
02300           for (i=0; i<NumPermuteIDs; i++) {
02301       FromRow = (int_type) A.GRID64(PermuteFromLIDs[i]);
02302       ToRow = (int_type) GRID64(PermuteToLIDs[i]);
02303       EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, maxNumEntries, NumEntries, values, Indices)); // Set pointers
02304       ierr = ReplaceOffsetValues(ToRow, NumEntries, values, Indexor->PermuteOffsets()[i]);
02305             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
02306           }
02307         }
02308         else {
02309           for (i=0; i<NumPermuteIDs; i++) {
02310       FromRow = (int_type) A.GRID64(PermuteFromLIDs[i]);
02311       ToRow = (int_type) GRID64(PermuteToLIDs[i]);
02312       EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, maxNumEntries, NumEntries, values, Indices)); // Set pointers
02313       ierr = ReplaceGlobalValues(ToRow, NumEntries, values, Indices);
02314             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
02315           }
02316         }
02317       }
02318       else {
02319         if(Indexor) {
02320           for (i=0; i<NumPermuteIDs; i++) {
02321       FromRow = (int_type) A.GRID64(PermuteFromLIDs[i]);
02322       ToRow = (int_type) GRID64(PermuteToLIDs[i]);
02323       EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, maxNumEntries, NumEntries, values, Indices)); // Set pointers
02324       ierr = InsertOffsetValues(ToRow, NumEntries, values, Indexor->PermuteOffsets()[i]);
02325       if (ierr<0) EPETRA_CHK_ERR(ierr);
02326           }
02327         }
02328         else {
02329           for (i=0; i<NumPermuteIDs; i++) {
02330       FromRow = (int_type) A.GRID64(PermuteFromLIDs[i]);
02331       ToRow = (int_type) GRID64(PermuteToLIDs[i]);
02332       EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, maxNumEntries, NumEntries, values, Indices)); // Set pointers
02333       ierr = InsertGlobalValues(ToRow, NumEntries, values, Indices);
02334       if (ierr<0) EPETRA_CHK_ERR(ierr);
02335           }
02336         }
02337       }
02338     }
02339     else { // A.IndicesAreGlobal()
02340       if (StaticGraph() || IndicesAreLocal()) {
02341         if(Indexor) {
02342           for (i=0; i<NumPermuteIDs; i++) {
02343       FromRow = (int_type) A.GRID64(PermuteFromLIDs[i]);
02344       ToRow = (int_type) GRID64(PermuteToLIDs[i]);
02345       EPETRA_CHK_ERR(A.ExtractGlobalRowView(FromRow, NumEntries, values, Indices)); // Set pointers
02346       ierr = ReplaceOffsetValues(ToRow, NumEntries, values, Indexor->PermuteOffsets()[i]);
02347       if (ierr<0) EPETRA_CHK_ERR(ierr);
02348           }
02349         }
02350         else {
02351           for (i=0; i<NumPermuteIDs; i++) {
02352       FromRow = (int_type) A.GRID64(PermuteFromLIDs[i]);
02353       ToRow = (int_type) GRID64(PermuteToLIDs[i]);
02354       EPETRA_CHK_ERR(A.ExtractGlobalRowView(FromRow, NumEntries, values, Indices)); // Set pointers
02355       ierr = ReplaceGlobalValues(ToRow, NumEntries, values, Indices);
02356       if (ierr<0) EPETRA_CHK_ERR(ierr);
02357           }
02358         }
02359       }
02360       else {
02361         if(Indexor) {
02362           for (i=0; i<NumPermuteIDs; i++) {
02363       FromRow = (int_type) A.GRID64(PermuteFromLIDs[i]);
02364       ToRow = (int_type) GRID64(PermuteToLIDs[i]);
02365       EPETRA_CHK_ERR(A.ExtractGlobalRowView(FromRow, NumEntries, values, Indices)); // Set pointers
02366       ierr = InsertOffsetValues(ToRow, NumEntries, values, Indexor->PermuteOffsets()[i]);
02367       if (ierr<0) EPETRA_CHK_ERR(ierr);
02368           }
02369         }
02370         else {
02371           for (i=0; i<NumPermuteIDs; i++) {
02372       FromRow = (int_type) A.GRID64(PermuteFromLIDs[i]);
02373       ToRow = (int_type) GRID64(PermuteToLIDs[i]);
02374       EPETRA_CHK_ERR(A.ExtractGlobalRowView(FromRow, NumEntries, values, Indices)); // Set pointers
02375       ierr = InsertGlobalValues(ToRow, NumEntries, values, Indices);
02376       if (ierr<0) EPETRA_CHK_ERR(ierr);
02377           }
02378         }
02379       }
02380     }
02381   }
02382 
02383   if (maxNumEntries>0 && A.IndicesAreLocal() ) { // Delete Temp Space
02384     delete [] values;
02385     delete [] Indices;
02386   }
02387 
02388   return(0);
02389 }
02390 
02391 int Epetra_CrsMatrix::CopyAndPermuteCrsMatrix(const Epetra_CrsMatrix & A,
02392                                               int NumSameIDs,
02393                 int NumPermuteIDs,
02394                                               int * PermuteToLIDs,
02395                 int *PermuteFromLIDs,
02396                                               const Epetra_OffsetIndex * Indexor)
02397 {
02398   if(!A.RowMap().GlobalIndicesTypeMatch(RowMap()))
02399     throw ReportError("Epetra_CrsMatrix::CopyAndPermuteCrsMatrix: Incoming global index type does not match the one for *this",-1);
02400 
02401   if(A.RowMap().GlobalIndicesInt())
02402 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
02403     return TCopyAndPermuteCrsMatrix<int>(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor);
02404 #else
02405     throw ReportError("Epetra_CrsMatrix::CopyAndPermuteCrsMatrix: ERROR, GlobalIndicesInt but no API for it.",-1);
02406 #endif
02407 
02408   if(A.RowMap().GlobalIndicesLongLong())
02409 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
02410     return TCopyAndPermuteCrsMatrix<long long>(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor);
02411 #else
02412     throw ReportError("Epetra_CrsMatrix::CopyAndPermuteCrsMatrix: ERROR, GlobalIndicesLongLong but no API for it.",-1);
02413 #endif
02414 
02415   throw ReportError("Epetra_CrsMatrix::CopyAndPermuteCrsMatrix: Internal error, unable to determine global index type of maps", -1);
02416 }
02417 
02418 //=========================================================================
02419 template<typename int_type>
02420 int Epetra_CrsMatrix::TCopyAndPermuteRowMatrix(const Epetra_RowMatrix & A,
02421                                               int NumSameIDs,
02422                 int NumPermuteIDs,
02423                                               int * PermuteToLIDs,
02424                 int *PermuteFromLIDs,
02425                                               const Epetra_OffsetIndex * Indexor ) {
02426   int i, j, ierr;
02427 
02428   int_type Row, ToRow;
02429   int NumEntries;
02430   int FromRow;
02431   int maxNumEntries = A.MaxNumEntries();
02432   int_type * gen_Indices = 0; // gen = general (int or long long)
02433   int * int_Indices = 0;
02434   double * values = 0;
02435 
02436   if (maxNumEntries>0) {
02437     if(StaticGraph() || IndicesAreLocal() || Indexor) {
02438       int_Indices = new int[maxNumEntries];
02439     }
02440     else {
02441       gen_Indices = new int_type[maxNumEntries];
02442       int_Indices = reinterpret_cast<int*>(gen_Indices);
02443     }
02444 
02445     values = new double[maxNumEntries]; // Must extract values even though we discard them
02446   }
02447 
02448   const Epetra_Map & rowMap = A.RowMatrixRowMap();
02449   const Epetra_Map & colMap = A.RowMatrixColMap();
02450 
02451   // Do copy first
02452   if (NumSameIDs>0) {
02453     if (StaticGraph() || IndicesAreLocal()) {
02454       if( Indexor ) {
02455         for (i=0; i<NumSameIDs; i++) {
02456           Row = (int_type) GRID64(i);
02457           int AlocalRow = rowMap.LID(Row);
02458           EPETRA_CHK_ERR(A.ExtractMyRowCopy(AlocalRow, maxNumEntries, NumEntries, values, int_Indices));
02459     ierr = ReplaceOffsetValues(Row, NumEntries, values, Indexor->SameOffsets()[i]);
02460           if (ierr<0) EPETRA_CHK_ERR(ierr);
02461         }
02462       }
02463       else {
02464         for (i=0; i<NumSameIDs; i++) {
02465           Row = (int_type) GRID64(i);
02466           int AlocalRow = rowMap.LID(Row);
02467           EPETRA_CHK_ERR(A.ExtractMyRowCopy(AlocalRow, maxNumEntries, NumEntries, values, int_Indices));
02468           for(j=0; j<NumEntries; ++j) {
02469             int_Indices[j] = LCID((int_type) colMap.GID64(int_Indices[j]));
02470           }
02471     ierr = ReplaceMyValues(i, NumEntries, values, int_Indices);
02472           if (ierr<0) EPETRA_CHK_ERR(ierr);
02473         }
02474       }
02475     }
02476     else {
02477       if( Indexor ) {
02478         for (i=0; i<NumSameIDs; i++) {
02479           EPETRA_CHK_ERR(A.ExtractMyRowCopy(i, maxNumEntries, NumEntries, values, int_Indices));
02480           Row = (int_type) GRID64(i);
02481     ierr = InsertOffsetValues(Row, NumEntries, values, Indexor->SameOffsets()[i]);
02482           if (ierr<0) EPETRA_CHK_ERR(ierr);
02483         }
02484       }
02485       else {
02486         for (i=0; i<NumSameIDs; i++) {
02487           EPETRA_CHK_ERR(A.ExtractMyRowCopy(i, maxNumEntries, NumEntries, values, int_Indices));
02488           Row = (int_type) GRID64(i);
02489 
02490           // convert to GIDs, start from right.
02491           for(j = NumEntries; j > 0;) {
02492             --j;
02493            gen_Indices[j] = (int_type) colMap.GID64(int_Indices[j]);
02494           }
02495     ierr = InsertGlobalValues(Row, NumEntries, values, gen_Indices);
02496           if (ierr<0) EPETRA_CHK_ERR(ierr);
02497         }
02498       }
02499     }
02500   }
02501 
02502   // Do local permutation next
02503   if (NumPermuteIDs>0) {
02504     if (StaticGraph() || IndicesAreLocal()) {
02505       if( Indexor ) {
02506         for (i=0; i<NumPermuteIDs; i++) {
02507           FromRow = PermuteFromLIDs[i];
02508           EPETRA_CHK_ERR(A.ExtractMyRowCopy(FromRow, maxNumEntries, NumEntries, values, int_Indices));
02509           ToRow = (int_type) GRID64(PermuteToLIDs[i]);
02510     ierr = ReplaceOffsetValues(ToRow, NumEntries, values, Indexor->PermuteOffsets()[i]);
02511           if (ierr<0) EPETRA_CHK_ERR(ierr);
02512         }
02513       }
02514       else {
02515         for (i=0; i<NumPermuteIDs; i++) {
02516           FromRow = PermuteFromLIDs[i];
02517           EPETRA_CHK_ERR(A.ExtractMyRowCopy(FromRow, maxNumEntries, NumEntries, values, int_Indices));
02518           ToRow = (int_type) GRID64(PermuteToLIDs[i]);
02519           for(j=0; j<NumEntries; ++j) {
02520             int_Indices[j] = LCID((int_type) colMap.GID64(int_Indices[j]));
02521           }
02522     ierr = ReplaceMyValues((int) ToRow, NumEntries, values, int_Indices);
02523           if (ierr<0) EPETRA_CHK_ERR(ierr);
02524         }
02525       }
02526     }
02527     else {
02528       if( Indexor ) {
02529         for (i=0; i<NumPermuteIDs; i++) {
02530           FromRow = PermuteFromLIDs[i];
02531           EPETRA_CHK_ERR(A.ExtractMyRowCopy(FromRow, maxNumEntries, NumEntries, values, int_Indices));
02532           ToRow = (int_type) GRID64(PermuteToLIDs[i]);
02533     ierr = InsertOffsetValues(ToRow, NumEntries, values, Indexor->PermuteOffsets()[i]);
02534           if (ierr<0) EPETRA_CHK_ERR(ierr);
02535         }
02536       }
02537       else {
02538         for (i=0; i<NumPermuteIDs; i++) {
02539           FromRow = PermuteFromLIDs[i];
02540           EPETRA_CHK_ERR(A.ExtractMyRowCopy(FromRow, maxNumEntries, NumEntries, values, int_Indices));
02541 
02542           // convert to GIDs, start from right.
02543           for(j = NumEntries; j > 0;) {
02544             --j;
02545            gen_Indices[j] = (int_type) colMap.GID64(int_Indices[j]);
02546           }
02547 
02548           ToRow = (int_type) GRID64(PermuteToLIDs[i]);
02549     ierr = InsertGlobalValues(ToRow, NumEntries, values, gen_Indices);
02550           if (ierr<0) EPETRA_CHK_ERR(ierr);
02551         }
02552       }
02553     }
02554   }
02555 
02556   if (maxNumEntries>0) {
02557     delete [] values;
02558     if(StaticGraph() || IndicesAreLocal() || Indexor) {
02559       delete [] int_Indices;
02560     }
02561     else {
02562       delete [] gen_Indices;
02563     }
02564   }
02565   return(0);
02566 }
02567 
02568 int Epetra_CrsMatrix::CopyAndPermuteRowMatrix(const Epetra_RowMatrix & A,
02569                                               int NumSameIDs,
02570                 int NumPermuteIDs,
02571                                               int * PermuteToLIDs,
02572                 int *PermuteFromLIDs,
02573                                               const Epetra_OffsetIndex * Indexor )
02574 {
02575   if(!A.Map().GlobalIndicesTypeMatch(RowMap()))
02576     throw ReportError("Epetra_CrsMatrix::CopyAndPermuteRowMatrix: Incoming global index type does not match the one for *this",-1);
02577 
02578   if(A.RowMatrixRowMap().GlobalIndicesInt())
02579 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
02580     return TCopyAndPermuteRowMatrix<int>(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor);
02581 #else
02582     throw ReportError("Epetra_CrsMatrix::CopyAndPermuteRowMatrix: ERROR, GlobalIndicesInt but no API for it.",-1);
02583 #endif
02584 
02585   if(A.RowMatrixRowMap().GlobalIndicesLongLong())
02586 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
02587     return TCopyAndPermuteRowMatrix<long long>(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor);
02588 #else
02589     throw ReportError("Epetra_CrsMatrix::CopyAndPermuteRowMatrix: ERROR, GlobalIndicesLongLong but no API for it.",-1);
02590 #endif
02591 
02592   throw ReportError("Epetra_CrsMatrix::CopyAndPermuteRowMatrix: Internal error, unable to determine global index type of maps", -1);
02593 }
02594 
02595 //=========================================================================
02596 int Epetra_CrsMatrix::PackAndPrepare(const Epetra_SrcDistObject & Source,
02597              int NumExportIDs,
02598                                      int * ExportLIDs,
02599              int & LenExports,
02600                                      char *& Exports,
02601              int & SizeOfPacket,
02602                                      int * Sizes,
02603                                      bool & VarSizes,
02604                                      Epetra_Distributor & Distor)
02605 {
02606   if(!Source.Map().GlobalIndicesTypeMatch(RowMap()))
02607     throw ReportError("Epetra_CrsMatrix::PackAndPrepare: Incoming global index type does not match the one for *this",-1);
02608 
02609   (void)Distor;
02610   // Rest of work can be done using RowMatrix only
02611   const Epetra_RowMatrix & A = dynamic_cast<const Epetra_RowMatrix &>(Source);
02612 
02613   VarSizes = true; //enable variable block size data comm
02614 
02615   int TotalSendLength = 0;
02616   int * IntSizes = 0;
02617   if( NumExportIDs>0 ) IntSizes = new int[NumExportIDs];
02618 
02619   int SizeofIntType = -1;
02620   if(Source.Map().GlobalIndicesInt())
02621     SizeofIntType = (int)sizeof(int);
02622   else if(Source.Map().GlobalIndicesLongLong())
02623     SizeofIntType = (int)sizeof(long long);
02624   else
02625     throw ReportError("Epetra_CrsMatrix::PackAndPrepare: Unable to determine source global index type",-1);
02626 
02627   for( int i = 0; i < NumExportIDs; ++i )
02628   {
02629     int NumEntries;
02630     A.NumMyRowEntries( ExportLIDs[i], NumEntries );
02631     // Will have NumEntries doubles, NumEntries +2 ints, pack them interleaved     Sizes[i] = NumEntries;
02632     Sizes[i] = NumEntries;
02633     IntSizes[i] = 1 + (((NumEntries+2)*SizeofIntType)/(int)sizeof(double));
02634     TotalSendLength += (Sizes[i]+IntSizes[i]);
02635   }
02636 
02637   double * DoubleExports = 0;
02638   SizeOfPacket = (int)sizeof(double);
02639 
02640   //setup buffer locally for memory management by this object
02641   if( TotalSendLength*SizeOfPacket > LenExports )
02642   {
02643     if( LenExports > 0 ) delete [] Exports;
02644     LenExports = TotalSendLength*SizeOfPacket;
02645     DoubleExports = new double[TotalSendLength];
02646     for( int i = 0; i < TotalSendLength; ++i ) DoubleExports[i] = 0.0;
02647     Exports = (char *) DoubleExports;
02648   }
02649 
02650   int NumEntries;
02651   double * values;
02652   double * valptr, * dintptr;
02653 
02654   // Each segment of Exports will be filled by a packed row of information for each row as follows:
02655   // 1st int: GRID of row where GRID is the global row ID for the source matrix
02656   // next int:  NumEntries, Number of indices in row.
02657   // next NumEntries: The actual indices for the row.
02658 
02659   const Epetra_Map & rowMap = A.RowMatrixRowMap();
02660   const Epetra_Map & colMap = A.RowMatrixColMap();
02661 
02662   if( NumExportIDs > 0 )
02663   {
02664     if(Source.Map().GlobalIndicesInt()) {
02665       int * Indices;
02666       int FromRow;
02667       int * intptr;
02668 
02669       int maxNumEntries = A.MaxNumEntries();
02670       dintptr = (double *) Exports;
02671       valptr = dintptr + IntSizes[0];
02672       intptr = (int *) dintptr;
02673       for (int i=0; i<NumExportIDs; i++)
02674       {
02675         FromRow = (int) rowMap.GID64(ExportLIDs[i]);
02676         intptr[0] = FromRow;
02677         values = valptr;
02678         Indices = intptr + 2;
02679         EPETRA_CHK_ERR(A.ExtractMyRowCopy(ExportLIDs[i], maxNumEntries, NumEntries, values, Indices));
02680         for (int j=0; j<NumEntries; j++) Indices[j] = (int) colMap.GID64(Indices[j]); // convert to GIDs
02681         intptr[1] = NumEntries; // Load second slot of segment
02682         if( i < (NumExportIDs-1) )
02683         {
02684           dintptr += (IntSizes[i]+Sizes[i]);
02685           valptr = dintptr + IntSizes[i+1];
02686           intptr = (int *) dintptr;
02687         }
02688       }
02689     }
02690     else if(Source.Map().GlobalIndicesLongLong()) {
02691       long long * LL_Indices;
02692       long long FromRow;
02693       long long * LLptr;
02694 
02695       int maxNumEntries = A.MaxNumEntries();
02696       dintptr = (double *) Exports;
02697       valptr = dintptr + IntSizes[0];
02698       LLptr = (long long *) dintptr;
02699       for (int i=0; i<NumExportIDs; i++)
02700       {
02701         FromRow = rowMap.GID64(ExportLIDs[i]);
02702         LLptr[0] = FromRow;
02703         values = valptr;
02704         LL_Indices = LLptr + 2;
02705         int * int_indices = reinterpret_cast<int*>(LL_Indices);
02706         EPETRA_CHK_ERR(A.ExtractMyRowCopy(ExportLIDs[i], maxNumEntries, NumEntries, values, int_indices));
02707 
02708         // convert to GIDs, start from right.
02709         for(int j = NumEntries; j > 0;) {
02710            --j;
02711            LL_Indices[j] = colMap.GID64(int_indices[j]);
02712         }
02713 
02714         LLptr[1] = NumEntries; // Load second slot of segment
02715         if( i < (NumExportIDs-1) )
02716         {
02717           dintptr += (IntSizes[i]+Sizes[i]);
02718           valptr = dintptr + IntSizes[i+1];
02719           LLptr = (long long *) dintptr;
02720         }
02721       }
02722     }
02723 
02724     for( int i = 0; i < NumExportIDs; ++i )
02725       Sizes[i] += IntSizes[i];
02726   }
02727 
02728   if( IntSizes ) delete [] IntSizes;
02729 
02730   return(0);
02731 }
02732 
02733 //=========================================================================
02734 template<typename int_type>
02735 int Epetra_CrsMatrix::TUnpackAndCombine(const Epetra_SrcDistObject & Source,
02736                int NumImportIDs,
02737                                        int * ImportLIDs,
02738                                        int LenImports,
02739                char * Imports,
02740                                        int & SizeOfPacket,
02741                Epetra_Distributor & Distor,
02742                Epetra_CombineMode CombineMode,
02743                                        const Epetra_OffsetIndex * Indexor )
02744 {
02745   (void)Source;
02746   (void)LenImports;
02747   (void)SizeOfPacket;
02748   (void)Distor;
02749   if (NumImportIDs<=0) return(0);
02750 
02751   if (   CombineMode != Add
02752    && CombineMode != Insert
02753    && CombineMode != Zero )
02754     EPETRA_CHK_ERR(-1); //Unsupported CombineMode, defaults to Zero
02755 
02756   int NumEntries;
02757   int_type * Indices;
02758   double * values;
02759   int_type ToRow;
02760   int i, ierr;
02761   int IntSize;
02762 
02763   double * valptr, *dintptr;
02764   int_type * intptr;
02765 
02766   // Each segment of Exports will be filled by a packed row of information for each row as follows:
02767   // 1st int: GRID of row where GRID is the global row ID for the source matrix
02768   // next int:  NumEntries, Number of indices in row.
02769   // next NumEntries: The actual indices for the row.
02770 
02771   dintptr = (double *) Imports;
02772   intptr = (int_type *) dintptr;
02773   NumEntries = (int) intptr[1];
02774   IntSize = 1 + (((NumEntries+2)*(int)sizeof(int_type))/(int)sizeof(double));
02775   valptr = dintptr + IntSize;
02776 
02777   for (i=0; i<NumImportIDs; i++)
02778   {
02779     ToRow = (int_type) GRID64(ImportLIDs[i]);
02780     assert((intptr[0])==ToRow); // Sanity check
02781     values = valptr;
02782     Indices = intptr + 2;
02783 
02784     if (CombineMode==Add) {
02785       if (StaticGraph() || IndicesAreLocal()) {
02786         if( Indexor )
02787           ierr = SumIntoOffsetValues(ToRow, NumEntries, values, Indexor->RemoteOffsets()[i]);
02788         else
02789           ierr = SumIntoGlobalValues(ToRow, NumEntries, values, Indices);
02790       }
02791       else {
02792         if( Indexor )
02793           ierr = InsertOffsetValues(ToRow, NumEntries, values, Indexor->RemoteOffsets()[i]);
02794         else
02795           ierr = InsertGlobalValues(ToRow, NumEntries, values, Indices);
02796       }
02797       if (ierr<0) EPETRA_CHK_ERR(ierr);
02798     }
02799     else if (CombineMode==Insert) {
02800       if (StaticGraph() || IndicesAreLocal()) {
02801         if( Indexor )
02802           ierr = ReplaceOffsetValues(ToRow, NumEntries, values, Indexor->RemoteOffsets()[i]);
02803         else
02804           ierr = ReplaceGlobalValues(ToRow, NumEntries, values, Indices);
02805       }
02806       else {
02807         if( Indexor )
02808           ierr = InsertOffsetValues(ToRow, NumEntries, values, Indexor->RemoteOffsets()[i]);
02809         else
02810           ierr = InsertGlobalValues(ToRow, NumEntries, values, Indices);
02811       }
02812       if (ierr<0) EPETRA_CHK_ERR(ierr);
02813     }
02814 
02815     if( i < (NumImportIDs-1) )
02816     {
02817       dintptr += IntSize + NumEntries;
02818       intptr = (int_type *) dintptr;
02819       NumEntries = (int) intptr[1];
02820       IntSize = 1 + (((NumEntries+2)*(int)sizeof(int_type))/(int)sizeof(double));
02821       valptr = dintptr + IntSize;
02822     }
02823   }
02824 
02825   return(0);
02826 }
02827 
02828 int Epetra_CrsMatrix::UnpackAndCombine(const Epetra_SrcDistObject & Source,
02829                int NumImportIDs,
02830                                        int * ImportLIDs,
02831                                        int LenImports,
02832                char * Imports,
02833                                        int & SizeOfPacket,
02834                Epetra_Distributor & Distor,
02835                Epetra_CombineMode CombineMode,
02836                                        const Epetra_OffsetIndex * Indexor )
02837 {
02838   if(!Source.Map().GlobalIndicesTypeMatch(RowMap()))
02839     throw ReportError("Epetra_CrsMatrix::UnpackAndCombine: Incoming global index type does not match the one for *this",-1);
02840 
02841   if(Source.Map().GlobalIndicesInt())
02842 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
02843     return TUnpackAndCombine<int>(Source, NumImportIDs, ImportLIDs, LenImports,
02844                Imports, SizeOfPacket, Distor, CombineMode, Indexor);
02845 #else
02846     throw ReportError("Epetra_CrsMatrix::UnpackAndCombine: ERROR, GlobalIndicesInt but no API for it.",-1);
02847 #endif
02848 
02849   if(Source.Map().GlobalIndicesLongLong())
02850 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
02851     return TUnpackAndCombine<long long>(Source, NumImportIDs, ImportLIDs, LenImports,
02852                Imports, SizeOfPacket, Distor, CombineMode, Indexor);
02853 #else
02854     throw ReportError("Epetra_CrsMatrix::UnpackAndCombine: ERROR, GlobalIndicesLongLong but no API for it.",-1);
02855 #endif
02856 
02857   throw ReportError("Epetra_CrsMatrix::UnpackAndCombine: Internal error, unable to determine global index type of maps", -1);
02858 }
02859 
02860 //=========================================================================
02861 
02862 void Epetra_CrsMatrix::Print(std::ostream& os) const {
02863   int MyPID = RowMap().Comm().MyPID();
02864   int NumProc = RowMap().Comm().NumProc();
02865 
02866   for (int iproc=0; iproc < NumProc; iproc++) {
02867     if (MyPID==iproc) {
02868       /*      const Epetra_fmtflags olda = os.setf(ios::right,ios::adjustfield);
02869         const Epetra_fmtflags oldf = os.setf(ios::scientific,ios::floatfield);
02870         const int             oldp = os.precision(12); */
02871       if (MyPID==0) {
02872   os <<  "\nNumber of Global Rows        = "; os << NumGlobalRows64(); os << std::endl;
02873   os <<    "Number of Global Cols        = "; os << NumGlobalCols64(); os << std::endl;
02874   os <<    "Number of Global Diagonals   = "; os << NumGlobalDiagonals64(); os << std::endl;
02875   os <<    "Number of Global Nonzeros    = "; os << NumGlobalNonzeros64(); os << std::endl;
02876   os <<    "Global Maximum Num Entries   = "; os << GlobalMaxNumEntries(); os << std::endl;
02877   if (LowerTriangular()) os <<    " ** Matrix is Lower Triangular **"; os << std::endl;
02878   if (UpperTriangular()) os <<    " ** Matrix is Upper Triangular **"; os << std::endl;
02879   if (NoDiagonal())      os <<    " ** Matrix has no diagonal     **"; os << std::endl; os << std::endl;
02880       }
02881 
02882       os <<  "\nNumber of My Rows        = "; os << NumMyRows(); os << std::endl;
02883       os <<    "Number of My Cols        = "; os << NumMyCols(); os << std::endl;
02884       os <<    "Number of My Diagonals   = "; os << NumMyDiagonals(); os << std::endl;
02885       os <<    "Number of My Nonzeros    = "; os << NumMyNonzeros(); os << std::endl;
02886       os <<    "My Maximum Num Entries   = "; os << MaxNumEntries(); os << std::endl; os << std::endl;
02887 
02888       os << std::flush;
02889 
02890       // Reset os flags
02891 
02892       /*      os.setf(olda,ios::adjustfield);
02893         os.setf(oldf,ios::floatfield);
02894         os.precision(oldp); */
02895     }
02896     // Do a few global ops to give I/O a chance to complete
02897     Comm().Barrier();
02898     Comm().Barrier();
02899     Comm().Barrier();
02900   }
02901 
02902   {for (int iproc=0; iproc < NumProc; iproc++) {
02903     if (MyPID==iproc) {
02904       int NumMyRows1 = NumMyRows();
02905       int MaxNumIndices = MaxNumEntries();
02906 
02907       int * Indices_int = 0;
02908       long long * Indices_LL = 0;
02909       if(RowMap().GlobalIndicesInt()) {
02910          Indices_int = new int[MaxNumIndices];
02911       }
02912       else if(RowMap().GlobalIndicesLongLong()) {
02913          Indices_LL = new long long[MaxNumIndices];
02914       }
02915       else
02916          throw ReportError("Epetra_CrsGraph::Print: Unable to determine source global index type",-1);
02917 
02918       double * values  = new double[MaxNumIndices];
02919 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
02920       int NumIndices;
02921       int j;
02922 #endif
02923       int i;
02924 
02925       if (MyPID==0) {
02926   os.width(8);
02927   os <<  "   Processor ";
02928   os.width(10);
02929   os <<  "   Row Index ";
02930   os.width(10);
02931   os <<  "   Col Index ";
02932   os.width(20);
02933   os <<  "   Value     ";
02934   os << std::endl;
02935       }
02936       for (i=0; i<NumMyRows1; i++) {
02937         if(RowMap().GlobalIndicesInt()) {
02938 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
02939            int Row = (int) GRID64(i); // Get global row number
02940            ExtractGlobalRowCopy(Row, MaxNumIndices, NumIndices, values, Indices_int);
02941 
02942            for (j = 0; j < NumIndices ; j++) {
02943               os.width(8);
02944               os <<  MyPID ; os << "    ";
02945               os.width(10);
02946               os <<  Row ; os << "    ";
02947               os.width(10);
02948               os <<  Indices_int[j]; os << "    ";
02949               os.width(20);
02950               os <<  values[j]; os << "    ";
02951               os << std::endl;
02952            }
02953 #else
02954     throw ReportError("Epetra_CrsMatrix::Print: ERROR, GlobalIndicesInt but no API for it.",-1);
02955 #endif
02956         }
02957         else if(RowMap().GlobalIndicesLongLong()) {
02958 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
02959            long long Row = GRID64(i); // Get global row number
02960            ExtractGlobalRowCopy(Row, MaxNumIndices, NumIndices, values, Indices_LL);
02961 
02962            for (j = 0; j < NumIndices ; j++) {
02963               os.width(8);
02964               os <<  MyPID ; os << "    ";
02965               os.width(10);
02966               os <<  Row ; os << "    ";
02967               os.width(10);
02968               os <<  Indices_LL[j]; os << "    ";
02969               os.width(20);
02970               os <<  values[j]; os << "    ";
02971               os << std::endl;
02972            }
02973 #else
02974     throw ReportError("Epetra_CrsMatrix::Print: ERROR, GlobalIndicesLongLong but no API for it.",-1);
02975 #endif
02976         }
02977       }
02978 
02979       if(RowMap().GlobalIndicesInt()) {
02980          delete [] Indices_int;
02981       }
02982       else if(RowMap().GlobalIndicesLongLong()) {
02983          delete [] Indices_LL;
02984       }
02985       delete [] values;
02986 
02987       os << std::flush;
02988 
02989     }
02990     // Do a few global ops to give I/O a chance to complete
02991     RowMap().Comm().Barrier();
02992     RowMap().Comm().Barrier();
02993     RowMap().Comm().Barrier();
02994   }}
02995 
02996   return;
02997 }
02998 //=============================================================================
02999 int Epetra_CrsMatrix::Multiply(bool TransA, const Epetra_Vector& x, Epetra_Vector& y) const {
03000 
03001 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS
03002   TEUCHOS_FUNC_TIME_MONITOR("Epetra_CrsMatrix::Multiply(TransA,x,y)");
03003 #endif
03004 
03005   //
03006   // This function forms the product y = A * x or y = A' * x
03007   //
03008 
03009   if(!Filled())
03010     EPETRA_CHK_ERR(-1); // Matrix must be filled.
03011 
03012   double* xp = (double*) x.Values();
03013   double* yp = (double*) y.Values();
03014 
03015   Epetra_Vector * xcopy = 0;
03016   if (&x==&y && Importer()==0 && Exporter()==0) {
03017     xcopy = new Epetra_Vector(x);
03018     xp = (double *) xcopy->Values();
03019   }
03020   UpdateImportVector(1); // Refresh import and output vectors if needed
03021   UpdateExportVector(1);
03022 
03023   if(!TransA) {
03024 
03025     // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
03026     if(Importer() != 0) {
03027       EPETRA_CHK_ERR(ImportVector_->Import(x, *Importer(), Insert));
03028       xp = (double*) ImportVector_->Values();
03029     }
03030 
03031     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
03032     if(Exporter() != 0)  yp = (double*) ExportVector_->Values();
03033 
03034     // Do actual computation
03035     GeneralMV(xp, yp);
03036 
03037     if(Exporter() != 0) {
03038       y.PutScalar(0.0); // Make sure target is zero
03039       EPETRA_CHK_ERR(y.Export(*ExportVector_, *Exporter(), Add)); // Fill y with Values from export vector
03040     }
03041     // Handle case of rangemap being a local replicated map
03042     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(y.Reduce());
03043   }
03044 
03045   else { // Transpose operation
03046 
03047     // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
03048     if(Exporter() != 0) {
03049       EPETRA_CHK_ERR(ExportVector_->Import(x, *Exporter(), Insert));
03050       xp = (double*) ExportVector_->Values();
03051     }
03052 
03053     // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
03054     if(Importer() != 0) yp = (double*) ImportVector_->Values();
03055 
03056     // Do actual computation
03057     GeneralMTV(xp, yp);
03058 
03059     if(Importer() != 0) {
03060       y.PutScalar(0.0); // Make sure target is zero
03061       EPETRA_CHK_ERR(y.Export(*ImportVector_, *Importer(), Add)); // Fill y with Values from export vector
03062     }
03063     // Handle case of rangemap being a local replicated map
03064     if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(y.Reduce());
03065   }
03066 
03067 
03068   UpdateFlops(2 * NumGlobalNonzeros64());
03069   if (xcopy!=0) {
03070     delete xcopy;
03071     EPETRA_CHK_ERR(1); // Return positive code to alert the user about needing extra copy of x
03072     return(1);
03073   }
03074   return(0);
03075 }
03076 
03077 //=============================================================================
03078 int Epetra_CrsMatrix::Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
03079 
03080 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS
03081   TEUCHOS_FUNC_TIME_MONITOR("Epetra_CrsMatrix::Multiply(TransA,X,Y)");
03082 #endif
03083 
03084 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
03085   Teuchos::RCP<Teuchos::FancyOStream>
03086     out = Teuchos::VerboseObjectBase::getDefaultOStream();
03087   Teuchos::OSTab tab(out);
03088   if(Epetra_CrsMatrixTraceDumpMultiply) {
03089     *out << std::boolalpha;
03090     *out << "\nEntering Epetra_CrsMatrix::Multipy("<<TransA<<",X,Y) ...\n";
03091     if(!TransA) {
03092       *out << "\nDomainMap =\n";
03093       this->DomainMap().Print(Teuchos::OSTab(out).o());
03094     }
03095     else {
03096       *out << "\nRangeMap =\n";
03097       this->RangeMap().Print(Teuchos::OSTab(out).o());
03098     }
03099     *out << "\nInitial input X with " << ( TransA ? "RangeMap" : "DomainMap" ) << " =\n\n";
03100     X.Print(Teuchos::OSTab(out).o());
03101   }
03102 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
03103 
03104   //
03105   // This function forms the product Y = A * Y or Y = A' * X
03106   //
03107   if(!Filled()) {
03108     EPETRA_CHK_ERR(-1); // Matrix must be filled.
03109   }
03110 
03111   int NumVectors = X.NumVectors();
03112   if (NumVectors!=Y.NumVectors()) {
03113     EPETRA_CHK_ERR(-2); // Need same number of vectors in each MV
03114   }
03115 
03116   double** Xp = (double**) X.Pointers();
03117   double** Yp = (double**) Y.Pointers();
03118 
03119   int LDX = X.ConstantStride() ? X.Stride() : 0;
03120   int LDY = Y.ConstantStride() ? Y.Stride() : 0;
03121 
03122   Epetra_MultiVector * Xcopy = 0;
03123   if (&X==&Y && Importer()==0 && Exporter()==0) {
03124     Xcopy = new Epetra_MultiVector(X);
03125     Xp = (double **) Xcopy->Pointers();
03126     LDX = Xcopy->ConstantStride() ? Xcopy->Stride() : 0;
03127   }
03128   UpdateImportVector(NumVectors); // Make sure Import and Export Vectors are compatible
03129   UpdateExportVector(NumVectors);
03130 
03131   if (!TransA) {
03132 
03133     // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
03134     if (Importer()!=0) {
03135       EPETRA_CHK_ERR(ImportVector_->Import(X, *Importer(), Insert));
03136       Xp = (double**)ImportVector_->Pointers();
03137       LDX = ImportVector_->ConstantStride() ? ImportVector_->Stride() : 0;
03138 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
03139       if(Epetra_CrsMatrixTraceDumpMultiply) {
03140         *out << "\nColMap =\n";
03141         this->ColMap().Print(Teuchos::OSTab(out).o());
03142         *out << "\nX after import from DomainMap to ColMap =\n\n";
03143         ImportVector_->Print(Teuchos::OSTab(out).o());
03144       }
03145 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
03146     }
03147 
03148     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
03149     if (Exporter()!=0) {
03150       Yp = (double**)ExportVector_->Pointers();
03151       LDY = ExportVector_->ConstantStride() ? ExportVector_->Stride() : 0;
03152     }
03153 
03154     // Do actual computation
03155     if (NumVectors==1)
03156       GeneralMV(*Xp, *Yp);
03157     else
03158       GeneralMM(Xp, LDX, Yp, LDY, NumVectors);
03159 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
03160     if(Epetra_CrsMatrixTraceDumpMultiply) {
03161       *out << "\nRowMap =\n";
03162       this->RowMap().Print(Teuchos::OSTab(out).o());
03163       *out << "\nY after local mat-vec where Y has RowMap =\n\n";
03164       if(Exporter()!=0)
03165         ExportVector_->Print(Teuchos::OSTab(out).o());
03166       else
03167         Y.Print(Teuchos::OSTab(out).o());
03168     }
03169 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
03170     if (Exporter()!=0) {
03171       Y.PutScalar(0.0);  // Make sure target is zero
03172       Y.Export(*ExportVector_, *Exporter(), Add); // Fill Y with Values from export vector
03173 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
03174       if(Epetra_CrsMatrixTraceDumpMultiply) {
03175         *out << "\nRangeMap =\n";
03176         this->RangeMap().Print(Teuchos::OSTab(out).o());
03177         *out << "\nY after export from RowMap to RangeMap = \n\n";
03178         Y.Print(Teuchos::OSTab(out).o());
03179       }
03180 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
03181     }
03182     // Handle case of rangemap being a local replicated map
03183     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Y.Reduce());
03184   }
03185   else { // Transpose operation
03186 
03187     // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
03188 
03189     if (Exporter()!=0) {
03190       EPETRA_CHK_ERR(ExportVector_->Import(X, *Exporter(), Insert));
03191       Xp = (double**)ExportVector_->Pointers();
03192       LDX = ExportVector_->ConstantStride() ? ExportVector_->Stride() : 0;
03193 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
03194       if(Epetra_CrsMatrixTraceDumpMultiply) {
03195         *out << "\nRowMap =\n";
03196         this->RowMap().Print(Teuchos::OSTab(out).o());
03197         *out << "\nX after import from RangeMap to RowMap =\n\n";
03198         ExportVector_->Print(Teuchos::OSTab(out).o());
03199       }
03200 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
03201     }
03202 
03203     // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
03204     if (Importer()!=0) {
03205       Yp = (double**)ImportVector_->Pointers();
03206       LDY = ImportVector_->ConstantStride() ? ImportVector_->Stride() : 0;
03207     }
03208 
03209     // Do actual computation
03210     if (NumVectors==1)
03211       GeneralMTV(*Xp, *Yp);
03212     else
03213       GeneralMTM(Xp, LDX, Yp, LDY, NumVectors);
03214 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
03215     if(Epetra_CrsMatrixTraceDumpMultiply) {
03216       *out << "\nColMap =\n";
03217       this->ColMap().Print(Teuchos::OSTab(out).o());
03218       *out << "\nY after local transpose mat-vec where Y has ColMap =\n\n";
03219       if(Importer()!=0)
03220         ImportVector_->Print(Teuchos::OSTab(out).o());
03221       else
03222         Y.Print(Teuchos::OSTab(out).o());
03223     }
03224 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
03225     if (Importer()!=0) {
03226       Y.PutScalar(0.0);  // Make sure target is zero
03227       EPETRA_CHK_ERR(Y.Export(*ImportVector_, *Importer(), Add)); // Fill Y with Values from export vector
03228 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
03229       if(Epetra_CrsMatrixTraceDumpMultiply) {
03230         *out << "\nDomainMap =\n";
03231         this->DomainMap().Print(Teuchos::OSTab(out).o());
03232         *out << "\nY after export from ColMap to DomainMap =\n\n";
03233         Y.Print(Teuchos::OSTab(out).o());
03234       }
03235 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
03236     }
03237     // Handle case of rangemap being a local replicated map
03238     if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1)  EPETRA_CHK_ERR(Y.Reduce());
03239   }
03240 
03241   UpdateFlops(2*NumVectors*NumGlobalNonzeros64());
03242   if (Xcopy!=0) {
03243     delete Xcopy;
03244     EPETRA_CHK_ERR(1); // Return positive code to alert the user about needing extra copy of X
03245     return(1);
03246   }
03247 
03248 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
03249   if(Epetra_CrsMatrixTraceDumpMultiply) {
03250     *out << "\nFinal output Y is the last Y printed above!\n";
03251     *out << "\nLeaving Epetra_CrsMatrix::Multipy("<<TransA<<",X,Y) ...\n";
03252   }
03253 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
03254 
03255   return(0);
03256 }
03257 //=======================================================================================================
03258 void Epetra_CrsMatrix::UpdateImportVector(int NumVectors) const {
03259   if(Importer() != 0) {
03260     if(ImportVector_ != 0) {
03261       if(ImportVector_->NumVectors() != NumVectors) {
03262   delete ImportVector_;
03263   ImportVector_= 0;
03264       }
03265     }
03266     if(ImportVector_ == 0)
03267       ImportVector_ = new Epetra_MultiVector(ColMap(),NumVectors); // Create import vector if needed
03268   }
03269   return;
03270 }
03271 //=======================================================================================================
03272 void Epetra_CrsMatrix::UpdateExportVector(int NumVectors) const {
03273   if(Exporter() != 0) {
03274     if(ExportVector_ != 0) {
03275       if(ExportVector_->NumVectors() != NumVectors) {
03276   delete ExportVector_;
03277   ExportVector_= 0;
03278       }
03279     }
03280     if(ExportVector_ == 0)
03281       ExportVector_ = new Epetra_MultiVector(RowMap(),NumVectors); // Create Export vector if needed
03282   }
03283   return;
03284 }
03285 //=======================================================================================================
03286 void Epetra_CrsMatrix::GeneralMV(double * x, double * y)  const {
03287 
03288 if (StorageOptimized() && Graph().StorageOptimized()) {
03289 
03290   double * values = All_Values();
03291   int * Indices = Graph().All_Indices();
03292   int * IndexOffset = Graph().IndexOffset();
03293 #if defined(Epetra_ENABLE_MKL_SPARSE)
03294     char transa = 'n';
03295     int m = NumMyRows_;
03296     int NumCols = NumMyCols();
03297     double alpha = 1, beta = 0;
03298     // MKL should ignore '/'. G = General, C = 0-based indexing.
03299     char matdescra[6] = "G//C/";
03300     mkl_dcsrmv(&transa, &m, &NumCols, &alpha, matdescra, values, Indices, IndexOffset, IndexOffset + 1, x, &beta, y);
03301 #elif defined(EPETRA_HAVE_OMP)
03302   const int numMyRows = NumMyRows_;
03303 #pragma omp parallel for default(none) shared(IndexOffset,values,Indices,y,x)
03304      for (int row=0; row<numMyRows; ++row)
03305         {
03306      const int curOffset = IndexOffset[row];
03307           const double *val_ptr    = values+curOffset;
03308           const int    *colnum_ptr = Indices+curOffset;
03309      double s = 0.;
03310      const double *const val_end_of_row = &values[IndexOffset[row+1]];
03311      while (val_ptr != val_end_of_row)
03312        s += *val_ptr++ * x[*colnum_ptr++];
03313      y[row] = s;
03314   }
03315 #elif defined(Epetra_ENABLE_CASK)
03316        cask_csr_dax_new(NumMyRows_, IndexOffset, Indices,
03317                         values, x, y, cask);
03318 #else
03319        const double *val_ptr    = values;
03320        const int    *colnum_ptr = Indices;
03321        double       * dst_ptr = y;
03322        for (int row=0; row<NumMyRows_; ++row)
03323    {
03324      double s = 0.;
03325      const double *const val_end_of_row = &values[IndexOffset[row+1]];
03326      while (val_ptr != val_end_of_row)
03327        s += *val_ptr++ * x[*colnum_ptr++];
03328      *dst_ptr++ = s;
03329    }
03330 #endif // Epetra_ENABLE_CASK or EPETRA_HAVE_OMP or Epetra_ENABLE_MKL_SPARSE
03331 
03332     return;
03333   }
03334   else if (!StorageOptimized() && !Graph().StorageOptimized()) {
03335 
03336 
03337     int* NumEntriesPerRow = Graph().NumIndicesPerRow();
03338     int** Indices = Graph().Indices();
03339     double** srcValues = Values();
03340         const int numMyRows = NumMyRows_;
03341 
03342     // Do actual computation
03343 #ifdef EPETRA_HAVE_OMP
03344 #pragma omp parallel for default(none) shared(NumEntriesPerRow,Indices,srcValues,y,x)
03345 #endif
03346     for(int i = 0; i < numMyRows; i++) {
03347       int     NumEntries = NumEntriesPerRow[i];
03348       int*    RowIndices = Indices[i];
03349       double* RowValues  = srcValues[i];
03350       double sum = 0.0;
03351       for(int j = 0; j < NumEntries; j++)
03352   sum += *RowValues++ * x[*RowIndices++];
03353 
03354       y[i] = sum;
03355 
03356     }
03357   }
03358   else { // Case where StorageOptimized is incompatible:  Use general accessors.
03359 
03360     const int numMyRows = NumMyRows_;
03361 
03362     // Do actual computation
03363 #ifdef EPETRA_HAVE_OMP
03364 #pragma omp parallel for default(none) shared(x,y)
03365 #endif
03366     for(int i = 0; i < numMyRows; i++) {
03367       int     NumEntries = NumMyEntries(i);
03368       int*    RowIndices = Graph().Indices(i);
03369       double* RowValues  = Values(i);
03370       double sum = 0.0;
03371       for(int j = 0; j < NumEntries; j++)
03372   sum += *RowValues++ * x[*RowIndices++];
03373 
03374       y[i] = sum;
03375 
03376     }
03377   }
03378   return;
03379 }
03380 //=======================================================================================================
03381 void Epetra_CrsMatrix::GeneralMTV(double * x, double * y) const {
03382 
03383 #if !defined(FORTRAN_DISABLED) || defined(Epetra_ENABLE_CASK) || defined(Epetra_ENABLE_MKL_SPARSE)
03384   if (StorageOptimized() && Graph().StorageOptimized()) {
03385     double * values = All_Values_;
03386     int * Indices = Graph().All_Indices();
03387     int * IndexOffset = Graph().IndexOffset();
03388     int NumCols = NumMyCols();
03389 #if defined(Epetra_ENABLE_MKL_SPARSE)
03390     char transa = 't';
03391     int m = NumMyRows_;
03392     double alpha = 1, beta = 0;
03393     // MKL should ignore '/'. G = General, C = 0-based indexing.
03394     char matdescra[6] = "G//C/";
03395     mkl_dcsrmv(&transa, &m, &NumCols, &alpha, matdescra, values, Indices, IndexOffset, IndexOffset + 1, x, &beta, y);
03396 #elif defined(Epetra_ENABLE_CASK)
03397    cask_csr_datx( NumMyRows_, NumCols, IndexOffset,  Indices,  values,x ,y );
03398 #else
03399    int ione = 1;
03400    EPETRA_DCRSMV_F77(&ione, &NumMyRows_, &NumCols, values, Indices, IndexOffset, x, y);
03401 #endif
03402 
03403     return;
03404   }
03405 #endif // FORTRAN_DISABLED etc
03406   int NumCols = NumMyCols();
03407   for(int i = 0; i < NumCols; i++)
03408     y[i] = 0.0; // Initialize y for transpose multiply
03409 
03410   if (StorageOptimized() && Graph().StorageOptimized()) {
03411     double * values = All_Values_;
03412     int * Indices = Graph().All_Indices();
03413     int * IndexOffset = Graph().IndexOffset();
03414     for(int i = 0; i < NumMyRows_; ++i) {
03415       int prevOffset = *IndexOffset++;
03416       int NumEntries = *IndexOffset - prevOffset;
03417       double xi = x[i];
03418       for(int j = 0; j < NumEntries; j++)
03419   y[*Indices++] += *values++ * xi;
03420     }
03421   }
03422   else if (!StorageOptimized() && !Graph().StorageOptimized()) {
03423 
03424     int* NumEntriesPerRow = Graph().NumIndicesPerRow();
03425     int** Indices = Graph().Indices();
03426     double** srcValues = Values();
03427 
03428     for(int i = 0; i < NumMyRows_; i++) {
03429       int     NumEntries = *NumEntriesPerRow++;
03430       int*    RowIndices = *Indices++;
03431       double* RowValues  = *srcValues++;
03432       double xi = x[i];
03433       for(int j = 0; j < NumEntries; j++)
03434   y[*RowIndices++] += *RowValues++ * xi;
03435     }
03436   }
03437   else { // Case where StorageOptimized is incompatible:  Use general accessors.
03438 
03439     for(int i = 0; i < NumMyRows_; i++) {
03440       int     NumEntries = NumMyEntries(i);
03441       int*    RowIndices = Graph().Indices(i);
03442       double* RowValues  = Values(i);
03443       double xi = x[i];
03444       for(int j = 0; j < NumEntries; j++)
03445   y[*RowIndices++] += *RowValues++ * xi;
03446     }
03447   }
03448 
03449   return;
03450 }
03451 //=======================================================================================================
03452 void Epetra_CrsMatrix::GeneralMM(double ** X, int LDX, double ** Y, int LDY, int NumVectors) const {
03453 
03454 #if !defined(FORTRAN_DISABLED) || defined(Epetra_ENABLE_CASK) || (defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM))
03455   if (StorageOptimized() && Graph().StorageOptimized()) {
03456     double * values = All_Values_;
03457     int * Indices = Graph().All_Indices();
03458     int * IndexOffset = Graph().IndexOffset();
03459 
03460     if (LDX!=0 && LDY!=0) {
03461 #if defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM)
03462     int * IndicesPlus1 = Graph().All_IndicesPlus1();
03463     char transa = 'n';
03464     int NumRows = NumMyRows_;
03465     int NumCols = NumMyCols();
03466     double alpha = 1, beta = 0;
03467     // MKL should ignore '/'. G = General, C = 0-based indexing, F = 1-based.
03468     // 1-based is used because X and Y are column-oriented and MKL does not
03469     // do 0-based and column-oriented.
03470     char matdescra[6] = "G//F/";
03471     mkl_dcsrmm(&transa, &NumRows, &NumVectors, &NumCols, &alpha, matdescra, values, IndicesPlus1, IndexOffset, IndexOffset + 1, *X, &LDX, &beta, *Y, &LDY);
03472 #elif defined(Epetra_ENABLE_CASK)
03473     cask_csr_dgesmm_new(0, 1.0, NumMyRows_, NumMyRows_,  NumVectors,
03474                     IndexOffset, Indices, values, *X, LDX, 0.0,  *Y, LDY,cask);
03475 #else
03476     int izero = 0;
03477     EPETRA_DCRSMM_F77(&izero, &NumMyRows_, &NumMyRows_, values, Indices, IndexOffset, *X, &LDX, *Y, &LDY, &NumVectors);
03478 #endif
03479     return;
03480     }
03481 
03482     double ** const xp = X;
03483     double ** const yp = Y;
03484     const int numMyRows = NumMyRows_;
03485 #ifdef EPETRA_HAVE_OMP
03486 #pragma omp parallel for default(none) shared(IndexOffset,Indices,values,NumVectors)
03487 #endif
03488     for (int i=0; i < numMyRows; i++) {
03489       int prevOffset = IndexOffset[i];
03490       int NumEntries = IndexOffset[i+1] - prevOffset;
03491       int *    RowIndices = Indices+prevOffset;
03492       double * RowValues  = values+prevOffset;
03493       for (int k=0; k<NumVectors; k++) {
03494   double sum = 0.0;
03495   const double * const x = xp[k];
03496   double * const y = yp[k];
03497   for (int j=0; j < NumEntries; j++) sum += RowValues[j] * x[RowIndices[j]];
03498   y[i] = sum;
03499       }
03500     }
03501   }
03502   else if (!StorageOptimized() && !Graph().StorageOptimized()) {
03503 #else
03504   if (!StorageOptimized() && !Graph().StorageOptimized()) {
03505 #endif // FORTRAN_DISABLED etc
03506 
03507     int* NumEntriesPerRow = Graph().NumIndicesPerRow();
03508     int** Indices = Graph().Indices();
03509     double** srcValues = Values();
03510     double ** const xp = X;
03511     double ** const yp = Y;
03512     const int numMyRows = NumMyRows_;
03513 
03514 #ifdef EPETRA_HAVE_OMP
03515 #pragma omp parallel for default(none) shared(NumEntriesPerRow,Indices,srcValues,NumVectors)
03516 #endif
03517     for (int i=0; i < numMyRows; i++) {
03518       int      NumEntries = NumEntriesPerRow[i];
03519       int *    RowIndices = Indices[i];
03520       double * RowValues  = srcValues[i];
03521       for (int k=0; k<NumVectors; k++) {
03522   double sum = 0.0;
03523   const double * const x = xp[k];
03524   double * const y = yp[k];
03525   for (int j=0; j < NumEntries; j++) sum += RowValues[j] * x[RowIndices[j]];
03526   y[i] = sum;
03527       }
03528     }
03529   }
03530   else {
03531 
03532     double ** const xp = X;
03533     double ** const yp = Y;
03534     const int numMyRows = NumMyRows_;
03535 #ifdef EPETRA_HAVE_OMP
03536 #pragma omp parallel for default(none) shared(NumVectors)
03537 #endif
03538     for (int i=0; i < numMyRows; i++) {
03539       int     NumEntries = NumMyEntries(i);
03540       int*    RowIndices = Graph().Indices(i);
03541       double* RowValues  = Values(i);
03542       for (int k=0; k<NumVectors; k++) {
03543   double sum = 0.0;
03544   double * x = xp[k];
03545   for (int j=0; j < NumEntries; j++) sum += RowValues[j] * x[RowIndices[j]];
03546   yp[k][i] = sum;
03547       }
03548     }
03549   }
03550   return;
03551 }
03552 //=======================================================================================================
03553 void Epetra_CrsMatrix::GeneralMTM(double ** X, int LDX, double ** Y, int LDY, int NumVectors)  const{
03554 
03555   int NumCols = NumMyCols();
03556 #if !defined(FORTRAN_DISABLED) || defined(Epetra_ENABLE_CASK) || (defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM))
03557   if (StorageOptimized() && Graph().StorageOptimized()) {
03558     if (LDX!=0 && LDY!=0) {
03559       double * values = All_Values_;
03560       int * Indices = Graph().All_Indices();
03561       int * IndexOffset = Graph().IndexOffset();
03562 #if defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM)
03563       int * IndicesPlus1 = Graph().All_IndicesPlus1();
03564       char transa = 't';
03565       int NumRows = NumMyRows_;
03566       int NumCols = NumMyCols();
03567       double alpha = 1, beta = 0;
03568       // MKL should ignore '/'. G = General, C = 0-based indexing, F = 1-based.
03569       // 1-based is used because X and Y are column-oriented and MKL does not
03570       // do 0-based and column-oriented.
03571       char matdescra[6] = "G//F/";
03572       mkl_dcsrmm(&transa, &NumRows, &NumVectors, &NumCols, &alpha, matdescra, values, IndicesPlus1, IndexOffset, IndexOffset + 1, *X, &LDX, &beta, *Y, &LDY);
03573 #elif defined(Epetra_ENABLE_CASK)
03574       cask_csr_dgesmm_new(1, 1.0, NumMyRows_, NumCols,  NumVectors,
03575                           IndexOffset, Indices, values, *X, LDX, 0.0,
03576                           *Y, LDY, cask);
03577 #else
03578       int ione = 1;
03579       EPETRA_DCRSMM_F77(&ione, &NumMyRows_, &NumCols, values, Indices, IndexOffset, *X, &LDX, *Y, &LDY, &NumVectors);
03580 #endif
03581       return;
03582     }
03583   }
03584 #endif // FORTRAN_DISABLED etc
03585   for (int k=0; k<NumVectors; k++)
03586     for (int i=0; i < NumCols; i++)
03587       Y[k][i] = 0.0; // Initialize y for transpose multiply
03588 
03589   if (StorageOptimized() && Graph().StorageOptimized()) {
03590     double * values = All_Values_;
03591     int * Indices = Graph().All_Indices();
03592     int * IndexOffset = Graph().IndexOffset();
03593     for (int i=0; i < NumMyRows_; i++) {
03594       int prevOffset = *IndexOffset++;
03595       int NumEntries = *IndexOffset - prevOffset;
03596       int *    RowIndices = Indices+prevOffset;
03597       double * RowValues  = values+prevOffset;
03598 
03599       for (int k=0; k<NumVectors; k++) {
03600   double * y = Y[k];
03601   double * x = X[k];
03602   for (int j=0; j < NumEntries; j++)
03603     y[RowIndices[j]] += RowValues[j] * x[i];
03604       }
03605     }
03606   }
03607   else if (!StorageOptimized() && !Graph().StorageOptimized()) {
03608 
03609     int* NumEntriesPerRow = Graph().NumIndicesPerRow();
03610     int** Indices = Graph().Indices();
03611     double** srcValues = Values();
03612 
03613     for (int i=0; i < NumMyRows_; i++) {
03614       int      NumEntries = *NumEntriesPerRow++;
03615       int *    RowIndices = *Indices++;
03616       double * RowValues  = *srcValues++;
03617       for (int k=0; k<NumVectors; k++) {
03618   double * y = Y[k];
03619   double * x = X[k];
03620   for (int j=0; j < NumEntries; j++)
03621     y[RowIndices[j]] += RowValues[j] * x[i];
03622       }
03623     }
03624   }
03625   else { // Case where StorageOptimized is incompatible:  Use general accessors.
03626 
03627     for (int i=0; i < NumMyRows_; i++) {
03628       int     NumEntries = NumMyEntries(i);
03629       int*    RowIndices = Graph().Indices(i);
03630       double* RowValues  = Values(i);
03631       for (int k=0; k<NumVectors; k++) {
03632   double * y = Y[k];
03633   double * x = X[k];
03634   for (int j=0; j < NumEntries; j++)
03635     y[RowIndices[j]] += RowValues[j] * x[i];
03636       }
03637     }
03638   }
03639   return;
03640 }
03641 //=======================================================================================================
03642 void Epetra_CrsMatrix::GeneralSV(bool Upper, bool Trans, bool UnitDiagonal, double * xp, double * yp)  const {
03643 
03644 
03645   int i, j, j0;
03646 
03647 #if !defined(FORTRAN_DISABLED) || defined(Epetra_ENABLE_CASK) || defined(Epetra_ENABLE_MKL_SPARSE)
03648   if (StorageOptimized() && Graph().StorageOptimized() && ((UnitDiagonal && NoDiagonal())|| (!UnitDiagonal && !NoDiagonal()))) {
03649     double * values = All_Values();
03650     int * Indices = Graph().All_Indices();
03651     int * IndexOffset = Graph().IndexOffset();
03652 
03653 #if !defined(Epetra_ENABLE_MKL_SPARSE)
03654     int iupper = Upper ? 1:0;
03655     int itrans = Trans ? 1:0;
03656     int udiag =  UnitDiagonal ? 1:0;
03657     int nodiag = NoDiagonal() ? 1:0;
03658     int xysame = (xp==yp) ? 1:0;
03659 #endif
03660 
03661 #if defined(Epetra_ENABLE_MKL_SPARSE)
03662     char transa = Trans ? 't' : 'n';
03663     int m = NumMyRows_;
03664     double alpha = 1;
03665     // T = Triangular, C = 0-based indexing. '/' should be ignored by MKL
03666     char matdescra[6] = {'T', Upper ? 'U' : 'L', UnitDiagonal ? 'U' : 'N', 'C', '/', '\0'};
03667     mkl_dcsrsv(&transa, &m, &alpha, matdescra, values, Indices, IndexOffset, IndexOffset + 1, xp, yp);
03668 #elif defined(Epetra_ENABLE_CASK)
03669     cask_csr_dtrsv_new( iupper, itrans, udiag, nodiag, 0, xysame, NumMyRows_,
03670                     NumMyRows_, IndexOffset, Indices, values, xp, yp, cask);
03671 #else
03672     EPETRA_DCRSSV_F77( &iupper, &itrans, &udiag, &nodiag, &NumMyRows_, &NumMyRows_, values, Indices, IndexOffset, xp, yp, &xysame);
03673 #endif
03674     return;
03675   }
03676   //=================================================================
03677   else { // !StorageOptimized()
03678   //=================================================================
03679 #endif
03680     if (!Trans) {
03681 
03682       if (Upper) {
03683 
03684   j0 = 1;
03685   if (NoDiagonal())
03686     j0--; // Include first term if no diagonal
03687   for (i=NumMyRows_-1; i >=0; i--) {
03688     int      NumEntries = NumMyEntries(i);
03689     int *    RowIndices = Graph().Indices(i);
03690     double * RowValues  = Values(i);
03691     double sum = 0.0;
03692     for (j=j0; j < NumEntries; j++)
03693       sum += RowValues[j] * yp[RowIndices[j]];
03694 
03695     if (UnitDiagonal)
03696       yp[i] = xp[i] - sum;
03697     else
03698       yp[i] = (xp[i] - sum)/RowValues[0];
03699 
03700   }
03701       }
03702       else {
03703   j0 = 1;
03704   if (NoDiagonal())
03705     j0--; // Include first term if no diagonal
03706   for (i=0; i < NumMyRows_; i++) {
03707     int      NumEntries = NumMyEntries(i) - j0;
03708     int *    RowIndices = Graph().Indices(i);
03709     double * RowValues  = Values(i);
03710     double sum = 0.0;
03711     for (j=0; j < NumEntries; j++)
03712       sum += RowValues[j] * yp[RowIndices[j]];
03713 
03714     if (UnitDiagonal)
03715       yp[i] = xp[i] - sum;
03716     else
03717       yp[i] = (xp[i] - sum)/RowValues[NumEntries];
03718 
03719   }
03720       }
03721     }
03722 
03723     // ***********  Transpose case *******************************
03724 
03725     else {
03726 
03727       if (xp!=yp)
03728   for (i=0; i < NumMyRows_; i++)
03729     yp[i] = xp[i]; // Initialize y for transpose solve
03730 
03731       if (Upper) {
03732 
03733   j0 = 1;
03734   if (NoDiagonal())
03735     j0--; // Include first term if no diagonal
03736 
03737   for (i=0; i < NumMyRows_; i++) {
03738     int      NumEntries = NumMyEntries(i);
03739     int *    RowIndices = Graph().Indices(i);
03740     double * RowValues  = Values(i);
03741     if (!UnitDiagonal)
03742       yp[i] = yp[i]/RowValues[0];
03743     double ytmp = yp[i];
03744     for (j=j0; j < NumEntries; j++)
03745       yp[RowIndices[j]] -= RowValues[j] * ytmp;
03746   }
03747       }
03748       else {
03749 
03750   j0 = 1;
03751   if (NoDiagonal())
03752     j0--; // Include first term if no diagonal
03753 
03754   for (i=NumMyRows_-1; i >= 0; i--) {
03755     int      NumEntries = NumMyEntries(i) - j0;
03756     int *    RowIndices = Graph().Indices(i);
03757     double * RowValues  = Values(i);
03758     if (!UnitDiagonal)
03759       yp[i] = yp[i]/RowValues[NumEntries];
03760     double ytmp = yp[i];
03761     for (j=0; j < NumEntries; j++)
03762       yp[RowIndices[j]] -= RowValues[j] * ytmp;
03763   }
03764       }
03765 
03766     }
03767 #if !defined(FORTRAN_DISABLED) || defined(Epetra_ENABLE_CASK) || defined(Epetra_ENABLE_MKL_SPARSE)
03768   }
03769 #endif
03770   return;
03771 }
03772 //=======================================================================================================
03773 void Epetra_CrsMatrix::GeneralSM(bool Upper, bool Trans, bool UnitDiagonal, double ** Xp, int LDX, double ** Yp, int LDY, int NumVectors)  const{
03774 
03775   int i, j, j0, k;
03776   double diag = 0.0;
03777 
03778   if (StorageOptimized() && Graph().StorageOptimized()) {
03779     double * values = All_Values();
03780     int * Indices = Graph().All_Indices();
03781     int * IndexOffset = Graph().IndexOffset();
03782 #if !defined(FORTRAN_DISABLED) || defined(Epetra_ENABLE_CASK) || (defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM))
03783     if (LDX!=0 && LDY!=0 && ((UnitDiagonal && NoDiagonal()) || (!UnitDiagonal && !NoDiagonal()))) {
03784 
03785 #if !defined(Epetra_ENABLE_MKL_SPARSE) || defined(Epetra_DISABLE_MKL_SPARSE_MM)
03786       int iupper = Upper ? 1:0;
03787       int itrans = Trans ? 1:0;
03788       int udiag =  UnitDiagonal ? 1:0;
03789       int nodiag = NoDiagonal() ? 1:0;
03790       int xysame = (Xp==Yp) ? 1:0;
03791 #endif
03792 
03793 #if defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM)
03794       int * IndicesPlus1 = Graph().All_IndicesPlus1();
03795       char transa = Trans ? 't' : 'n';
03796       int NumRows = NumMyRows_;
03797       double alpha = 1;
03798       // T = Triangular, C = 0-based indexing, F = 1-based. '/' should be ignored by MKL.
03799       // 1-based is used because X and Y are column-oriented and MKL does not
03800       // do 0-based and column-oriented.
03801       char matdescra[6] = {'T', Upper ? 'U' : 'L', UnitDiagonal ? 'U' : 'N', 'F', '/', '\0'};
03802       mkl_dcsrsm(&transa, &NumRows, &NumVectors, &alpha, matdescra, values, IndicesPlus1, IndexOffset, IndexOffset + 1, *Xp, &LDX, *Yp, &LDY);
03803 #elif defined(Epetra_ENABLE_CASK)
03804       cask_csr_dtrsm( iupper, itrans, udiag, nodiag, 0, xysame,  NumMyRows_,
03805                       NumMyRows_, NumVectors, IndexOffset, Indices, values,
03806                       *Xp, LDX, *Yp, LDY);
03807 #else
03808       EPETRA_DCRSSM_F77( &iupper, &itrans, &udiag, &nodiag, &NumMyRows_, &NumMyRows_, values, Indices, IndexOffset,
03809        *Xp, &LDX, *Yp, &LDY, &xysame, &NumVectors);
03810 #endif
03811       return;
03812     }
03813 #endif // FORTRAN_DISABLED etc
03814     if(!Trans) {
03815       if(Upper) {
03816   j0 = 1;
03817   if(NoDiagonal())
03818     j0--; // Include first term if no diagonal
03819   for(i = NumMyRows_ - 1; i >= 0; i--) {
03820     int Offset = IndexOffset[i];
03821     int      NumEntries = IndexOffset[i+1]-Offset;
03822     int *    RowIndices = Indices+Offset;
03823     double * RowValues  = values+Offset;
03824     if(!UnitDiagonal)
03825       diag = 1.0/RowValues[0]; // Take inverse of diagonal once for later use
03826     for(k = 0; k < NumVectors; k++) {
03827       double sum = 0.0;
03828       for(j = j0; j < NumEntries; j++)
03829         sum += RowValues[j] * Yp[k][RowIndices[j]];
03830 
03831       if(UnitDiagonal)
03832         Yp[k][i] = Xp[k][i] - sum;
03833       else
03834         Yp[k][i] = (Xp[k][i] - sum) * diag;
03835     }
03836   }
03837       }
03838       else {
03839   j0 = 1;
03840   if(NoDiagonal())
03841     j0--; // Include first term if no diagonal
03842   for(i = 0; i < NumMyRows_; i++) {
03843     int Offset = IndexOffset[i];
03844     int      NumEntries = IndexOffset[i+1]-Offset - j0;
03845     int *    RowIndices = Indices+Offset;
03846     double * RowValues  = values+Offset;
03847     if(!UnitDiagonal)
03848       diag = 1.0/RowValues[NumEntries]; // Take inverse of diagonal once for later use
03849     for(k = 0; k < NumVectors; k++) {
03850       double sum = 0.0;
03851       for(j = 0; j < NumEntries; j++)
03852         sum += RowValues[j] * Yp[k][RowIndices[j]];
03853 
03854       if(UnitDiagonal)
03855         Yp[k][i] = Xp[k][i] - sum;
03856       else
03857         Yp[k][i] = (Xp[k][i] - sum)*diag;
03858     }
03859   }
03860       }
03861     }
03862     // ***********  Transpose case *******************************
03863 
03864     else {
03865       for(k = 0; k < NumVectors; k++)
03866   if(Yp[k] != Xp[k])
03867     for(i = 0; i < NumMyRows_; i++)
03868       Yp[k][i] = Xp[k][i]; // Initialize y for transpose multiply
03869 
03870       if(Upper) {
03871   j0 = 1;
03872   if(NoDiagonal())
03873     j0--; // Include first term if no diagonal
03874 
03875   for(i = 0; i < NumMyRows_; i++) {
03876     int Offset = IndexOffset[i];
03877     int      NumEntries = IndexOffset[i+1]-Offset;
03878     int *    RowIndices = Indices+Offset;
03879     double * RowValues  = values+Offset;
03880     if(!UnitDiagonal)
03881       diag = 1.0/RowValues[0]; // Take inverse of diagonal once for later use
03882     for(k = 0; k < NumVectors; k++) {
03883       if(!UnitDiagonal)
03884         Yp[k][i] = Yp[k][i]*diag;
03885       double ytmp = Yp[k][i];
03886       for(j = j0; j < NumEntries; j++)
03887         Yp[k][RowIndices[j]] -= RowValues[j] * ytmp;
03888     }
03889   }
03890       }
03891       else {
03892   j0 = 1;
03893   if(NoDiagonal())
03894     j0--; // Include first term if no diagonal
03895   for(i = NumMyRows_ - 1; i >= 0; i--) {
03896     int Offset = IndexOffset[i];
03897     int      NumEntries = IndexOffset[i+1]-Offset - j0;
03898     int *    RowIndices = Indices+Offset;
03899     double * RowValues  = values+Offset;
03900     if(!UnitDiagonal)
03901       diag = 1.0/RowValues[NumEntries]; // Take inverse of diagonal once for later use
03902     for(k = 0; k < NumVectors; k++) {
03903       if(!UnitDiagonal)
03904         Yp[k][i] = Yp[k][i]*diag;
03905       double ytmp = Yp[k][i];
03906       for(j = 0; j < NumEntries; j++)
03907         Yp[k][RowIndices[j]] -= RowValues[j] * ytmp;
03908     }
03909   }
03910       }
03911     }
03912   }
03913     // ========================================================
03914   else { // !StorageOptimized()
03915     // ========================================================
03916 
03917     if(!Trans) {
03918       if(Upper) {
03919   j0 = 1;
03920   if(NoDiagonal())
03921     j0--; // Include first term if no diagonal
03922   for(i = NumMyRows_ - 1; i >= 0; i--) {
03923     int     NumEntries = NumMyEntries(i);
03924     int*    RowIndices = Graph().Indices(i);
03925     double* RowValues  = Values(i);
03926     if(!UnitDiagonal)
03927       diag = 1.0/RowValues[0]; // Take inverse of diagonal once for later use
03928     for(k = 0; k < NumVectors; k++) {
03929       double sum = 0.0;
03930       for(j = j0; j < NumEntries; j++)
03931         sum += RowValues[j] * Yp[k][RowIndices[j]];
03932 
03933       if(UnitDiagonal)
03934         Yp[k][i] = Xp[k][i] - sum;
03935       else
03936         Yp[k][i] = (Xp[k][i] - sum) * diag;
03937     }
03938   }
03939       }
03940       else {
03941   j0 = 1;
03942   if(NoDiagonal())
03943     j0--; // Include first term if no diagonal
03944   for(i = 0; i < NumMyRows_; i++) {
03945     int     NumEntries = NumMyEntries(i) - j0;
03946     int*    RowIndices = Graph().Indices(i);
03947     double* RowValues  = Values(i);
03948     if(!UnitDiagonal)
03949       diag = 1.0/RowValues[NumEntries]; // Take inverse of diagonal once for later use
03950     for(k = 0; k < NumVectors; k++) {
03951       double sum = 0.0;
03952       for(j = 0; j < NumEntries; j++)
03953         sum += RowValues[j] * Yp[k][RowIndices[j]];
03954 
03955       if(UnitDiagonal)
03956         Yp[k][i] = Xp[k][i] - sum;
03957       else
03958         Yp[k][i] = (Xp[k][i] - sum)*diag;
03959     }
03960   }
03961       }
03962     }
03963     // ***********  Transpose case *******************************
03964 
03965     else {
03966       for(k = 0; k < NumVectors; k++)
03967   if(Yp[k] != Xp[k])
03968     for(i = 0; i < NumMyRows_; i++)
03969       Yp[k][i] = Xp[k][i]; // Initialize y for transpose multiply
03970 
03971       if(Upper) {
03972   j0 = 1;
03973   if(NoDiagonal())
03974     j0--; // Include first term if no diagonal
03975 
03976   for(i = 0; i < NumMyRows_; i++) {
03977     int     NumEntries = NumMyEntries(i);
03978     int*    RowIndices = Graph().Indices(i);
03979     double* RowValues  = Values(i);
03980     if(!UnitDiagonal)
03981       diag = 1.0/RowValues[0]; // Take inverse of diagonal once for later use
03982     for(k = 0; k < NumVectors; k++) {
03983       if(!UnitDiagonal)
03984         Yp[k][i] = Yp[k][i]*diag;
03985       double ytmp = Yp[k][i];
03986       for(j = j0; j < NumEntries; j++)
03987         Yp[k][RowIndices[j]] -= RowValues[j] * ytmp;
03988     }
03989   }
03990       }
03991       else {
03992   j0 = 1;
03993   if(NoDiagonal())
03994     j0--; // Include first term if no diagonal
03995   for(i = NumMyRows_ - 1; i >= 0; i--) {
03996     int     NumEntries = NumMyEntries(i) - j0;
03997     int*    RowIndices = Graph().Indices(i);
03998     double* RowValues  = Values(i);
03999     if(!UnitDiagonal)
04000       diag = 1.0/RowValues[NumEntries]; // Take inverse of diagonal once for later use
04001     for(k = 0; k < NumVectors; k++) {
04002       if(!UnitDiagonal)
04003         Yp[k][i] = Yp[k][i]*diag;
04004       double ytmp = Yp[k][i];
04005       for(j = 0; j < NumEntries; j++)
04006         Yp[k][RowIndices[j]] -= RowValues[j] * ytmp;
04007     }
04008   }
04009       }
04010     }
04011   }
04012 
04013   return;
04014 }
04015 //=============================================================================
04016 int Epetra_CrsMatrix::Multiply1(bool TransA, const Epetra_Vector& x, Epetra_Vector& y) const {
04017 
04018 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS
04019   TEUCHOS_FUNC_TIME_MONITOR("Epetra_CrsMatrix::Multiply1(TransA,x,y)");
04020 #endif
04021 
04022   //
04023   // This function forms the product y = A * x or y = A' * x
04024   //
04025 
04026   if(!Filled())
04027     EPETRA_CHK_ERR(-1); // Matrix must be filled.
04028 
04029   int i, j;
04030   double* xp = (double*) x.Values();
04031   double* yp = (double*) y.Values();
04032   int NumMyCols_ = NumMyCols();
04033 
04034   if(!TransA) {
04035 
04036     // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
04037     if(Importer() != 0) {
04038       if(ImportVector_ != 0) {
04039   if(ImportVector_->NumVectors() != 1) {
04040     delete ImportVector_;
04041     ImportVector_= 0;
04042   }
04043       }
04044       if(ImportVector_ == 0)
04045   ImportVector_ = new Epetra_MultiVector(ColMap(),1); // Create import vector if needed
04046       EPETRA_CHK_ERR(ImportVector_->Import(x, *Importer(), Insert));
04047       xp = (double*) ImportVector_->Values();
04048     }
04049 
04050     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
04051     if(Exporter() != 0) {
04052       if(ExportVector_ != 0) {
04053   if(ExportVector_->NumVectors() != 1) {
04054     delete ExportVector_;
04055     ExportVector_= 0;
04056   }
04057       }
04058       if(ExportVector_ == 0)
04059   ExportVector_ = new Epetra_MultiVector(RowMap(),1); // Create Export vector if needed
04060       yp = (double*) ExportVector_->Values();
04061     }
04062 
04063     // Do actual computation
04064     for(i = 0; i < NumMyRows_; i++) {
04065       int     NumEntries = NumMyEntries(i);
04066       int*    RowIndices = Graph().Indices(i);
04067       double* RowValues  = Values(i);
04068       double sum = 0.0;
04069       for(j = 0; j < NumEntries; j++)
04070   sum += RowValues[j] * xp[RowIndices[j]];
04071 
04072       yp[i] = sum;
04073 
04074     }
04075     if(Exporter() != 0) {
04076       y.PutScalar(0.0); // Make sure target is zero
04077       EPETRA_CHK_ERR(y.Export(*ExportVector_, *Exporter(), Add)); // Fill y with Values from export vector
04078     }
04079     // Handle case of rangemap being a local replicated map
04080     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(y.Reduce());
04081   }
04082 
04083   else { // Transpose operation
04084 
04085     // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
04086     if(Exporter() != 0) {
04087       if(ExportVector_ != 0) {
04088   if(ExportVector_->NumVectors() != 1) {
04089     delete ExportVector_;
04090     ExportVector_= 0;
04091   }
04092       }
04093       if(ExportVector_ == 0)
04094   ExportVector_ = new Epetra_MultiVector(RowMap(),1); // Create Export vector if needed
04095       EPETRA_CHK_ERR(ExportVector_->Import(x, *Exporter(), Insert));
04096       xp = (double*) ExportVector_->Values();
04097     }
04098 
04099     // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
04100     if(Importer() != 0) {
04101       if(ImportVector_ != 0) {
04102   if(ImportVector_->NumVectors() != 1) {
04103     delete ImportVector_;
04104     ImportVector_= 0;
04105   }
04106       }
04107       if(ImportVector_ == 0)
04108   ImportVector_ = new Epetra_MultiVector(ColMap(),1); // Create import vector if needed
04109       yp = (double*) ImportVector_->Values();
04110     }
04111 
04112     // Do actual computation
04113     for(i = 0; i < NumMyCols_; i++)
04114       yp[i] = 0.0; // Initialize y for transpose multiply
04115 
04116     for(i = 0; i < NumMyRows_; i++) {
04117       int     NumEntries = NumMyEntries(i);
04118       int*    RowIndices = Graph().Indices(i);
04119       double* RowValues  = Values(i);
04120       for(j = 0; j < NumEntries; j++)
04121   yp[RowIndices[j]] += RowValues[j] * xp[i];
04122     }
04123     if(Importer() != 0) {
04124       y.PutScalar(0.0); // Make sure target is zero
04125       EPETRA_CHK_ERR(y.Export(*ImportVector_, *Importer(), Add)); // Fill y with Values from export vector
04126     }
04127     // Handle case of rangemap being a local replicated map
04128     if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(y.Reduce());
04129   }
04130 
04131   UpdateFlops(2 * NumGlobalNonzeros64());
04132   return(0);
04133 }
04134 //=============================================================================
04135 int Epetra_CrsMatrix::Multiply1(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
04136 
04137 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS
04138   TEUCHOS_FUNC_TIME_MONITOR("Epetra_CrsMatrix::Multiply1(TransA,X,Y)");
04139 #endif
04140 
04141   //
04142   // This function forms the product Y = A * Y or Y = A' * X
04143   //
04144   if((X.NumVectors() == 1) && (Y.NumVectors() == 1)) {
04145     double* xp = (double*) X[0];
04146     double* yp = (double*) Y[0];
04147     Epetra_Vector x(View, X.Map(), xp);
04148     Epetra_Vector y(View, Y.Map(), yp);
04149     EPETRA_CHK_ERR(Multiply1(TransA, x, y));
04150     return(0);
04151   }
04152   if(!Filled()) {
04153     EPETRA_CHK_ERR(-1); // Matrix must be filled.
04154   }
04155 
04156   int i, j, k;
04157 
04158   double** Xp = (double**) X.Pointers();
04159   double** Yp = (double**) Y.Pointers();
04160 
04161   int NumVectors = X.NumVectors();
04162   int NumMyCols_ = NumMyCols();
04163 
04164 
04165   // Need to better manage the Import and Export vectors:
04166   // - Need accessor functions
04167   // - Need to make the NumVector match (use a View to do this)
04168   // - Need to look at RightScale and ColSum routines too.
04169 
04170   if (!TransA) {
04171 
04172     // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
04173     if (Importer()!=0) {
04174       if (ImportVector_!=0) {
04175   if (ImportVector_->NumVectors()!=NumVectors) {
04176     delete ImportVector_; ImportVector_= 0;}
04177       }
04178       if (ImportVector_==0)
04179   ImportVector_ = new Epetra_MultiVector(ColMap(),NumVectors); // Create import vector if needed
04180       EPETRA_CHK_ERR(ImportVector_->Import(X, *Importer(), Insert));
04181       Xp = (double**)ImportVector_->Pointers();
04182     }
04183 
04184     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
04185     if (Exporter()!=0) {
04186       if (ExportVector_!=0) {
04187   if (ExportVector_->NumVectors()!=NumVectors) {
04188     delete ExportVector_; ExportVector_= 0;}
04189       }
04190       if (ExportVector_==0)
04191   ExportVector_ = new Epetra_MultiVector(RowMap(),NumVectors); // Create Export vector if needed
04192       Yp = (double**)ExportVector_->Pointers();
04193     }
04194 
04195     // Do actual computation
04196 
04197     for (i=0; i < NumMyRows_; i++) {
04198       int      NumEntries = NumMyEntries(i);
04199       int *    RowIndices = Graph().Indices(i);
04200       double * RowValues  = Values(i);
04201       for (k=0; k<NumVectors; k++) {
04202   double sum = 0.0;
04203   for (j=0; j < NumEntries; j++) sum += RowValues[j] * Xp[k][RowIndices[j]];
04204   Yp[k][i] = sum;
04205       }
04206     }
04207     if (Exporter()!=0) {
04208       Y.PutScalar(0.0); // Make sure target is zero
04209       Y.Export(*ExportVector_, *Exporter(), Add); // Fill Y with Values from export vector
04210     }
04211     // Handle case of rangemap being a local replicated map
04212     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Y.Reduce());
04213   }
04214   else { // Transpose operation
04215 
04216 
04217     // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
04218 
04219     if (Exporter()!=0) {
04220       if (ExportVector_!=0) {
04221   if (ExportVector_->NumVectors()!=NumVectors) {
04222     delete ExportVector_; ExportVector_= 0;}
04223       }
04224       if (ExportVector_==0)
04225   ExportVector_ = new Epetra_MultiVector(RowMap(),NumVectors); // Create Export vector if needed
04226       EPETRA_CHK_ERR(ExportVector_->Import(X, *Exporter(), Insert));
04227       Xp = (double**)ExportVector_->Pointers();
04228     }
04229 
04230     // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
04231     if (Importer()!=0) {
04232       if (ImportVector_!=0) {
04233   if (ImportVector_->NumVectors()!=NumVectors) {
04234     delete ImportVector_; ImportVector_= 0;}
04235       }
04236       if (ImportVector_==0)
04237   ImportVector_ = new Epetra_MultiVector(ColMap(),NumVectors); // Create import vector if needed
04238       Yp = (double**)ImportVector_->Pointers();
04239     }
04240 
04241     // Do actual computation
04242 
04243 
04244 
04245     for (k=0; k<NumVectors; k++)
04246       for (i=0; i < NumMyCols_; i++)
04247   Yp[k][i] = 0.0; // Initialize y for transpose multiply
04248 
04249     for (i=0; i < NumMyRows_; i++) {
04250       int      NumEntries = NumMyEntries(i);
04251       int *    RowIndices = Graph().Indices(i);
04252       double * RowValues  = Values(i);
04253       for (k=0; k<NumVectors; k++) {
04254   for (j=0; j < NumEntries; j++)
04255     Yp[k][RowIndices[j]] += RowValues[j] * Xp[k][i];
04256       }
04257     }
04258     if (Importer()!=0) {
04259       Y.PutScalar(0.0); // Make sure target is zero
04260       EPETRA_CHK_ERR(Y.Export(*ImportVector_, *Importer(), Add)); // Fill Y with Values from export vector
04261     }
04262     // Handle case of rangemap being a local replicated map
04263     if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1)  EPETRA_CHK_ERR(Y.Reduce());
04264   }
04265 
04266   UpdateFlops(2*NumVectors*NumGlobalNonzeros64());
04267   return(0);
04268 }
04269 
04270 //=============================================================================
04271 int Epetra_CrsMatrix::Solve1(bool Upper, bool Trans, bool UnitDiagonal,
04272           const Epetra_Vector& x, Epetra_Vector& y) const
04273 {
04274 
04275 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS
04276   TEUCHOS_FUNC_TIME_MONITOR("Epetra_CrsMatrix::Solve1(Upper,Trans,UnitDiag,x,y)");
04277 #endif
04278 
04279   //
04280   // This function finds y such that Ly = x or Uy = x or the transpose cases.
04281   //
04282 
04283   if (!Filled()) {
04284     EPETRA_CHK_ERR(-1); // Matrix must be filled.
04285   }
04286 
04287   if ((Upper) && (!UpperTriangular()))
04288     EPETRA_CHK_ERR(-2);
04289   if ((!Upper) && (!LowerTriangular()))
04290     EPETRA_CHK_ERR(-3);
04291   if ((!UnitDiagonal) && (NoDiagonal()))
04292     EPETRA_CHK_ERR(-4); // If matrix has no diagonal, we must use UnitDiagonal
04293   if ((!UnitDiagonal) && (NumMyDiagonals()<NumMyRows_))
04294     EPETRA_CHK_ERR(-5); // Need each row to have a diagonal
04295 
04296 
04297   int i, j, j0;
04298   int * NumEntriesPerRow = Graph().NumIndicesPerRow();
04299   int ** Indices = Graph().Indices();
04300   double ** Vals = Values();
04301   int NumMyCols_ = NumMyCols();
04302 
04303   // If upper, point to last row
04304   if ((Upper && !Trans) || (!Upper && Trans)) {
04305     NumEntriesPerRow += NumMyRows_-1;
04306     Indices += NumMyRows_-1;
04307     Vals += NumMyRows_-1;
04308   }
04309 
04310   double *xp = (double*)x.Values();
04311   double *yp = (double*)y.Values();
04312 
04313   if (!Trans) {
04314 
04315     if (Upper) {
04316 
04317       j0 = 1;
04318       if (NoDiagonal())
04319   j0--; // Include first term if no diagonal
04320       for (i=NumMyRows_-1; i >=0; i--) {
04321   int      NumEntries = *NumEntriesPerRow--;
04322   int *    RowIndices = *Indices--;
04323   double * RowValues  = *Vals--;
04324   double sum = 0.0;
04325   for (j=j0; j < NumEntries; j++)
04326     sum += RowValues[j] * yp[RowIndices[j]];
04327 
04328   if (UnitDiagonal)
04329     yp[i] = xp[i] - sum;
04330   else
04331     yp[i] = (xp[i] - sum)/RowValues[0];
04332 
04333       }
04334     }
04335     else {
04336       j0 = 1;
04337       if (NoDiagonal())
04338   j0--; // Include first term if no diagonal
04339       for (i=0; i < NumMyRows_; i++) {
04340   int      NumEntries = *NumEntriesPerRow++ - j0;
04341   int *    RowIndices = *Indices++;
04342   double * RowValues  = *Vals++;
04343   double sum = 0.0;
04344   for (j=0; j < NumEntries; j++)
04345     sum += RowValues[j] * yp[RowIndices[j]];
04346 
04347   if (UnitDiagonal)
04348     yp[i] = xp[i] - sum;
04349   else
04350     yp[i] = (xp[i] - sum)/RowValues[NumEntries];
04351 
04352       }
04353     }
04354   }
04355 
04356   // ***********  Transpose case *******************************
04357 
04358   else {
04359 
04360     if (xp!=yp)
04361       for (i=0; i < NumMyCols_; i++)
04362   yp[i] = xp[i]; // Initialize y for transpose solve
04363 
04364     if (Upper) {
04365 
04366       j0 = 1;
04367       if (NoDiagonal())
04368   j0--; // Include first term if no diagonal
04369 
04370       for (i=0; i < NumMyRows_; i++) {
04371   int      NumEntries = *NumEntriesPerRow++;
04372   int *    RowIndices = *Indices++;
04373   double * RowValues  = *Vals++;
04374   if (!UnitDiagonal)
04375     yp[i] = yp[i]/RowValues[0];
04376      double ytmp = yp[i];
04377   for (j=j0; j < NumEntries; j++)
04378     yp[RowIndices[j]] -= RowValues[j] * ytmp;
04379       }
04380     }
04381     else {
04382 
04383       j0 = 1;
04384       if (NoDiagonal())
04385   j0--; // Include first term if no diagonal
04386 
04387       for (i=NumMyRows_-1; i >= 0; i--) {
04388   int      NumEntries = *NumEntriesPerRow-- - j0;
04389   int *    RowIndices = *Indices--;
04390   double * RowValues  = *Vals--;
04391   if (!UnitDiagonal)
04392     yp[i] = yp[i]/RowValues[NumEntries];
04393      double ytmp = yp[i];
04394   for (j=0; j < NumEntries; j++)
04395     yp[RowIndices[j]] -= RowValues[j] * ytmp;
04396       }
04397     }
04398 
04399   }
04400   UpdateFlops(2*NumGlobalNonzeros64());
04401   return(0);
04402 }
04403 
04404 //=============================================================================
04405 int Epetra_CrsMatrix::Solve1(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
04406 
04407 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS
04408   TEUCHOS_FUNC_TIME_MONITOR("Epetra_CrsMatrix::Solve(Upper,Trans,UnitDiag,X,Y)");
04409 #endif
04410 
04411   //
04412   // This function find Y such that LY = X or UY = X or the transpose cases.
04413   //
04414   if((X.NumVectors() == 1) && (Y.NumVectors() == 1)) {
04415     double* xp = (double*) X[0];
04416     double* yp = (double*) Y[0];
04417     Epetra_Vector x(View, X.Map(), xp);
04418     Epetra_Vector y(View, Y.Map(), yp);
04419     EPETRA_CHK_ERR(Solve1(Upper, Trans, UnitDiagonal, x, y));
04420     return(0);
04421   }
04422   if(!Filled())
04423     EPETRA_CHK_ERR(-1); // Matrix must be filled.
04424 
04425   if((Upper) && (!UpperTriangular()))
04426     EPETRA_CHK_ERR(-2);
04427   if((!Upper) && (!LowerTriangular()))
04428     EPETRA_CHK_ERR(-3);
04429   if((!UnitDiagonal) && (NoDiagonal()))
04430     EPETRA_CHK_ERR(-4); // If matrix has no diagonal, we must use UnitDiagonal
04431   if((!UnitDiagonal) && (NumMyDiagonals()<NumMyRows_))
04432     EPETRA_CHK_ERR(-5); // Need each row to have a diagonal
04433 
04434   int i, j, j0, k;
04435   int* NumEntriesPerRow = Graph().NumIndicesPerRow();
04436   int** Indices = Graph().Indices();
04437   double** Vals = Values();
04438   double diag = 0.0;
04439 
04440   // If upper, point to last row
04441   if((Upper && !Trans) || (!Upper && Trans)) {
04442     NumEntriesPerRow += NumMyRows_-1;
04443     Indices += NumMyRows_-1;
04444     Vals += NumMyRows_-1;
04445   }
04446 
04447   double** Xp = (double**) X.Pointers();
04448   double** Yp = (double**) Y.Pointers();
04449 
04450   int NumVectors = X.NumVectors();
04451 
04452   if(!Trans) {
04453     if(Upper) {
04454       j0 = 1;
04455       if(NoDiagonal())
04456   j0--; // Include first term if no diagonal
04457       for(i = NumMyRows_ - 1; i >= 0; i--) {
04458   int     NumEntries = *NumEntriesPerRow--;
04459   int*    RowIndices = *Indices--;
04460   double* RowValues  = *Vals--;
04461   if(!UnitDiagonal)
04462     diag = 1.0/RowValues[0]; // Take inverse of diagonal once for later use
04463   for(k = 0; k < NumVectors; k++) {
04464     double sum = 0.0;
04465     for(j = j0; j < NumEntries; j++)
04466       sum += RowValues[j] * Yp[k][RowIndices[j]];
04467 
04468     if(UnitDiagonal)
04469       Yp[k][i] = Xp[k][i] - sum;
04470     else
04471       Yp[k][i] = (Xp[k][i] - sum) * diag;
04472   }
04473       }
04474     }
04475     else {
04476       j0 = 1;
04477       if(NoDiagonal())
04478   j0--; // Include first term if no diagonal
04479       for(i = 0; i < NumMyRows_; i++) {
04480   int     NumEntries = *NumEntriesPerRow++ - j0;
04481   int*    RowIndices = *Indices++;
04482   double* RowValues  = *Vals++;
04483   if(!UnitDiagonal)
04484     diag = 1.0/RowValues[NumEntries]; // Take inverse of diagonal once for later use
04485   for(k = 0; k < NumVectors; k++) {
04486     double sum = 0.0;
04487     for(j = 0; j < NumEntries; j++)
04488       sum += RowValues[j] * Yp[k][RowIndices[j]];
04489 
04490     if(UnitDiagonal)
04491       Yp[k][i] = Xp[k][i] - sum;
04492     else
04493       Yp[k][i] = (Xp[k][i] - sum)*diag;
04494   }
04495       }
04496     }
04497   }
04498   // ***********  Transpose case *******************************
04499 
04500   else {
04501     for(k = 0; k < NumVectors; k++)
04502       if(Yp[k] != Xp[k])
04503   for(i = 0; i < NumMyRows_; i++)
04504     Yp[k][i] = Xp[k][i]; // Initialize y for transpose multiply
04505 
04506     if(Upper) {
04507       j0 = 1;
04508       if(NoDiagonal())
04509   j0--; // Include first term if no diagonal
04510 
04511       for(i = 0; i < NumMyRows_; i++) {
04512   int     NumEntries = *NumEntriesPerRow++;
04513   int*    RowIndices = *Indices++;
04514   double* RowValues  = *Vals++;
04515   if(!UnitDiagonal)
04516     diag = 1.0/RowValues[0]; // Take inverse of diagonal once for later use
04517   for(k = 0; k < NumVectors; k++) {
04518     if(!UnitDiagonal)
04519       Yp[k][i] = Yp[k][i]*diag;
04520        double ytmp = Yp[k][i];
04521     for(j = j0; j < NumEntries; j++)
04522       Yp[k][RowIndices[j]] -= RowValues[j] * ytmp;
04523   }
04524       }
04525     }
04526     else {
04527       j0 = 1;
04528       if(NoDiagonal())
04529   j0--; // Include first term if no diagonal
04530       for(i = NumMyRows_ - 1; i >= 0; i--) {
04531   int     NumEntries = *NumEntriesPerRow-- - j0;
04532   int*    RowIndices = *Indices--;
04533   double* RowValues  = *Vals--;
04534      if (!UnitDiagonal)
04535        diag = 1.0/RowValues[NumEntries]; // Take inverse of diagonal once for later use
04536   for(k = 0; k < NumVectors; k++) {
04537     if(!UnitDiagonal)
04538       Yp[k][i] = Yp[k][i]*diag;
04539        double ytmp = Yp[k][i];
04540     for(j = 0; j < NumEntries; j++)
04541       Yp[k][RowIndices[j]] -= RowValues[j] * ytmp;
04542         }
04543       }
04544     }
04545   }
04546 
04547   UpdateFlops(2 * NumVectors * NumGlobalNonzeros64());
04548   return(0);
04549 }
04550 
04551 
04552 
04553 //=============================================================================
04554 Epetra_IntSerialDenseVector& Epetra_CrsMatrix::ExpertExtractIndexOffset(){
04555    return Graph_.CrsGraphData_->IndexOffset_;
04556  }
04557 
04558 //=============================================================================
04559 Epetra_IntSerialDenseVector& Epetra_CrsMatrix::ExpertExtractIndices() {
04560   return Graph_.CrsGraphData_->data->All_Indices_;
04561  }
04562 
04563 //=============================================================================
04564 int Epetra_CrsMatrix::ExpertMakeUniqueCrsGraphData(){
04565    if(Graph_.CrsGraphData_->ReferenceCount() > 1){
04566      Graph_.CrsGraphData_->DecrementReferenceCount();
04567      Graph_.CrsGraphData_=new Epetra_CrsGraphData(Copy, RowMap(),ColMap(),true); // true is for StaticProfile
04568    }
04569    return (0);
04570  }
04571 
04572 //=============================================================================
04573 int Epetra_CrsMatrix::ExpertStaticFillComplete(const Epetra_Map & DomainMap,const Epetra_Map & RangeMap, const Epetra_Import * Importer, const Epetra_Export * Exporter, int NumMyDiagonals){
04574 
04575   Epetra_CrsGraphData& D=*Graph_.CrsGraphData_;
04576   int m=D.RowMap_.NumMyElements();
04577 
04578   bool UseLL=false;
04579   if(D.RowMap_.GlobalIndicesLongLong()) UseLL=true;
04580 
04581   // This only works for constant element size matrices w/ size=1
04582   if(!D.RowMap_.ConstantElementSize() || D.RowMap_.ElementSize()!=1 ||
04583      !D.ColMap_.ConstantElementSize() || D.ColMap_.ElementSize()!=1)
04584     EPETRA_CHK_ERR(-1);
04585 
04586   // Maps
04587   D.DomainMap_ = DomainMap;
04588   D.RangeMap_  = RangeMap;
04589 
04590   // Create import, if needed
04591   if (!D.ColMap_.SameAs(D.DomainMap_)) {
04592     if (D.Importer_ != 0) {
04593       delete D.Importer_;
04594       D.Importer_ = 0;
04595     }
04596     if(Importer && Importer->SourceMap().SameAs(D.DomainMap_) && Importer->TargetMap().SameAs(D.ColMap_)){
04597       D.Importer_=Importer;
04598     }
04599     else {
04600       delete Importer;
04601       D.Importer_ = new Epetra_Import(D.ColMap_, D.DomainMap_);
04602     }
04603   }
04604 
04605   // Create export, if needed
04606   if (!D.RowMap_.SameAs(D.RangeMap_)) {
04607     if (D.Exporter_ != 0) {
04608       delete D.Exporter_;
04609       D.Exporter_ = 0;
04610     }
04611     if(Exporter && Exporter->SourceMap().SameAs(D.RowMap_) && Exporter->TargetMap().SameAs(D.RangeMap_)){
04612       D.Exporter_=Exporter;
04613     }
04614     else {
04615       delete Exporter;
04616       D.Exporter_ = new Epetra_Export(D.RowMap_,D.RangeMap_);       
04617     }
04618   }
04619 
04620 
04621   // Matrix constants
04622   Allocated_                  = true;
04623   StaticGraph_                = true;
04624   UseTranspose_               = false;
04625   constructedWithFilledGraph_ = true;
04626   matrixFillCompleteCalled_   = true;
04627   StorageOptimized_           = true;
04628   squareFillCompleteCalled_   = false; // Always false for the full FillComplete
04629   // Note: Entries in D.Indices_ array point to memory we don't need to dealloc, but the D.Indices_
04630   // array itself needs deallocation.
04631 
04632   // Cleanup existing data
04633   for(int i=0;i<m;i++){
04634     if(Values_)            delete [] Values_[i];
04635   }
04636   D.data->SortedEntries_.clear();
04637 
04638   delete [] Values_;                 Values_=0;
04639   delete [] Values_alloc_lengths_;   Values_alloc_lengths_=0;
04640   delete [] D.data->Indices_;        D.data->Indices_=0;
04641   D.NumAllocatedIndicesPerRow_.Resize(0);
04642 
04643 
04644   // GraphData constants
04645   D.Filled_                                = true;
04646   D.Allocated_                             = true;
04647   D.Sorted_                                = false;
04648   D.StorageOptimized_                      = true;
04649   D.NoRedundancies_                        = true;
04650   D.IndicesAreGlobal_                      = false;
04651   D.IndicesAreLocal_                       = true;
04652   D.IndicesAreContiguous_                  = true;
04653   D.GlobalConstantsComputed_               = true;
04654   D.StaticProfile_                         = true;
04655   D.SortGhostsAssociatedWithEachProcessor_ = true;
04656   D.HaveColMap_                            = true;
04657 
04658   // Don't change the index base
04659 
04660   // Local quantities (no computing)
04661   int nnz=D.IndexOffset_[m]-D.IndexOffset_[0];
04662   D.NumMyRows_       = D.NumMyBlockRows_ = m;
04663   D.NumMyCols_       = D.NumMyBlockCols_ = D.ColMap_.NumMyElements();
04664   D.NumMyNonzeros_   = D.NumMyEntries_   = nnz;
04665   D.MaxRowDim_       = 1;
04666   D.MaxColDim_       = 1;
04667 
04668   // A priori globals
04669   D.GlobalMaxRowDim_ = 1;
04670   D.GlobalMaxColDim_ = 1;
04671 
04672   // Compute max NZ per row, indices, diagonals, triangular
04673   D.MaxNumIndices_=0;
04674   for(int i=0; i<m; i++){
04675     int NumIndices=D.IndexOffset_[i+1]-D.IndexOffset_[i];
04676     D.MaxNumIndices_=EPETRA_MAX(D.MaxNumIndices_,NumIndices);
04677 
04678     // Code copied from Epetra_CrsGraph.Determine_Triangular()
04679     if(NumIndices > 0) {
04680       const int* col_indices = & D.data->All_Indices_[D.IndexOffset_[i]];
04681 
04682       int jl_0 = col_indices[0];
04683       int jl_n = col_indices[NumIndices-1];
04684 
04685       if(jl_n > i) D.LowerTriangular_ = false;
04686       if(jl_0 < i) D.UpperTriangular_ = false;
04687 
04688 
04689       if(NumMyDiagonals == -1) {
04690   // Binary search in GID space
04691   // NOTE: This turns out to be noticibly faster than doing a single LID lookup and then binary searching
04692   // on the LIDs directly.
04693   int insertPoint = -1;
04694   if(UseLL)  {
04695 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
04696     long long jg = D.RowMap_.GID64(i);
04697     if (Epetra_Util_binary_search_aux(jg, col_indices, D.ColMap_.MyGlobalElements64(), NumIndices, insertPoint)>-1) {
04698       D.NumMyBlockDiagonals_++;
04699       D.NumMyDiagonals_ ++;
04700     }
04701 #endif
04702   }
04703   else {
04704 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
04705     int jg = D.RowMap_.GID(i);
04706     if (Epetra_Util_binary_search_aux(jg, col_indices, D.ColMap_.MyGlobalElements(), NumIndices, insertPoint)>-1) {
04707       D.NumMyBlockDiagonals_++;
04708       D.NumMyDiagonals_ ++;
04709     }
04710 #endif
04711   }
04712       }
04713     }
04714   } // end DetermineTriangular
04715 
04716   if(NumMyDiagonals > -1) D.NumMyDiagonals_ = D.NumMyBlockDiagonals_ = NumMyDiagonals;
04717 
04718   D.MaxNumNonzeros_=D.MaxNumIndices_;
04719 
04720   // Compute global constants - Code copied from Epetra_CrsGraph.ComputeGlobalConstants()
04721   Epetra_IntSerialDenseVector tempvec(8); // Temp space
04722   tempvec[0] = D.NumMyEntries_;
04723   tempvec[1] = D.NumMyBlockDiagonals_;
04724   tempvec[2] = D.NumMyDiagonals_;
04725   tempvec[3] = D.NumMyNonzeros_;
04726 
04727   Comm().SumAll(&tempvec[0], &tempvec[4], 4);
04728 
04729   D.NumGlobalEntries_        = tempvec[4];
04730   D.NumGlobalBlockDiagonals_ = tempvec[5];
04731   D.NumGlobalDiagonals_      = tempvec[6];
04732   D.NumGlobalNonzeros_       = tempvec[7];
04733 
04734   tempvec[0] = D.MaxNumIndices_;
04735   tempvec[1] = D.MaxNumNonzeros_;
04736 
04737   Comm().MaxAll(&tempvec[0], &tempvec[2], 2);
04738 
04739   D.GlobalMaxNumIndices_  = tempvec[2];
04740   D.GlobalMaxNumNonzeros_ = tempvec[3];
04741   //end  Epetra_CrsGraph.ComputeGlobalConstants()
04742 
04743 
04744   // Global Info
04745   if(UseLL) {
04746     D.NumGlobalRows_ = D.RangeMap_.NumGlobalPoints64();
04747     D.NumGlobalCols_ = D.DomainMap_.NumGlobalPoints64();
04748   }
04749   else {
04750     D.NumGlobalRows_ = (int) D.RangeMap_.NumGlobalPoints64();
04751     D.NumGlobalCols_ = (int) D.DomainMap_.NumGlobalPoints64();
04752   }
04753   return 0;
04754 }
04755 
04756 // ===================================================================
04757 template<class TransferType>
04758   void Epetra_CrsMatrix::FusedTransfer(const Epetra_CrsMatrix & SourceMatrix, const TransferType & RowTransfer,const Epetra_Map * DomainMap, const Epetra_Map * RangeMap,bool RestrictCommunicator)
04759 {
04760   // Fused constructor, import & FillComplete
04761   int rv;
04762   int N = NumMyRows();
04763 
04764   bool communication_needed = RowTransfer.SourceMap().DistributedGlobal();
04765 
04766   bool UseLL=false;
04767   if(SourceMatrix.RowMap().GlobalIndicesInt()) {
04768     UseLL=false;
04769   }
04770   else if(SourceMatrix.RowMap().GlobalIndicesLongLong()) {
04771     UseLL=true;
04772   }
04773   else
04774     throw ReportError("Epetra_CrsMatrix: Fused import/export constructor unable to determine source global index type",-1);
04775 
04776 
04777   // The basic algorithm here is:
04778   // 1) Call Distor.Do to handle the import.
04779   // 2) Copy all the Imported and Copy/Permuted data into the raw Epetra_CrsMatrix / Epetra_CrsGraphData pointers, still using GIDs.
04780   // 3) Call an optimized version of MakeColMap that avoids the Directory lookups (since the importer knows who owns all the gids) AND
04781   //    reindexes to LIDs.
04782   // 4) Call ExpertStaticFillComplete()
04783 
04784   // Sanity Check
04785   if (!SourceMatrix.RowMap().SameAs(RowTransfer.SourceMap()))
04786     throw ReportError("Epetra_CrsMatrix: Fused import/export constructor requires RowTransfer.SourceMap() to match SourceMatrix.RowMap()",-2);
04787 
04788   // Get information from the Importer
04789   int NumSameIDs             = RowTransfer.NumSameIDs();
04790   int NumPermuteIDs          = RowTransfer.NumPermuteIDs();
04791   int NumRemoteIDs           = RowTransfer.NumRemoteIDs();
04792   int NumExportIDs           = RowTransfer.NumExportIDs();
04793   int* ExportLIDs            = RowTransfer.ExportLIDs();
04794   int* RemoteLIDs            = RowTransfer.RemoteLIDs();
04795   int* PermuteToLIDs         = RowTransfer.PermuteToLIDs();
04796   int* PermuteFromLIDs       = RowTransfer.PermuteFromLIDs();
04797   Epetra_Distributor& Distor = RowTransfer.Distributor();
04798 
04799   // Pull and pointers & allocate memory (rowptr should be the right size already)
04800   Epetra_IntSerialDenseVector & CSR_rowptr = ExpertExtractIndexOffset();
04801   Epetra_IntSerialDenseVector & CSR_colind = ExpertExtractIndices();
04802   double *&                     CSR_vals   = ExpertExtractValues();
04803   CSR_rowptr.Resize(N+1);
04804 
04805   // Unused if we're not in LL mode
04806   std::vector<long long> CSR_colind_LL;
04807 
04808   // Owning PIDs
04809   std::vector<int> SourcePids;
04810   std::vector<int> TargetPids;
04811   int MyPID = Comm().MyPID();
04812 
04813   // The new Domain & Range maps
04814   const Epetra_Map* MyRowMap        = &RowMap();
04815   const Epetra_Map* MyDomainMap     = DomainMap ? DomainMap : &SourceMatrix.DomainMap();
04816   const Epetra_Map* MyRangeMap      = RangeMap  ? RangeMap  : &RowMap();
04817   const Epetra_Map* BaseRowMap      = &RowMap();
04818   const Epetra_Map* BaseDomainMap   = MyDomainMap;
04819 
04820   // Temp variables for sub-communicators
04821   Epetra_Map *ReducedRowMap=0, *ReducedDomainMap=0, *ReducedRangeMap=0;
04822   const Epetra_Comm * ReducedComm = 0;
04823 
04824   /***************************************************/
04825   /***** 1) First communicator restriction phase ****/
04826   /***************************************************/
04827   if(RestrictCommunicator) {
04828     ReducedRowMap = BaseRowMap->RemoveEmptyProcesses();
04829     ReducedComm   = ReducedRowMap ? &ReducedRowMap->Comm() : 0;
04830 
04831     // Now we have to strip down the communicators for the Domain & Range Maps.  Check first if we're lucky
04832     ReducedDomainMap = RowMap().DataPtr() == MyDomainMap->DataPtr() ? ReducedRowMap : MyDomainMap->ReplaceCommWithSubset(ReducedComm);
04833     ReducedRangeMap  = RowMap().DataPtr() == MyRangeMap->DataPtr()  ? ReducedRowMap : MyRangeMap->ReplaceCommWithSubset(ReducedComm);
04834 
04835     // Reset the "my" maps
04836     MyRowMap    = ReducedRowMap;
04837     MyDomainMap = ReducedDomainMap;
04838     MyRangeMap  = ReducedRangeMap;
04839 
04840     // Update my PID, if we've restricted the communicator
04841     if(ReducedComm) MyPID = ReducedComm->MyPID();
04842     else MyPID=-2; // For Debugging
04843   }
04844 
04845   /***************************************************/
04846   /***** 2) From Epetra_DistObject::DoTransfer() *****/
04847   /***************************************************/
04848   rv=CheckSizes(SourceMatrix);
04849   if(rv) throw ReportError("Epetra_CrsMatrix: Fused import/export constructor failed in CheckSizes()",-3);
04850 
04851   int SizeOfPacket;
04852   bool VarSizes = false;
04853   if( NumExportIDs > 0) {
04854     delete [] Sizes_;
04855     Sizes_ = new int[NumExportIDs];
04856   }
04857 
04858   // Get the owning PIDs
04859   const Epetra_Import *MyImporter= SourceMatrix.Importer();
04860 
04861   if(!RestrictCommunicator && MyImporter && BaseDomainMap->SameAs(SourceMatrix.DomainMap())) {
04862     // Same domain map as source matrix
04863     // NOTE: This won't work for restrictComm (because the Importer doesn't know the restricted PIDs), though
04864     // writing and optimized version for that case would be easy (Import an IntVector of the new PIDs).  Might
04865     // want to add this later.
04866     Epetra_Util::GetPids(*MyImporter,SourcePids,false);
04867   }
04868   else if(!MyImporter && BaseDomainMap->SameAs(SourceMatrix.DomainMap())) {
04869     // Matrix has no off-processor entries
04870     SourcePids.resize(SourceMatrix.ColMap().NumMyElements());
04871     SourcePids.assign(SourceMatrix.ColMap().NumMyElements(),MyPID);
04872   }
04873   else if(BaseDomainMap->SameAs(*BaseRowMap) && SourceMatrix.DomainMap().SameAs(SourceMatrix.RowMap())){
04874     // We can use the RowTransfer + SourceMatrix' importer to find out who owns what.
04875     Epetra_IntVector TargetRow_pids(*DomainMap,true);
04876     Epetra_IntVector SourceRow_pids(SourceMatrix.RowMap());
04877     SourcePids.resize(SourceMatrix.ColMap().NumMyElements(),0);
04878     Epetra_IntVector SourceCol_pids(View,SourceMatrix.ColMap(),&SourcePids[0]);
04879 
04880     TargetRow_pids.PutValue(MyPID);
04881     if(typeid(TransferType)==typeid(Epetra_Import))
04882       SourceRow_pids.Export(TargetRow_pids,RowTransfer,Insert);
04883     else if(typeid(TransferType)==typeid(Epetra_Export))
04884       SourceRow_pids.Import(TargetRow_pids,RowTransfer,Insert);
04885     else
04886       throw ReportError("Epetra_CrsMatrix: Fused import/export constructor TransferType must be Epetra_Import or Epetra_Export",-31);
04887     SourceCol_pids.Import(SourceRow_pids,*MyImporter,Insert);
04888   }
04889   else
04890     throw ReportError("Epetra_CrsMatrix: Fused import/export constructor only supports *DomainMap==SourceMatrix.DomainMap() || *DomainMap==RowTransfer.TargetMap() && SourceMatrix.DomainMap() == SourceMatrix.RowMap().",-4);
04891 
04892   // Pack & Prepare w/ owning PIDs
04893   rv=Epetra_Import_Util::PackAndPrepareWithOwningPIDs(SourceMatrix,
04894                   NumExportIDs,ExportLIDs,
04895                   LenExports_,Exports_,SizeOfPacket,
04896                   Sizes_,VarSizes,SourcePids);
04897   if(rv) throw ReportError("Epetra_CrsMatrix: Fused import/export constructor failed in PackAndPrepareWithOwningPIDs()",-5);
04898 
04899   if (communication_needed) {
04900     // Do the exchange of remote data
04901     if( VarSizes )
04902       rv=Distor.Do(Exports_, SizeOfPacket, Sizes_, LenImports_, Imports_);
04903     else
04904       rv=Distor.Do(Exports_, SizeOfPacket, LenImports_, Imports_);
04905   }
04906   if(rv) throw ReportError("Epetra_CrsMatrix: Fused import/export constructor failed in Distor.Do",-6);
04907 
04908 
04909   /*********************************************************************/
04910   /**** 3) Copy all of the Same/Permute/Remote data into CSR_arrays ****/
04911   /*********************************************************************/
04912   // Count nnz
04913   int mynnz = Epetra_Import_Util::UnpackWithOwningPIDsCount(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,PermuteToLIDs,PermuteFromLIDs,LenImports_,Imports_);
04914   // Presume the RowPtr is the right size
04915   // Allocate CSR_colind & CSR_values arrays
04916   CSR_colind.Resize(mynnz);
04917   if(UseLL) CSR_colind_LL.resize(mynnz);
04918   delete [] CSR_vals; // should be a noop.
04919   CSR_vals = new double[mynnz];
04920 
04921   // Unpack into arrays
04922   if(UseLL)
04923     Epetra_Import_Util::UnpackAndCombineIntoCrsArrays(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,PermuteToLIDs,PermuteFromLIDs,LenImports_,Imports_,NumMyRows(),mynnz,MyPID,CSR_rowptr.Values(),CSR_colind_LL.size()?&CSR_colind_LL[0]:0,CSR_vals,SourcePids,TargetPids);
04924   else
04925     Epetra_Import_Util::UnpackAndCombineIntoCrsArrays(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,PermuteToLIDs,PermuteFromLIDs,LenImports_,Imports_,NumMyRows(),mynnz,MyPID,CSR_rowptr.Values(),CSR_colind.Values(),CSR_vals,SourcePids,TargetPids);
04926 
04927   /***************************************************/
04928   /**** 4) Second communicator restriction phase  ****/
04929   /***************************************************/
04930   if(RestrictCommunicator) {
04931     // Dangerous: If we're not in the new communicator, call it quits here.  The user is then responsible
04932     // for not breaking anything on the return.  Thankfully, the dummy RowMap should report no local unknowns,
04933     // so the user can at least test for this particular case.
04934     RemoveEmptyProcessesInPlace(ReducedRowMap);
04935 
04936     if(ReducedComm) {
04937       // Replace the RowMap
04938       Graph_.CrsGraphData_->RowMap_ = *MyRowMap;
04939       Comm_ = &Graph_.CrsGraphData_->RowMap_.Comm();
04940     }
04941     else {
04942       // Replace all the maps with dummy maps with SerialComm, then quit
04943 #if defined(EPETRA_NO_32BIT_GLOBAL_INDICES) && defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
04944       Epetra_Map DummyMap;
04945 #else
04946       Epetra_SerialComm SComm;
04947       Epetra_Map DummyMap(0,0,SComm);
04948 #endif
04949       Graph_.CrsGraphData_->RowMap_    = DummyMap;
04950       Graph_.CrsGraphData_->ColMap_    = DummyMap;
04951       Graph_.CrsGraphData_->RangeMap_  = DummyMap;
04952       Graph_.CrsGraphData_->DomainMap_ = DummyMap;
04953       Comm_ = &Graph_.CrsGraphData_->RowMap_.Comm();
04954       return;
04955     }
04956   }
04957 
04958   /**************************************************************/
04959   /**** 5) Call Optimized MakeColMap w/ no Directory Lookups ****/
04960   /**************************************************************/
04961   //Call an optimized version of MakeColMap that avoids the Directory lookups (since the importer knows who owns all the gids).
04962   std::vector<int> RemotePIDs;
04963   int * pids_ptr = TargetPids.size() ? &TargetPids[0] : 0;
04964   if(UseLL) {
04965 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
04966     long long * CSR_colind_LL_ptr = CSR_colind_LL.size() ? &CSR_colind_LL[0] : 0;
04967     Epetra_Import_Util::LowCommunicationMakeColMapAndReindex(N,CSR_rowptr.Values(),CSR_colind.Values(),CSR_colind_LL_ptr,
04968                    *MyDomainMap,pids_ptr,
04969                    Graph_.CrsGraphData_->SortGhostsAssociatedWithEachProcessor_,RemotePIDs,
04970                    Graph_.CrsGraphData_->ColMap_);
04971     Graph_.CrsGraphData_->HaveColMap_ = true;
04972 #endif
04973   }
04974   else {
04975 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
04976     Epetra_Import_Util::LowCommunicationMakeColMapAndReindex(N,CSR_rowptr.Values(),CSR_colind.Values(),*MyDomainMap,pids_ptr,
04977                    Graph_.CrsGraphData_->SortGhostsAssociatedWithEachProcessor_,RemotePIDs,
04978                    Graph_.CrsGraphData_->ColMap_);
04979     Graph_.CrsGraphData_->HaveColMap_ = true;
04980 #endif
04981   }
04982 
04983   /***************************************************/
04984   /**** 6) Sort (and merge if needed)            ****/
04985   /***************************************************/
04986   // Sort the entries
04987   const Epetra_Import* xferAsImport = dynamic_cast<const Epetra_Import*> (&RowTransfer);
04988   if(xferAsImport)
04989     Epetra_Util::SortCrsEntries(N, CSR_rowptr.Values(), CSR_colind.Values(), CSR_vals);
04990   else 
04991     Epetra_Util::SortAndMergeCrsEntries(N, CSR_rowptr.Values(), CSR_colind.Values(), CSR_vals);
04992 
04993   /***************************************************/
04994   /**** 7) Build Importer & Call ESFC             ****/
04995   /***************************************************/
04996   // Pre-build the importer using the existing PIDs
04997   Epetra_Import * MyImport=0;
04998   int NumRemotePIDs = RemotePIDs.size();
04999   int *RemotePIDs_ptr = RemotePIDs.size() ? &RemotePIDs[0] : 0;
05000   //  if(!RestrictCommunicator && !MyDomainMap->SameAs(ColMap()))
05001   if(!MyDomainMap->SameAs(ColMap()))
05002     MyImport = new Epetra_Import(ColMap(),*MyDomainMap,NumRemotePIDs,RemotePIDs_ptr);
05003 
05004   // Note: At the moment, the RemotePIDs_ptr won't work with the restricted communicator.
05005   // This should be fixed.
05006   ExpertStaticFillComplete(*MyDomainMap,*MyRangeMap,MyImport);
05007 
05008   // Note: ExpertStaticFillComplete assumes ownership of the importer, if we made one...
05009   // We are not to deallocate that here.
05010 
05011   // Cleanup
05012   if(ReducedDomainMap!=ReducedRowMap) delete ReducedDomainMap;
05013   if(ReducedRangeMap !=ReducedRowMap) delete ReducedRangeMap;
05014   delete ReducedRowMap;
05015 }// end FusedTransfer
05016 
05017 
05018 
05019 // ===================================================================
05020 Epetra_CrsMatrix::Epetra_CrsMatrix(const Epetra_CrsMatrix & SourceMatrix, const Epetra_Import & RowImporter,const Epetra_Map * DomainMap, const Epetra_Map * RangeMap, bool RestrictCommunicator)
05021   : Epetra_DistObject(RowImporter.TargetMap(), "Epetra::CrsMatrix"),
05022   Epetra_CompObject(),
05023   Epetra_BLAS(),
05024   Graph_(Copy, RowImporter.TargetMap(), 0, false),
05025   Values_(0),
05026   Values_alloc_lengths_(0),
05027   All_Values_(0),
05028   NumMyRows_(RowImporter.TargetMap().NumMyPoints()),
05029   ImportVector_(0),
05030   ExportVector_(0),
05031   CV_(Copy)
05032 {
05033   FusedTransfer<Epetra_Import>(SourceMatrix,RowImporter,DomainMap,RangeMap,RestrictCommunicator);
05034 }// end fused import constructor
05035 
05036 
05037 // ===================================================================
05038 Epetra_CrsMatrix::Epetra_CrsMatrix(const Epetra_CrsMatrix & SourceMatrix, const Epetra_Export & RowExporter,const Epetra_Map * DomainMap, const Epetra_Map * RangeMap, bool RestrictCommunicator)
05039    : Epetra_DistObject(RowExporter.TargetMap(), "Epetra::CrsMatrix"),
05040    Epetra_CompObject(),
05041    Epetra_BLAS(),
05042    Graph_(Copy, RowExporter.TargetMap(), 0, false),
05043    Values_(0),
05044    Values_alloc_lengths_(0),
05045    All_Values_(0),
05046    NumMyRows_(RowExporter.TargetMap().NumMyPoints()),
05047    ImportVector_(0),
05048    ExportVector_(0),
05049    CV_(Copy)
05050 {
05051   FusedTransfer<Epetra_Export>(SourceMatrix,RowExporter,DomainMap,RangeMap,RestrictCommunicator);
05052 } // end fused export constructor
05053 
05054 
05055 // ===================================================================
05056 void Epetra_CrsMatrix::FusedImport(const Epetra_CrsMatrix & SourceMatrix,
05057            const Epetra_Import & RowImporter,
05058            const Epetra_Map * DomainMap,
05059            const Epetra_Map * RangeMap,
05060            bool RestrictCommunicator) {
05061   FusedTransfer<Epetra_Import>(SourceMatrix,RowImporter,DomainMap,RangeMap,RestrictCommunicator);
05062 }  // end fused import non-constructor
05063 
05064 // ===================================================================
05065 void Epetra_CrsMatrix::FusedExport(const Epetra_CrsMatrix & SourceMatrix,
05066            const Epetra_Export & RowExporter,
05067            const Epetra_Map * DomainMap,
05068            const Epetra_Map * RangeMap,
05069            bool RestrictCommunicator) {
05070   FusedTransfer<Epetra_Export>(SourceMatrix,RowExporter,DomainMap,RangeMap,RestrictCommunicator);
05071 } // end fused export non-constructor
05072 
05073 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines