Epetra Package Browser (Single Doxygen Collection) Development
Epetra_FastCrsMatrix.cpp
Go to the documentation of this file.
00001 /*
00002 //@HEADER
00003 // ************************************************************************
00004 // 
00005 //               Epetra: Linear Algebra Services Package 
00006 //                 Copyright 2011 Sandia Corporation
00007 // 
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00039 // 
00040 // ************************************************************************
00041 //@HEADER
00042 */
00043 
00044 #include "Epetra_FastCrsMatrix.h"
00045 #include "Epetra_Map.h"
00046 #include "Epetra_Import.h"
00047 #include "Epetra_Export.h"
00048 #include "Epetra_Vector.h"
00049 #include "Epetra_MultiVector.h"
00050 #include "Epetra_Comm.h"
00051 #include "Epetra_Distributor.h"
00052 #include <limits.h>
00053 
00054 //==============================================================================
00055 Epetra_FastCrsMatrix::Epetra_FastCrsMatrix(const Epetra_CrsMatrix & Matrix, bool UseFloats) 
00056   : CrsMatrix_(Matrix),
00057     Values_(0),
00058     NumMyRows_(Matrix.NumMyRows()),
00059     NumMyNonzeros(Matrix.NumMyNonzeros()),
00060     ImportVector_(0),
00061     ExportVector_(0),
00062     CV_(Copy)
00063 {
00064   if (!CrsMatrix_.Filled()) throw CrsMatrix_.ReportError("Input matrix must have called FillComplete()", -1);
00065   Allocate(UseFloats);
00066 }
00067 
00068 //==============================================================================
00069 Epetra_FastCrsMatrix::~Epetra_FastCrsMatrix(){
00070 
00071   if (Values_!=0 && ValuesAllocated_) delete [] Values_;
00072   if (FloatValues_!=0) delete [] FloatValues_;
00073   if (Indices_!=0) delete [] Indices_;
00074   if (ShortIndices_!=0) delete [] ShortIndices_;
00075 
00076   if (ImportVector_!=0) delete ImportVector_;
00077   ImportVector_=0;
00078   if (ExportVector_!=0) delete ExportVector_;
00079   ExportVector_=0;
00080 }
00081 
00082 //==============================================================================
00083 int Epetra_FastCrsMatrix::UpdateValues(const Epetra_CrsMatrix & Matrix) {
00084 }
00085 
00086 //==============================================================================
00087 int Epetra_FastCrsMatrix::Allocate(bool UseFloats) {
00088 
00089   int i, j;
00090   
00091   // First determine how to handle values.  Three possibilities:
00092   // 1) Values are contiguously stored in user matrix, UseFloats is false so we point to the values in
00093   //    the user matrix.
00094   // 2) Values are not contiguously stored, UseFloats is false, so we copy values into a contiguous double array.
00095   // 3) UseFloats is true so we create a single precision copy of matrix values.
00096 
00097   if (CrsMatrix.IndicesAreContiguous() && !UseFloats) {
00098     ValuesAllocated_ = false;
00099     Values_ = CrsMatrix_[0]; // First value in the user's matrix
00100   }
00101   else if (!UseFloats) {
00102     // Allocate Values array
00103     Values_ = new double[NumMyNonzeros_];
00104     double * ptr = Values_;
00105     for (i=0; i<NumMyRows_; i++) {
00106       int NumEntries = CrsMatrix_.NumMyEntries(i);
00107       double * rowi = CrsMatrix_[i];
00108       for (j=0; j< NumEntries; j++) *ptr++ = *rowi++; // Copy values to contiguous storage
00109     }
00110   }
00111   else { // UseFloats == true
00112     FloatValues_ = new float[NumMyNonzeros_];
00113     float * ptr = FloatValues_;
00114     for (i=0; i<NumMyRows_; i++) {
00115       int NumEntries = CrsMatrix_.NumMyEntries(i);
00116       double * rowi = CrsMatrix_[i];
00117       for (j=0; j< NumEntries; j++) *ptr++ = (float) *rowi++; // convert and copy values
00118     }
00119   }
00120 
00121   // Next determine how to handle integers.  Two possibilities:
00122   // 1) Local column range is within the range of unsigned short ints, so we copy the indices to an array of unsigned shorts.
00123   // 2) Local column range is outside range of unsigned shorts and we copy to array of ints.
00124   // In both cases we embed the nonzero count per row into the index array.
00125 
00126   if (CrsMatrix_.NumMyCols()>USHRT_MAX) {
00127     // Allocate Values array
00128     Indices_ = new int[NumMyNonzeros_+NumMyRows_];
00129     int * ptr = Indices_;
00130     for (i=0; i<NumMyRows_; i++) {
00131       int NumEntries = CrsMatrix_.NumMyEntries(i);
00132       int * rowi = CrsMatrix_.Graph()[i];
00133       *ptr++ = NumEntries; // Pack this value first
00134       for (j=0; j< NumEntries; j++) *ptr++ = *rowi++; // Copy values to contiguous storage
00135     }
00136   }
00137   else { // CrsMatrix_.NumMyCols()<=USHRT_MAX
00138     ShortIndices_ = new unsigned short[NumMyNonzeros_+NumMyRows_];
00139     unsigned short * ptr = ShortIndices_;
00140     for (i=0; i<NumMyRows_; i++) {
00141       int NumEntries = CrsMatrix_.NumMyEntries(i);
00142       unsigned short * rowi = CrsMatrix_Graph()[i];
00143       *ptr++ = NumEntries; // Pack this value first
00144       for (j=0; j< NumEntries; j++) *ptr++ = (unsigned short) *rowi++; // convert and copy indices
00145     }
00146   }
00147 
00148   return(0);
00149 }
00150 //=============================================================================
00151 int Epetra_FastCrsMatrix::Multiply(bool TransA, const Epetra_Vector& x, Epetra_Vector& y) const {
00152   //
00153   // This function forms the product y = A * x or y = A' * x
00154   //
00155 
00156   int i, j;
00157   double * xp = (double*)x.Values();
00158   double *yp = (double*)y.Values();
00159   int NumMyCols_ = NumMyCols();
00160 
00161 
00162   if (!TransA) {
00163 
00164     // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
00165     if (Importer()!=0) {
00166       if (ImportVector_!=0) {
00167   if (ImportVector_->NumVectors()!=1) { delete ImportVector_; ImportVector_= 0;}
00168       }
00169       if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),1); // Create import vector if needed
00170       ImportVector_->Import(x, *Importer(), Insert);
00171       xp = (double*)ImportVector_->Values();
00172     }
00173 
00174     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
00175     if (Exporter()!=0) {
00176       if (ExportVector_!=0) {
00177   if (ExportVector_->NumVectors()!=1) { delete ExportVector_; ExportVector_= 0;}
00178       }
00179       if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),1); // Create Export vector if needed
00180       yp = (double*)ExportVector_->Values();
00181     }
00182 
00183     // Do actual computation
00184 
00185     for (i=0; i < NumMyRows_; i++) {
00186       int      NumEntries = *NumEntriesPerRow++;
00187       int *    RowIndices = *Indices++;
00188       double * RowValues  = *Values++;
00189       double sum = 0.0;
00190       for (j=0; j < NumEntries; j++) sum += RowValues[j] * xp[RowIndices[j]];
00191 
00192       yp[i] = sum;
00193 
00194     }
00195     if (Exporter()!=0) y.Export(*ExportVector_, *Exporter(), Add); // Fill y with Values from export vector
00196   }
00197 
00198   else { // Transpose operation
00199 
00200 
00201     // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
00202 
00203     if (Exporter()!=0) {
00204       if (ExportVector_!=0) {
00205   if (ExportVector_->NumVectors()!=1) { delete ExportVector_; ExportVector_= 0;}
00206       }
00207       if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),1); // Create Export vector if needed
00208       ExportVector_->Import(x, *Exporter(), Insert);
00209       xp = (double*)ExportVector_->Values();
00210     }
00211 
00212     // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
00213     if (Importer()!=0) {
00214       if (ImportVector_!=0) {
00215   if (ImportVector_->NumVectors()!=1) { delete ImportVector_; ImportVector_= 0;}
00216       }
00217       if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),1); // Create import vector if needed
00218       yp = (double*)ImportVector_->Values();
00219     }
00220 
00221     // Do actual computation
00222 
00223     for (i=0; i < NumMyCols_; i++) yp[i] = 0.0; // Initialize y for transpose multiply
00224         
00225     for (i=0; i < NumMyRows_; i++) {
00226       int      NumEntries = *NumEntriesPerRow++;
00227       int *    RowIndices = *Indices++;
00228       double * RowValues  = *Values++;
00229       for (j=0; j < NumEntries; j++) yp[RowIndices[j]] += RowValues[j] * xp[i];
00230     }
00231     if (Importer()!=0) y.Export(*ImportVector_, *Importer(), Add); // Fill y with Values from export vector
00232   }
00233 
00234   UpdateFlops(2*NumGlobalNonzeros64());
00235   return(0);
00236 }
00237 
00238 //=============================================================================
00239 int Epetra_FastCrsMatrix::Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
00240   //
00241   // This function forms the product Y = A * Y or Y = A' * X
00242   //
00243   if (X.NumVectors()==1 && Y.NumVectors()==1) {
00244     double * xp = (double *) X[0];
00245     double * yp = (double *) Y[0];
00246     Epetra_Vector x(View, X.Map(), xp);
00247     Epetra_Vector y(View, Y.Map(), yp);
00248     return(Multiply(TransA, x, y));
00249   }
00250   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
00251 
00252   int i, j, k;
00253   int * NumEntriesPerRow = NumEntriesPerRow_;
00254   int ** Indices = Indices_;
00255   double ** Values = Values_;
00256 
00257   double **Xp = (double**)X.Pointers();
00258   double **Yp = (double**)Y.Pointers();
00259 
00260   int NumVectors = X.NumVectors();
00261   int NumMyCols_ = NumMyCols();
00262 
00263 
00264   // Need to better manage the Import and Export vectors:
00265   // - Need accessor functions
00266   // - Need to make the NumVector match (use a View to do this)
00267   // - Need to look at RightScale and ColSum routines too.
00268 
00269   if (!TransA) {
00270 
00271     // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
00272     if (Importer()!=0) {
00273       if (ImportVector_!=0) {
00274   if (ImportVector_->NumVectors()!=NumVectors) { delete ImportVector_; ImportVector_= 0;}
00275       }
00276       if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),NumVectors); // Create import vector if needed
00277       ImportVector_->Import(X, *Importer(), Insert);
00278       Xp = (double**)ImportVector_->Pointers();
00279     }
00280 
00281     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
00282     if (Exporter()!=0) {
00283       if (ExportVector_!=0) {
00284   if (ExportVector_->NumVectors()!=NumVectors) { delete ExportVector_; ExportVector_= 0;}
00285       }
00286       if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),NumVectors); // Create Export vector if needed
00287       Yp = (double**)ExportVector_->Pointers();
00288     }
00289 
00290     // Do actual computation
00291 
00292     for (i=0; i < NumMyRows_; i++) {
00293       int      NumEntries = *NumEntriesPerRow++;
00294       int *    RowIndices = *Indices++;
00295       double * RowValues  = *Values++;
00296       for (k=0; k<NumVectors; k++) {
00297   double sum = 0.0;
00298   for (j=0; j < NumEntries; j++) sum += RowValues[j] * Xp[k][RowIndices[j]];
00299   Yp[k][i] = sum;
00300       }
00301     }
00302     if (Exporter()!=0) Y.Export(*ExportVector_, *Exporter(), Add); // Fill Y with Values from export vector
00303   }
00304   else { // Transpose operation
00305 
00306 
00307     // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
00308 
00309     if (Exporter()!=0) {
00310       if (ExportVector_!=0) {
00311   if (ExportVector_->NumVectors()!=NumVectors) { delete ExportVector_; ExportVector_= 0;}
00312       }
00313       if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),NumVectors); // Create Export vector if needed
00314       ExportVector_->Import(X, *Exporter(), Insert);
00315       Xp = (double**)ExportVector_->Pointers();
00316     }
00317 
00318     // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
00319     if (Importer()!=0) {
00320       if (ImportVector_!=0) {
00321   if (ImportVector_->NumVectors()!=NumVectors) { delete ImportVector_; ImportVector_= 0;}
00322       }
00323       if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),NumVectors); // Create import vector if needed
00324       Yp = (double**)ImportVector_->Pointers();
00325     }
00326 
00327     // Do actual computation
00328 
00329 
00330 
00331     for (k=0; k<NumVectors; k++) 
00332       for (i=0; i < NumMyCols_; i++) Yp[k][i] = 0.0; // Initialize y for transpose multiply
00333     
00334     for (i=0; i < NumMyRows_; i++) {
00335       int      NumEntries = *NumEntriesPerRow++;
00336       int *    RowIndices = *Indices++;
00337       double * RowValues  = *Values++;
00338       for (k=0; k<NumVectors; k++) {
00339   for (j=0; j < NumEntries; j++) Yp[k][RowIndices[j]] += RowValues[j] * Xp[k][i];
00340       }
00341     }
00342     if (Importer()!=0) Y.Export(*ImportVector_, *Importer(), Add); // Fill Y with Values from export vector
00343   }
00344 
00345   UpdateFlops(2*NumVectors*NumGlobalNonzeros64());
00346   return(0);
00347 }
00348 
00349 //=============================================================================
00350 int Epetra_FastCrsMatrix::Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_Vector& x, Epetra_Vector& y) const {
00351   //
00352   // This function find y such that Ly = x or Uy = x or the transpose cases.
00353   //
00354 
00355   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
00356 
00357   if ((Upper) && (!UpperTriangular())) EPETRA_CHK_ERR(-2);
00358   if ((!Upper) && (!LowerTriangular())) EPETRA_CHK_ERR(-3);
00359   if ((!UnitDiagonal) && (NoDiagonal())) EPETRA_CHK_ERR(-4); // If matrix has no diagonal, we must use UnitDiagonal
00360   if ((!UnitDiagonal) && (NumMyDiagonals()<NumMyRows_)) EPETRA_CHK_ERR(-5); // Need each row to have a diagonal
00361       
00362 
00363   int i, j, j0;
00364   int * NumEntriesPerRow = NumEntriesPerRow_;
00365   int ** Indices = Indices_;
00366   double ** Values = Values_;
00367   int NumMyCols_ = NumMyCols();
00368 
00369   // If upper, point to last row
00370   if ((Upper && !Trans) || (!Upper && Trans)) {
00371     NumEntriesPerRow += NumMyRows_-1;
00372     Indices += NumMyRows_-1;
00373     Values += NumMyRows_-1;
00374   }
00375     
00376   double *xp = (double*)x.Values();
00377   double *yp = (double*)y.Values();
00378 
00379   if (!Trans) {
00380 
00381     if (Upper) {
00382 
00383       j0 = 1;
00384       if (NoDiagonal()) j0--; // Include first term if no diagonal
00385       for (i=NumMyRows_-1; i >=0; i--) {
00386   int      NumEntries = *NumEntriesPerRow--;
00387   int *    RowIndices = *Indices--;
00388   double * RowValues  = *Values--;
00389   double sum = 0.0;
00390   for (j=j0; j < NumEntries; j++) sum += RowValues[j] * yp[RowIndices[j]];
00391   
00392   if (UnitDiagonal) yp[i] = xp[i] - sum;
00393   else yp[i] = (xp[i] - sum)/RowValues[0];
00394 
00395       }
00396     }
00397     else {
00398       j0 = 1;
00399       if (NoDiagonal()) j0--; // Include first term if no diagonal
00400       for (i=0; i < NumMyRows_; i++) {
00401   int      NumEntries = *NumEntriesPerRow++ - j0;
00402   int *    RowIndices = *Indices++;
00403   double * RowValues  = *Values++;
00404   double sum = 0.0;
00405   for (j=0; j < NumEntries; j++) sum += RowValues[j] * yp[RowIndices[j]];
00406   
00407   if (UnitDiagonal) yp[i] = xp[i] - sum;
00408   else yp[i] = (xp[i] - sum)/RowValues[NumEntries];
00409 
00410       }
00411     }
00412   }
00413 
00414   // ***********  Transpose case *******************************
00415 
00416   else {
00417 
00418     if (xp!=yp) for (i=0; i < NumMyCols_; i++) yp[i] = xp[i]; // Initialize y for transpose solve
00419     
00420     if (Upper) {
00421 
00422       j0 = 1;
00423       if (NoDiagonal()) j0--; // Include first term if no diagonal
00424     
00425       for (i=0; i < NumMyRows_; i++) {
00426   int      NumEntries = *NumEntriesPerRow++;
00427   int *    RowIndices = *Indices++;
00428   double * RowValues  = *Values++;
00429   if (!UnitDiagonal) yp[i] = yp[i]/RowValues[0];
00430   for (j=j0; j < NumEntries; j++) yp[RowIndices[j]] -= RowValues[j] * yp[i];
00431       }
00432     }
00433     else {
00434 
00435       j0 = 1;
00436       if (NoDiagonal()) j0--; // Include first term if no diagonal
00437     
00438       for (i=NumMyRows_-1; i >= 0; i--) {
00439   int      NumEntries = *NumEntriesPerRow-- - j0;
00440   int *    RowIndices = *Indices--;
00441   double * RowValues  = *Values--;
00442   if (!UnitDiagonal)  yp[i] = yp[i]/RowValues[NumEntries];
00443   for (j=0; j < NumEntries; j++) yp[RowIndices[j]] -= RowValues[j] * yp[i];
00444       }
00445     }
00446 
00447   }
00448   UpdateFlops(2*NumGlobalNonzeros64());
00449   return(0);
00450 }
00451 
00452 //=============================================================================
00453 int Epetra_FastCrsMatrix::Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
00454   //
00455   // This function find Y such that LY = X or UY = X or the transpose cases.
00456   //
00457   if (X.NumVectors()==1 && Y.NumVectors()==1) {
00458     double * xp = (double *) X[0];
00459     double * yp = (double *) Y[0];
00460     Epetra_Vector x(View, X.Map(), xp);
00461     Epetra_Vector y(View, Y.Map(), yp);
00462     return(Solve(Upper, Trans, UnitDiagonal, x, y));
00463   }
00464   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
00465 
00466   if ((Upper) && (!UpperTriangular())) EPETRA_CHK_ERR(-2);
00467   if ((!Upper) && (!LowerTriangular())) EPETRA_CHK_ERR(-3);
00468   if ((!UnitDiagonal) && (NoDiagonal())) EPETRA_CHK_ERR(-4); // If matrix has no diagonal, we must use UnitDiagonal
00469   if ((!UnitDiagonal) && (NumMyDiagonals()<NumMyRows_)) EPETRA_CHK_ERR(-5); // Need each row to have a diagonal
00470 
00471   int i, j, j0, k;
00472   int * NumEntriesPerRow = NumEntriesPerRow_;
00473   int ** Indices = Indices_;
00474   double ** Values = Values_;
00475   double diag;
00476 
00477   // If upper, point to last row
00478   if ((Upper && !Trans) || (!Upper && Trans)) {
00479     NumEntriesPerRow += NumMyRows_-1;
00480     Indices += NumMyRows_-1;
00481     Values += NumMyRows_-1;
00482   }
00483 
00484   double **Xp = (double**)X.Pointers();
00485   double **Yp = (double**)Y.Pointers();
00486 
00487   int NumVectors = X.NumVectors();
00488 
00489   if (!Trans) {
00490     
00491 
00492     if (Upper) {
00493       
00494       j0 = 1;
00495       if (NoDiagonal()) j0--; // Include first term if no diagonal
00496       for (i=NumMyRows_-1; i >=0; i--) {
00497   int      NumEntries = *NumEntriesPerRow--;
00498   int *    RowIndices = *Indices--;
00499   double * RowValues  = *Values--;
00500   if (!UnitDiagonal) diag = 1.0/RowValues[0]; // Take inverse of diagonal once for later use
00501   for (k=0; k<NumVectors; k++) {
00502     double sum = 0.0;
00503     for (j=j0; j < NumEntries; j++) sum += RowValues[j] * Yp[k][RowIndices[j]];
00504     
00505     if (UnitDiagonal) Yp[k][i] = Xp[k][i] - sum;
00506     else Yp[k][i] = (Xp[k][i] - sum)*diag;
00507   }
00508       }
00509     }
00510     else {
00511       j0 = 1;
00512       if (NoDiagonal()) j0--; // Include first term if no diagonal
00513       for (i=0; i < NumMyRows_; i++) {
00514   int      NumEntries = *NumEntriesPerRow++ - j0;
00515   int *    RowIndices = *Indices++;
00516   double * RowValues  = *Values++;
00517   if (!UnitDiagonal) diag = 1.0/RowValues[NumEntries]; // Take inverse of diagonal once for later use
00518   for (k=0; k<NumVectors; k++) {
00519     double sum = 0.0;
00520     for (j=0; j < NumEntries; j++) sum += RowValues[j] * Yp[k][RowIndices[j]];
00521     
00522     if (UnitDiagonal) Yp[k][i] = Xp[k][i] - sum;
00523     else Yp[k][i] = (Xp[k][i] - sum)*diag;
00524   }
00525       }
00526     }
00527   }
00528   // ***********  Transpose case *******************************
00529 
00530   else {
00531 
00532     for (k=0; k<NumVectors; k++) 
00533       if (Yp[k]!=Xp[k]) for (i=0; i < NumMyRows_; i++) Yp[k][i] = Xp[k][i]; // Initialize y for transpose multiply
00534     
00535     if (Upper) {
00536       
00537       j0 = 1;
00538       if (NoDiagonal()) j0--; // Include first term if no diagonal
00539       
00540       for (i=0; i < NumMyRows_; i++) {
00541   int      NumEntries = *NumEntriesPerRow++;
00542   int *    RowIndices = *Indices++;
00543   double * RowValues  = *Values++;
00544   if (!UnitDiagonal) diag = 1.0/RowValues[j0]; // Take inverse of diagonal once for later use
00545   for (k=0; k<NumVectors; k++) {
00546     if (!UnitDiagonal) Yp[k][i] = Yp[k][i]*diag;
00547     for (j=j0; j < NumEntries; j++) Yp[k][RowIndices[j]] -= RowValues[j] * Yp[k][i];
00548   }
00549       }
00550     }
00551     else {
00552       
00553       j0 = 1;
00554       if (NoDiagonal()) j0--; // Include first term if no diagonal
00555       
00556       for (i=NumMyRows_-1; i>=0; i--) {
00557   int      NumEntries = *NumEntriesPerRow-- - j0;
00558   int *    RowIndices = *Indices--;
00559   double * RowValues  = *Values--;
00560   for (k=0; k<NumVectors; k++) {
00561     if (!UnitDiagonal)  Yp[k][i] = Yp[k][i]/Xp[k][i];
00562     for (j=0; j < NumEntries; j++) Yp[k][RowIndices[j]] -= RowValues[j] * Yp[k][i];
00563         }
00564       }
00565     }
00566   }
00567   
00568   UpdateFlops(2*NumVectors*NumGlobalNonzeros64());
00569   return(0);
00570 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines