Epetra_FastCrsMatrix.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_FastCrsMatrix.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 <limits.h>
00039 
00040 //==============================================================================
00041 Epetra_FastCrsMatrix::Epetra_FastCrsMatrix(const Epetra_CrsMatrix & Matrix, bool UseFloats) 
00042   : CrsMatrix_(Matrix),
00043     Values_(0),
00044     NumMyRows_(Matrix.NumMyRows()),
00045     NumMyNonzeros(Matrix.NumMyNonzeros()),
00046     ImportVector_(0),
00047     ExportVector_(0),
00048     CV_(Copy)
00049 {
00050   if (!CrsMatrix_.Filled()) throw CrsMatrix_.ReportError("Input matrix must have called FillComplete()", -1);
00051   Allocate(UseFloats);
00052 }
00053 
00054 //==============================================================================
00055 Epetra_FastCrsMatrix::~Epetra_FastCrsMatrix(){
00056 
00057   if (Values_!=0 && ValuesAllocated_) delete [] Values_;
00058   if (FloatValues_!=0) delete [] FloatValues_;
00059   if (Indices_!=0) delete [] Indices_;
00060   if (ShortIndices_!=0) delete [] ShortIndices_;
00061 
00062   if (ImportVector_!=0) delete ImportVector_;
00063   ImportVector_=0;
00064   if (ExportVector_!=0) delete ExportVector_;
00065   ExportVector_=0;
00066 }
00067 
00068 //==============================================================================
00069 int Epetra_FastCrsMatrix::UpdateValues(const Epetra_CrsMatrix & Matrix) {
00070 }
00071 
00072 //==============================================================================
00073 int Epetra_FastCrsMatrix::Allocate(bool UseFloats) {
00074 
00075   int i, j;
00076   
00077   // First determine how to handle values.  Three possibilities:
00078   // 1) Values are contiguously stored in user matrix, UseFloats is false so we point to the values in
00079   //    the user matrix.
00080   // 2) Values are not contiguously stored, UseFloats is false, so we copy values into a contiguous double array.
00081   // 3) UseFloats is true so we create a single precision copy of matrix values.
00082 
00083   if (CrsMatrix.IndicesAreContiguous() && !UseFloats) {
00084     ValuesAllocated_ = false;
00085     Values_ = CrsMatrix_[0]; // First value in the user's matrix
00086   }
00087   else if (!UseFloats) {
00088     // Allocate Values array
00089     Values_ = new double[NumMyNonzeros_];
00090     double * ptr = Values_;
00091     for (i=0; i<NumMyRows_; i++) {
00092       int NumEntries = CrsMatrix_.NumMyEntries(i);
00093       double * rowi = CrsMatrix_[i];
00094       for (j=0; j< NumEntries; j++) *ptr++ = *rowi++; // Copy values to contiguous storage
00095     }
00096   }
00097   else { // UseFloats == true
00098     FloatValues_ = new float[NumMyNonzeros_];
00099     float * ptr = FloatValues_;
00100     for (i=0; i<NumMyRows_; i++) {
00101       int NumEntries = CrsMatrix_.NumMyEntries(i);
00102       double * rowi = CrsMatrix_[i];
00103       for (j=0; j< NumEntries; j++) *ptr++ = (float) *rowi++; // convert and copy values
00104     }
00105   }
00106 
00107   // Next determine how to handle integers.  Two possibilities:
00108   // 1) Local column range is within the range of unsigned short ints, so we copy the indices to an array of unsigned shorts.
00109   // 2) Local column range is outside range of unsigned shorts and we copy to array of ints.
00110   // In both cases we embed the nonzero count per row into the index array.
00111 
00112   if (CrsMatrix_.NumMyCols()>USHRT_MAX) {
00113     // Allocate Values array
00114     Indices_ = new int[NumMyNonzeros_+NumMyRows_];
00115     int * ptr = Indices_;
00116     for (i=0; i<NumMyRows_; i++) {
00117       int NumEntries = CrsMatrix_.NumMyEntries(i);
00118       int * rowi = CrsMatrix_.Graph()[i];
00119       *ptr++ = NumEntries; // Pack this value first
00120       for (j=0; j< NumEntries; j++) *ptr++ = *rowi++; // Copy values to contiguous storage
00121     }
00122   }
00123   else { // CrsMatrix_.NumMyCols()<=USHRT_MAX
00124     ShortIndices_ = new unsigned short[NumMyNonzeros_+NumMyRows_];
00125     unsigned short * ptr = ShortIndices_;
00126     for (i=0; i<NumMyRows_; i++) {
00127       int NumEntries = CrsMatrix_.NumMyEntries(i);
00128       unsigned short * rowi = CrsMatrix_Graph()[i];
00129       *ptr++ = NumEntries; // Pack this value first
00130       for (j=0; j< NumEntries; j++) *ptr++ = (unsigned short) *rowi++; // convert and copy indices
00131     }
00132   }
00133 
00134   return(0);
00135 }
00136 //=============================================================================
00137 int Epetra_FastCrsMatrix::Multiply(bool TransA, const Epetra_Vector& x, Epetra_Vector& y) const {
00138   //
00139   // This function forms the product y = A * x or y = A' * x
00140   //
00141 
00142   int i, j;
00143   double * xp = (double*)x.Values();
00144   double *yp = (double*)y.Values();
00145   int NumMyCols_ = NumMyCols();
00146 
00147 
00148   if (!TransA) {
00149 
00150     // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
00151     if (Importer()!=0) {
00152       if (ImportVector_!=0) {
00153   if (ImportVector_->NumVectors()!=1) { delete ImportVector_; ImportVector_= 0;}
00154       }
00155       if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),1); // Create import vector if needed
00156       ImportVector_->Import(x, *Importer(), Insert);
00157       xp = (double*)ImportVector_->Values();
00158     }
00159 
00160     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
00161     if (Exporter()!=0) {
00162       if (ExportVector_!=0) {
00163   if (ExportVector_->NumVectors()!=1) { delete ExportVector_; ExportVector_= 0;}
00164       }
00165       if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),1); // Create Export vector if needed
00166       yp = (double*)ExportVector_->Values();
00167     }
00168 
00169     // Do actual computation
00170 
00171     for (i=0; i < NumMyRows_; i++) {
00172       int      NumEntries = *NumEntriesPerRow++;
00173       int *    RowIndices = *Indices++;
00174       double * RowValues  = *Values++;
00175       double sum = 0.0;
00176       for (j=0; j < NumEntries; j++) sum += RowValues[j] * xp[RowIndices[j]];
00177 
00178       yp[i] = sum;
00179 
00180     }
00181     if (Exporter()!=0) y.Export(*ExportVector_, *Exporter(), Add); // Fill y with Values from export vector
00182   }
00183 
00184   else { // Transpose operation
00185 
00186 
00187     // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
00188 
00189     if (Exporter()!=0) {
00190       if (ExportVector_!=0) {
00191   if (ExportVector_->NumVectors()!=1) { delete ExportVector_; ExportVector_= 0;}
00192       }
00193       if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),1); // Create Export vector if needed
00194       ExportVector_->Import(x, *Exporter(), Insert);
00195       xp = (double*)ExportVector_->Values();
00196     }
00197 
00198     // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
00199     if (Importer()!=0) {
00200       if (ImportVector_!=0) {
00201   if (ImportVector_->NumVectors()!=1) { delete ImportVector_; ImportVector_= 0;}
00202       }
00203       if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),1); // Create import vector if needed
00204       yp = (double*)ImportVector_->Values();
00205     }
00206 
00207     // Do actual computation
00208 
00209     for (i=0; i < NumMyCols_; i++) yp[i] = 0.0; // Initialize y for transpose multiply
00210         
00211     for (i=0; i < NumMyRows_; i++) {
00212       int      NumEntries = *NumEntriesPerRow++;
00213       int *    RowIndices = *Indices++;
00214       double * RowValues  = *Values++;
00215       for (j=0; j < NumEntries; j++) yp[RowIndices[j]] += RowValues[j] * xp[i];
00216     }
00217     if (Importer()!=0) y.Export(*ImportVector_, *Importer(), Add); // Fill y with Values from export vector
00218   }
00219 
00220   UpdateFlops(2*NumGlobalNonzeros());
00221   return(0);
00222 }
00223 
00224 //=============================================================================
00225 int Epetra_FastCrsMatrix::Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
00226   //
00227   // This function forms the product Y = A * Y or Y = A' * X
00228   //
00229   if (X.NumVectors()==1 && Y.NumVectors()==1) {
00230     double * xp = (double *) X[0];
00231     double * yp = (double *) Y[0];
00232     Epetra_Vector x(View, X.Map(), xp);
00233     Epetra_Vector y(View, Y.Map(), yp);
00234     return(Multiply(TransA, x, y));
00235   }
00236   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
00237 
00238   int i, j, k;
00239   int * NumEntriesPerRow = NumEntriesPerRow_;
00240   int ** Indices = Indices_;
00241   double ** Values = Values_;
00242 
00243   double **Xp = (double**)X.Pointers();
00244   double **Yp = (double**)Y.Pointers();
00245 
00246   int NumVectors = X.NumVectors();
00247   int NumMyCols_ = NumMyCols();
00248 
00249 
00250   // Need to better manage the Import and Export vectors:
00251   // - Need accessor functions
00252   // - Need to make the NumVector match (use a View to do this)
00253   // - Need to look at RightScale and ColSum routines too.
00254 
00255   if (!TransA) {
00256 
00257     // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
00258     if (Importer()!=0) {
00259       if (ImportVector_!=0) {
00260   if (ImportVector_->NumVectors()!=NumVectors) { delete ImportVector_; ImportVector_= 0;}
00261       }
00262       if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),NumVectors); // Create import vector if needed
00263       ImportVector_->Import(X, *Importer(), Insert);
00264       Xp = (double**)ImportVector_->Pointers();
00265     }
00266 
00267     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
00268     if (Exporter()!=0) {
00269       if (ExportVector_!=0) {
00270   if (ExportVector_->NumVectors()!=NumVectors) { delete ExportVector_; ExportVector_= 0;}
00271       }
00272       if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),NumVectors); // Create Export vector if needed
00273       Yp = (double**)ExportVector_->Pointers();
00274     }
00275 
00276     // Do actual computation
00277 
00278     for (i=0; i < NumMyRows_; i++) {
00279       int      NumEntries = *NumEntriesPerRow++;
00280       int *    RowIndices = *Indices++;
00281       double * RowValues  = *Values++;
00282       for (k=0; k<NumVectors; k++) {
00283   double sum = 0.0;
00284   for (j=0; j < NumEntries; j++) sum += RowValues[j] * Xp[k][RowIndices[j]];
00285   Yp[k][i] = sum;
00286       }
00287     }
00288     if (Exporter()!=0) Y.Export(*ExportVector_, *Exporter(), Add); // Fill Y with Values from export vector
00289   }
00290   else { // Transpose operation
00291 
00292 
00293     // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
00294 
00295     if (Exporter()!=0) {
00296       if (ExportVector_!=0) {
00297   if (ExportVector_->NumVectors()!=NumVectors) { delete ExportVector_; ExportVector_= 0;}
00298       }
00299       if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),NumVectors); // Create Export vector if needed
00300       ExportVector_->Import(X, *Exporter(), Insert);
00301       Xp = (double**)ExportVector_->Pointers();
00302     }
00303 
00304     // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
00305     if (Importer()!=0) {
00306       if (ImportVector_!=0) {
00307   if (ImportVector_->NumVectors()!=NumVectors) { delete ImportVector_; ImportVector_= 0;}
00308       }
00309       if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),NumVectors); // Create import vector if needed
00310       Yp = (double**)ImportVector_->Pointers();
00311     }
00312 
00313     // Do actual computation
00314 
00315 
00316 
00317     for (k=0; k<NumVectors; k++) 
00318       for (i=0; i < NumMyCols_; i++) Yp[k][i] = 0.0; // Initialize y for transpose multiply
00319     
00320     for (i=0; i < NumMyRows_; i++) {
00321       int      NumEntries = *NumEntriesPerRow++;
00322       int *    RowIndices = *Indices++;
00323       double * RowValues  = *Values++;
00324       for (k=0; k<NumVectors; k++) {
00325   for (j=0; j < NumEntries; j++) Yp[k][RowIndices[j]] += RowValues[j] * Xp[k][i];
00326       }
00327     }
00328     if (Importer()!=0) Y.Export(*ImportVector_, *Importer(), Add); // Fill Y with Values from export vector
00329   }
00330 
00331   UpdateFlops(2*NumVectors*NumGlobalNonzeros());
00332   return(0);
00333 }
00334 
00335 //=============================================================================
00336 int Epetra_FastCrsMatrix::Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_Vector& x, Epetra_Vector& y) const {
00337   //
00338   // This function find y such that Ly = x or Uy = x or the transpose cases.
00339   //
00340 
00341   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
00342 
00343   if ((Upper) && (!UpperTriangular())) EPETRA_CHK_ERR(-2);
00344   if ((!Upper) && (!LowerTriangular())) EPETRA_CHK_ERR(-3);
00345   if ((!UnitDiagonal) && (NoDiagonal())) EPETRA_CHK_ERR(-4); // If matrix has no diagonal, we must use UnitDiagonal
00346   if ((!UnitDiagonal) && (NumMyDiagonals()<NumMyRows_)) EPETRA_CHK_ERR(-5); // Need each row to have a diagonal
00347       
00348 
00349   int i, j, j0;
00350   int * NumEntriesPerRow = NumEntriesPerRow_;
00351   int ** Indices = Indices_;
00352   double ** Values = Values_;
00353   int NumMyCols_ = NumMyCols();
00354 
00355   // If upper, point to last row
00356   if ((Upper && !Trans) || (!Upper && Trans)) {
00357     NumEntriesPerRow += NumMyRows_-1;
00358     Indices += NumMyRows_-1;
00359     Values += NumMyRows_-1;
00360   }
00361     
00362   double *xp = (double*)x.Values();
00363   double *yp = (double*)y.Values();
00364 
00365   if (!Trans) {
00366 
00367     if (Upper) {
00368 
00369       j0 = 1;
00370       if (NoDiagonal()) j0--; // Include first term if no diagonal
00371       for (i=NumMyRows_-1; i >=0; i--) {
00372   int      NumEntries = *NumEntriesPerRow--;
00373   int *    RowIndices = *Indices--;
00374   double * RowValues  = *Values--;
00375   double sum = 0.0;
00376   for (j=j0; j < NumEntries; j++) sum += RowValues[j] * yp[RowIndices[j]];
00377   
00378   if (UnitDiagonal) yp[i] = xp[i] - sum;
00379   else yp[i] = (xp[i] - sum)/RowValues[0];
00380 
00381       }
00382     }
00383     else {
00384       j0 = 1;
00385       if (NoDiagonal()) j0--; // Include first term if no diagonal
00386       for (i=0; i < NumMyRows_; i++) {
00387   int      NumEntries = *NumEntriesPerRow++ - j0;
00388   int *    RowIndices = *Indices++;
00389   double * RowValues  = *Values++;
00390   double sum = 0.0;
00391   for (j=0; j < NumEntries; j++) sum += RowValues[j] * yp[RowIndices[j]];
00392   
00393   if (UnitDiagonal) yp[i] = xp[i] - sum;
00394   else yp[i] = (xp[i] - sum)/RowValues[NumEntries];
00395 
00396       }
00397     }
00398   }
00399 
00400   // ***********  Transpose case *******************************
00401 
00402   else {
00403 
00404     if (xp!=yp) for (i=0; i < NumMyCols_; i++) yp[i] = xp[i]; // Initialize y for transpose solve
00405     
00406     if (Upper) {
00407 
00408       j0 = 1;
00409       if (NoDiagonal()) j0--; // Include first term if no diagonal
00410     
00411       for (i=0; i < NumMyRows_; i++) {
00412   int      NumEntries = *NumEntriesPerRow++;
00413   int *    RowIndices = *Indices++;
00414   double * RowValues  = *Values++;
00415   if (!UnitDiagonal) yp[i] = yp[i]/RowValues[0];
00416   for (j=j0; j < NumEntries; j++) yp[RowIndices[j]] -= RowValues[j] * yp[i];
00417       }
00418     }
00419     else {
00420 
00421       j0 = 1;
00422       if (NoDiagonal()) j0--; // Include first term if no diagonal
00423     
00424       for (i=NumMyRows_-1; i >= 0; i--) {
00425   int      NumEntries = *NumEntriesPerRow-- - j0;
00426   int *    RowIndices = *Indices--;
00427   double * RowValues  = *Values--;
00428   if (!UnitDiagonal)  yp[i] = yp[i]/RowValues[NumEntries];
00429   for (j=0; j < NumEntries; j++) yp[RowIndices[j]] -= RowValues[j] * yp[i];
00430       }
00431     }
00432 
00433   }
00434   UpdateFlops(2*NumGlobalNonzeros());
00435   return(0);
00436 }
00437 
00438 //=============================================================================
00439 int Epetra_FastCrsMatrix::Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
00440   //
00441   // This function find Y such that LY = X or UY = X or the transpose cases.
00442   //
00443   if (X.NumVectors()==1 && Y.NumVectors()==1) {
00444     double * xp = (double *) X[0];
00445     double * yp = (double *) Y[0];
00446     Epetra_Vector x(View, X.Map(), xp);
00447     Epetra_Vector y(View, Y.Map(), yp);
00448     return(Solve(Upper, Trans, UnitDiagonal, x, y));
00449   }
00450   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
00451 
00452   if ((Upper) && (!UpperTriangular())) EPETRA_CHK_ERR(-2);
00453   if ((!Upper) && (!LowerTriangular())) EPETRA_CHK_ERR(-3);
00454   if ((!UnitDiagonal) && (NoDiagonal())) EPETRA_CHK_ERR(-4); // If matrix has no diagonal, we must use UnitDiagonal
00455   if ((!UnitDiagonal) && (NumMyDiagonals()<NumMyRows_)) EPETRA_CHK_ERR(-5); // Need each row to have a diagonal
00456 
00457   int i, j, j0, k;
00458   int * NumEntriesPerRow = NumEntriesPerRow_;
00459   int ** Indices = Indices_;
00460   double ** Values = Values_;
00461   double diag;
00462 
00463   // If upper, point to last row
00464   if ((Upper && !Trans) || (!Upper && Trans)) {
00465     NumEntriesPerRow += NumMyRows_-1;
00466     Indices += NumMyRows_-1;
00467     Values += NumMyRows_-1;
00468   }
00469 
00470   double **Xp = (double**)X.Pointers();
00471   double **Yp = (double**)Y.Pointers();
00472 
00473   int NumVectors = X.NumVectors();
00474 
00475   if (!Trans) {
00476     
00477 
00478     if (Upper) {
00479       
00480       j0 = 1;
00481       if (NoDiagonal()) j0--; // Include first term if no diagonal
00482       for (i=NumMyRows_-1; i >=0; i--) {
00483   int      NumEntries = *NumEntriesPerRow--;
00484   int *    RowIndices = *Indices--;
00485   double * RowValues  = *Values--;
00486   if (!UnitDiagonal) diag = 1.0/RowValues[0]; // Take inverse of diagonal once for later use
00487   for (k=0; k<NumVectors; k++) {
00488     double sum = 0.0;
00489     for (j=j0; j < NumEntries; j++) sum += RowValues[j] * Yp[k][RowIndices[j]];
00490     
00491     if (UnitDiagonal) Yp[k][i] = Xp[k][i] - sum;
00492     else Yp[k][i] = (Xp[k][i] - sum)*diag;
00493   }
00494       }
00495     }
00496     else {
00497       j0 = 1;
00498       if (NoDiagonal()) j0--; // Include first term if no diagonal
00499       for (i=0; i < NumMyRows_; i++) {
00500   int      NumEntries = *NumEntriesPerRow++ - j0;
00501   int *    RowIndices = *Indices++;
00502   double * RowValues  = *Values++;
00503   if (!UnitDiagonal) diag = 1.0/RowValues[NumEntries]; // Take inverse of diagonal once for later use
00504   for (k=0; k<NumVectors; k++) {
00505     double sum = 0.0;
00506     for (j=0; j < NumEntries; j++) sum += RowValues[j] * Yp[k][RowIndices[j]];
00507     
00508     if (UnitDiagonal) Yp[k][i] = Xp[k][i] - sum;
00509     else Yp[k][i] = (Xp[k][i] - sum)*diag;
00510   }
00511       }
00512     }
00513   }
00514   // ***********  Transpose case *******************************
00515 
00516   else {
00517 
00518     for (k=0; k<NumVectors; k++) 
00519       if (Yp[k]!=Xp[k]) for (i=0; i < NumMyRows_; i++) Yp[k][i] = Xp[k][i]; // Initialize y for transpose multiply
00520     
00521     if (Upper) {
00522       
00523       j0 = 1;
00524       if (NoDiagonal()) j0--; // Include first term if no diagonal
00525       
00526       for (i=0; i < NumMyRows_; i++) {
00527   int      NumEntries = *NumEntriesPerRow++;
00528   int *    RowIndices = *Indices++;
00529   double * RowValues  = *Values++;
00530   if (!UnitDiagonal) diag = 1.0/RowValues[j0]; // Take inverse of diagonal once for later use
00531   for (k=0; k<NumVectors; k++) {
00532     if (!UnitDiagonal) Yp[k][i] = Yp[k][i]*diag;
00533     for (j=j0; j < NumEntries; j++) Yp[k][RowIndices[j]] -= RowValues[j] * Yp[k][i];
00534   }
00535       }
00536     }
00537     else {
00538       
00539       j0 = 1;
00540       if (NoDiagonal()) j0--; // Include first term if no diagonal
00541       
00542       for (i=NumMyRows_-1; i>=0; i--) {
00543   int      NumEntries = *NumEntriesPerRow-- - j0;
00544   int *    RowIndices = *Indices--;
00545   double * RowValues  = *Values--;
00546   for (k=0; k<NumVectors; k++) {
00547     if (!UnitDiagonal)  Yp[k][i] = Yp[k][i]/Xp[k][i];
00548     for (j=0; j < NumEntries; j++) Yp[k][RowIndices[j]] -= RowValues[j] * Yp[k][i];
00549         }
00550       }
00551     }
00552   }
00553   
00554   UpdateFlops(2*NumVectors*NumGlobalNonzeros());
00555   return(0);
00556 }

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