Epetra_SerialDenseMatrix.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_SerialDenseMatrix.h"
00031 #include "Epetra_SerialSymDenseMatrix.h"
00032 #include "Epetra_SerialSymDenseMatrix.h"
00033 #include "Epetra_Util.h"
00034 
00035 //=============================================================================
00036 Epetra_SerialDenseMatrix::Epetra_SerialDenseMatrix(bool set_object_label)
00037   : Epetra_CompObject(),
00038     Epetra_Object(-1, false),
00039     M_(0),
00040     N_(0),
00041     A_Copied_(false),
00042     CV_(Copy),
00043     LDA_(0),
00044     A_(0),
00045     UseTranspose_(false)
00046 {
00047   if (set_object_label) {
00048     SetLabel("Epetra::SerialDenseMatrix");
00049   }
00050 }
00051 
00052 //=============================================================================
00053 Epetra_SerialDenseMatrix::Epetra_SerialDenseMatrix(int NumRows, int NumCols,
00054                                                    bool set_object_label)
00055   : Epetra_CompObject(),
00056     Epetra_Object(-1, false),
00057     M_(0),
00058     N_(0),
00059     A_Copied_(false),
00060     CV_(Copy),
00061     LDA_(0),
00062     A_(0),
00063     UseTranspose_(false)
00064 {
00065   if (set_object_label) {
00066     SetLabel("Epetra::SerialDenseMatrix");
00067   }
00068   if(NumRows < 0)
00069   throw ReportError("NumRows = " + toString(NumRows) + ". Should be >= 0", -1);
00070   if(NumCols < 0)
00071   throw ReportError("NumCols = " + toString(NumCols) + ". Should be >= 0", -1);
00072 
00073   int errorcode = Shape(NumRows, NumCols);
00074   if(errorcode != 0)
00075     throw ReportError("Shape returned non-zero value", errorcode);
00076 }
00077 
00078 //=============================================================================
00079 Epetra_SerialDenseMatrix::Epetra_SerialDenseMatrix(Epetra_DataAccess CV_in, double* A_in, int LDA_in,
00080                                                    int NumRows, int NumCols,
00081                                                    bool set_object_label)
00082   : Epetra_CompObject(),
00083     Epetra_Object(-1, false),
00084     M_(NumRows),
00085     N_(NumCols),
00086     A_Copied_(false),    
00087     CV_(CV_in),
00088     LDA_(LDA_in),
00089     A_(A_in),
00090     UseTranspose_(false)
00091 {
00092   if (set_object_label) {
00093     SetLabel("Epetra::SerialDenseMatrix");
00094   }
00095   if(A_in == 0)
00096   throw ReportError("Null pointer passed as A parameter.", -3);
00097   if(NumRows < 0)
00098     throw ReportError("NumRows = " + toString(NumRows) + ". Should be >= 0", -1);
00099   if(NumCols < 0)
00100   throw ReportError("NumCols = " + toString(NumCols) + ". Should be >= 0", -1);
00101   if(LDA_in < 0)
00102   throw ReportError("LDA = " + toString(LDA_in) + ". Should be >= 0", -1);
00103 
00104   if (CV_in == Copy) {
00105     LDA_ = M_;
00106     const int newsize = LDA_ * N_;
00107     if (newsize > 0) {
00108       A_ = new double[newsize];
00109       CopyMat(A_in, LDA_in, M_, N_, A_, LDA_);
00110       A_Copied_ = true;
00111     }
00112     else {
00113       A_ = 0;
00114     }
00115   }
00116 }
00117 //=============================================================================
00118 Epetra_SerialDenseMatrix::Epetra_SerialDenseMatrix(const Epetra_SerialDenseMatrix& Source)
00119   : Epetra_CompObject(Source),  
00120     M_(Source.M_),
00121     N_(Source.N_),
00122     A_Copied_(false),
00123     CV_(Source.CV_),
00124     LDA_(Source.LDA_),
00125     A_(Source.A_),
00126     UseTranspose_(false)
00127 {
00128   SetLabel(Source.Label());
00129   if(CV_ == Copy) {
00130     LDA_ = M_;
00131     const int newsize = LDA_ * N_;
00132     if(newsize > 0) {
00133       A_ = new double[newsize];
00134       CopyMat(Source.A_, Source.LDA_, M_, N_, A_, LDA_);
00135       A_Copied_ = true;
00136     }
00137     else {
00138       A_ = 0;
00139     }
00140   }
00141 }
00142 
00143 //=============================================================================
00144 int Epetra_SerialDenseMatrix::Reshape(int NumRows, int NumCols) {
00145   if(NumRows < 0 || NumCols < 0)
00146     return(-1);
00147 
00148   double* A_tmp = 0;
00149   const int newsize = NumRows * NumCols;
00150 
00151   if(newsize > 0) {
00152     // Allocate space for new matrix
00153     A_tmp = new double[newsize];
00154     for (int k = 0; k < newsize; k++) 
00155       A_tmp[k] = 0.0; // Zero out values
00156     int M_tmp = EPETRA_MIN(M_, NumRows);
00157     int N_tmp = EPETRA_MIN(N_, NumCols);
00158     if (A_ != 0) 
00159       CopyMat(A_, LDA_, M_tmp, N_tmp, A_tmp, NumRows); // Copy principal submatrix of A to new A
00160   }
00161   CleanupData(); // Get rid of anything that might be already allocated  
00162   M_ = NumRows;
00163   N_ = NumCols;
00164   LDA_ = M_;
00165   if(newsize > 0) {
00166     A_ = A_tmp; // Set pointer to new A
00167     A_Copied_ = true;
00168   }
00169 
00170   return(0);
00171 }
00172 //=============================================================================
00173 int Epetra_SerialDenseMatrix::Shape(int NumRows, int NumCols) {
00174   if(NumRows < 0 || NumCols < 0)
00175     return(-1);
00176 
00177   CleanupData(); // Get rid of anything that might be already allocated
00178   M_ = NumRows;
00179   N_ = NumCols;
00180   LDA_ = M_;
00181   const int newsize = LDA_ * N_;
00182   if(newsize > 0) {
00183     A_ = new double[newsize];
00184     for (int k = 0; k < newsize; k++) 
00185       A_[k] = 0.0; // Zero out values
00186     A_Copied_ = true;
00187   }
00188 
00189   return(0);
00190 }
00191 //=============================================================================
00192 Epetra_SerialDenseMatrix::~Epetra_SerialDenseMatrix()
00193 {
00194   CleanupData();
00195 }
00196 //=============================================================================
00197 void Epetra_SerialDenseMatrix::CleanupData()
00198 {
00199   if (A_Copied_)
00200     delete [] A_;
00201   A_ = 0;
00202   A_Copied_ = false;
00203   M_ = 0;
00204   N_ = 0;
00205   LDA_ = 0;
00206 }
00207 //=============================================================================
00208 Epetra_SerialDenseMatrix& Epetra_SerialDenseMatrix::operator = (const Epetra_SerialDenseMatrix& Source) {
00209   if(this == &Source)
00210     return(*this); // Special case of source same as target
00211   if((CV_ == View) && (Source.CV_ == View) && (A_ == Source.A_))
00212     return(*this); // Special case of both are views to same data.
00213 
00214   if(strcmp(Label(), Source.Label()) != 0)
00215     throw ReportError("operator= type mismatch (lhs = " + string(Label()) + 
00216                       ", rhs = " + string(Source.Label()) + ").", -5);
00217   
00218   if(Source.CV_ == View) {
00219     if(CV_ == Copy) { // C->V only
00220       CleanupData();
00221       CV_ = View;
00222     }
00223     M_ = Source.M_; // C->V and V->V
00224     N_ = Source.N_;
00225     LDA_ = Source.LDA_;
00226     A_ = Source.A_;
00227   }
00228   else {
00229     if(CV_ == View) { // V->C
00230       CV_ = Copy;
00231       M_ = Source.M_;
00232       N_ = Source.N_;
00233       LDA_ = Source.M_;
00234       const int newsize = LDA_ * N_;
00235       if(newsize > 0) {
00236         A_ = new double[newsize];
00237         A_Copied_ = true;
00238       }
00239       else {
00240         A_ = 0;
00241         A_Copied_ = false;
00242       }
00243     }
00244     else { // C->C
00245       if((Source.M_ <= LDA_) && (Source.N_ == N_)) { // we don't need to reallocate
00246         M_ = Source.M_;
00247         N_ = Source.N_;
00248       }
00249       else { // we need to allocate more space (or less space)
00250         CleanupData();
00251         M_ = Source.M_;
00252         N_ = Source.N_;
00253         LDA_ = Source.M_;
00254         const int newsize = LDA_ * N_;
00255         if(newsize > 0) {
00256           A_ = new double[newsize];
00257           A_Copied_ = true;
00258         }
00259       }
00260     }
00261     CopyMat(Source.A_, Source.LDA_, M_, N_, A_, LDA_); // V->C and C->C
00262   }
00263   
00264   return(*this);
00265 }
00266 
00267 
00268 //=============================================================================
00269 bool Epetra_SerialDenseMatrix::operator==(const Epetra_SerialDenseMatrix& rhs) const
00270 {
00271   if (M_ != rhs.M_ || N_ != rhs.N_) return(false);
00272 
00273   const double* A = A_;
00274   const double* rhsA = rhs.A_;
00275 
00276   for(int j=0; j<N_; ++j) {
00277     int offset = j*LDA_;
00278     int rhsOffset = j*rhs.LDA_;
00279     for(int i=0; i<M_; ++i) {
00280       if (std::abs(A[offset+i] - rhsA[rhsOffset+i]) > Epetra_MinDouble) {
00281   return(false);
00282       }
00283     }
00284   }
00285 
00286   return(true);
00287 }
00288 
00289 //=============================================================================
00290 Epetra_SerialDenseMatrix& Epetra_SerialDenseMatrix::operator+= ( const Epetra_SerialDenseMatrix & Source) {
00291   if (M() != Source.M()) 
00292     throw ReportError("Row dimension of source = " + toString(Source.M()) +
00293                       " is different than  row dimension of target = " + toString(LDA()), -1);
00294   if (N() != Source.N()) 
00295     throw ReportError("Column dimension of source = " + toString(Source.N()) +
00296                       " is different than column dimension of target = " + toString(N()), -2);
00297 
00298   CopyMat(Source.A(), Source.LDA(), Source.M(), Source.N(), A(), LDA(), true);
00299   return(*this);
00300 }
00301 //=============================================================================
00302 void Epetra_SerialDenseMatrix::CopyMat(double* Source,
00303                                        int Source_LDA,
00304                                        int NumRows,
00305                                        int NumCols, 
00306                                        double* Target,
00307                                        int Target_LDA,
00308                                        bool add)
00309 {
00310   int i;
00311   double* tptr = Target;
00312   double* sptr = Source;
00313   if (add) {
00314     for(int j=0; j<NumCols; ++j) {
00315       for(i=0; i<NumRows; ++i) {
00316         tptr[i] += sptr[i];
00317       }
00318 
00319       tptr += Target_LDA;
00320       sptr += Source_LDA;
00321     }
00322   }
00323   else {
00324     for(int j=0; j<NumCols; ++j) {
00325       for(i=0; i<NumRows; ++i) {
00326         tptr[i] = sptr[i];
00327       }
00328 
00329       tptr += Target_LDA;
00330       sptr += Source_LDA;
00331     }
00332   }
00333 /*
00334   double* targetPtr = Target;
00335   double* sourcePtr = Source;
00336   double* targetEnd = 0;
00337  
00338   for (j = 0; j < NumCols; j++) {
00339     targetPtr = Target + j * Target_LDA;
00340     targetEnd = targetPtr+NumRows;
00341     sourcePtr = Source + j * Source_LDA;
00342     if (add) {
00343       while(targetPtr != targetEnd) 
00344     *targetPtr++ += *sourcePtr++;
00345     }
00346     else {
00347       while(targetPtr != targetEnd) 
00348     *targetPtr++ = *sourcePtr++;
00349     }
00350   }
00351 */
00352   return;
00353 }
00354 //=============================================================================
00355 double Epetra_SerialDenseMatrix::NormOne() const {
00356 
00357   int i, j;
00358 
00359     double anorm = 0.0;
00360     double * ptr;
00361     for (j=0; j<N_; j++) {
00362       double sum=0.0;
00363       ptr = A_ + j*LDA_;
00364       for (i=0; i<M_; i++) sum += std::abs(*ptr++);
00365       anorm = EPETRA_MAX(anorm, sum);
00366     }
00367     UpdateFlops((double)N_*(double)N_);
00368     return(anorm);
00369 }
00370 //=============================================================================
00371 double Epetra_SerialDenseMatrix::NormInf() const {
00372 
00373   int i, j;
00374 
00375     double anorm = 0.0;
00376     double * ptr;
00377 
00378     // Loop across columns in inner loop.  Most expensive memory access, but 
00379     // requires no extra storage.
00380     for (i=0; i<M_; i++) {
00381       double sum=0.0;
00382       ptr = A_ + i;
00383       for (j=0; j<N_; j++) {
00384   sum += std::abs(*ptr);
00385   ptr += LDA_;
00386       }
00387       anorm = EPETRA_MAX(anorm, sum);
00388     }
00389     UpdateFlops((double)N_*(double)N_);
00390     return(anorm);
00391 }
00392 //=============================================================================
00393 int Epetra_SerialDenseMatrix::Scale(double ScalarA) {
00394 
00395   int i, j;
00396   
00397   double * ptr;
00398   for (j=0; j<N_; j++) {
00399     ptr = A_ + j*LDA_;
00400     for (i=0; i<M_; i++) { *ptr = ScalarA * (*ptr); ptr++; }
00401   }
00402   UpdateFlops((double)N_*(double)N_);
00403   return(0);
00404   
00405 }
00406 //=========================================================================
00407 int Epetra_SerialDenseMatrix::Multiply (char TransA, char TransB, double ScalarAB, 
00408           const Epetra_SerialDenseMatrix& A, 
00409           const Epetra_SerialDenseMatrix& B,
00410           double ScalarThis ) {
00411   // Check for compatible dimensions
00412 
00413   if (TransA!='T' && TransA!='N') EPETRA_CHK_ERR(-2); // Return error
00414   if (TransB!='T' && TransB!='N') EPETRA_CHK_ERR(-3);
00415   
00416   int A_nrows = (TransA=='T') ? A.N() : A.M();
00417   int A_ncols = (TransA=='T') ? A.M() : A.N();
00418   int B_nrows = (TransB=='T') ? B.N() : B.M();
00419   int B_ncols = (TransB=='T') ? B.M() : B.N();
00420   
00421   if (M_        != A_nrows     ||
00422       A_ncols   != B_nrows     ||
00423       N_        != B_ncols  ) EPETRA_CHK_ERR(-1); // Return error
00424 
00425     
00426   // Call GEMM function
00427   GEMM(TransA, TransB, M_, N_, A_ncols, ScalarAB, A.A(), A.LDA(), 
00428        B.A(), B.LDA(), ScalarThis, A_, LDA_);
00429   long int nflops = 2*M_;
00430   nflops *= N_;
00431   nflops *= A_ncols;
00432   if (ScalarAB != 1.0) nflops += M_*N_;
00433   if (ScalarThis != 0.0) nflops += M_*N_;
00434   UpdateFlops((double)nflops);
00435 
00436   return(0);
00437 }
00438 
00439 //=========================================================================
00440 int  Epetra_SerialDenseMatrix::Multiply (bool transA,
00441                                          const Epetra_SerialDenseMatrix& x,
00442                                          Epetra_SerialDenseMatrix& y)
00443 {
00444   int A_nrows = M();
00445   int x_nrows = x.M();
00446   int y_nrows = y.M();
00447   int A_ncols = N();
00448   int x_ncols = x.N();
00449   int y_ncols = y.N();
00450 
00451   if (transA) {
00452     if (x_nrows != A_nrows) EPETRA_CHK_ERR(-1);
00453     if (y_ncols != x_ncols || y_nrows != A_ncols) y.Reshape(A_ncols, x_ncols);
00454   }
00455   else {
00456     if (x_nrows != A_ncols) EPETRA_CHK_ERR(-1);
00457     if (y_ncols != x_ncols || y_nrows != A_nrows) y.Reshape(A_nrows, x_ncols);
00458   }
00459 
00460   double scalar0 = 0.0;
00461   double scalar1 = 1.0;
00462 
00463   int err = 0;
00464   if (transA) {
00465     err = y.Multiply('T', 'N', scalar1, *this, x, scalar0);
00466   }
00467   else {
00468     err = y.Multiply('N', 'N', scalar1, *this, x, scalar0);
00469   }
00470 
00471   return(0);
00472 }
00473 
00474 //=========================================================================
00475 int  Epetra_SerialDenseMatrix::Multiply (char SideA, double ScalarAB, 
00476               const Epetra_SerialSymDenseMatrix& A, 
00477               const Epetra_SerialDenseMatrix& B,
00478               double ScalarThis ) {
00479   // Check for compatible dimensions
00480   
00481   if (SideA=='R') {
00482     if (M_ != B.M() || 
00483   N_ != A.N() ||
00484   B.N() != A.M() ) EPETRA_CHK_ERR(-1); // Return error
00485   }
00486   else if (SideA=='L') {
00487     if (M_ != A.M() || 
00488   N_ != B.N() ||
00489   A.N() != B.M() ) EPETRA_CHK_ERR(-1); // Return error
00490   }
00491   else {
00492     EPETRA_CHK_ERR(-2); // Return error, incorrect value for SideA
00493   }
00494   
00495   // Call SYMM function
00496   SYMM(SideA, A.UPLO(), M_, N_, ScalarAB, A.A(), A.LDA(), 
00497        B.A(), B.LDA(), ScalarThis, A_, LDA_);
00498   long int nflops = 2*M_;
00499   nflops *= N_;
00500   nflops *= A.N();
00501   if (ScalarAB != 1.0) nflops += M_*N_;
00502   if (ScalarThis != 0.0) nflops += M_*N_;
00503   UpdateFlops((double)nflops);
00504 
00505   return(0);
00506 }
00507 //=========================================================================
00508 void Epetra_SerialDenseMatrix::Print(ostream& os) const {
00509   os << endl;
00510   if(CV_ == Copy)
00511     os << "Data access mode: Copy" << endl;
00512   else
00513     os << "Data access mode: View" << endl;
00514   if(A_Copied_)
00515     os << "A_Copied: yes" << endl;
00516   else
00517     os << "A_Copied: no" << endl;
00518   os << "Rows(M): " << M_ << endl;
00519   os << "Columns(N): " << N_ << endl;
00520   os << "LDA: " << LDA_ << endl;
00521   if(M_ == 0 || N_ == 0)
00522     os << "(matrix is empty, no values to display)" << endl;
00523   else
00524     for(int i = 0; i < M_; i++) {
00525       for(int j = 0; j < N_; j++){
00526         os << (*this)(i,j) << " ";
00527       }
00528       os << endl;
00529     }
00530 }
00531 //=========================================================================
00532 int Epetra_SerialDenseMatrix::Random() {
00533 
00534   Epetra_Util util;
00535 
00536   for(int j = 0; j < N_; j++) {
00537     double* arrayPtr = A_ + (j * LDA_);
00538     for(int i = 0; i < M_; i++) {
00539       *arrayPtr++ = util.RandomDouble();
00540     }
00541   }
00542   
00543   return(0);
00544 }
00545   
00546 //=========================================================================
00547 int Epetra_SerialDenseMatrix::Apply(const Epetra_SerialDenseMatrix& X,
00548             Epetra_SerialDenseMatrix& Y) {
00549   return Multiply(UseTranspose(), X, Y);
00550 }

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