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