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_CrsMatrix.h"
00045 #include "Epetra_Map.h"
00046 #include "Epetra_Import.h"
00047 #include "Epetra_Export.h"
00048 #include "Epetra_Vector.h"
00049 #include "Epetra_MultiVector.h"
00050 #include "Epetra_Comm.h"
00051 #include "Epetra_Distributor.h"
00052 #include "Epetra_OffsetIndex.h"
00053 #include "Epetra_BLAS_wrappers.h"
00054 
00055 
00056 #include <cstdlib>
00057 
00058 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
00059 # include "Teuchos_VerboseObject.hpp"
00060 bool Epetra_CrsMatrixTraceDumpMultiply = false;
00061 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
00062 
00063 #ifdef HAVE_EPETRA_TEUCHOS
00064 // Define this macro to see some timers for some of these functions
00065 # define EPETRA_CRSMATRIX_TEUCHOS_TIMERS
00066 #endif
00067 
00068 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS
00069 #  include "Teuchos_TimeMonitor.hpp"
00070 #endif
00071 
00072 //==============================================================================
00073 Epetra_CrsMatrix::Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_Map& rowMap, const int* NumEntriesPerRow, bool StaticProfile) 
00074   : Epetra_DistObject(rowMap, "Epetra::CrsMatrix"),
00075     Epetra_CompObject(),
00076     Epetra_BLAS(),
00077     Graph_(CV, rowMap, NumEntriesPerRow, StaticProfile),
00078     Allocated_(false),
00079     StaticGraph_(false),
00080     UseTranspose_(false),
00081     constructedWithFilledGraph_(false),
00082     matrixFillCompleteCalled_(false),
00083     StorageOptimized_(false),
00084     Values_(0),
00085     Values_alloc_lengths_(0),
00086     All_Values_(0),
00087     NormInf_(0.0),
00088     NormOne_(0.0),
00089     NormFrob_(0.0),
00090     NumMyRows_(rowMap.NumMyPoints()),
00091     ImportVector_(0),
00092     ExportVector_(0),
00093     CV_(CV),
00094     squareFillCompleteCalled_(false)
00095 {
00096   InitializeDefaults();
00097   Allocate();
00098 }
00099 
00100 //==============================================================================
00101 Epetra_CrsMatrix::Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_Map& rowMap, int NumEntriesPerRow, bool StaticProfile) 
00102   : Epetra_DistObject(rowMap, "Epetra::CrsMatrix"),
00103     Epetra_CompObject(),
00104     Epetra_BLAS(),
00105     Graph_(CV, rowMap, NumEntriesPerRow, StaticProfile),
00106     Allocated_(false),
00107     StaticGraph_(false),
00108     UseTranspose_(false),
00109     constructedWithFilledGraph_(false),
00110     matrixFillCompleteCalled_(false),
00111     StorageOptimized_(false),
00112     Values_(0),
00113     Values_alloc_lengths_(0),
00114     All_Values_(0),
00115     NormInf_(0.0),
00116     NormOne_(0.0),
00117     NormFrob_(0.0),
00118     NumMyRows_(rowMap.NumMyPoints()),
00119     ImportVector_(0),
00120     ExportVector_(0),
00121     CV_(CV),
00122     squareFillCompleteCalled_(false)
00123 {
00124   InitializeDefaults();
00125   Allocate();
00126 }
00127 //==============================================================================
00128 Epetra_CrsMatrix::Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_Map& rowMap, 
00129            const Epetra_Map& colMap, const int* NumEntriesPerRow, bool StaticProfile) 
00130   : Epetra_DistObject(rowMap, "Epetra::CrsMatrix"),
00131     Epetra_CompObject(),
00132     Epetra_BLAS(),
00133     Graph_(CV, rowMap, colMap, NumEntriesPerRow, StaticProfile),
00134     Allocated_(false),
00135     StaticGraph_(false),
00136     UseTranspose_(false),
00137     constructedWithFilledGraph_(false),
00138     matrixFillCompleteCalled_(false),
00139     StorageOptimized_(false),
00140     Values_(0),
00141     Values_alloc_lengths_(0),
00142     All_Values_(0),
00143     NormInf_(0.0),
00144     NormOne_(0.0),
00145     NormFrob_(0.0),
00146     NumMyRows_(rowMap.NumMyPoints()),
00147     ImportVector_(0),
00148     ExportVector_(0),
00149     CV_(CV),
00150     squareFillCompleteCalled_(false)
00151 {
00152   InitializeDefaults();
00153   Allocate();
00154 }
00155 
00156 //==============================================================================
00157 Epetra_CrsMatrix::Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_Map& rowMap, 
00158            const Epetra_Map& colMap, int NumEntriesPerRow, bool StaticProfile) 
00159   : Epetra_DistObject(rowMap, "Epetra::CrsMatrix"),
00160     Epetra_CompObject(),
00161     Epetra_BLAS(),
00162     Graph_(CV, rowMap, colMap,  NumEntriesPerRow, StaticProfile),
00163     Allocated_(false),
00164     StaticGraph_(false),
00165     UseTranspose_(false),
00166     constructedWithFilledGraph_(false),
00167     matrixFillCompleteCalled_(false),
00168     StorageOptimized_(false),
00169     Values_(0),
00170     Values_alloc_lengths_(0),
00171     All_Values_(0),
00172     NormInf_(0.0),
00173     NormOne_(0.0),
00174     NormFrob_(0.0),
00175     NumMyRows_(rowMap.NumMyPoints()),
00176     ImportVector_(0),
00177     ExportVector_(0),
00178     CV_(CV),
00179     squareFillCompleteCalled_(false)
00180 {
00181   InitializeDefaults();
00182   Allocate();
00183 }
00184 //==============================================================================
00185 Epetra_CrsMatrix::Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_CrsGraph& graph) 
00186   : Epetra_DistObject(graph.Map(), "Epetra::CrsMatrix"),
00187     Epetra_CompObject(),
00188     Epetra_BLAS(),
00189     Graph_(graph),
00190     Allocated_(false),
00191     StaticGraph_(true),
00192     UseTranspose_(false),
00193     constructedWithFilledGraph_(false),
00194     matrixFillCompleteCalled_(false),
00195     StorageOptimized_(false),
00196     Values_(0),
00197     Values_alloc_lengths_(0),
00198     All_Values_(0),
00199     NormInf_(0.0),
00200     NormOne_(0.0),
00201     NormFrob_(0.0),
00202     NumMyRows_(graph.NumMyRows()),
00203     ImportVector_(0),
00204     ExportVector_(0),
00205     CV_(CV),
00206     squareFillCompleteCalled_(false)
00207 {
00208   constructedWithFilledGraph_ = graph.Filled();
00209   InitializeDefaults();
00210   Allocate();
00211 }
00212 
00213 //==============================================================================
00214 Epetra_CrsMatrix::Epetra_CrsMatrix(const Epetra_CrsMatrix& Matrix) 
00215   : Epetra_DistObject(Matrix),
00216     Epetra_CompObject(Matrix),
00217     Epetra_BLAS(),
00218     Graph_(Matrix.Graph()),
00219     Allocated_(false),
00220     StaticGraph_(true),
00221     UseTranspose_(Matrix.UseTranspose_),
00222     constructedWithFilledGraph_(false),
00223     matrixFillCompleteCalled_(false),
00224     StorageOptimized_(false),
00225     Values_(0),
00226     Values_alloc_lengths_(0),
00227     All_Values_(0),
00228     NormInf_(0.0),
00229     NormOne_(0.0),
00230     NormFrob_(0.0),
00231     NumMyRows_(Matrix.NumMyRows()),
00232     ImportVector_(0),
00233     ExportVector_(0),
00234     CV_(Copy),
00235     squareFillCompleteCalled_(false)
00236 {
00237   InitializeDefaults();
00238   operator=(Matrix);
00239 }
00240 
00241 //==============================================================================
00242 Epetra_CrsMatrix& Epetra_CrsMatrix::operator=(const Epetra_CrsMatrix& src)
00243 {
00244   if (this == &src) {
00245     return( *this );
00246   }
00247 
00248   if (!src.Filled()) throw ReportError("Copying an Epetra_CrsMatrix requires source matrix to have Filled()==true", -1);
00249 
00250   Graph_ = src.Graph_; // Copy graph
00251 
00252   DeleteMemory();
00253 
00254   StaticGraph_ = true;
00255   UseTranspose_ = src.UseTranspose_;
00256   constructedWithFilledGraph_ = src.constructedWithFilledGraph_;
00257   matrixFillCompleteCalled_ = src.matrixFillCompleteCalled_;
00258   Values_ = 0;
00259   Values_alloc_lengths_ = 0;
00260   All_Values_ = 0;
00261   NormInf_ = -1.0;
00262   NormOne_ = -1.0;
00263   NormFrob_ = -1.0;
00264   NumMyRows_ = src.NumMyRows_;
00265   ImportVector_ = 0;
00266   ExportVector_ = 0;
00267 
00268   CV_ = Copy;
00269 
00270   StorageOptimized_ = src.StorageOptimized_;
00271   if (src.StorageOptimized()) { // Special copy for case where storage is optimized
00272 
00273     int numMyNonzeros = Graph().NumMyEntries();
00274     if (numMyNonzeros>0) All_Values_ = new double[numMyNonzeros];
00275     double * srcValues = src.All_Values();
00276     for (int i=0; i<numMyNonzeros; ++i) All_Values_[i] = srcValues[i];
00277     Allocated_ = true;
00278 #ifdef Epetra_ENABLE_CASK
00279     if( matrixFillCompleteCalled_  ) {
00280       cask = cask_handler_copy(src.cask);
00281     }
00282 #endif
00283   }
00284   else { // copy for non-optimized storage
00285     
00286     Allocate();
00287     for (int i=0; i<NumMyRows_; i++) {
00288       int NumEntries = src.NumMyEntries(i);
00289       double * const srcValues = src.Values(i);
00290       double * targValues = Values(i);
00291       for (int j=0; j< NumEntries; j++) targValues[j] = srcValues[j];
00292     }
00293   }
00294 
00295   return( *this );
00296 }
00297 
00298 //==============================================================================
00299 void Epetra_CrsMatrix::InitializeDefaults() { // Initialize all attributes that have trivial default values
00300 
00301   UseTranspose_ = false;
00302   Values_ = 0;
00303   Values_alloc_lengths_ = 0;
00304   All_Values_ = 0;
00305   NormInf_ = -1.0;
00306   NormOne_ = -1.0;
00307   NormFrob_ = -1.0;
00308   ImportVector_ = 0;
00309   ExportVector_ = 0;
00310 
00311   return;
00312 }
00313 
00314 //==============================================================================
00315 int Epetra_CrsMatrix::Allocate() {
00316 
00317   int i, j;
00318 
00319   // Allocate Values array
00320   Values_ = NumMyRows_ > 0 ? new double*[NumMyRows_] : NULL;
00321   Values_alloc_lengths_ = NumMyRows_ > 0 ? new int[NumMyRows_] : NULL;
00322   if (NumMyRows_ > 0) {
00323     for(j=0; j<NumMyRows_; ++j) Values_alloc_lengths_[j] = 0;
00324   }
00325 
00326   // Allocate and initialize entries if we are copying data
00327   if (CV_==Copy) {
00328     if (Graph().StaticProfile() || Graph().StorageOptimized()) {
00329       int numMyNonzeros = Graph().NumMyEntries();
00330       if (numMyNonzeros>0) All_Values_ = new double[numMyNonzeros];
00331       if(Graph().StorageOptimized()){
00332         StorageOptimized_ = true;
00333       }
00334     }
00335     double * all_values = All_Values_;
00336     for (i=0; i<NumMyRows_; i++) {
00337       int NumAllocatedEntries = Graph().NumAllocatedMyIndices(i);
00338 
00339       if (NumAllocatedEntries > 0) {
00340         if (Graph().StaticProfile() || Graph().StorageOptimized()) {
00341           Values_[i] = all_values;
00342           all_values += NumAllocatedEntries;
00343         }
00344         else {
00345           Values_[i] = new double[NumAllocatedEntries];
00346           Values_alloc_lengths_[i] = NumAllocatedEntries;
00347         }
00348       }
00349       else 
00350         Values_[i] = 0;
00351 
00352       for(j=0; j< NumAllocatedEntries; j++) 
00353         Values_[i][j] = 0.0; // Fill values with zero
00354     }
00355   }  
00356   else {
00357     for (i=0; i<NumMyRows_; i++) {
00358       Values_[i] = 0;
00359     }
00360   }
00361   SetAllocated(true);
00362 #ifdef Epetra_ENABLE_CASK
00363   cask=NULL;
00364 #endif
00365   return(0);
00366 }
00367 //==============================================================================
00368 Epetra_CrsMatrix::~Epetra_CrsMatrix()
00369 {
00370   DeleteMemory();
00371 }
00372 
00373 //==============================================================================
00374 void Epetra_CrsMatrix::DeleteMemory()
00375 {
00376   int i;
00377 
00378   if (CV_==Copy) {
00379     if (All_Values_!=0)
00380       delete [] All_Values_;
00381     else if (Values_!=0)
00382       for (i=0; i<NumMyRows_; i++) 
00383         if (Graph().NumAllocatedMyIndices(i) >0) 
00384           delete [] Values_[i];
00385   }
00386 
00387   if (ImportVector_!=0) 
00388     delete ImportVector_;
00389   ImportVector_=0;
00390 
00391   if (ExportVector_!=0)
00392     delete ExportVector_;
00393   ExportVector_=0;
00394 
00395   delete [] Values_;
00396   Values_ = 0;
00397 
00398   delete [] Values_alloc_lengths_;
00399   Values_alloc_lengths_ = 0;
00400 
00401   NumMyRows_ = 0;
00402 
00403 #ifdef Epetra_ENABLE_CASK
00404   if( StorageOptimized_  )
00405   {
00406     if( cask != NULL ) 
00407       cask_handler_destroy(cask);
00408 
00409     cask = NULL;
00410   }
00411 #endif
00412 
00413   Allocated_ = false;
00414 }
00415 
00416 //==============================================================================
00417 int Epetra_CrsMatrix::ReplaceRowMap(const Epetra_BlockMap& newmap)
00418 {
00419   int err = Graph_.ReplaceRowMap(newmap);
00420   if (err == 0) {
00421     //update export vector.
00422 
00423     if (ExportVector_ != 0) {
00424       delete ExportVector_; 
00425       ExportVector_= 0;
00426     }
00427 
00428     ExportVector_ = new Epetra_MultiVector(RowMap(),1);
00429   }
00430   return(err);
00431 }
00432 
00433 //==============================================================================
00434 int Epetra_CrsMatrix::ReplaceColMap(const Epetra_BlockMap& newmap)
00435 {
00436   int err = Graph_.ReplaceColMap(newmap);
00437   if (err == 0) {
00438     //update import vector.
00439 
00440     if (ImportVector_ != 0) {
00441       delete ImportVector_;
00442       ImportVector_= 0;
00443     }
00444 
00445     ImportVector_ = new Epetra_MultiVector(ColMap(),1);
00446   }
00447   return(err);
00448 }
00449 
00450 //==============================================================================
00451 int Epetra_CrsMatrix::PutScalar(double ScalarConstant) 
00452 {
00453   if (StorageOptimized()) {
00454     int length = NumMyNonzeros();
00455     for (int i=0; i<length; ++i) All_Values_[i] = ScalarConstant;
00456   }
00457   else {
00458     for(int i=0; i<NumMyRows_; i++) {
00459       int NumEntries = Graph().NumMyIndices(i);
00460       double * targValues = Values(i);
00461       for(int j=0; j< NumEntries; j++) 
00462   targValues[j] = ScalarConstant;
00463     }
00464   }
00465   return(0);
00466 }
00467 //==============================================================================
00468 int Epetra_CrsMatrix::Scale(double ScalarConstant) 
00469 {
00470   if (StorageOptimized()) {
00471     int length = NumMyNonzeros();
00472     for (int i=0; i<length; ++i) All_Values_[i] *= ScalarConstant;
00473   }
00474   else {
00475     for(int i=0; i<NumMyRows_; i++) {
00476       int NumEntries = Graph().NumMyIndices(i);
00477       double * targValues = Values(i);
00478       for(int j=0; j< NumEntries; j++) 
00479   targValues[j] *= ScalarConstant;
00480     }
00481   }
00482   return(0);
00483 }
00484 
00485 //==========================================================================
00486 int Epetra_CrsMatrix::InsertGlobalValues(int Row, int NumEntries,
00487            const double* values,
00488            const int* Indices)
00489 {
00490   if(IndicesAreLocal()) 
00491     EPETRA_CHK_ERR(-2); // Cannot insert global values into local graph
00492   if(IndicesAreContiguous()) 
00493     EPETRA_CHK_ERR(-3); // Indices cannot be individually deleted and newed
00494   Graph_.SetIndicesAreGlobal(true);
00495   Row = Graph_.LRID(Row); // Find local row number for this global row index
00496 
00497   EPETRA_CHK_ERR( InsertValues(Row, NumEntries, values, Indices) );
00498 
00499   return(0);
00500 }
00501 
00502 //==========================================================================
00503 int Epetra_CrsMatrix::InsertGlobalValues(int Row, int NumEntries,
00504            double* values,
00505            int* Indices)
00506 {
00507   if(IndicesAreLocal()) 
00508     EPETRA_CHK_ERR(-2); // Cannot insert global values into local graph
00509   if(IndicesAreContiguous()) 
00510     EPETRA_CHK_ERR(-3); // Indices cannot be individually deleted and newed
00511   Graph_.SetIndicesAreGlobal(true);
00512   Row = Graph_.LRID(Row); // Find local row number for this global row index
00513 
00514   EPETRA_CHK_ERR( InsertValues(Row, NumEntries, values, Indices) );
00515 
00516   return(0);
00517 }
00518 
00519 //==========================================================================
00520 int Epetra_CrsMatrix::InsertMyValues(int Row, int NumEntries,
00521              const double* values,
00522              const int* Indices)
00523 {
00524   if(IndicesAreGlobal()) 
00525     EPETRA_CHK_ERR(-2); // Cannot insert global values into filled graph
00526   if(IndicesAreContiguous() && CV_==Copy) 
00527     EPETRA_CHK_ERR(-3); // Indices cannot be individually deleted and new
00528   Graph_.SetIndicesAreLocal(true);
00529 
00530   EPETRA_CHK_ERR( InsertValues(Row, NumEntries, values, Indices) );
00531 
00532   return(0);
00533 
00534 }
00535 
00536 //==========================================================================
00537 int Epetra_CrsMatrix::InsertMyValues(int Row, int NumEntries,
00538              double* values,
00539              int* Indices)
00540 {
00541   if(IndicesAreGlobal()) 
00542     EPETRA_CHK_ERR(-2); // Cannot insert global values into filled graph
00543   if(IndicesAreContiguous() && CV_==Copy) 
00544     EPETRA_CHK_ERR(-3); // Indices cannot be individually deleted and new
00545   Graph_.SetIndicesAreLocal(true);
00546 
00547   EPETRA_CHK_ERR( InsertValues(Row, NumEntries, values, Indices) );
00548 
00549   return(0);
00550 
00551 }
00552 
00553 //==========================================================================
00554 int Epetra_CrsMatrix::InsertValues(int Row, int NumEntries,
00555            const double* values,
00556            const int* Indices)
00557 {
00558   if(CV_ == View){
00559     //cannot allow View mode with const pointers
00560     EPETRA_CHK_ERR(-4);
00561   }
00562   else{
00563     //to avoid code duplication I am using a cheap tactic of removing the constness of
00564     //values and Indices. Since this is only called in copy mode the passed in variables
00565     //will not be modified.
00566     return(InsertValues(Row, NumEntries, const_cast<double*>(values), const_cast<int*>(Indices)));
00567   }
00568   return 0;
00569 }
00570 
00571 //==========================================================================
00572 int Epetra_CrsMatrix::InsertValues(int Row, int NumEntries,
00573            double* values,
00574            int* Indices)
00575 {
00576   int j;
00577   int ierr = 0;
00578 
00579   if(Row < 0 || Row >= NumMyRows_) 
00580     EPETRA_CHK_ERR(-1); // Not in Row range
00581 
00582   if(CV_ == View) {
00583     //test indices in static graph
00584     if(StaticGraph()) {
00585       int testNumEntries;
00586       int* testIndices;
00587       int testRow = Row;
00588       if(IndicesAreGlobal()) 
00589         testRow = Graph_.LRID( Row );
00590       EPETRA_CHK_ERR(Graph_.ExtractMyRowView(testRow, testNumEntries, testIndices));
00591 
00592       bool match = true;
00593       if(NumEntries != testNumEntries) 
00594         match = false;
00595       for(int i = 0; i < NumEntries; ++i)
00596         match = match && (Indices[i]==testIndices[i]);
00597 
00598       if(!match)
00599         ierr = -3;
00600     }
00601 
00602     if(Values_[Row] != 0) 
00603       ierr = 2; // This row has been defined already.  Issue warning.
00604     Values_[Row] = values;
00605   }
00606   else {    
00607     if(StaticGraph()) 
00608       EPETRA_CHK_ERR(-2); // If the matrix graph is fully constructed, we cannot insert new values
00609 
00610     int tmpNumEntries = NumEntries;
00611 
00612     if(Graph_.HaveColMap()) { //must insert only valid indices, values
00613       const double* tmpValues = values;
00614       values = new double[NumEntries];
00615       int loc = 0;
00616       if(IndicesAreLocal()) {
00617         for(int i = 0; i < NumEntries; ++i)
00618           if(Graph_.ColMap().MyLID(Indices[i])) 
00619             values[loc++] = tmpValues[i];
00620       }
00621       else {
00622         for(int i = 0; i < NumEntries; ++i)
00623           if(Graph_.ColMap().MyGID(Indices[i])) 
00624             values[loc++] = tmpValues[i];
00625       }
00626       if(NumEntries != loc) 
00627         ierr = 2; //Some columns excluded
00628       NumEntries = loc;
00629     } 
00630 
00631     int start = Graph().NumMyIndices(Row);
00632     int stop = start + NumEntries;
00633     int NumAllocatedEntries = Values_alloc_lengths_[Row];
00634     if(stop > NumAllocatedEntries) {
00635       if (Graph().StaticProfile() && stop > Graph().NumAllocatedMyIndices(Row)) {
00636         EPETRA_CHK_ERR(-2); // Cannot expand graph storage if graph created using StaticProfile
00637       }
00638       if(NumAllocatedEntries == 0) {
00639         Values_[Row] = new double[NumEntries]; // Row was never allocated, so do it
00640         Values_alloc_lengths_[Row] = NumEntries;
00641       }
00642       else {
00643         ierr = 1; // Out of room.  Must delete and allocate more space...
00644         double* tmp_Values = new double[stop];
00645         for(j = 0; j < start; j++) 
00646           tmp_Values[j] = Values_[Row][j]; // Copy existing entries
00647         delete[] Values_[Row]; // Delete old storage
00648         Values_[Row] = tmp_Values; // Set pointer to new storage
00649         Values_alloc_lengths_[Row] = stop;
00650       }
00651     }
00652 
00653     for(j = start; j < stop; j++) 
00654       Values_[Row][j] = values[j-start];
00655 
00656 
00657     NumEntries = tmpNumEntries;
00658     if(Graph_.HaveColMap()) 
00659       delete[] values;
00660   }
00661 
00662   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
00663   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
00664   NormFrob_ = -1.0;
00665 
00666   if(!StaticGraph()) {
00667     EPETRA_CHK_ERR(Graph_.InsertIndices(Row, NumEntries, Indices));
00668   }
00669 
00670   EPETRA_CHK_ERR(ierr);
00671   return(0);
00672 }
00673 
00674 //==========================================================================
00675 int Epetra_CrsMatrix::InsertOffsetValues(int Row, int NumEntries,
00676            double* values,
00677            int* Indices)
00678 {
00679   return ReplaceOffsetValues(Row, NumEntries, values, Indices);
00680 }
00681 
00682 //==========================================================================
00683 int Epetra_CrsMatrix::ReplaceGlobalValues(int Row, int NumEntries, const double * srcValues, const int *Indices) {
00684 
00685   int j;
00686   int ierr = 0;
00687   int Loc;
00688 
00689   Row = Graph_.LRID(Row); // Normalize row range
00690     
00691   if (Row < 0 || Row >= NumMyRows_) {
00692     EPETRA_CHK_ERR(-1); // Not in Row range
00693   }
00694   double * targValues = Values(Row);
00695   for (j=0; j<NumEntries; j++) {
00696     int Index = Indices[j];
00697     if (Graph_.FindGlobalIndexLoc(Row,Index,j,Loc)) 
00698       targValues[Loc] = srcValues[j];
00699     else 
00700       ierr = 2; // Value Excluded
00701   }
00702 
00703   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
00704   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
00705   NormFrob_ = -1.0;
00706 
00707   EPETRA_CHK_ERR(ierr);
00708   return(0);
00709 }
00710 
00711 //==========================================================================
00712 int Epetra_CrsMatrix::ReplaceMyValues(int Row, int NumEntries, const double * srcValues, const int *Indices) {
00713 
00714   if (!IndicesAreLocal()) 
00715     EPETRA_CHK_ERR(-4); // Indices must be local.
00716 
00717   int j;
00718   int ierr = 0;
00719   int Loc;
00720 
00721   if (Row < 0 || Row >= NumMyRows_) {
00722     EPETRA_CHK_ERR(-1); // Not in Row range
00723   }
00724 
00725   double* RowValues = Values(Row); 
00726   for (j=0; j<NumEntries; j++) {
00727     int Index = Indices[j];
00728     if (Graph_.FindMyIndexLoc(Row,Index,j,Loc)) 
00729       RowValues[Loc] = srcValues[j];
00730     else 
00731       ierr = 2; // Value Excluded
00732   }
00733 
00734   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
00735   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
00736   NormFrob_ = -1.0;
00737 
00738   EPETRA_CHK_ERR(ierr);
00739   return(0);
00740 }
00741 
00742 //==========================================================================
00743 int Epetra_CrsMatrix::ReplaceOffsetValues(int Row, int NumEntries,
00744             const double * srcValues, const int *Offsets)
00745 {
00746   int j;
00747   int ierr = 0;
00748 
00749   Row = Graph_.LRID(Row); // Normalize row range
00750     
00751   if (Row < 0 || Row >= NumMyRows_) {
00752     EPETRA_CHK_ERR(-1); // Not in Row range
00753   }
00754 
00755   double* RowValues = Values(Row); 
00756   for(j=0; j<NumEntries; j++) {
00757     if( Offsets[j] != -1 )
00758       RowValues[Offsets[j]] = srcValues[j];
00759   }
00760 
00761   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
00762   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
00763   NormFrob_ = -1.0;
00764 
00765   EPETRA_CHK_ERR(ierr);
00766   return(0);
00767 }
00768 
00769 //==========================================================================
00770 int Epetra_CrsMatrix::SumIntoGlobalValues(int Row,
00771             int NumEntries,
00772             const double * srcValues,
00773             const int *Indices)
00774 {
00775   int j;
00776   int ierr = 0;
00777   int Loc = 0;
00778 
00779   Row = Graph_.LRID(Row); // Normalize row range
00780     
00781   if (Row < 0 || Row >= NumMyRows_) {
00782     EPETRA_CHK_ERR(-1); // Not in Row range
00783   }
00784 
00785   if (StaticGraph() && !Graph_.HaveColMap()) {
00786     EPETRA_CHK_ERR(-1);
00787   }
00788 
00789   double * RowValues = Values(Row);
00790 
00791   if (!StaticGraph()) {
00792     for (j=0; j<NumEntries; j++) {
00793       int Index = Indices[j];
00794       if (Graph_.FindGlobalIndexLoc(Row,Index,j,Loc))
00795         RowValues[Loc] += srcValues[j];
00796       else
00797         ierr = 2; // Value Excluded
00798     }
00799   }
00800   else {
00801     const Epetra_BlockMap& colmap = Graph_.ColMap();
00802     int NumColIndices = Graph_.NumMyIndices(Row);
00803     const int* ColIndices = Graph_.Indices(Row);
00804 
00805     if (Graph_.Sorted()) {
00806       int insertPoint;
00807       for (j=0; j<NumEntries; j++) {
00808         int Index = colmap.LID(Indices[j]);
00809 
00810         // Check whether the next added element is the subsequent element in
00811         // the graph indices, then we can skip the binary search
00812         if (Loc < NumColIndices && Index == ColIndices[Loc])
00813           RowValues[Loc] += srcValues[j];
00814         else {
00815           Loc = Epetra_Util_binary_search(Index, ColIndices, NumColIndices, insertPoint);
00816           if (Loc > -1)
00817             RowValues[Loc] += srcValues[j];
00818           else 
00819             ierr = 2; // Value Excluded
00820         }
00821         ++Loc;
00822       }
00823     }
00824     else
00825       for (j=0; j<NumEntries; j++) {
00826         int Index = colmap.LID(Indices[j]);
00827         if (Graph_.FindMyIndexLoc(NumColIndices,ColIndices,Index,j,Loc)) 
00828           RowValues[Loc] += srcValues[j];
00829         else 
00830           ierr = 2; // Value Excluded
00831       }
00832   }
00833 
00834   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
00835   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
00836   NormFrob_ = -1.0;
00837 
00838   EPETRA_CHK_ERR(ierr);
00839 
00840   return(0);
00841 }
00842 
00843 //==========================================================================
00844 int Epetra_CrsMatrix::SumIntoMyValues(int Row, int NumEntries, const double * srcValues, const int *Indices) {
00845 
00846   if (!IndicesAreLocal()) 
00847     EPETRA_CHK_ERR(-4); // Indices must be local.
00848 
00849   int j;
00850   int ierr = 0;
00851   int Loc = 0;
00852   int insertPoint;
00853 
00854   if (Row < 0 || Row >= NumMyRows_) {
00855     EPETRA_CHK_ERR(-1); // Not in Row range
00856   }
00857 
00858   double* RowValues = Values(Row);
00859   int NumColIndices = Graph_.NumMyIndices(Row);
00860   const int* ColIndices = Graph_.Indices(Row);
00861   if (Graph_.Sorted()) {
00862     for (j=0; j<NumEntries; j++) {
00863       int Index = Indices[j];
00864 
00865       // Check whether the next added element is the subsequent element in
00866       // the graph indices.
00867       if (Loc < NumColIndices && Index == ColIndices[Loc])
00868         RowValues[Loc] += srcValues[j];
00869       else {
00870         Loc = Epetra_Util_binary_search(Index, ColIndices, NumColIndices, insertPoint);
00871         if (Loc > -1)
00872           RowValues[Loc] += srcValues[j];
00873         else 
00874           ierr = 2; // Value Excluded
00875       }
00876       ++Loc;
00877     }
00878   }
00879   else {
00880     for (j=0; j<NumEntries; j++) {
00881       int Index = Indices[j];
00882       if (Graph_.FindMyIndexLoc(Row,Index,j,Loc)) 
00883         RowValues[Loc] += srcValues[j];
00884       else 
00885         ierr = 2; // Value Excluded
00886     }
00887   }
00888 
00889   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
00890   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
00891   NormFrob_ = -1.0;
00892 
00893   EPETRA_CHK_ERR(ierr);
00894 
00895   return(0);
00896 }
00897 
00898 //==========================================================================
00899 int Epetra_CrsMatrix::SumIntoOffsetValues(int Row, int NumEntries, const double * srcValues, const int *Offsets) {
00900 
00901   int j;
00902   int ierr = 0;
00903 
00904   Row = Graph_.LRID(Row); // Normalize row range
00905     
00906   if (Row < 0 || Row >= NumMyRows_) {
00907     EPETRA_CHK_ERR(-1); // Not in Row range
00908   }
00909 
00910   double* RowValues = Values(Row);
00911   for (j=0; j<NumEntries; j++) {
00912     if( Offsets[j] != -1 )
00913       RowValues[Offsets[j]] += srcValues[j];
00914   }
00915 
00916   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
00917   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
00918   NormFrob_ = -1.0;
00919 
00920   EPETRA_CHK_ERR(ierr);
00921 
00922   return(0);
00923 }
00924 
00925 //==========================================================================
00926 int Epetra_CrsMatrix::FillComplete(bool OptimizeDataStorage) {
00927   squareFillCompleteCalled_ = true;
00928   EPETRA_CHK_ERR(FillComplete(RowMap(), RowMap(), OptimizeDataStorage));
00929   return(0);
00930 }
00931 
00932 //==========================================================================
00933 int Epetra_CrsMatrix::FillComplete(const Epetra_Map& domain_map,
00934            const Epetra_Map& range_map, bool OptimizeDataStorage)
00935 {
00936   int returnValue = 0;
00937 
00938   if (Graph_.Filled()) {
00939     if (!constructedWithFilledGraph_ && !matrixFillCompleteCalled_) {
00940       returnValue = 2;
00941     }
00942   }
00943 
00944   if (!StaticGraph()) {
00945     if (Graph_.MakeIndicesLocal(domain_map, range_map) < 0) {
00946       return(-1);
00947     }
00948   }
00949   SortEntries();  // Sort column entries from smallest to largest
00950   MergeRedundantEntries(); // Get rid of any redundant index values
00951   if (!StaticGraph()) {
00952     if (Graph_.FillComplete(domain_map, range_map) < 0) {
00953       return(-2);
00954     }
00955   }
00956 
00957   matrixFillCompleteCalled_ = true;
00958 
00959   if (squareFillCompleteCalled_) {
00960     if (DomainMap().NumGlobalElements() != RangeMap().NumGlobalElements()) {
00961       returnValue = 3;
00962     }
00963     squareFillCompleteCalled_ = false;
00964     EPETRA_CHK_ERR(returnValue);
00965   }
00966 
00967   if (OptimizeDataStorage) { EPETRA_CHK_ERR(OptimizeStorage()); }
00968 
00969   return(returnValue);
00970 }
00971 
00972 //==========================================================================
00973 int Epetra_CrsMatrix::TransformToLocal() {
00974   EPETRA_CHK_ERR(FillComplete());
00975   return(0);
00976 }
00977 
00978 //==========================================================================
00979 int Epetra_CrsMatrix::TransformToLocal(const Epetra_Map* domainMap, const Epetra_Map* rangeMap) {
00980   EPETRA_CHK_ERR(FillComplete(*domainMap, *rangeMap));
00981   return(0);
00982 }
00983 
00984 //==========================================================================
00985 int Epetra_CrsMatrix::SortEntries() {
00986 
00987   if(!IndicesAreLocal()) 
00988     EPETRA_CHK_ERR(-1);
00989   if(Sorted())
00990     return(0);
00991 
00992   // For each row, sort column entries from smallest to largest.
00993   // Use shell sort. Stable sort so it is fast if indices are already sorted.
00994 
00995   
00996   for(int i = 0; i < NumMyRows_; i++){
00997 
00998     double* locValues = Values(i);
00999     int NumEntries = Graph().NumMyIndices(i);
01000     int* locIndices = Graph().Indices(i);
01001     
01002     int n = NumEntries;
01003     int m = n/2;
01004     
01005     while(m > 0) {
01006       int max = n - m;
01007       for(int j = 0; j < max; j++) {
01008   for(int k = j; k >= 0; k-=m) {
01009     if(locIndices[k+m] >= locIndices[k])
01010       break;
01011     double dtemp = locValues[k+m];
01012     locValues[k+m] = locValues[k];
01013     locValues[k] = dtemp;
01014     int itemp = locIndices[k+m];
01015     locIndices[k+m] = locIndices[k];
01016     locIndices[k] = itemp;
01017   }
01018       }
01019       m = m/2;
01020     }
01021   }
01022   Graph_.SetSorted(true); // This also sorted the graph
01023   return(0);
01024 }
01025 
01026 //==========================================================================
01027 int Epetra_CrsMatrix::MergeRedundantEntries() {
01028 
01029   int i;
01030 
01031   if(NoRedundancies()) 
01032     return(0);
01033   if(!Sorted()) 
01034     EPETRA_CHK_ERR(-1);  // Must have sorted entries
01035 
01036   // For each row, remove column indices that are repeated.
01037   // Also, determine if matrix is upper or lower triangular or has no diagonal (Done in graph)
01038   // Note:  This function assumes that SortEntries was already called.
01039 
01040   for(i = 0; i<NumMyRows_; i++) {
01041     int NumEntries = Graph().NumMyIndices(i);
01042     if(NumEntries > 1) {
01043       double* const locValues = Values(i);
01044       int* const locIndices = Graph().Indices(i);   
01045       int curEntry =0;
01046       double curValue = locValues[0];
01047       for(int k = 1; k < NumEntries; k++) {
01048   if(locIndices[k] == locIndices[k-1]) 
01049     curValue += locValues[k];
01050   else {
01051     locValues[curEntry++] = curValue;
01052     curValue = locValues[k];
01053   }
01054       }
01055       locValues[curEntry] = curValue;
01056       
01057     }
01058   }
01059   
01060   EPETRA_CHK_ERR(Graph_.RemoveRedundantIndices()); // Remove redundant indices and then return
01061   return(0);
01062 }
01063 
01064 //==========================================================================
01065 int Epetra_CrsMatrix::OptimizeStorage() {
01066 
01067 
01068   if (StorageOptimized()) 
01069     return(0); // Have we been here before?
01070   if (!Filled()) EPETRA_CHK_ERR(-1); // Cannot optimize storage before calling FillComplete()
01071 
01072 
01073   int ierr = Graph_.OptimizeStorage();
01074   if (ierr!=0) EPETRA_CHK_ERR(ierr);  // In order for OptimizeStorage to make sense for the matrix, it must work on the graph.
01075 
01076   bool Contiguous = true; // Assume contiguous is true
01077   for (int i=1; i<NumMyRows_; i++){
01078     int NumEntries = Graph().NumMyIndices(i-1);
01079     
01080     // check if end of beginning of current row starts immediately after end of previous row.
01081     if (Values_[i]!=Values_[i-1]+NumEntries) {
01082       Contiguous = false;
01083       break;
01084     }
01085   }
01086 
01087   // NOTE:  At the end of the above loop set, there is a possibility that NumEntries and NumAllocatedEntries
01088   //        for the last row could be different, but I don't think it matters.
01089 
01090 
01091   if ((CV_==View) && !Contiguous) 
01092     EPETRA_CHK_ERR(-1);  // This is user data, it's not contiguous and we can't make it so.
01093 
01094   
01095   if(!Contiguous) { // Must pack indices if not already contiguous.  Pack values into All_values_
01096     
01097     if (All_Values_==0) {
01098       // Compute Number of Nonzero entries (Done in FillComplete, but we may not have been there yet.)
01099       int numMyNonzeros = Graph_.NumMyNonzeros();
01100     
01101       // Allocate one big array for all values
01102       All_Values_ = new double[numMyNonzeros];
01103       if(All_Values_ == 0) throw ReportError("Error with All_Values_ allocation.", -99);
01104 
01105       const int *const  IndexOffset = Graph().IndexOffset();
01106       const int numMyRows = NumMyRows_;
01107       double ** Values_s = Values_;
01108       double * All_Values_s = All_Values_;
01109 #ifdef EPETRA_HAVE_OMP
01110 #pragma omp parallel for default(none) shared(Values_s,All_Values_s)
01111 #endif
01112       for (int i=0; i<numMyRows; i++) {
01113         int NumEntries = Graph().NumMyIndices(i);
01114         int curOffset = IndexOffset[i];
01115         double * values = Values_s[i];
01116         double * newValues = All_Values_s+curOffset;
01117   for (int j=0; j<NumEntries; j++) newValues[j] = values[j];
01118       }
01119     }
01120     else { // Static Profile, so just pack into existing storage (can't be threaded)
01121       double * tmp = All_Values_;
01122       for (int i=0; i<NumMyRows_; i++) {
01123         int NumEntries = Graph().NumMyIndices(i);
01124         double * values = Values_[i];
01125         if (tmp!=values) // Copy values if not pointing to same location
01126           for (int j=0; j<NumEntries; j++) tmp[j] = values[j];
01127         tmp += NumEntries;
01128       }
01129     }
01130 
01131     // Free Values_ arrays
01132     for (int i=0;i<NumMyRows_; ++i) {
01133       if (Values_alloc_lengths_[i] != 0) delete [] Values_[i];
01134     }
01135     
01136     delete [] Values_alloc_lengths_; Values_alloc_lengths_ = 0;
01137   } // End of !Contiguous section
01138   else {
01139     //if already contiguous, we'll simply set All_Values_ to be
01140     //a copy of Values_[0].
01141     All_Values_ = NumMyRows_ > 0 ? Values_[0] : NULL;
01142   }
01143   
01144   // Delete unneeded storage
01145   delete [] Values_; Values_=0;
01146 
01147 #ifdef Epetra_ENABLE_CASK 
01148   if (cask == NULL  && Graph().StorageOptimized() )  {
01149      int * Indices = Graph().All_Indices();
01150      int * IndexOffset = Graph().IndexOffset();
01151      int NumMyCols_ = NumMyCols();
01152      cask_handler_initialize(&cask);
01153      cask_csr_analysis(NumMyRows_, NumMyCols_, IndexOffset, Indices,
01154                        NumGlobalNonzeros(),cask);
01155   }
01156 #endif
01157 
01158 
01159   StorageOptimized_ = true;
01160 
01161   
01162   return(0);
01163 }
01164 //==========================================================================
01165 int Epetra_CrsMatrix::ExtractGlobalRowCopy(int Row, int Length, int & NumEntries, double * values,
01166              int * Indices) const 
01167 {
01168 
01169   int ierr = Graph_.ExtractGlobalRowCopy(Row, Length, NumEntries, Indices);
01170   if (ierr) 
01171     EPETRA_CHK_ERR(ierr);
01172 
01173   EPETRA_CHK_ERR(ExtractGlobalRowCopy(Row, Length, NumEntries, values));
01174   return(0);
01175 }
01176 //==========================================================================
01177 int Epetra_CrsMatrix::ExtractMyRowCopy(int Row, int Length, int & NumEntries, double * values,
01178                int * Indices) const 
01179 {
01180 
01181   int ierr = Graph_.ExtractMyRowCopy(Row, Length, NumEntries, Indices);
01182   if (ierr) 
01183     EPETRA_CHK_ERR(ierr);
01184 
01185   EPETRA_CHK_ERR(ExtractMyRowCopy(Row, Length, NumEntries, values));
01186   return(0);
01187 }
01188 //==========================================================================
01189 int Epetra_CrsMatrix::NumMyRowEntries(int Row, int & NumEntries) const 
01190 {
01191 
01192   if (!MyLRID(Row)) 
01193     EPETRA_CHK_ERR(-1); // Not in the range of local rows
01194   NumEntries = NumMyEntries(Row);
01195   return(0);
01196 }
01197 //==========================================================================
01198 int Epetra_CrsMatrix::ExtractGlobalRowCopy(int Row, int Length, int & NumEntries, double * values) const 
01199 {
01200 
01201   int Row0 = Graph_.RowMap().LID(Row); // Normalize row range
01202 
01203   EPETRA_CHK_ERR(ExtractMyRowCopy(Row0, Length, NumEntries, values));
01204   return(0);
01205 }
01206 
01207 
01208 //==========================================================================
01209 int Epetra_CrsMatrix::ExtractMyRowCopy(int Row, int Length, int & NumEntries, double * targValues) const 
01210 {
01211   int j;
01212 
01213   if (Row < 0 || Row >= NumMyRows_) 
01214     EPETRA_CHK_ERR(-1); // Not in Row range
01215 
01216   NumEntries = Graph().NumMyIndices(Row);
01217   if (Length < NumEntries) 
01218     EPETRA_CHK_ERR(-2); // Not enough space for copy. Needed size is passed back in NumEntries
01219 
01220   double * srcValues = Values(Row);
01221 
01222   for(j=0; j<NumEntries; j++)
01223     targValues[j] = srcValues[j];
01224   
01225   return(0);
01226 }
01227 
01228 
01229 //==============================================================================
01230 int Epetra_CrsMatrix::ExtractDiagonalCopy(Epetra_Vector & Diagonal) const {
01231   
01232   if(!Filled()) 
01233     EPETRA_CHK_ERR(-1); // Can't get diagonal unless matrix is filled (and in local index space)
01234   if(!RowMap().SameAs(Diagonal.Map())) 
01235     EPETRA_CHK_ERR(-2); // Maps must be the same
01236 
01237   for(int i = 0; i < NumMyRows_; i++) {
01238     int ii = GRID(i);
01239     int NumEntries = Graph().NumMyIndices(i);
01240     int* Indices = Graph().Indices(i);
01241     double * srcValues = Values(i);
01242     
01243     Diagonal[i] = 0.0;
01244     for(int j = 0; j < NumEntries; j++) {
01245       if(ii == GCID(Indices[j])) {
01246   Diagonal[i] = srcValues[j];
01247   break;
01248       }
01249     }
01250   }
01251   return(0);
01252 }
01253 //==============================================================================
01254 int Epetra_CrsMatrix::ReplaceDiagonalValues(const Epetra_Vector & Diagonal) {
01255   
01256   if(!Filled())
01257     EPETRA_CHK_ERR(-1); // Can't replace diagonal unless matrix is filled (and in local index space)
01258   if(!RowMap().SameAs(Diagonal.Map())) 
01259     EPETRA_CHK_ERR(-2); // Maps must be the same
01260 
01261   int ierr = 0;
01262   for(int i = 0; i < NumMyRows_; i++) {
01263     int ii = GRID(i);
01264     int NumEntries = Graph().NumMyIndices(i);
01265     int* Indices = Graph().Indices(i);
01266     double * targValues = Values(i);
01267     bool DiagMissing = true;
01268     for(int j = 0; j < NumEntries; j++) {
01269       if(ii == GCID(Indices[j])) {
01270   targValues[j] = Diagonal[i];
01271   DiagMissing = false;
01272   break;
01273       }
01274     }
01275     if(DiagMissing) 
01276       ierr = 1; // flag a warning error
01277   }
01278 
01279   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
01280   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
01281   NormFrob_ = -1.0;
01282 
01283   EPETRA_CHK_ERR(ierr);
01284 
01285   return(0);
01286 }
01287 //==========================================================================
01288 int Epetra_CrsMatrix::ExtractGlobalRowView(int Row, int & NumEntries, double *& values, int *& Indices) const 
01289 {
01290 
01291   int ierr = Graph_.ExtractGlobalRowView(Row, NumEntries, Indices);
01292   if (ierr) 
01293     EPETRA_CHK_ERR(ierr);
01294 
01295   EPETRA_CHK_ERR(ExtractGlobalRowView(Row, NumEntries, values));
01296   return(0);
01297 }
01298 //==========================================================================
01299 int Epetra_CrsMatrix::ExtractMyRowView(int Row, int & NumEntries, double *& values, int *& Indices) const 
01300 {
01301   int ierr = Graph_.ExtractMyRowView(Row, NumEntries, Indices);
01302   if (ierr) 
01303     EPETRA_CHK_ERR(ierr);
01304 
01305   EPETRA_CHK_ERR(ExtractMyRowView(Row, NumEntries, values));
01306   return(0);
01307 }
01308 //==========================================================================
01309 int Epetra_CrsMatrix::ExtractGlobalRowView(int Row, int & NumEntries, double *& values) const 
01310 {
01311 
01312   int Row0 = Graph_.RowMap().LID(Row); // Normalize row range
01313 
01314   EPETRA_CHK_ERR(ExtractMyRowView(Row0, NumEntries, values));
01315   return(0);
01316 }
01317 
01318 //==========================================================================
01319 int Epetra_CrsMatrix::ExtractMyRowView(int Row, int & NumEntries, double *& targValues) const 
01320 {
01321 
01322   if (Row < 0 || Row >= NumMyRows_) 
01323     EPETRA_CHK_ERR(-1); // Not in Row range
01324 
01325   NumEntries = Graph().NumMyIndices(Row);
01326 
01327   targValues = Values(Row);
01328   
01329   return(0);
01330 }
01331 
01332 //=============================================================================
01333 int Epetra_CrsMatrix::Solve(bool Upper, bool Trans, bool UnitDiagonal,
01334           const Epetra_Vector& x, Epetra_Vector& y) const
01335 {
01336 
01337 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS
01338   TEUCHOS_FUNC_TIME_MONITOR("Epetra_CrsMatrix::Solve(Upper,Trans,UnitDiag,x,y)");
01339 #endif
01340 
01341   //
01342   // This function finds y such that Ly = x or Uy = x or the transpose cases.
01343   //
01344 
01345   // First short-circuit to the pre-5.0 version if no storage optimization was performed
01346   if (!StorageOptimized() && !Graph().StorageOptimized()) {
01347     EPETRA_CHK_ERR(Solve1(Upper, Trans, UnitDiagonal, x, y));
01348     return(0);
01349   }
01350 
01351   if (!Filled()) {
01352     EPETRA_CHK_ERR(-1); // Matrix must be filled.
01353   }
01354 
01355   if ((Upper) && (!UpperTriangular())) 
01356     EPETRA_CHK_ERR(-2);
01357   if ((!Upper) && (!LowerTriangular())) 
01358     EPETRA_CHK_ERR(-3);
01359   if ((!UnitDiagonal) && (NoDiagonal())) 
01360     EPETRA_CHK_ERR(-4); // If matrix has no diagonal, we must use UnitDiagonal
01361   if ((!UnitDiagonal) && (NumMyDiagonals()<NumMyRows_)) 
01362     EPETRA_CHK_ERR(-5); // Need each row to have a diagonal
01363 
01364   double *xp = (double*)x.Values();
01365   double *yp = (double*)y.Values();
01366 
01367       
01368   GeneralSV(Upper, Trans, UnitDiagonal, xp, yp);
01369 
01370   UpdateFlops(2*NumGlobalNonzeros());
01371   return(0);
01372 }
01373 
01374 //=============================================================================
01375 int Epetra_CrsMatrix::Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
01376 
01377 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS
01378   TEUCHOS_FUNC_TIME_MONITOR("Epetra_CrsMatrix::Solve(Upper,Trans,UnitDiag,X,Y)");
01379 #endif
01380 
01381   //
01382   // This function find Y such that LY = X or UY = X or the transpose cases.
01383   //
01384 
01385   // First short-circuit to the pre-5.0 version if no storage optimization was performed
01386   if (!StorageOptimized() && !Graph().StorageOptimized()) {
01387     EPETRA_CHK_ERR(Solve1(Upper, Trans, UnitDiagonal, X, Y));
01388     return(0);
01389   }
01390   if(!Filled()) 
01391     EPETRA_CHK_ERR(-1); // Matrix must be filled.
01392 
01393   if((Upper) && (!UpperTriangular()))
01394     EPETRA_CHK_ERR(-2);
01395   if((!Upper) && (!LowerTriangular()))
01396     EPETRA_CHK_ERR(-3);
01397   if((!UnitDiagonal) && (NoDiagonal()))
01398     EPETRA_CHK_ERR(-4); // If matrix has no diagonal, we must use UnitDiagonal
01399   if((!UnitDiagonal) && (NumMyDiagonals()<NumMyRows_))
01400     EPETRA_CHK_ERR(-5); // Need each row to have a diagonal
01401 
01402   double** Xp = (double**) X.Pointers();
01403   double** Yp = (double**) Y.Pointers();
01404   int LDX = X.ConstantStride() ? X.Stride() : 0;
01405   int LDY = Y.ConstantStride() ? Y.Stride() : 0;
01406   int NumVectors = X.NumVectors();
01407 
01408 
01409     // Do actual computation
01410   if (NumVectors==1)
01411     GeneralSV(Upper, Trans, UnitDiagonal, *Xp, *Yp);
01412   else
01413     GeneralSM(Upper, Trans, UnitDiagonal, Xp, LDX, Yp, LDY, NumVectors);
01414 
01415   UpdateFlops(2 * NumVectors * NumGlobalNonzeros());
01416   return(0);
01417 }
01418 
01419 //=============================================================================
01420 int Epetra_CrsMatrix::InvRowSums(Epetra_Vector& x) const {
01421   //
01422   // Put inverse of the sum of absolute values of the ith row of A in x[i].
01423   //
01424 
01425   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
01426   int ierr = 0;
01427   int i, j;
01428   x.PutScalar(0.0); // Make sure we sum into a vector of zeros.
01429   double * xp = (double*)x.Values();
01430   if (Graph().RangeMap().SameAs(x.Map()) && Exporter() != 0) {
01431     Epetra_Vector x_tmp(RowMap());
01432     x_tmp.PutScalar(0.0);
01433     double * x_tmp_p = (double*)x_tmp.Values();
01434     for (i=0; i < NumMyRows_; i++) {
01435       int      NumEntries = NumMyEntries(i);
01436       double * RowValues  = Values(i);
01437       for (j=0; j < NumEntries; j++)  x_tmp_p[i] += std::abs(RowValues[j]);
01438     }
01439     EPETRA_CHK_ERR(x.Export(x_tmp, *Exporter(), Add)); //Export partial row sums to x.
01440     int myLength = x.MyLength();
01441     for (i=0; i<myLength; i++) { 
01442       if (xp[i]<Epetra_MinDouble) {
01443         if (xp[i]==0.0) ierr = 1; // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)
01444         else if (ierr!=1) ierr = 2;
01445         xp[i] = Epetra_MaxDouble;
01446       }
01447       else
01448         xp[i] = 1.0/xp[i];
01449     }
01450   }
01451   else if (Graph().RowMap().SameAs(x.Map())) {
01452     for (i=0; i < NumMyRows_; i++) {
01453       int      NumEntries = NumMyEntries(i);
01454       double * RowValues  = Values(i);
01455       double scale = 0.0;
01456       for (j=0; j < NumEntries; j++) scale += std::abs(RowValues[j]);
01457       if (scale<Epetra_MinDouble) {
01458         if (scale==0.0) ierr = 1; // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)
01459         else if (ierr!=1) ierr = 2;
01460         xp[i] = Epetra_MaxDouble;
01461       }
01462       else
01463         xp[i] = 1.0/scale;
01464     }
01465   }
01466   else { // x.Map different than both Graph().RowMap() and Graph().RangeMap()
01467     EPETRA_CHK_ERR(-2); // The map of x must be the RowMap or RangeMap of A.
01468   }
01469   UpdateFlops(NumGlobalNonzeros());
01470   EPETRA_CHK_ERR(ierr);
01471   return(0);
01472 }
01473 
01474 //=============================================================================
01475 int Epetra_CrsMatrix::InvRowMaxs(Epetra_Vector& x) const {
01476   //
01477   // Put inverse of the max of absolute values of the ith row of A in x[i].
01478   //
01479 
01480   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
01481   int ierr = 0;
01482   int i, j;
01483   bool needExport = false;
01484   double * xp = (double*)x.Values();
01485   Epetra_Vector* x_tmp = 0;
01486   if (Graph().RangeMap().SameAs(x.Map())) {
01487     if (Exporter() != 0) {
01488       needExport = true; //Having this information later avoids a .SameAs
01489       x_tmp = new Epetra_Vector(RowMap()); // Create import vector if needed
01490       xp = (double*)x_tmp->Values();
01491     }
01492   }
01493   else if (!Graph().RowMap().SameAs(x.Map())) {
01494     EPETRA_CHK_ERR(-2); // The map of x must be the RowMap or RangeMap of A.
01495   }
01496   for (i=0; i < NumMyRows_; i++) {
01497     int      NumEntries = NumMyEntries(i);
01498     double * RowValues  = Values(i);
01499     double scale = 0.0;
01500     for (j=0; j < NumEntries; j++) scale = EPETRA_MAX(std::abs(RowValues[j]),scale);
01501     if (scale<Epetra_MinDouble) {
01502       if (scale==0.0) ierr = 1; // Set error to 1 to signal that zero rowmax found (supercedes ierr = 2)
01503       else if (ierr!=1) ierr = 2;
01504       xp[i] = Epetra_MaxDouble;
01505     }
01506     else
01507       xp[i] = 1.0/scale;
01508   }
01509   if(needExport) {
01510     x.PutScalar(0.0);
01511     EPETRA_CHK_ERR(x.Export(*x_tmp, *Exporter(), Insert)); // Fill x with values from temp vector
01512     delete x_tmp;
01513   }
01514   UpdateFlops(NumGlobalNonzeros());
01515   EPETRA_CHK_ERR(ierr);
01516   return(0);
01517 }
01518 
01519 //=============================================================================
01520 int Epetra_CrsMatrix::InvColSums(Epetra_Vector& x) const {
01521   //
01522   // Put inverse of the sum of absolute values of the jth column of A in x[j].
01523   //
01524 
01525   if(!Filled())  EPETRA_CHK_ERR(-1); // Matrix must be filled.
01526   int ierr = 0;
01527   int i, j;
01528   int MapNumMyElements = x.Map().NumMyElements();
01529   x.PutScalar(0.0); // Make sure we sum into a vector of zeros.
01530   double* xp = (double*)x.Values();
01531   if(Graph().DomainMap().SameAs(x.Map()) && Importer() != 0) {
01532     Epetra_Vector x_tmp(ColMap());
01533     x_tmp.PutScalar(0.0);
01534     double * x_tmp_p = (double*)x_tmp.Values();
01535     for(i = 0; i < NumMyRows_; i++) {
01536       int     NumEntries = NumMyEntries(i);
01537       int*    ColIndices = Graph().Indices(i);
01538       double* RowValues  = Values(i);
01539       for(j = 0; j < NumEntries; j++) 
01540         x_tmp_p[ColIndices[j]] += std::abs(RowValues[j]);
01541     }
01542     EPETRA_CHK_ERR(x.Export(x_tmp, *Importer(), Add)); // Fill x with partial column sums
01543   }
01544   else if(Graph().ColMap().SameAs(x.Map())) {
01545     for(i = 0; i < NumMyRows_; i++) {
01546       int     NumEntries = NumMyEntries(i);
01547       int*    ColIndices = Graph().Indices(i);
01548       double* RowValues  = Values(i);
01549       for(j = 0; j < NumEntries; j++) 
01550         xp[ColIndices[j]] += std::abs(RowValues[j]);
01551     }
01552   }
01553   else { //x.Map different than both Graph().ColMap() and Graph().DomainMap()
01554     EPETRA_CHK_ERR(-2); // x must have the same distribution as the domain of A
01555   }
01556 
01557   // Invert values, don't allow them to get too large
01558   for(i = 0; i < MapNumMyElements; i++) {
01559     double scale = xp[i];
01560     if(scale < Epetra_MinDouble) {
01561       if(scale == 0.0) 
01562   ierr = 1; // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)
01563       else if(ierr != 1) 
01564   ierr = 2;
01565       xp[i] = Epetra_MaxDouble;
01566     }
01567     else
01568       xp[i] = 1.0 / scale;
01569   }
01570 
01571   UpdateFlops(NumGlobalNonzeros());
01572   EPETRA_CHK_ERR(ierr);
01573   return(0);
01574 }
01575 
01576 //=============================================================================
01577 int Epetra_CrsMatrix::InvColMaxs(Epetra_Vector& x) const {
01578   //
01579   // Put inverse of the max of absolute values of the jth column of A in x[j].
01580   //
01581 
01582   if(!Filled())  EPETRA_CHK_ERR(-1); // Matrix must be filled.
01583   int ierr = 0;
01584   int i, j;
01585   int MapNumMyElements = x.Map().NumMyElements();
01586   x.PutScalar(0.0); // Make sure we sum into a vector of zeros.
01587   double* xp = (double*)x.Values();
01588   if(Graph().DomainMap().SameAs(x.Map()) && Importer() != 0) {
01589     Epetra_Vector x_tmp(ColMap());
01590     x_tmp.PutScalar(0.0);
01591     double * x_tmp_p = (double*)x_tmp.Values();
01592     for(i = 0; i < NumMyRows_; i++) {
01593       int     NumEntries = NumMyEntries(i);
01594       int*    ColIndices = Graph().Indices(i);
01595       double* RowValues  = Values(i);
01596       for(j = 0; j < NumEntries; j++) 
01597         x_tmp_p[ColIndices[j]] = EPETRA_MAX(std::abs(RowValues[j]),x_tmp_p[ColIndices[j]]);
01598     }
01599     EPETRA_CHK_ERR(x.Export(x_tmp, *Importer(), AbsMax)); // Fill x with partial column sums
01600   }
01601   else if(Graph().ColMap().SameAs(x.Map())) {
01602     for(i = 0; i < NumMyRows_; i++) {
01603       int     NumEntries = NumMyEntries(i);
01604       int*    ColIndices = Graph().Indices(i);
01605       double* RowValues  = Values(i);
01606       for(j = 0; j < NumEntries; j++) 
01607         xp[ColIndices[j]] = EPETRA_MAX(std::abs(RowValues[j]),xp[ColIndices[j]]);
01608     }
01609   }
01610   else { //x.Map different than both Graph().ColMap() and Graph().DomainMap()
01611     EPETRA_CHK_ERR(-2); // x must have the same distribution as the domain of A
01612   }
01613 
01614   // Invert values, don't allow them to get too large
01615   for(i = 0; i < MapNumMyElements; i++) {
01616     double scale = xp[i];
01617     if(scale < Epetra_MinDouble) {
01618       if(scale == 0.0) 
01619   ierr = 1; // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)
01620       else if(ierr != 1) 
01621   ierr = 2;
01622       xp[i] = Epetra_MaxDouble;
01623     }
01624     else
01625       xp[i] = 1.0 / scale;
01626   }
01627 
01628   UpdateFlops(NumGlobalNonzeros());
01629   EPETRA_CHK_ERR(ierr);
01630   return(0);
01631 }
01632 
01633 //=============================================================================
01634 int Epetra_CrsMatrix::LeftScale(const Epetra_Vector& x) {
01635   //
01636   // This function scales the ith row of A by x[i].
01637   //
01638 
01639   if(!Filled()) 
01640     EPETRA_CHK_ERR(-1); // Matrix must be filled.
01641   double* xp = 0;
01642   if(Graph().RangeMap().SameAs(x.Map()))  
01643     // If we have a non-trivial exporter, we must import elements that are 
01644     // permuted or are on other processors.  (We will use the exporter to
01645     // perform the import.)
01646     if(Exporter() != 0) {
01647       UpdateExportVector(1);
01648       EPETRA_CHK_ERR(ExportVector_->Import(x,*Exporter(), Insert));
01649       xp = (double*) ExportVector_->Values();
01650     }
01651     else
01652       xp = (double*)x.Values();
01653   else if (Graph().RowMap().SameAs(x.Map()))
01654     xp = (double*)x.Values();
01655   else {
01656     EPETRA_CHK_ERR(-2); // The Map of x must be the RowMap or RangeMap of A.
01657   }
01658   int i, j;
01659 
01660   for(i = 0; i < NumMyRows_; i++) {
01661     int      NumEntries = NumMyEntries(i);
01662     double* RowValues  = Values(i);
01663     double scale = xp[i];
01664     for(j = 0; j < NumEntries; j++)  
01665       RowValues[j] *= scale;
01666   }
01667   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
01668   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
01669   NormFrob_ = -1.0;
01670 
01671   UpdateFlops(NumGlobalNonzeros());
01672 
01673   return(0);
01674 }
01675 
01676 //=============================================================================
01677 int Epetra_CrsMatrix::RightScale(const Epetra_Vector& x) {
01678   //
01679   // This function scales the jth column of A by x[j].
01680   //
01681 
01682   if(!Filled()) 
01683     EPETRA_CHK_ERR(-1); // Matrix must be filled.
01684   double* xp = 0;
01685   if(Graph().DomainMap().SameAs(x.Map())) 
01686     // If we have a non-trivial exporter, we must import elements that are 
01687     // permuted or are on other processors.
01688     if(Importer() != 0) {
01689       UpdateImportVector(1);
01690       EPETRA_CHK_ERR(ImportVector_->Import(x, *Importer(), Insert));
01691       xp = (double*) ImportVector_->Values();
01692     }
01693     else
01694       xp = (double*)x.Values();
01695   else if(Graph().ColMap().SameAs(x.Map()))
01696     xp = (double*)x.Values(); 
01697   else
01698     EPETRA_CHK_ERR(-2); // The Map of x must be the RowMap or RangeMap of A.
01699   int i, j;
01700 
01701   for(i = 0; i < NumMyRows_; i++) {
01702     int     NumEntries = NumMyEntries(i);
01703     int*    ColIndices = Graph().Indices(i);
01704     double* RowValues  = Values(i);
01705     for(j = 0; j < NumEntries; j++)  
01706       RowValues[j] *=  xp[ColIndices[j]];
01707   }
01708   NormOne_ = -1.0; // Reset Norm so it will be recomputed.
01709   NormInf_ = -1.0; // Reset Norm so it will be recomputed.
01710   NormFrob_ = -1.0;
01711 
01712   UpdateFlops(NumGlobalNonzeros());
01713   return(0);
01714 }
01715 
01716 //=============================================================================
01717 double Epetra_CrsMatrix::NormInf() const {
01718 
01719 #if 0
01720   //
01721   //  Commenting this section out disables caching, ie. 
01722   //  causes the norm to be computed each time NormInf is called.
01723   //  See bug #1151 for a full discussion.  
01724   //
01725   double MinNorm ; 
01726   Comm().MinAll( &NormInf_, &MinNorm, 1 ) ; 
01727 
01728   if( MinNorm >= 0.0) 
01729     return(NormInf_);
01730 #endif
01731 
01732   if(!Filled()) 
01733     EPETRA_CHK_ERR(-1); // Matrix must be filled.
01734 
01735   Epetra_Vector x(RangeMap()); // Need temp vector for row sums
01736   double* xp = (double*)x.Values();
01737   Epetra_MultiVector* x_tmp = 0;
01738 
01739   // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
01740   if(Exporter() != 0) {
01741     x_tmp = new Epetra_Vector(RowMap()); // Create temporary import vector if needed
01742     xp = (double*)x_tmp->Values();
01743   }
01744   int i, j;
01745 
01746   for(i = 0; i < NumMyRows_; i++) {
01747     xp[i] = 0.0;
01748     int     NumEntries = NumMyEntries(i);
01749     double* RowValues  = Values(i);
01750     for(j = 0; j < NumEntries; j++) 
01751       xp[i] += std::abs(RowValues[j]);
01752   }
01753   if(Exporter() != 0) {
01754     x.PutScalar(0.0);
01755     EPETRA_CHK_ERR(x.Export(*x_tmp, *Exporter(), Add)); // Fill x with Values from temp vector
01756   }
01757   x.MaxValue(&NormInf_); // Find max
01758   if(x_tmp != 0) 
01759     delete x_tmp;
01760   UpdateFlops(NumGlobalNonzeros());
01761   return(NormInf_);
01762 }
01763 //=============================================================================
01764 double Epetra_CrsMatrix::NormOne() const {
01765 
01766 #if 0
01767   //
01768   //  Commenting this section out disables caching, ie. 
01769   //  causes the norm to be computed each time NormOne is called.  
01770   //  See bug #1151 for a full discussion.  
01771   //
01772   double MinNorm ; 
01773   Comm().MinAll( &NormOne_, &MinNorm, 1 ) ; 
01774 
01775   if( MinNorm >= 0.0) 
01776     return(NormOne_);
01777 #endif
01778 
01779   if(!Filled()) 
01780     EPETRA_CHK_ERR(-1); // Matrix must be filled.
01781 
01782   Epetra_Vector x(DomainMap()); // Need temp vector for column sums
01783   
01784   double* xp = (double*)x.Values();
01785   Epetra_MultiVector* x_tmp = 0;
01786   int NumCols = NumMyCols();
01787   
01788 
01789   // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
01790   if(Importer() != 0) {
01791     x_tmp = new Epetra_Vector(ColMap()); // Create temporary import vector if needed
01792     xp = (double*)x_tmp->Values();
01793   }
01794   int i, j;
01795 
01796   for(i = 0; i < NumCols; i++) 
01797     xp[i] = 0.0;
01798 
01799   for(i = 0; i < NumMyRows_; i++) {
01800     int     NumEntries = NumMyEntries(i);
01801     int*    ColIndices = Graph().Indices(i);
01802     double* RowValues  = Values(i);
01803     for(j = 0; j < NumEntries; j++) 
01804       xp[ColIndices[j]] += std::abs(RowValues[j]);
01805   }
01806   if(Importer() != 0) {
01807     x.PutScalar(0.0);
01808     EPETRA_CHK_ERR(x.Export(*x_tmp, *Importer(), Add)); // Fill x with Values from temp vector
01809   }
01810   x.MaxValue(&NormOne_); // Find max
01811   if(x_tmp != 0) 
01812     delete x_tmp;
01813   UpdateFlops(NumGlobalNonzeros());
01814   return(NormOne_);
01815 }
01816 //=============================================================================
01817 double Epetra_CrsMatrix::NormFrobenius() const {
01818 
01819 #if 0
01820   //
01821   //  Commenting this section out disables caching, ie. 
01822   //  causes the norm to be computed each time NormFrobenius is called.  
01823   //  See bug #1151 for a full discussion.  
01824   //
01825   double MinNorm ; 
01826   Comm().MinAll( &NormFrob_, &MinNorm, 1 ) ; 
01827 
01828   if( MinNorm >= 0.0) 
01829     return(NormFrob_);
01830 #endif
01831 
01832   if(!Filled()) 
01833     EPETRA_CHK_ERR(-1); // Matrix must be filled.
01834 
01835   double local_sum = 0.0;
01836 
01837   for(int i = 0; i < NumMyRows_; i++) {
01838     int     NumEntries = NumMyEntries(i);
01839     double* RowValues  = Values(i);
01840     for(int j = 0; j < NumEntries; j++) {
01841       local_sum += RowValues[j]*RowValues[j];
01842     }
01843   }
01844 
01845   double global_sum = 0.0;
01846   Comm().SumAll(&local_sum, &global_sum, 1);
01847 
01848   NormFrob_ = std::sqrt(global_sum);
01849 
01850   UpdateFlops(NumGlobalNonzeros());
01851 
01852   return(NormFrob_);
01853 }
01854 //=========================================================================
01855 int Epetra_CrsMatrix::CheckSizes(const Epetra_SrcDistObject & Source) {
01856   try {
01857     const Epetra_CrsMatrix & A = dynamic_cast<const Epetra_CrsMatrix &>(Source);
01858     if (!A.Graph().GlobalConstantsComputed()) EPETRA_CHK_ERR(-1); // Must have global constants to proceed
01859   }
01860   catch (...) {
01861     return(0); // No error at this point, object could be a RowMatrix
01862   }
01863   return(0);
01864 }
01865 //=========================================================================
01866 int Epetra_CrsMatrix::CopyAndPermute(const Epetra_SrcDistObject & Source,
01867              int NumSameIDs, 
01868              int NumPermuteIDs,
01869                                      int * PermuteToLIDs,
01870              int *PermuteFromLIDs,
01871                                      const Epetra_OffsetIndex * Indexor ) {
01872  
01873   try {
01874     const Epetra_CrsMatrix & A = dynamic_cast<const Epetra_CrsMatrix &>(Source);
01875     EPETRA_CHK_ERR(CopyAndPermuteCrsMatrix(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs,
01876              PermuteFromLIDs,Indexor));
01877   }
01878   catch (...) {
01879     try {
01880       const Epetra_RowMatrix & A = dynamic_cast<const Epetra_RowMatrix &>(Source);
01881       EPETRA_CHK_ERR(CopyAndPermuteRowMatrix(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs,
01882                PermuteFromLIDs,Indexor));
01883     }
01884     catch (...) {
01885       EPETRA_CHK_ERR(-1); // Incompatible SrcDistObject
01886     }
01887   }
01888   
01889   return(0);
01890 }
01891 
01892 //=========================================================================
01893 int Epetra_CrsMatrix::CopyAndPermuteCrsMatrix(const Epetra_CrsMatrix & A,
01894                                               int NumSameIDs, 
01895                 int NumPermuteIDs,
01896                                               int * PermuteToLIDs,
01897                 int *PermuteFromLIDs,
01898                                               const Epetra_OffsetIndex * Indexor) {
01899   
01900   int i, ierr;
01901   
01902   int Row, NumEntries;
01903   int maxNumEntries = A.MaxNumEntries();
01904   int * Indices = 0;
01905   double * values = 0;
01906 
01907   if (maxNumEntries>0 && A.IndicesAreLocal() ) { //Need Temp Space
01908     Indices = new int[maxNumEntries];
01909     values = new double[maxNumEntries];
01910   }
01911   
01912   // Do copy first
01913   if (NumSameIDs>0) {
01914     if (A.IndicesAreLocal()) {
01915       if (StaticGraph() || IndicesAreLocal()) {
01916         if(Indexor) {
01917           for (i=0; i<NumSameIDs; i++) {
01918       Row = GRID(i);
01919       EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(Row, maxNumEntries, NumEntries, values, Indices)); // Set pointers
01920             ierr = ReplaceOffsetValues(Row, NumEntries, values, Indexor->SameOffsets()[i]);
01921             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
01922           }
01923         }
01924         else {
01925           for (i=0; i<NumSameIDs; i++) {
01926       Row = GRID(i);
01927       EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(Row, maxNumEntries, NumEntries, values, Indices)); // Set pointers
01928             ierr = ReplaceGlobalValues(Row, NumEntries, values, Indices);
01929             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
01930           }
01931         }
01932       }
01933       else {
01934         if(Indexor) {
01935           for (i=0; i<NumSameIDs; i++) {
01936       Row = GRID(i);
01937       EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(Row, maxNumEntries, NumEntries, values, Indices)); // Set pointers
01938             ierr = InsertOffsetValues(Row, NumEntries, values, Indexor->SameOffsets()[i]);
01939             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
01940           }
01941         }
01942         else {
01943           for (i=0; i<NumSameIDs; i++) {
01944       Row = GRID(i);
01945       EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(Row, maxNumEntries, NumEntries, values, Indices)); // Set pointers
01946             ierr = InsertGlobalValues(Row, NumEntries, values, Indices);
01947             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
01948           }
01949         }
01950       } 
01951     }
01952     else { // A.IndicesAreGlobal()
01953       if (StaticGraph() || IndicesAreLocal()) {
01954         if(Indexor) {
01955           for (i=0; i<NumSameIDs; i++) {
01956       Row = GRID(i);
01957       EPETRA_CHK_ERR(A.ExtractGlobalRowView(Row, NumEntries, values, Indices)); // Set pointers
01958             ierr = ReplaceOffsetValues(Row, NumEntries, values, Indexor->SameOffsets()[i]);
01959             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
01960           }
01961         }
01962         else {
01963           for (i=0; i<NumSameIDs; i++) {
01964       Row = GRID(i);
01965       EPETRA_CHK_ERR(A.ExtractGlobalRowView(Row, NumEntries, values, Indices)); // Set pointers
01966             ierr = ReplaceGlobalValues(Row, NumEntries, values, Indices);
01967             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
01968           }
01969         }
01970       }
01971       else {
01972         if(Indexor) {
01973           for (i=0; i<NumSameIDs; i++) {
01974       Row = GRID(i);
01975       EPETRA_CHK_ERR(A.ExtractGlobalRowView(Row, NumEntries, values, Indices)); // Set pointers
01976             ierr = InsertOffsetValues(Row, NumEntries, values, Indexor->SameOffsets()[i]);
01977             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
01978           }
01979         }
01980         else {
01981           for (i=0; i<NumSameIDs; i++) {
01982       Row = GRID(i);
01983       EPETRA_CHK_ERR(A.ExtractGlobalRowView(Row, NumEntries, values, Indices)); // Set pointers
01984             ierr = InsertGlobalValues(Row, NumEntries, values, Indices);
01985             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
01986           }
01987         }
01988       } 
01989     }
01990   }
01991   
01992   // Do local permutation next
01993   int FromRow, ToRow;
01994   if (NumPermuteIDs>0) {
01995     if (A.IndicesAreLocal()) {
01996       if (StaticGraph() || IndicesAreLocal()) {
01997         if(Indexor) {
01998           for (i=0; i<NumPermuteIDs; i++) {
01999       FromRow = A.GRID(PermuteFromLIDs[i]);
02000       ToRow = GRID(PermuteToLIDs[i]);
02001       EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, maxNumEntries, NumEntries, values, Indices)); // Set pointers
02002       ierr = ReplaceOffsetValues(ToRow, NumEntries, values, Indexor->PermuteOffsets()[i]);
02003             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
02004           }
02005         }
02006         else {
02007           for (i=0; i<NumPermuteIDs; i++) {
02008       FromRow = A.GRID(PermuteFromLIDs[i]);
02009       ToRow = GRID(PermuteToLIDs[i]);
02010       EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, maxNumEntries, NumEntries, values, Indices)); // Set pointers
02011       ierr = ReplaceGlobalValues(ToRow, NumEntries, values, Indices);
02012             if( ierr<0 ) EPETRA_CHK_ERR(ierr);
02013           }
02014         }
02015       }
02016       else {
02017         if(Indexor) {
02018           for (i=0; i<NumPermuteIDs; i++) {
02019       FromRow = A.GRID(PermuteFromLIDs[i]);
02020       ToRow = GRID(PermuteToLIDs[i]);
02021       EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, maxNumEntries, NumEntries, values, Indices)); // Set pointers
02022       ierr = InsertOffsetValues(ToRow, NumEntries, values, Indexor->PermuteOffsets()[i]);
02023       if (ierr<0) EPETRA_CHK_ERR(ierr);
02024           }
02025         }
02026         else {
02027           for (i=0; i<NumPermuteIDs; i++) {
02028       FromRow = A.GRID(PermuteFromLIDs[i]);
02029       ToRow = GRID(PermuteToLIDs[i]);
02030       EPETRA_CHK_ERR(A.ExtractGlobalRowCopy(FromRow, maxNumEntries, NumEntries, values, Indices)); // Set pointers
02031       ierr = InsertGlobalValues(ToRow, NumEntries, values, Indices);
02032       if (ierr<0) EPETRA_CHK_ERR(ierr);
02033           }
02034         }
02035       }
02036     }
02037     else { // A.IndicesAreGlobal()
02038       if (StaticGraph() || IndicesAreLocal()) {
02039         if(Indexor) {
02040           for (i=0; i<NumPermuteIDs; i++) {
02041       FromRow = A.GRID(PermuteFromLIDs[i]);
02042       ToRow = GRID(PermuteToLIDs[i]);
02043       EPETRA_CHK_ERR(A.ExtractGlobalRowView(FromRow, NumEntries, values, Indices)); // Set pointers
02044       ierr = ReplaceOffsetValues(ToRow, NumEntries, values, Indexor->PermuteOffsets()[i]);
02045       if (ierr<0) EPETRA_CHK_ERR(ierr);
02046           }
02047         }
02048         else {
02049           for (i=0; i<NumPermuteIDs; i++) {
02050       FromRow = A.GRID(PermuteFromLIDs[i]);
02051       ToRow = GRID(PermuteToLIDs[i]);
02052       EPETRA_CHK_ERR(A.ExtractGlobalRowView(FromRow, NumEntries, values, Indices)); // Set pointers
02053       ierr = ReplaceGlobalValues(ToRow, NumEntries, values, Indices);
02054       if (ierr<0) EPETRA_CHK_ERR(ierr);
02055           }
02056         }
02057       }
02058       else {
02059         if(Indexor) {
02060           for (i=0; i<NumPermuteIDs; i++) {
02061       FromRow = A.GRID(PermuteFromLIDs[i]);
02062       ToRow = GRID(PermuteToLIDs[i]);
02063       EPETRA_CHK_ERR(A.ExtractGlobalRowView(FromRow, NumEntries, values, Indices)); // Set pointers
02064       ierr = InsertOffsetValues(ToRow, NumEntries, values, Indexor->PermuteOffsets()[i]);
02065       if (ierr<0) EPETRA_CHK_ERR(ierr);
02066           }
02067         }
02068         else {
02069           for (i=0; i<NumPermuteIDs; i++) {
02070       FromRow = A.GRID(PermuteFromLIDs[i]);
02071       ToRow = GRID(PermuteToLIDs[i]);
02072       EPETRA_CHK_ERR(A.ExtractGlobalRowView(FromRow, NumEntries, values, Indices)); // Set pointers
02073       ierr = InsertGlobalValues(ToRow, NumEntries, values, Indices);
02074       if (ierr<0) EPETRA_CHK_ERR(ierr);
02075           }
02076         }
02077       }
02078     }
02079   }
02080 
02081   if (maxNumEntries>0 && A.IndicesAreLocal() ) { // Delete Temp Space
02082     delete [] values;
02083     delete [] Indices;
02084   }
02085   
02086   return(0);
02087 }
02088 
02089 //=========================================================================
02090 int Epetra_CrsMatrix::CopyAndPermuteRowMatrix(const Epetra_RowMatrix & A,
02091                                               int NumSameIDs, 
02092                 int NumPermuteIDs,
02093                                               int * PermuteToLIDs,
02094                 int *PermuteFromLIDs,
02095                                               const Epetra_OffsetIndex * Indexor ) {
02096   
02097   int i, j, ierr;
02098   
02099   int Row, NumEntries;
02100   int FromRow, ToRow;
02101   int maxNumEntries = A.MaxNumEntries();
02102   int * Indices = 0;
02103   double * values = 0;
02104 
02105   if (maxNumEntries>0) {
02106     Indices = new int[maxNumEntries];
02107     values = new double[maxNumEntries]; // Must extract values even though we discard them
02108   }
02109   
02110   const Epetra_Map & rowMap = A.RowMatrixRowMap();
02111   const Epetra_Map & colMap = A.RowMatrixColMap();
02112 
02113   // Do copy first
02114   if (NumSameIDs>0) {
02115     if (StaticGraph() || IndicesAreLocal()) {
02116       if( Indexor ) {
02117         for (i=0; i<NumSameIDs; i++) {
02118           Row = GRID(i);
02119           int AlocalRow = rowMap.LID(Row);
02120           EPETRA_CHK_ERR(A.ExtractMyRowCopy(AlocalRow, maxNumEntries, NumEntries, values, Indices));
02121     ierr = ReplaceOffsetValues(Row, NumEntries, values, Indexor->SameOffsets()[i]);
02122           if (ierr<0) EPETRA_CHK_ERR(ierr);
02123         }
02124       }
02125       else {
02126         for (i=0; i<NumSameIDs; i++) {
02127           Row = GRID(i);
02128           int AlocalRow = rowMap.LID(Row);
02129           EPETRA_CHK_ERR(A.ExtractMyRowCopy(AlocalRow, maxNumEntries, NumEntries, values, Indices));
02130           for(j=0; j<NumEntries; ++j) {
02131             Indices[j] = LCID(colMap.GID(Indices[j]));
02132           }
02133     ierr = ReplaceMyValues(i, NumEntries, values, Indices);
02134           if (ierr<0) EPETRA_CHK_ERR(ierr);
02135         }
02136       }
02137     }
02138     else {
02139       if( Indexor ) {
02140         for (i=0; i<NumSameIDs; i++) {
02141           EPETRA_CHK_ERR(A.ExtractMyRowCopy(i, maxNumEntries, NumEntries, values, Indices));
02142           Row = GRID(i);
02143     ierr = InsertOffsetValues(Row, NumEntries, values, Indexor->SameOffsets()[i]); 
02144           if (ierr<0) EPETRA_CHK_ERR(ierr);
02145         }
02146       }
02147       else {
02148         for (i=0; i<NumSameIDs; i++) {
02149           EPETRA_CHK_ERR(A.ExtractMyRowCopy(i, maxNumEntries, NumEntries, values, Indices));
02150           Row = GRID(i);
02151           for( j=0; j<NumEntries; ++j ) Indices[j] = colMap.GID(Indices[j]); //convert to GIDs
02152     ierr = InsertGlobalValues(Row, NumEntries, values, Indices); 
02153           if (ierr<0) EPETRA_CHK_ERR(ierr);
02154         }
02155       }
02156     }
02157   }
02158   
02159   // Do local permutation next
02160   if (NumPermuteIDs>0) {
02161     if (StaticGraph() || IndicesAreLocal()) {
02162       if( Indexor ) {
02163         for (i=0; i<NumPermuteIDs; i++) {
02164           FromRow = PermuteFromLIDs[i];
02165           EPETRA_CHK_ERR(A.ExtractMyRowCopy(FromRow, maxNumEntries, NumEntries, values, Indices));
02166           ToRow = GRID(PermuteToLIDs[i]);
02167     ierr = ReplaceOffsetValues(ToRow, NumEntries, values, Indexor->PermuteOffsets()[i]);
02168           if (ierr<0) EPETRA_CHK_ERR(ierr);
02169         }
02170       }
02171       else {
02172         for (i=0; i<NumPermuteIDs; i++) {
02173           FromRow = PermuteFromLIDs[i];
02174           EPETRA_CHK_ERR(A.ExtractMyRowCopy(FromRow, maxNumEntries, NumEntries, values, Indices));
02175           ToRow = GRID(PermuteToLIDs[i]);
02176           for(j=0; j<NumEntries; ++j) {
02177             Indices[j] = LCID(colMap.GID(Indices[j]));
02178           }
02179     ierr = ReplaceMyValues(ToRow, NumEntries, values, Indices);
02180           if (ierr<0) EPETRA_CHK_ERR(ierr);
02181         }
02182       }
02183     }
02184     else {
02185       if( Indexor ) {
02186         for (i=0; i<NumPermuteIDs; i++) {
02187           FromRow = PermuteFromLIDs[i];
02188           EPETRA_CHK_ERR(A.ExtractMyRowCopy(FromRow, maxNumEntries, NumEntries, values, Indices));
02189           ToRow = GRID(PermuteToLIDs[i]);
02190     ierr = InsertOffsetValues(ToRow, NumEntries, values, Indexor->PermuteOffsets()[i]); 
02191           if (ierr<0) EPETRA_CHK_ERR(ierr);
02192         }
02193       }
02194       else {
02195         for (i=0; i<NumPermuteIDs; i++) {
02196           FromRow = PermuteFromLIDs[i];
02197           EPETRA_CHK_ERR(A.ExtractMyRowCopy(FromRow, maxNumEntries, NumEntries, values, Indices));
02198     for (j=0; j<NumEntries; j++) Indices[j] = colMap.GID(Indices[j]); // convert to GIDs
02199           ToRow = GRID(PermuteToLIDs[i]);
02200     ierr = InsertGlobalValues(ToRow, NumEntries, values, Indices); 
02201           if (ierr<0) EPETRA_CHK_ERR(ierr);
02202         }
02203       }
02204     }
02205   } 
02206 
02207   if (maxNumEntries>0) {
02208     delete [] values;
02209     delete [] Indices;
02210   }
02211   
02212   return(0);
02213 }
02214 
02215 //=========================================================================
02216 int Epetra_CrsMatrix::PackAndPrepare(const Epetra_SrcDistObject & Source, 
02217              int NumExportIDs,
02218                                      int * ExportLIDs,
02219              int & LenExports,
02220                                      char *& Exports,
02221              int & SizeOfPacket,
02222                                      int * Sizes,
02223                                      bool & VarSizes,
02224                                      Epetra_Distributor & Distor)
02225 {
02226   (void)Distor; 
02227   // Rest of work can be done using RowMatrix only  
02228   const Epetra_RowMatrix & A = dynamic_cast<const Epetra_RowMatrix &>(Source);
02229 
02230   VarSizes = true; //enable variable block size data comm
02231 
02232   int TotalSendLength = 0;
02233   int * IntSizes = 0; 
02234   if( NumExportIDs>0 ) IntSizes = new int[NumExportIDs];
02235 
02236   for( int i = 0; i < NumExportIDs; ++i )
02237   {    
02238     int NumEntries;
02239     A.NumMyRowEntries( ExportLIDs[i], NumEntries );
02240     // Will have NumEntries doubles, NumEntries +2 ints, pack them interleaved     Sizes[i] = NumEntries;
02241     Sizes[i] = NumEntries;
02242     IntSizes[i] = 1 + (((NumEntries+2)*(int)sizeof(int))/(int)sizeof(double));
02243     TotalSendLength += (Sizes[i]+IntSizes[i]);
02244   }    
02245          
02246   double * DoubleExports = 0; 
02247   SizeOfPacket = (int)sizeof(double);
02248        
02249   //setup buffer locally for memory management by this object
02250   if( TotalSendLength*SizeOfPacket > LenExports )
02251   {
02252     if( LenExports > 0 ) delete [] Exports;
02253     LenExports = TotalSendLength*SizeOfPacket;
02254     DoubleExports = new double[TotalSendLength];
02255     for( int i = 0; i < TotalSendLength; ++i ) DoubleExports[i] = 0.0;
02256     Exports = (char *) DoubleExports;
02257   } 
02258  
02259   int NumEntries;
02260   int * Indices;
02261   double * values;
02262   int FromRow; 
02263   double * valptr, * dintptr; 
02264   int * intptr;                         
02265  
02266   // Each segment of Exports will be filled by a packed row of information for each row as follows:
02267   // 1st int: GRID of row where GRID is the global row ID for the source matrix
02268   // next int:  NumEntries, Number of indices in row.
02269   // next NumEntries: The actual indices for the row.
02270  
02271   const Epetra_Map & rowMap = A.RowMatrixRowMap();
02272   const Epetra_Map & colMap = A.RowMatrixColMap();
02273  
02274   if( NumExportIDs > 0 )
02275   {
02276     int maxNumEntries = A.MaxNumEntries();
02277     dintptr = (double *) Exports;
02278     valptr = dintptr + IntSizes[0];
02279     intptr = (int *) dintptr;
02280     for (int i=0; i<NumExportIDs; i++)
02281     {
02282       FromRow = rowMap.GID(ExportLIDs[i]);
02283       intptr[0] = FromRow;
02284       values = valptr;
02285       Indices = intptr + 2;
02286       EPETRA_CHK_ERR(A.ExtractMyRowCopy(ExportLIDs[i], maxNumEntries, NumEntries, values, Indices));
02287       for (int j=0; j<NumEntries; j++) Indices[j] = colMap.GID(Indices[j]); // convert to GIDs
02288       intptr[1] = NumEntries; // Load second slot of segment
02289       if( i < (NumExportIDs-1) )
02290       {
02291         dintptr += (IntSizes[i]+Sizes[i]);
02292         valptr = dintptr + IntSizes[i+1];
02293         intptr = (int *) dintptr;
02294       }
02295     }
02296  
02297     for( int i = 0; i < NumExportIDs; ++i )
02298       Sizes[i] += IntSizes[i];
02299   }
02300  
02301   if( IntSizes ) delete [] IntSizes;
02302  
02303   return(0);
02304 }
02305 
02306 //=========================================================================
02307 int Epetra_CrsMatrix::UnpackAndCombine(const Epetra_SrcDistObject & Source, 
02308                int NumImportIDs,
02309                                        int * ImportLIDs, 
02310                                        int LenImports,
02311                char * Imports,
02312                                        int & SizeOfPacket, 
02313                Epetra_Distributor & Distor, 
02314                Epetra_CombineMode CombineMode,
02315                                        const Epetra_OffsetIndex * Indexor )
02316 {
02317   (void)Source;
02318   (void)LenImports;
02319   (void)SizeOfPacket;
02320   (void)Distor;
02321   if (NumImportIDs<=0) return(0);
02322   
02323   if (   CombineMode != Add
02324    && CombineMode != Insert
02325    && CombineMode != Zero )
02326     EPETRA_CHK_ERR(-1); //Unsupported CombineMode, defaults to Zero
02327 
02328   int NumEntries;
02329   int * Indices;
02330   double * values;
02331   int ToRow;
02332   int i, ierr;
02333   int IntSize;
02334   
02335   double * valptr, *dintptr;
02336   int * intptr;
02337 
02338   // Each segment of Exports will be filled by a packed row of information for each row as follows:
02339   // 1st int: GRID of row where GRID is the global row ID for the source matrix
02340   // next int:  NumEntries, Number of indices in row.
02341   // next NumEntries: The actual indices for the row.
02342 
02343   dintptr = (double *) Imports;
02344   intptr = (int *) dintptr;
02345   NumEntries = intptr[1];
02346   IntSize = 1 + (((NumEntries+2)*(int)sizeof(int))/(int)sizeof(double));
02347   valptr = dintptr + IntSize;
02348  
02349   for (i=0; i<NumImportIDs; i++)
02350   {
02351     ToRow = GRID(ImportLIDs[i]);
02352     assert((intptr[0])==ToRow); // Sanity check
02353     values = valptr;
02354     Indices = intptr + 2;
02355  
02356     if (CombineMode==Add) {
02357       if (StaticGraph() || IndicesAreLocal()) {
02358         if( Indexor )
02359           ierr = SumIntoOffsetValues(ToRow, NumEntries, values, Indexor->RemoteOffsets()[i]);
02360         else
02361           ierr = SumIntoGlobalValues(ToRow, NumEntries, values, Indices);
02362       }
02363       else {
02364         if( Indexor )
02365           ierr = InsertOffsetValues(ToRow, NumEntries, values, Indexor->RemoteOffsets()[i]);
02366         else
02367           ierr = InsertGlobalValues(ToRow, NumEntries, values, Indices);
02368       }
02369       if (ierr<0) EPETRA_CHK_ERR(ierr);
02370     }
02371     else if (CombineMode==Insert) {
02372       if (StaticGraph() || IndicesAreLocal()) {
02373         if( Indexor )
02374           ierr = ReplaceOffsetValues(ToRow, NumEntries, values, Indexor->RemoteOffsets()[i]);
02375         else
02376           ierr = ReplaceGlobalValues(ToRow, NumEntries, values, Indices);
02377       }
02378       else {
02379         if( Indexor )
02380           ierr = InsertOffsetValues(ToRow, NumEntries, values, Indexor->RemoteOffsets()[i]);
02381         else
02382           ierr = InsertGlobalValues(ToRow, NumEntries, values, Indices);
02383       }
02384       if (ierr<0) EPETRA_CHK_ERR(ierr);
02385     }
02386  
02387     if( i < (NumImportIDs-1) )
02388     {
02389       dintptr += IntSize + NumEntries;
02390       intptr = (int *) dintptr;
02391       NumEntries = intptr[1];
02392       IntSize = 1 + (((NumEntries+2)*(int)sizeof(int))/(int)sizeof(double));
02393       valptr = dintptr + IntSize;
02394     }
02395   }
02396 
02397   return(0);
02398 }
02399 
02400 //=========================================================================
02401 
02402 void Epetra_CrsMatrix::Print(ostream& os) const {
02403   int MyPID = RowMap().Comm().MyPID();
02404   int NumProc = RowMap().Comm().NumProc();
02405 
02406   for (int iproc=0; iproc < NumProc; iproc++) {
02407     if (MyPID==iproc) {
02408       /*      const Epetra_fmtflags olda = os.setf(ios::right,ios::adjustfield);
02409         const Epetra_fmtflags oldf = os.setf(ios::scientific,ios::floatfield);
02410         const int             oldp = os.precision(12); */
02411       if (MyPID==0) {
02412   os <<  "\nNumber of Global Rows        = "; os << NumGlobalRows(); os << endl;
02413   os <<    "Number of Global Cols        = "; os << NumGlobalCols(); os << endl;
02414   os <<    "Number of Global Diagonals   = "; os << NumGlobalDiagonals(); os << endl;
02415   os <<    "Number of Global Nonzeros    = "; os << NumGlobalNonzeros(); os << endl;
02416   os <<    "Global Maximum Num Entries   = "; os << GlobalMaxNumEntries(); os << endl;
02417   if (LowerTriangular()) os <<    " ** Matrix is Lower Triangular **"; os << endl;
02418   if (UpperTriangular()) os <<    " ** Matrix is Upper Triangular **"; os << endl;
02419   if (NoDiagonal())      os <<    " ** Matrix has no diagonal     **"; os << endl; os << endl;
02420       }
02421       
02422       os <<  "\nNumber of My Rows        = "; os << NumMyRows(); os << endl;
02423       os <<    "Number of My Cols        = "; os << NumMyCols(); os << endl;
02424       os <<    "Number of My Diagonals   = "; os << NumMyDiagonals(); os << endl;
02425       os <<    "Number of My Nonzeros    = "; os << NumMyNonzeros(); os << endl;
02426       os <<    "My Maximum Num Entries   = "; os << MaxNumEntries(); os << endl; os << endl;
02427 
02428       os << flush;
02429       
02430       // Reset os flags
02431       
02432       /*      os.setf(olda,ios::adjustfield);
02433         os.setf(oldf,ios::floatfield);
02434         os.precision(oldp); */
02435     }
02436     // Do a few global ops to give I/O a chance to complete
02437     Comm().Barrier();
02438     Comm().Barrier();
02439     Comm().Barrier();
02440   }
02441   
02442   {for (int iproc=0; iproc < NumProc; iproc++) {
02443     if (MyPID==iproc) {
02444       int NumMyRows1 = NumMyRows();
02445       int MaxNumIndices = MaxNumEntries();
02446       int * Indices  = new int[MaxNumIndices];
02447       double * values  = new double[MaxNumIndices];
02448       int NumIndices;
02449       int i, j;
02450       
02451       if (MyPID==0) {
02452   os.width(8);
02453   os <<  "   Processor ";
02454   os.width(10);
02455   os <<  "   Row Index ";
02456   os.width(10);
02457   os <<  "   Col Index ";
02458   os.width(20);
02459   os <<  "   Value     ";
02460   os << endl;
02461       }
02462       for (i=0; i<NumMyRows1; i++) {
02463   int Row = GRID(i); // Get global row number
02464   ExtractGlobalRowCopy(Row, MaxNumIndices, NumIndices, values, Indices);
02465         
02466   for (j = 0; j < NumIndices ; j++) {   
02467     os.width(8);
02468     os <<  MyPID ; os << "    ";  
02469     os.width(10);
02470     os <<  Row ; os << "    ";  
02471     os.width(10);
02472     os <<  Indices[j]; os << "    ";
02473     os.width(20);
02474     os <<  values[j]; os << "    ";
02475     os << endl;
02476   }
02477       }
02478       
02479       delete [] Indices;
02480       delete [] values;
02481       
02482       os << flush;
02483       
02484     }
02485     // Do a few global ops to give I/O a chance to complete
02486     RowMap().Comm().Barrier();
02487     RowMap().Comm().Barrier();
02488     RowMap().Comm().Barrier();
02489   }}
02490   
02491   return;
02492 }
02493 //=============================================================================
02494 int Epetra_CrsMatrix::Multiply(bool TransA, const Epetra_Vector& x, Epetra_Vector& y) const {
02495 
02496 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS
02497   TEUCHOS_FUNC_TIME_MONITOR("Epetra_CrsMatrix::Multiply(TransA,x,y)");
02498 #endif
02499 
02500   //
02501   // This function forms the product y = A * x or y = A' * x
02502   //
02503 
02504   if(!Filled()) 
02505     EPETRA_CHK_ERR(-1); // Matrix must be filled.
02506 
02507   double* xp = (double*) x.Values();
02508   double* yp = (double*) y.Values();
02509 
02510   Epetra_Vector * xcopy = 0;
02511   if (&x==&y && Importer()==0 && Exporter()==0) {
02512     xcopy = new Epetra_Vector(x);
02513     xp = (double *) xcopy->Values();
02514   }
02515   UpdateImportVector(1); // Refresh import and output vectors if needed
02516   UpdateExportVector(1);
02517 
02518   if(!TransA) {
02519 
02520     // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
02521     if(Importer() != 0) {
02522       EPETRA_CHK_ERR(ImportVector_->Import(x, *Importer(), Insert));
02523       xp = (double*) ImportVector_->Values();
02524     }
02525     
02526     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
02527     if(Exporter() != 0)  yp = (double*) ExportVector_->Values();
02528     
02529     // Do actual computation
02530     GeneralMV(xp, yp);
02531 
02532     if(Exporter() != 0) {
02533       y.PutScalar(0.0); // Make sure target is zero
02534       EPETRA_CHK_ERR(y.Export(*ExportVector_, *Exporter(), Add)); // Fill y with Values from export vector
02535     }
02536     // Handle case of rangemap being a local replicated map
02537     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(y.Reduce());
02538   }
02539   
02540   else { // Transpose operation
02541 
02542     // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
02543     if(Exporter() != 0) {
02544       EPETRA_CHK_ERR(ExportVector_->Import(x, *Exporter(), Insert));
02545       xp = (double*) ExportVector_->Values();
02546     }
02547 
02548     // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
02549     if(Importer() != 0) yp = (double*) ImportVector_->Values();
02550 
02551     // Do actual computation
02552     GeneralMTV(xp, yp);
02553 
02554     if(Importer() != 0) {
02555       y.PutScalar(0.0); // Make sure target is zero
02556       EPETRA_CHK_ERR(y.Export(*ImportVector_, *Importer(), Add)); // Fill y with Values from export vector
02557     }
02558     // Handle case of rangemap being a local replicated map
02559     if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(y.Reduce());
02560   }
02561 
02562 
02563   UpdateFlops(2 * NumGlobalNonzeros());
02564   if (xcopy!=0) {
02565     delete xcopy;
02566     EPETRA_CHK_ERR(1); // Return positive code to alert the user about needing extra copy of x
02567     return(1);
02568   }
02569   return(0);
02570 }
02571 
02572 //=============================================================================
02573 int Epetra_CrsMatrix::Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
02574 
02575 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS
02576   TEUCHOS_FUNC_TIME_MONITOR("Epetra_CrsMatrix::Multiply(TransA,X,Y)");
02577 #endif
02578 
02579 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
02580   Teuchos::RCP<Teuchos::FancyOStream>
02581     out = Teuchos::VerboseObjectBase::getDefaultOStream();
02582   Teuchos::OSTab tab(out);
02583   if(Epetra_CrsMatrixTraceDumpMultiply) {
02584     *out << std::boolalpha;
02585     *out << "\nEntering Epetra_CrsMatrix::Multipy("<<TransA<<",X,Y) ...\n";
02586     if(!TransA) {
02587       *out << "\nDomainMap =\n";
02588       this->DomainMap().Print(Teuchos::OSTab(out).o());
02589     }
02590     else {
02591       *out << "\nRangeMap =\n";
02592       this->RangeMap().Print(Teuchos::OSTab(out).o());
02593     }
02594     *out << "\nInitial input X with " << ( TransA ? "RangeMap" : "DomainMap" ) << " =\n\n";
02595     X.Print(Teuchos::OSTab(out).o());
02596   }
02597 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
02598 
02599   //
02600   // This function forms the product Y = A * Y or Y = A' * X
02601   //
02602   if(!Filled()) {
02603     EPETRA_CHK_ERR(-1); // Matrix must be filled.
02604   }
02605 
02606   int NumVectors = X.NumVectors();
02607   if (NumVectors!=Y.NumVectors()) {
02608     EPETRA_CHK_ERR(-2); // Need same number of vectors in each MV
02609   }
02610 
02611   double** Xp = (double**) X.Pointers();
02612   double** Yp = (double**) Y.Pointers();
02613 
02614   int LDX = X.ConstantStride() ? X.Stride() : 0;
02615   int LDY = Y.ConstantStride() ? Y.Stride() : 0;
02616 
02617   Epetra_MultiVector * Xcopy = 0;
02618   if (&X==&Y && Importer()==0 && Exporter()==0) {
02619     Xcopy = new Epetra_MultiVector(X);
02620     Xp = (double **) Xcopy->Pointers();
02621     LDX = Xcopy->ConstantStride() ? Xcopy->Stride() : 0;
02622   }
02623   UpdateImportVector(NumVectors); // Make sure Import and Export Vectors are compatible
02624   UpdateExportVector(NumVectors);
02625 
02626   if (!TransA) {
02627 
02628     // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
02629     if (Importer()!=0) {
02630       EPETRA_CHK_ERR(ImportVector_->Import(X, *Importer(), Insert));
02631       Xp = (double**)ImportVector_->Pointers();
02632       LDX = ImportVector_->ConstantStride() ? ImportVector_->Stride() : 0;
02633 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
02634       if(Epetra_CrsMatrixTraceDumpMultiply) {
02635         *out << "\nColMap =\n";
02636         this->ColMap().Print(Teuchos::OSTab(out).o());
02637         *out << "\nX after import from DomainMap to ColMap =\n\n";
02638         ImportVector_->Print(Teuchos::OSTab(out).o());
02639       }
02640 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
02641     }
02642 
02643     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
02644     if (Exporter()!=0) {
02645       Yp = (double**)ExportVector_->Pointers();
02646       LDY = ExportVector_->ConstantStride() ? ExportVector_->Stride() : 0;
02647     }
02648 
02649     // Do actual computation
02650     if (NumVectors==1)
02651       GeneralMV(*Xp, *Yp);
02652     else
02653       GeneralMM(Xp, LDX, Yp, LDY, NumVectors);
02654 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
02655     if(Epetra_CrsMatrixTraceDumpMultiply) {
02656       *out << "\nRowMap =\n";
02657       this->RowMap().Print(Teuchos::OSTab(out).o());
02658       *out << "\nY after local mat-vec where Y has RowMap =\n\n";
02659       if(Exporter()!=0)
02660         ExportVector_->Print(Teuchos::OSTab(out).o());
02661       else
02662         Y.Print(Teuchos::OSTab(out).o());
02663     }
02664 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
02665     if (Exporter()!=0) {
02666       Y.PutScalar(0.0);  // Make sure target is zero
02667       Y.Export(*ExportVector_, *Exporter(), Add); // Fill Y with Values from export vector
02668 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
02669       if(Epetra_CrsMatrixTraceDumpMultiply) {
02670         *out << "\nRangeMap =\n";
02671         this->RangeMap().Print(Teuchos::OSTab(out).o());
02672         *out << "\nY after export from RowMap to RangeMap = \n\n";
02673         Y.Print(Teuchos::OSTab(out).o());
02674       }
02675 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
02676     }
02677     // Handle case of rangemap being a local replicated map
02678     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Y.Reduce());
02679   }
02680   else { // Transpose operation
02681 
02682     // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
02683 
02684     if (Exporter()!=0) {
02685       EPETRA_CHK_ERR(ExportVector_->Import(X, *Exporter(), Insert));
02686       Xp = (double**)ExportVector_->Pointers();
02687       LDX = ExportVector_->ConstantStride() ? ExportVector_->Stride() : 0;
02688 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
02689       if(Epetra_CrsMatrixTraceDumpMultiply) {
02690         *out << "\nRowMap =\n";
02691         this->RowMap().Print(Teuchos::OSTab(out).o());
02692         *out << "\nX after import from RangeMap to RowMap =\n\n";
02693         ExportVector_->Print(Teuchos::OSTab(out).o());
02694       }
02695 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
02696     }
02697 
02698     // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
02699     if (Importer()!=0) {
02700       Yp = (double**)ImportVector_->Pointers();
02701       LDY = ImportVector_->ConstantStride() ? ImportVector_->Stride() : 0;
02702     }
02703 
02704     // Do actual computation
02705     if (NumVectors==1)
02706       GeneralMTV(*Xp, *Yp);
02707     else
02708       GeneralMTM(Xp, LDX, Yp, LDY, NumVectors);
02709 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
02710     if(Epetra_CrsMatrixTraceDumpMultiply) {
02711       *out << "\nColMap =\n";
02712       this->ColMap().Print(Teuchos::OSTab(out).o());
02713       *out << "\nY after local transpose mat-vec where Y has ColMap =\n\n";
02714       if(Importer()!=0)
02715         ImportVector_->Print(Teuchos::OSTab(out).o());
02716       else
02717         Y.Print(Teuchos::OSTab(out).o());
02718     }
02719 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
02720     if (Importer()!=0) {
02721       Y.PutScalar(0.0);  // Make sure target is zero
02722       EPETRA_CHK_ERR(Y.Export(*ImportVector_, *Importer(), Add)); // Fill Y with Values from export vector
02723 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
02724       if(Epetra_CrsMatrixTraceDumpMultiply) {
02725         *out << "\nDomainMap =\n";
02726         this->DomainMap().Print(Teuchos::OSTab(out).o());
02727         *out << "\nY after export from ColMap to DomainMap =\n\n";
02728         Y.Print(Teuchos::OSTab(out).o());
02729       }
02730 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
02731     }
02732     // Handle case of rangemap being a local replicated map
02733     if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1)  EPETRA_CHK_ERR(Y.Reduce());
02734   }
02735 
02736   UpdateFlops(2*NumVectors*NumGlobalNonzeros());
02737   if (Xcopy!=0) {
02738     delete Xcopy;
02739     EPETRA_CHK_ERR(1); // Return positive code to alert the user about needing extra copy of X
02740     return(1);
02741   }
02742 
02743 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
02744   if(Epetra_CrsMatrixTraceDumpMultiply) {
02745     *out << "\nFinal output Y is the last Y printed above!\n";
02746     *out << "\nLeaving Epetra_CrsMatrix::Multipy("<<TransA<<",X,Y) ...\n";
02747   }
02748 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
02749 
02750   return(0);
02751 }
02752 //=======================================================================================================
02753 void Epetra_CrsMatrix::UpdateImportVector(int NumVectors) const {    
02754   if(Importer() != 0) {
02755     if(ImportVector_ != 0) {
02756       if(ImportVector_->NumVectors() != NumVectors) { 
02757   delete ImportVector_; 
02758   ImportVector_= 0;
02759       }
02760     }
02761     if(ImportVector_ == 0) 
02762       ImportVector_ = new Epetra_MultiVector(ColMap(),NumVectors); // Create import vector if needed
02763   }
02764   return;
02765 }
02766 //=======================================================================================================
02767 void Epetra_CrsMatrix::UpdateExportVector(int NumVectors) const {    
02768   if(Exporter() != 0) {
02769     if(ExportVector_ != 0) {
02770       if(ExportVector_->NumVectors() != NumVectors) { 
02771   delete ExportVector_; 
02772   ExportVector_= 0;
02773       }
02774     }
02775     if(ExportVector_ == 0) 
02776       ExportVector_ = new Epetra_MultiVector(RowMap(),NumVectors); // Create Export vector if needed
02777   }
02778   return;
02779 }
02780 //=======================================================================================================
02781 void Epetra_CrsMatrix::GeneralMV(double * x, double * y)  const {
02782   
02783 if (StorageOptimized() && Graph().StorageOptimized()) {
02784 
02785   double * values = All_Values();
02786   int * Indices = Graph().All_Indices();
02787   int * IndexOffset = Graph().IndexOffset();
02788 #ifdef EPETRA_HAVE_OMP
02789   const int numMyRows = NumMyRows_;
02790 #pragma omp parallel for default(none) shared(IndexOffset,values,Indices,y,x)
02791      for (int row=0; row<numMyRows; ++row)
02792         {
02793     const int curOffset = IndexOffset[row];
02794           const double *val_ptr    = values+curOffset;
02795           const int    *colnum_ptr = Indices+curOffset;
02796     double s = 0.;
02797     const double *const val_end_of_row = &values[IndexOffset[row+1]];
02798     while (val_ptr != val_end_of_row)
02799       s += *val_ptr++ * x[*colnum_ptr++];
02800     y[row] = s;
02801   }
02802 #else
02803 #ifndef Epetra_ENABLE_CASK
02804        const double *val_ptr    = values;
02805        const int    *colnum_ptr = Indices;
02806        double       * dst_ptr = y;
02807        for (int row=0; row<NumMyRows_; ++row)
02808   {
02809     double s = 0.;
02810     const double *const val_end_of_row = &values[IndexOffset[row+1]];
02811     while (val_ptr != val_end_of_row)
02812       s += *val_ptr++ * x[*colnum_ptr++];
02813     *dst_ptr++ = s;
02814   }
02815 #else
02816        cask_csr_dax_new(NumMyRows_, IndexOffset, Indices,
02817                         values, x, y, cask);
02818 #endif // Epetra_ENABLE_CASK
02819 #endif // EPETRA_HAVE_OMP
02820 
02821     return;
02822   }
02823   else if (!StorageOptimized() && !Graph().StorageOptimized()) {
02824 
02825 
02826     int* NumEntriesPerRow = Graph().NumIndicesPerRow();
02827     int** Indices = Graph().Indices();
02828     double** srcValues = Values();
02829         const int numMyRows = NumMyRows_;
02830 
02831     // Do actual computation
02832 #ifdef EPETRA_HAVE_OMP
02833 #pragma omp parallel for default(none) shared(NumEntriesPerRow,Indices,srcValues,y,x)
02834 #endif
02835     for(int i = 0; i < numMyRows; i++) {
02836       int     NumEntries = NumEntriesPerRow[i];
02837       int*    RowIndices = Indices[i];
02838       double* RowValues  = srcValues[i];
02839       double sum = 0.0;
02840       for(int j = 0; j < NumEntries; j++) 
02841   sum += *RowValues++ * x[*RowIndices++];
02842       
02843       y[i] = sum;
02844       
02845     }
02846   } 
02847   else { // Case where StorageOptimized is incompatible:  Use general accessors.
02848 
02849     const int numMyRows = NumMyRows_;
02850 
02851     // Do actual computation
02852 #ifdef EPETRA_HAVE_OMP
02853 #pragma omp parallel for default(none) shared(x,y)
02854 #endif
02855     for(int i = 0; i < numMyRows; i++) {
02856       int     NumEntries = NumMyEntries(i);
02857       int*    RowIndices = Graph().Indices(i);
02858       double* RowValues  = Values(i);
02859       double sum = 0.0;
02860       for(int j = 0; j < NumEntries; j++) 
02861   sum += *RowValues++ * x[*RowIndices++];
02862       
02863       y[i] = sum;
02864       
02865     }
02866   } 
02867   return;
02868 }
02869 //=======================================================================================================
02870 void Epetra_CrsMatrix::GeneralMTV(double * x, double * y) const {
02871 
02872   int NumCols = NumMyCols();
02873 #if !defined(FORTRAN_DISABLED) || defined(Epetra_ENABLE_CASK)
02874   if (StorageOptimized() && Graph().StorageOptimized()) {
02875     double * values = All_Values_;
02876     int * Indices = Graph().All_Indices();
02877     int * IndexOffset = Graph().IndexOffset();
02878 #ifndef Epetra_ENABLE_CASK
02879    int ione = 1;
02880    EPETRA_DCRSMV_F77(&ione, &NumMyRows_, &NumCols, values, Indices, IndexOffset, x, y);
02881 #else
02882    cask_csr_datx( NumMyRows_, NumCols, IndexOffset,  Indices,  values,x ,y );
02883 #endif
02884 
02885     return;
02886   }
02887 #endif // FORTRAN_DISABLED
02888   for(int i = 0; i < NumCols; i++) 
02889     y[i] = 0.0; // Initialize y for transpose multiply
02890 
02891   if (StorageOptimized() && Graph().StorageOptimized()) {
02892     double * values = All_Values_;
02893     int * Indices = Graph().All_Indices();
02894     int * IndexOffset = Graph().IndexOffset();
02895     for(int i = 0; i < NumMyRows_; ++i) {
02896       int prevOffset = *IndexOffset++;
02897       int NumEntries = *IndexOffset - prevOffset;
02898       double xi = x[i];
02899       for(int j = 0; j < NumEntries; j++) 
02900   y[*Indices++] += *values++ * xi;
02901     }
02902   }
02903   else if (!StorageOptimized() && !Graph().StorageOptimized()) {
02904 
02905     int* NumEntriesPerRow = Graph().NumIndicesPerRow();
02906     int** Indices = Graph().Indices();
02907     double** srcValues = Values();
02908     
02909     for(int i = 0; i < NumMyRows_; i++) {
02910       int     NumEntries = *NumEntriesPerRow++;
02911       int*    RowIndices = *Indices++;
02912       double* RowValues  = *srcValues++;
02913       double xi = x[i];
02914       for(int j = 0; j < NumEntries; j++) 
02915   y[*RowIndices++] += *RowValues++ * xi;
02916     }
02917   }
02918   else { // Case where StorageOptimized is incompatible:  Use general accessors.
02919   
02920     for(int i = 0; i < NumMyRows_; i++) {
02921       int     NumEntries = NumMyEntries(i);
02922       int*    RowIndices = Graph().Indices(i);
02923       double* RowValues  = Values(i);
02924       double xi = x[i];
02925       for(int j = 0; j < NumEntries; j++) 
02926   y[*RowIndices++] += *RowValues++ * xi;
02927     }
02928   }
02929 
02930   return;
02931 }
02932 //=======================================================================================================
02933 void Epetra_CrsMatrix::GeneralMM(double ** X, int LDX, double ** Y, int LDY, int NumVectors) const {
02934 
02935 #if !defined(FORTRAN_DISABLED) || defined(Epetra_ENABLE_CASK)
02936   if (StorageOptimized() && Graph().StorageOptimized()) {
02937     double * values = All_Values_;
02938     int * Indices = Graph().All_Indices();
02939     int * IndexOffset = Graph().IndexOffset();
02940 
02941     if (LDX!=0 && LDY!=0) {
02942 #ifndef Epetra_ENABLE_CASK
02943     int izero = 0;
02944     EPETRA_DCRSMM_F77(&izero, &NumMyRows_, &NumMyRows_, values, Indices, IndexOffset, *X, &LDX, *Y, &LDY, &NumVectors);
02945 #else
02946     cask_csr_dgesmm_new(0, 1.0, NumMyRows_, NumMyRows_,  NumVectors,
02947                     IndexOffset, Indices, values, *X, LDX, 0.0,  *Y, LDY,cask);
02948 #endif
02949     return;
02950     }
02951 
02952     double ** const xp = X;
02953     double ** const yp = Y;
02954     const int numMyRows = NumMyRows_;
02955 #ifdef EPETRA_HAVE_OMP
02956 #pragma omp parallel for default(none) shared(IndexOffset,Indices,values,NumVectors)
02957 #endif
02958     for (int i=0; i < numMyRows; i++) {
02959       int prevOffset = IndexOffset[i];
02960       int NumEntries = IndexOffset[i+1] - prevOffset;
02961       int *    RowIndices = Indices+prevOffset;
02962       double * RowValues  = values+prevOffset;
02963       for (int k=0; k<NumVectors; k++) {
02964   double sum = 0.0;
02965   const double * const x = xp[k];
02966   double * const y = yp[k];
02967   for (int j=0; j < NumEntries; j++) sum += RowValues[j] * x[RowIndices[j]];
02968   y[i] = sum;
02969       }
02970     }
02971   }
02972   else if (!StorageOptimized() && !Graph().StorageOptimized()) {
02973 #else
02974   if (!StorageOptimized() && !Graph().StorageOptimized()) {
02975 #endif
02976 
02977     int* NumEntriesPerRow = Graph().NumIndicesPerRow();
02978     int** Indices = Graph().Indices();
02979     double** srcValues = Values();
02980     double ** const xp = X;
02981     double ** const yp = Y;
02982     const int numMyRows = NumMyRows_;
02983 
02984 #ifdef EPETRA_HAVE_OMP
02985 #pragma omp parallel for default(none) shared(NumEntriesPerRow,Indices,srcValues,NumVectors)
02986 #endif
02987     for (int i=0; i < numMyRows; i++) {
02988       int      NumEntries = NumEntriesPerRow[i];
02989       int *    RowIndices = Indices[i];
02990       double * RowValues  = srcValues[i];
02991       for (int k=0; k<NumVectors; k++) {
02992   double sum = 0.0;
02993   const double * const x = xp[k];
02994   double * const y = yp[k];
02995   for (int j=0; j < NumEntries; j++) sum += RowValues[j] * x[RowIndices[j]];
02996   y[i] = sum;
02997       }
02998     }
02999   }
03000   else {
03001 
03002     double ** const xp = X;
03003     double ** const yp = Y;
03004     const int numMyRows = NumMyRows_;
03005 #ifdef EPETRA_HAVE_OMP
03006 #pragma omp parallel for default(none) shared(NumVectors)
03007 #endif
03008     for (int i=0; i < numMyRows; i++) {
03009       int     NumEntries = NumMyEntries(i);
03010       int*    RowIndices = Graph().Indices(i);
03011       double* RowValues  = Values(i);
03012       for (int k=0; k<NumVectors; k++) {
03013   double sum = 0.0;
03014   double * x = xp[k];
03015   for (int j=0; j < NumEntries; j++) sum += RowValues[j] * x[RowIndices[j]];
03016   yp[k][i] = sum;
03017       }
03018     }
03019   }
03020   return;
03021 }
03022 //=======================================================================================================
03023 void Epetra_CrsMatrix::GeneralMTM(double ** X, int LDX, double ** Y, int LDY, int NumVectors)  const{
03024 
03025   int NumCols = NumMyCols();
03026 #if !defined(FORTRAN_DISABLED) || defined(Epetra_ENABLE_CASK)
03027   if (StorageOptimized() && Graph().StorageOptimized()) {
03028     if (LDX!=0 && LDY!=0) {
03029       double * values = All_Values_;
03030       int * Indices = Graph().All_Indices();
03031       int * IndexOffset = Graph().IndexOffset();
03032 
03033 #ifndef Epetra_ENABLE_CASK
03034       int ione = 1;
03035       EPETRA_DCRSMM_F77(&ione, &NumMyRows_, &NumCols, values, Indices, IndexOffset, *X, &LDX, *Y, &LDY, &NumVectors);
03036 #else
03037       cask_csr_dgesmm_new(1, 1.0, NumMyRows_, NumCols,  NumVectors,
03038                           IndexOffset, Indices, values, *X, LDX, 0.0,
03039                           *Y, LDY, cask);
03040 #endif
03041       return;
03042     }
03043   }
03044 #endif
03045   for (int k=0; k<NumVectors; k++) 
03046     for (int i=0; i < NumCols; i++) 
03047       Y[k][i] = 0.0; // Initialize y for transpose multiply
03048   
03049   if (StorageOptimized() && Graph().StorageOptimized()) {
03050     double * values = All_Values_;
03051     int * Indices = Graph().All_Indices();
03052     int * IndexOffset = Graph().IndexOffset();
03053     for (int i=0; i < NumMyRows_; i++) {
03054       int prevOffset = *IndexOffset++;
03055       int NumEntries = *IndexOffset - prevOffset;
03056       int *    RowIndices = Indices+prevOffset;
03057       double * RowValues  = values+prevOffset;
03058       
03059       for (int k=0; k<NumVectors; k++) {
03060   double * y = Y[k];
03061   double * x = X[k];
03062   for (int j=0; j < NumEntries; j++) 
03063     y[RowIndices[j]] += RowValues[j] * x[i];
03064       }
03065     }
03066   }
03067   else if (!StorageOptimized() && !Graph().StorageOptimized()) {
03068     
03069     int* NumEntriesPerRow = Graph().NumIndicesPerRow();
03070     int** Indices = Graph().Indices();
03071     double** srcValues = Values();
03072 
03073     for (int i=0; i < NumMyRows_; i++) {
03074       int      NumEntries = *NumEntriesPerRow++;
03075       int *    RowIndices = *Indices++;
03076       double * RowValues  = *srcValues++;
03077       for (int k=0; k<NumVectors; k++) {
03078   double * y = Y[k];
03079   double * x = X[k];
03080   for (int j=0; j < NumEntries; j++) 
03081     y[RowIndices[j]] += RowValues[j] * x[i];
03082       }
03083     }
03084   }
03085   else { // Case where StorageOptimized is incompatible:  Use general accessors.
03086     
03087     for (int i=0; i < NumMyRows_; i++) {
03088       int     NumEntries = NumMyEntries(i);
03089       int*    RowIndices = Graph().Indices(i);
03090       double* RowValues  = Values(i);
03091       for (int k=0; k<NumVectors; k++) {
03092   double * y = Y[k];
03093   double * x = X[k];
03094   for (int j=0; j < NumEntries; j++) 
03095     y[RowIndices[j]] += RowValues[j] * x[i];
03096       }
03097     }
03098   }
03099   return;
03100 }
03101 //=======================================================================================================
03102 void Epetra_CrsMatrix::GeneralSV(bool Upper, bool Trans, bool UnitDiagonal, double * xp, double * yp)  const {
03103 
03104 
03105   int i, j, j0;
03106 
03107 #if !defined(FORTRAN_DISABLED) || defined(Epetra_ENABLE_CASK)
03108   if (StorageOptimized() && Graph().StorageOptimized() && ((UnitDiagonal && NoDiagonal())|| (!UnitDiagonal && !NoDiagonal()))) {
03109     double * values = All_Values();
03110     int * Indices = Graph().All_Indices();
03111     int * IndexOffset = Graph().IndexOffset();
03112 
03113     int iupper = Upper ? 1:0;
03114     int itrans = Trans ? 1:0;
03115     int udiag =  UnitDiagonal ? 1:0;
03116     int nodiag = NoDiagonal() ? 1:0;
03117     int xysame = (xp==yp) ? 1:0;
03118 
03119 #ifndef Epetra_ENABLE_CASK
03120     EPETRA_DCRSSV_F77( &iupper, &itrans, &udiag, &nodiag, &NumMyRows_, &NumMyRows_, values, Indices, IndexOffset, xp, yp, &xysame);
03121 #else
03122     cask_csr_dtrsv_new( iupper, itrans, udiag, nodiag, 0, xysame, NumMyRows_,
03123                     NumMyRows_, IndexOffset, Indices, values, xp, yp, cask);
03124 #endif
03125     return;
03126   }
03127   //=================================================================
03128   else { // !StorageOptimized()
03129   //=================================================================
03130 #endif    
03131     if (!Trans) {
03132       
03133       if (Upper) {
03134   
03135   j0 = 1;
03136   if (NoDiagonal()) 
03137     j0--; // Include first term if no diagonal
03138   for (i=NumMyRows_-1; i >=0; i--) {
03139     int      NumEntries = NumMyEntries(i);
03140     int *    RowIndices = Graph().Indices(i);
03141     double * RowValues  = Values(i);
03142     double sum = 0.0;
03143     for (j=j0; j < NumEntries; j++) 
03144       sum += RowValues[j] * yp[RowIndices[j]];
03145     
03146     if (UnitDiagonal) 
03147       yp[i] = xp[i] - sum;
03148     else 
03149       yp[i] = (xp[i] - sum)/RowValues[0];
03150     
03151   }
03152       }
03153       else {
03154   j0 = 1;
03155   if (NoDiagonal())
03156     j0--; // Include first term if no diagonal
03157   for (i=0; i < NumMyRows_; i++) {
03158     int      NumEntries = NumMyEntries(i) - j0;
03159     int *    RowIndices = Graph().Indices(i);
03160     double * RowValues  = Values(i);
03161     double sum = 0.0;
03162     for (j=0; j < NumEntries; j++) 
03163       sum += RowValues[j] * yp[RowIndices[j]];
03164     
03165     if (UnitDiagonal) 
03166       yp[i] = xp[i] - sum;
03167     else 
03168       yp[i] = (xp[i] - sum)/RowValues[NumEntries];
03169     
03170   }
03171       }
03172     }
03173     
03174     // ***********  Transpose case *******************************
03175     
03176     else {
03177       
03178       if (xp!=yp) 
03179   for (i=0; i < NumMyRows_; i++) 
03180     yp[i] = xp[i]; // Initialize y for transpose solve
03181       
03182       if (Upper) {
03183   
03184   j0 = 1;
03185   if (NoDiagonal()) 
03186     j0--; // Include first term if no diagonal
03187   
03188   for (i=0; i < NumMyRows_; i++) {
03189     int      NumEntries = NumMyEntries(i);
03190     int *    RowIndices = Graph().Indices(i);
03191     double * RowValues  = Values(i);
03192     if (!UnitDiagonal) 
03193       yp[i] = yp[i]/RowValues[0];
03194     double ytmp = yp[i];
03195     for (j=j0; j < NumEntries; j++) 
03196       yp[RowIndices[j]] -= RowValues[j] * ytmp;
03197   }
03198       }
03199       else {
03200   
03201   j0 = 1;
03202   if (NoDiagonal()) 
03203     j0--; // Include first term if no diagonal
03204   
03205   for (i=NumMyRows_-1; i >= 0; i--) {
03206     int      NumEntries = NumMyEntries(i) - j0;
03207     int *    RowIndices = Graph().Indices(i);
03208     double * RowValues  = Values(i);
03209     if (!UnitDiagonal) 
03210       yp[i] = yp[i]/RowValues[NumEntries];
03211     double ytmp = yp[i];
03212     for (j=0; j < NumEntries; j++) 
03213       yp[RowIndices[j]] -= RowValues[j] * ytmp;
03214   }
03215       }
03216       
03217     }
03218 #if !defined(FORTRAN_DISABLED) || defined(Epetra_ENABLE_CASK)
03219   }
03220 #endif
03221   return;
03222 }
03223 //=======================================================================================================
03224 void Epetra_CrsMatrix::GeneralSM(bool Upper, bool Trans, bool UnitDiagonal, double ** Xp, int LDX, double ** Yp, int LDY, int NumVectors)  const{
03225 
03226   int i, j, j0, k;
03227   double diag = 0.0;
03228 
03229   if (StorageOptimized() && Graph().StorageOptimized()) {
03230     double * values = All_Values();
03231     int * Indices = Graph().All_Indices();
03232     int * IndexOffset = Graph().IndexOffset();
03233 #if !defined(FORTRAN_DISABLED) || defined(Epetra_ENABLE_CASK)
03234     if (LDX!=0 && LDY!=0 && ((UnitDiagonal && NoDiagonal()) || (!UnitDiagonal && !NoDiagonal()))) {
03235       int iupper = Upper ? 1:0;
03236       int itrans = Trans ? 1:0;
03237       int udiag =  UnitDiagonal ? 1:0;
03238       int nodiag = NoDiagonal() ? 1:0;
03239       int xysame = (Xp==Yp) ? 1:0;
03240 
03241 #ifndef Epetra_ENABLE_CASK
03242       EPETRA_DCRSSM_F77( &iupper, &itrans, &udiag, &nodiag, &NumMyRows_, &NumMyRows_, values, Indices, IndexOffset, 
03243        *Xp, &LDX, *Yp, &LDY, &xysame, &NumVectors);
03244 #else
03245       cask_csr_dtrsm( iupper, itrans, udiag, nodiag, 0, xysame,  NumMyRows_, 
03246                       NumMyRows_, NumVectors, IndexOffset, Indices, values, 
03247                       *Xp, LDX, *Yp, LDY);
03248 #endif
03249       return;
03250     }
03251 #endif
03252     if(!Trans) {   
03253       if(Upper) {   
03254   j0 = 1;
03255   if(NoDiagonal()) 
03256     j0--; // Include first term if no diagonal
03257   for(i = NumMyRows_ - 1; i >= 0; i--) {
03258     int Offset = IndexOffset[i];
03259     int      NumEntries = IndexOffset[i+1]-Offset;
03260     int *    RowIndices = Indices+Offset;
03261     double * RowValues  = values+Offset;
03262     if(!UnitDiagonal) 
03263       diag = 1.0/RowValues[0]; // Take inverse of diagonal once for later use
03264     for(k = 0; k < NumVectors; k++) {
03265       double sum = 0.0;
03266       for(j = j0; j < NumEntries; j++) 
03267         sum += RowValues[j] * Yp[k][RowIndices[j]];
03268           
03269       if(UnitDiagonal) 
03270         Yp[k][i] = Xp[k][i] - sum;
03271       else 
03272         Yp[k][i] = (Xp[k][i] - sum) * diag;
03273     }
03274   }
03275       }
03276       else {
03277   j0 = 1;
03278   if(NoDiagonal()) 
03279     j0--; // Include first term if no diagonal
03280   for(i = 0; i < NumMyRows_; i++) {
03281     int Offset = IndexOffset[i];
03282     int      NumEntries = IndexOffset[i+1]-Offset - j0;
03283     int *    RowIndices = Indices+Offset;
03284     double * RowValues  = values+Offset;
03285     if(!UnitDiagonal)
03286       diag = 1.0/RowValues[NumEntries]; // Take inverse of diagonal once for later use
03287     for(k = 0; k < NumVectors; k++) {
03288       double sum = 0.0;
03289       for(j = 0; j < NumEntries; j++) 
03290         sum += RowValues[j] * Yp[k][RowIndices[j]];
03291           
03292       if(UnitDiagonal) 
03293         Yp[k][i] = Xp[k][i] - sum;
03294       else 
03295         Yp[k][i] = (Xp[k][i] - sum)*diag;
03296     }
03297   }
03298       }
03299     }
03300     // ***********  Transpose case *******************************
03301 
03302     else {
03303       for(k = 0; k < NumVectors; k++) 
03304   if(Yp[k] != Xp[k]) 
03305     for(i = 0; i < NumMyRows_; i++)
03306       Yp[k][i] = Xp[k][i]; // Initialize y for transpose multiply
03307     
03308       if(Upper) {
03309   j0 = 1;
03310   if(NoDiagonal()) 
03311     j0--; // Include first term if no diagonal
03312       
03313   for(i = 0; i < NumMyRows_; i++) {
03314     int Offset = IndexOffset[i];
03315     int      NumEntries = IndexOffset[i+1]-Offset;
03316     int *    RowIndices = Indices+Offset;
03317     double * RowValues  = values+Offset;
03318     if(!UnitDiagonal) 
03319       diag = 1.0/RowValues[0]; // Take inverse of diagonal once for later use
03320     for(k = 0; k < NumVectors; k++) {
03321       if(!UnitDiagonal) 
03322         Yp[k][i] = Yp[k][i]*diag;
03323       double ytmp = Yp[k][i];
03324       for(j = j0; j < NumEntries; j++) 
03325         Yp[k][RowIndices[j]] -= RowValues[j] * ytmp;
03326     }
03327   }
03328       }
03329       else {
03330   j0 = 1;
03331   if(NoDiagonal()) 
03332     j0--; // Include first term if no diagonal  
03333   for(i = NumMyRows_ - 1; i >= 0; i--) {
03334     int Offset = IndexOffset[i];
03335     int      NumEntries = IndexOffset[i+1]-Offset - j0;
03336     int *    RowIndices = Indices+Offset;
03337     double * RowValues  = values+Offset;
03338     if(!UnitDiagonal) 
03339       diag = 1.0/RowValues[NumEntries]; // Take inverse of diagonal once for later use
03340     for(k = 0; k < NumVectors; k++) {
03341       if(!UnitDiagonal)  
03342         Yp[k][i] = Yp[k][i]*diag;
03343       double ytmp = Yp[k][i];
03344       for(j = 0; j < NumEntries; j++)
03345         Yp[k][RowIndices[j]] -= RowValues[j] * ytmp;
03346     }
03347   }
03348       }
03349     }
03350   }
03351     // ========================================================
03352   else { // !StorageOptimized()
03353     // ========================================================
03354 
03355     if(!Trans) {   
03356       if(Upper) {   
03357   j0 = 1;
03358   if(NoDiagonal()) 
03359     j0--; // Include first term if no diagonal
03360   for(i = NumMyRows_ - 1; i >= 0; i--) {
03361     int     NumEntries = NumMyEntries(i);
03362     int*    RowIndices = Graph().Indices(i);
03363     double* RowValues  = Values(i);
03364     if(!UnitDiagonal) 
03365       diag = 1.0/RowValues[0]; // Take inverse of diagonal once for later use
03366     for(k = 0; k < NumVectors; k++) {
03367       double sum = 0.0;
03368       for(j = j0; j < NumEntries; j++) 
03369         sum += RowValues[j] * Yp[k][RowIndices[j]];
03370           
03371       if(UnitDiagonal) 
03372         Yp[k][i] = Xp[k][i] - sum;
03373       else 
03374         Yp[k][i] = (Xp[k][i] - sum) * diag;
03375     }
03376   }
03377       }
03378       else {
03379   j0 = 1;
03380   if(NoDiagonal()) 
03381     j0--; // Include first term if no diagonal
03382   for(i = 0; i < NumMyRows_; i++) {
03383     int     NumEntries = NumMyEntries(i) - j0;
03384     int*    RowIndices = Graph().Indices(i);
03385     double* RowValues  = Values(i);
03386     if(!UnitDiagonal)
03387       diag = 1.0/RowValues[NumEntries]; // Take inverse of diagonal once for later use
03388     for(k = 0; k < NumVectors; k++) {
03389       double sum = 0.0;
03390       for(j = 0; j < NumEntries; j++) 
03391         sum += RowValues[j] * Yp[k][RowIndices[j]];
03392           
03393       if(UnitDiagonal) 
03394         Yp[k][i] = Xp[k][i] - sum;
03395       else 
03396         Yp[k][i] = (Xp[k][i] - sum)*diag;
03397     }
03398   }
03399       }
03400     }
03401     // ***********  Transpose case *******************************
03402 
03403     else {
03404       for(k = 0; k < NumVectors; k++) 
03405   if(Yp[k] != Xp[k]) 
03406     for(i = 0; i < NumMyRows_; i++)
03407       Yp[k][i] = Xp[k][i]; // Initialize y for transpose multiply
03408     
03409       if(Upper) {
03410   j0 = 1;
03411   if(NoDiagonal()) 
03412     j0--; // Include first term if no diagonal
03413       
03414   for(i = 0; i < NumMyRows_; i++) {
03415     int     NumEntries = NumMyEntries(i);
03416     int*    RowIndices = Graph().Indices(i);
03417     double* RowValues  = Values(i);
03418     if(!UnitDiagonal) 
03419       diag = 1.0/RowValues[0]; // Take inverse of diagonal once for later use
03420     for(k = 0; k < NumVectors; k++) {
03421       if(!UnitDiagonal) 
03422         Yp[k][i] = Yp[k][i]*diag;
03423       double ytmp = Yp[k][i];
03424       for(j = j0; j < NumEntries; j++) 
03425         Yp[k][RowIndices[j]] -= RowValues[j] * ytmp;
03426     }
03427   }
03428       }
03429       else {
03430   j0 = 1;
03431   if(NoDiagonal()) 
03432     j0--; // Include first term if no diagonal  
03433   for(i = NumMyRows_ - 1; i >= 0; i--) {
03434     int     NumEntries = NumMyEntries(i) - j0;
03435     int*    RowIndices = Graph().Indices(i);
03436     double* RowValues  = Values(i);
03437     if(!UnitDiagonal) 
03438       diag = 1.0/RowValues[NumEntries]; // Take inverse of diagonal once for later use
03439     for(k = 0; k < NumVectors; k++) {
03440       if(!UnitDiagonal)  
03441         Yp[k][i] = Yp[k][i]*diag;
03442       double ytmp = Yp[k][i];
03443       for(j = 0; j < NumEntries; j++)
03444         Yp[k][RowIndices[j]] -= RowValues[j] * ytmp;
03445     }
03446   }
03447       }
03448     }
03449   }
03450   
03451   return;
03452 }
03453 //=============================================================================
03454 int Epetra_CrsMatrix::Multiply1(bool TransA, const Epetra_Vector& x, Epetra_Vector& y) const {
03455 
03456 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS
03457   TEUCHOS_FUNC_TIME_MONITOR("Epetra_CrsMatrix::Multiply1(TransA,x,y)");
03458 #endif
03459 
03460   //
03461   // This function forms the product y = A * x or y = A' * x
03462   //
03463 
03464   if(!Filled()) 
03465     EPETRA_CHK_ERR(-1); // Matrix must be filled.
03466 
03467   int i, j;
03468   double* xp = (double*) x.Values();
03469   double* yp = (double*) y.Values();
03470   int NumMyCols_ = NumMyCols();
03471 
03472   if(!TransA) {
03473 
03474     // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
03475     if(Importer() != 0) {
03476       if(ImportVector_ != 0) {
03477   if(ImportVector_->NumVectors() != 1) { 
03478     delete ImportVector_; 
03479     ImportVector_= 0;
03480   }
03481       }
03482       if(ImportVector_ == 0) 
03483   ImportVector_ = new Epetra_MultiVector(ColMap(),1); // Create import vector if needed
03484       EPETRA_CHK_ERR(ImportVector_->Import(x, *Importer(), Insert));
03485       xp = (double*) ImportVector_->Values();
03486     }
03487     
03488     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
03489     if(Exporter() != 0) {
03490       if(ExportVector_ != 0) {
03491   if(ExportVector_->NumVectors() != 1) { 
03492     delete ExportVector_; 
03493     ExportVector_= 0;
03494   }
03495       }
03496       if(ExportVector_ == 0) 
03497   ExportVector_ = new Epetra_MultiVector(RowMap(),1); // Create Export vector if needed
03498       yp = (double*) ExportVector_->Values();
03499     }
03500     
03501     // Do actual computation
03502     for(i = 0; i < NumMyRows_; i++) {
03503       int     NumEntries = NumMyEntries(i);
03504       int*    RowIndices = Graph().Indices(i);
03505       double* RowValues  = Values(i);
03506       double sum = 0.0;
03507       for(j = 0; j < NumEntries; j++) 
03508   sum += RowValues[j] * xp[RowIndices[j]];
03509       
03510       yp[i] = sum;
03511       
03512     }
03513     if(Exporter() != 0) {
03514       y.PutScalar(0.0); // Make sure target is zero
03515       EPETRA_CHK_ERR(y.Export(*ExportVector_, *Exporter(), Add)); // Fill y with Values from export vector
03516     }
03517     // Handle case of rangemap being a local replicated map
03518     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(y.Reduce());
03519   }
03520   
03521   else { // Transpose operation
03522 
03523     // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
03524     if(Exporter() != 0) {
03525       if(ExportVector_ != 0) {
03526   if(ExportVector_->NumVectors() != 1) { 
03527     delete ExportVector_; 
03528     ExportVector_= 0;
03529   }
03530       }
03531       if(ExportVector_ == 0) 
03532   ExportVector_ = new Epetra_MultiVector(RowMap(),1); // Create Export vector if needed
03533       EPETRA_CHK_ERR(ExportVector_->Import(x, *Exporter(), Insert));
03534       xp = (double*) ExportVector_->Values();
03535     }
03536 
03537     // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
03538     if(Importer() != 0) {
03539       if(ImportVector_ != 0) {
03540   if(ImportVector_->NumVectors() != 1) { 
03541     delete ImportVector_; 
03542     ImportVector_= 0;
03543   }
03544       }
03545       if(ImportVector_ == 0) 
03546   ImportVector_ = new Epetra_MultiVector(ColMap(),1); // Create import vector if needed
03547       yp = (double*) ImportVector_->Values();
03548     }
03549 
03550     // Do actual computation
03551     for(i = 0; i < NumMyCols_; i++) 
03552       yp[i] = 0.0; // Initialize y for transpose multiply
03553         
03554     for(i = 0; i < NumMyRows_; i++) {
03555       int     NumEntries = NumMyEntries(i);
03556       int*    RowIndices = Graph().Indices(i);
03557       double* RowValues  = Values(i);
03558       for(j = 0; j < NumEntries; j++) 
03559   yp[RowIndices[j]] += RowValues[j] * xp[i];
03560     }
03561     if(Importer() != 0) {
03562       y.PutScalar(0.0); // Make sure target is zero
03563       EPETRA_CHK_ERR(y.Export(*ImportVector_, *Importer(), Add)); // Fill y with Values from export vector
03564     }
03565     // Handle case of rangemap being a local replicated map
03566     if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(y.Reduce());
03567   }
03568 
03569   UpdateFlops(2 * NumGlobalNonzeros());
03570   return(0);
03571 }
03572 //=============================================================================
03573 int Epetra_CrsMatrix::Multiply1(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
03574 
03575 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS
03576   TEUCHOS_FUNC_TIME_MONITOR("Epetra_CrsMatrix::Multiply1(TransA,X,Y)");
03577 #endif
03578 
03579   //
03580   // This function forms the product Y = A * Y or Y = A' * X
03581   //
03582   if((X.NumVectors() == 1) && (Y.NumVectors() == 1)) {
03583     double* xp = (double*) X[0];
03584     double* yp = (double*) Y[0];
03585     Epetra_Vector x(View, X.Map(), xp);
03586     Epetra_Vector y(View, Y.Map(), yp);
03587     EPETRA_CHK_ERR(Multiply1(TransA, x, y));
03588     return(0);
03589   }
03590   if(!Filled()) {
03591     EPETRA_CHK_ERR(-1); // Matrix must be filled.
03592   }
03593 
03594   int i, j, k;
03595 
03596   double** Xp = (double**) X.Pointers();
03597   double** Yp = (double**) Y.Pointers();
03598 
03599   int NumVectors = X.NumVectors();
03600   int NumMyCols_ = NumMyCols();
03601 
03602 
03603   // Need to better manage the Import and Export vectors:
03604   // - Need accessor functions
03605   // - Need to make the NumVector match (use a View to do this)
03606   // - Need to look at RightScale and ColSum routines too.
03607 
03608   if (!TransA) {
03609 
03610     // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
03611     if (Importer()!=0) {
03612       if (ImportVector_!=0) {
03613   if (ImportVector_->NumVectors()!=NumVectors) { 
03614     delete ImportVector_; ImportVector_= 0;}
03615       }
03616       if (ImportVector_==0) 
03617   ImportVector_ = new Epetra_MultiVector(ColMap(),NumVectors); // Create import vector if needed
03618       EPETRA_CHK_ERR(ImportVector_->Import(X, *Importer(), Insert));
03619       Xp = (double**)ImportVector_->Pointers();
03620     }
03621 
03622     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
03623     if (Exporter()!=0) {
03624       if (ExportVector_!=0) {
03625   if (ExportVector_->NumVectors()!=NumVectors) { 
03626     delete ExportVector_; ExportVector_= 0;}
03627       }
03628       if (ExportVector_==0) 
03629   ExportVector_ = new Epetra_MultiVector(RowMap(),NumVectors); // Create Export vector if needed
03630       Yp = (double**)ExportVector_->Pointers();
03631     }
03632 
03633     // Do actual computation
03634 
03635     for (i=0; i < NumMyRows_; i++) {
03636       int      NumEntries = NumMyEntries(i);
03637       int *    RowIndices = Graph().Indices(i);
03638       double * RowValues  = Values(i);
03639       for (k=0; k<NumVectors; k++) {
03640   double sum = 0.0;
03641   for (j=0; j < NumEntries; j++) sum += RowValues[j] * Xp[k][RowIndices[j]];
03642   Yp[k][i] = sum;
03643       }
03644     }
03645     if (Exporter()!=0) {
03646       Y.PutScalar(0.0); // Make sure target is zero
03647       Y.Export(*ExportVector_, *Exporter(), Add); // Fill Y with Values from export vector
03648     }
03649     // Handle case of rangemap being a local replicated map
03650     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Y.Reduce());
03651   }
03652   else { // Transpose operation
03653     
03654 
03655     // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
03656 
03657     if (Exporter()!=0) {
03658       if (ExportVector_!=0) {
03659   if (ExportVector_->NumVectors()!=NumVectors) { 
03660     delete ExportVector_; ExportVector_= 0;}
03661       }
03662       if (ExportVector_==0) 
03663   ExportVector_ = new Epetra_MultiVector(RowMap(),NumVectors); // Create Export vector if needed
03664       EPETRA_CHK_ERR(ExportVector_->Import(X, *Exporter(), Insert));
03665       Xp = (double**)ExportVector_->Pointers();
03666     }
03667 
03668     // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
03669     if (Importer()!=0) {
03670       if (ImportVector_!=0) {
03671   if (ImportVector_->NumVectors()!=NumVectors) { 
03672     delete ImportVector_; ImportVector_= 0;}
03673       }
03674       if (ImportVector_==0) 
03675   ImportVector_ = new Epetra_MultiVector(ColMap(),NumVectors); // Create import vector if needed
03676       Yp = (double**)ImportVector_->Pointers();
03677     }
03678 
03679     // Do actual computation
03680 
03681 
03682 
03683     for (k=0; k<NumVectors; k++) 
03684       for (i=0; i < NumMyCols_; i++) 
03685   Yp[k][i] = 0.0; // Initialize y for transpose multiply
03686     
03687     for (i=0; i < NumMyRows_; i++) {
03688       int      NumEntries = NumMyEntries(i);
03689       int *    RowIndices = Graph().Indices(i);
03690       double * RowValues  = Values(i);
03691       for (k=0; k<NumVectors; k++) {
03692   for (j=0; j < NumEntries; j++) 
03693     Yp[k][RowIndices[j]] += RowValues[j] * Xp[k][i];
03694       }
03695     }
03696     if (Importer()!=0) {
03697       Y.PutScalar(0.0); // Make sure target is zero
03698       EPETRA_CHK_ERR(Y.Export(*ImportVector_, *Importer(), Add)); // Fill Y with Values from export vector
03699     }
03700     // Handle case of rangemap being a local replicated map
03701     if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1)  EPETRA_CHK_ERR(Y.Reduce());
03702   }
03703 
03704   UpdateFlops(2*NumVectors*NumGlobalNonzeros());
03705   return(0);
03706 }
03707 
03708 //=============================================================================
03709 int Epetra_CrsMatrix::Solve1(bool Upper, bool Trans, bool UnitDiagonal,
03710           const Epetra_Vector& x, Epetra_Vector& y) const
03711 {
03712 
03713 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS
03714   TEUCHOS_FUNC_TIME_MONITOR("Epetra_CrsMatrix::Solve1(Upper,Trans,UnitDiag,x,y)");
03715 #endif
03716 
03717   //
03718   // This function finds y such that Ly = x or Uy = x or the transpose cases.
03719   //
03720 
03721   if (!Filled()) {
03722     EPETRA_CHK_ERR(-1); // Matrix must be filled.
03723   }
03724 
03725   if ((Upper) && (!UpperTriangular())) 
03726     EPETRA_CHK_ERR(-2);
03727   if ((!Upper) && (!LowerTriangular())) 
03728     EPETRA_CHK_ERR(-3);
03729   if ((!UnitDiagonal) && (NoDiagonal())) 
03730     EPETRA_CHK_ERR(-4); // If matrix has no diagonal, we must use UnitDiagonal
03731   if ((!UnitDiagonal) && (NumMyDiagonals()<NumMyRows_)) 
03732     EPETRA_CHK_ERR(-5); // Need each row to have a diagonal
03733       
03734 
03735   int i, j, j0;
03736   int * NumEntriesPerRow = Graph().NumIndicesPerRow();
03737   int ** Indices = Graph().Indices();
03738   double ** Vals = Values();
03739   int NumMyCols_ = NumMyCols();
03740 
03741   // If upper, point to last row
03742   if ((Upper && !Trans) || (!Upper && Trans)) {
03743     NumEntriesPerRow += NumMyRows_-1;
03744     Indices += NumMyRows_-1;
03745     Vals += NumMyRows_-1;
03746   }
03747     
03748   double *xp = (double*)x.Values();
03749   double *yp = (double*)y.Values();
03750 
03751   if (!Trans) {
03752 
03753     if (Upper) {
03754 
03755       j0 = 1;
03756       if (NoDiagonal()) 
03757   j0--; // Include first term if no diagonal
03758       for (i=NumMyRows_-1; i >=0; i--) {
03759   int      NumEntries = *NumEntriesPerRow--;
03760   int *    RowIndices = *Indices--;
03761   double * RowValues  = *Vals--;
03762   double sum = 0.0;
03763   for (j=j0; j < NumEntries; j++) 
03764     sum += RowValues[j] * yp[RowIndices[j]];
03765         
03766   if (UnitDiagonal) 
03767     yp[i] = xp[i] - sum;
03768   else 
03769     yp[i] = (xp[i] - sum)/RowValues[0];
03770         
03771       }
03772     }
03773     else {
03774       j0 = 1;
03775       if (NoDiagonal())
03776   j0--; // Include first term if no diagonal
03777       for (i=0; i < NumMyRows_; i++) {
03778   int      NumEntries = *NumEntriesPerRow++ - j0;
03779   int *    RowIndices = *Indices++;
03780   double * RowValues  = *Vals++;
03781   double sum = 0.0;
03782   for (j=0; j < NumEntries; j++) 
03783     sum += RowValues[j] * yp[RowIndices[j]];
03784         
03785   if (UnitDiagonal) 
03786     yp[i] = xp[i] - sum;
03787   else 
03788     yp[i] = (xp[i] - sum)/RowValues[NumEntries];
03789         
03790       }
03791     }
03792   }
03793   
03794   // ***********  Transpose case *******************************
03795   
03796   else {
03797 
03798     if (xp!=yp) 
03799       for (i=0; i < NumMyCols_; i++) 
03800   yp[i] = xp[i]; // Initialize y for transpose solve
03801     
03802     if (Upper) {
03803 
03804       j0 = 1;
03805       if (NoDiagonal()) 
03806   j0--; // Include first term if no diagonal
03807     
03808       for (i=0; i < NumMyRows_; i++) {
03809   int      NumEntries = *NumEntriesPerRow++;
03810   int *    RowIndices = *Indices++;
03811   double * RowValues  = *Vals++;
03812   if (!UnitDiagonal) 
03813     yp[i] = yp[i]/RowValues[0];
03814      double ytmp = yp[i];
03815   for (j=j0; j < NumEntries; j++) 
03816     yp[RowIndices[j]] -= RowValues[j] * ytmp;
03817       }
03818     }
03819     else {
03820       
03821       j0 = 1;
03822       if (NoDiagonal()) 
03823   j0--; // Include first term if no diagonal
03824     
03825       for (i=NumMyRows_-1; i >= 0; i--) {
03826   int      NumEntries = *NumEntriesPerRow-- - j0;
03827   int *    RowIndices = *Indices--;
03828   double * RowValues  = *Vals--;
03829   if (!UnitDiagonal) 
03830     yp[i] = yp[i]/RowValues[NumEntries];
03831      double ytmp = yp[i];
03832   for (j=0; j < NumEntries; j++) 
03833     yp[RowIndices[j]] -= RowValues[j] * ytmp;
03834       }
03835     }
03836     
03837   }
03838   UpdateFlops(2*NumGlobalNonzeros());
03839   return(0);
03840 }
03841 
03842 //=============================================================================
03843 int Epetra_CrsMatrix::Solve1(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
03844 
03845 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS
03846   TEUCHOS_FUNC_TIME_MONITOR("Epetra_CrsMatrix::Solve(Upper,Trans,UnitDiag,X,Y)");
03847 #endif
03848 
03849   //
03850   // This function find Y such that LY = X or UY = X or the transpose cases.
03851   //
03852   if((X.NumVectors() == 1) && (Y.NumVectors() == 1)) {
03853     double* xp = (double*) X[0];
03854     double* yp = (double*) Y[0];
03855     Epetra_Vector x(View, X.Map(), xp);
03856     Epetra_Vector y(View, Y.Map(), yp);
03857     EPETRA_CHK_ERR(Solve1(Upper, Trans, UnitDiagonal, x, y));
03858     return(0);
03859   }
03860   if(!Filled()) 
03861     EPETRA_CHK_ERR(-1); // Matrix must be filled.
03862 
03863   if((Upper) && (!UpperTriangular()))
03864     EPETRA_CHK_ERR(-2);
03865   if((!Upper) && (!LowerTriangular()))
03866     EPETRA_CHK_ERR(-3);
03867   if((!UnitDiagonal) && (NoDiagonal()))
03868     EPETRA_CHK_ERR(-4); // If matrix has no diagonal, we must use UnitDiagonal
03869   if((!UnitDiagonal) && (NumMyDiagonals()<NumMyRows_))
03870     EPETRA_CHK_ERR(-5); // Need each row to have a diagonal
03871 
03872   int i, j, j0, k;
03873   int* NumEntriesPerRow = Graph().NumIndicesPerRow();
03874   int** Indices = Graph().Indices();
03875   double** Vals = Values();
03876   double diag = 0.0;
03877 
03878   // If upper, point to last row
03879   if((Upper && !Trans) || (!Upper && Trans)) {
03880     NumEntriesPerRow += NumMyRows_-1;
03881     Indices += NumMyRows_-1;
03882     Vals += NumMyRows_-1;
03883   }
03884 
03885   double** Xp = (double**) X.Pointers();
03886   double** Yp = (double**) Y.Pointers();
03887 
03888   int NumVectors = X.NumVectors();
03889 
03890   if(!Trans) {   
03891     if(Upper) {   
03892       j0 = 1;
03893       if(NoDiagonal()) 
03894   j0--; // Include first term if no diagonal
03895       for(i = NumMyRows_ - 1; i >= 0; i--) {
03896   int     NumEntries = *NumEntriesPerRow--;
03897   int*    RowIndices = *Indices--;
03898   double* RowValues  = *Vals--;
03899   if(!UnitDiagonal) 
03900     diag = 1.0/RowValues[0]; // Take inverse of diagonal once for later use
03901   for(k = 0; k < NumVectors; k++) {
03902     double sum = 0.0;
03903     for(j = j0; j < NumEntries; j++) 
03904       sum += RowValues[j] * Yp[k][RowIndices[j]];
03905           
03906     if(UnitDiagonal) 
03907       Yp[k][i] = Xp[k][i] - sum;
03908     else 
03909       Yp[k][i] = (Xp[k][i] - sum) * diag;
03910   }
03911       }
03912     }
03913     else {
03914       j0 = 1;
03915       if(NoDiagonal()) 
03916   j0--; // Include first term if no diagonal
03917       for(i = 0; i < NumMyRows_; i++) {
03918   int     NumEntries = *NumEntriesPerRow++ - j0;
03919   int*    RowIndices = *Indices++;
03920   double* RowValues  = *Vals++;
03921   if(!UnitDiagonal)
03922     diag = 1.0/RowValues[NumEntries]; // Take inverse of diagonal once for later use
03923   for(k = 0; k < NumVectors; k++) {
03924     double sum = 0.0;
03925     for(j = 0; j < NumEntries; j++) 
03926       sum += RowValues[j] * Yp[k][RowIndices[j]];
03927           
03928     if(UnitDiagonal) 
03929       Yp[k][i] = Xp[k][i] - sum;
03930     else 
03931       Yp[k][i] = (Xp[k][i] - sum)*diag;
03932   }
03933       }
03934     }
03935   }
03936   // ***********  Transpose case *******************************
03937 
03938   else {
03939     for(k = 0; k < NumVectors; k++) 
03940       if(Yp[k] != Xp[k]) 
03941   for(i = 0; i < NumMyRows_; i++)
03942     Yp[k][i] = Xp[k][i]; // Initialize y for transpose multiply
03943     
03944     if(Upper) {
03945       j0 = 1;
03946       if(NoDiagonal()) 
03947   j0--; // Include first term if no diagonal
03948       
03949       for(i = 0; i < NumMyRows_; i++) {
03950   int     NumEntries = *NumEntriesPerRow++;
03951   int*    RowIndices = *Indices++;
03952   double* RowValues  = *Vals++;
03953   if(!UnitDiagonal) 
03954     diag = 1.0/RowValues[0]; // Take inverse of diagonal once for later use
03955   for(k = 0; k < NumVectors; k++) {
03956     if(!UnitDiagonal) 
03957       Yp[k][i] = Yp[k][i]*diag;
03958        double ytmp = Yp[k][i];
03959     for(j = j0; j < NumEntries; j++) 
03960       Yp[k][RowIndices[j]] -= RowValues[j] * ytmp;
03961   }
03962       }
03963     }
03964     else {
03965       j0 = 1;
03966       if(NoDiagonal()) 
03967   j0--; // Include first term if no diagonal  
03968       for(i = NumMyRows_ - 1; i >= 0; i--) {
03969   int     NumEntries = *NumEntriesPerRow-- - j0;
03970   int*    RowIndices = *Indices--;
03971   double* RowValues  = *Vals--;
03972      if (!UnitDiagonal)
03973        diag = 1.0/RowValues[NumEntries]; // Take inverse of diagonal once for later use
03974   for(k = 0; k < NumVectors; k++) {
03975     if(!UnitDiagonal)  
03976       Yp[k][i] = Yp[k][i]*diag;
03977        double ytmp = Yp[k][i];
03978     for(j = 0; j < NumEntries; j++)
03979       Yp[k][RowIndices[j]] -= RowValues[j] * ytmp;
03980         }
03981       }
03982     }
03983   }
03984   
03985   UpdateFlops(2 * NumVectors * NumGlobalNonzeros());
03986   return(0);
03987 }
03988 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines