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

Generated on Thu Sep 18 12:37:57 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1