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

Generated on Wed May 12 21:41:05 2010 for Epetra Package Browser (Single Doxygen Collection) by  doxygen 1.4.7