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, double* A, int LDA,
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),
00088     LDA_(LDA),
00089     A_(A),
00090     UseTranspose_(false)
00091 {
00092   if (set_object_label) {
00093     SetLabel("Epetra::SerialDenseMatrix");
00094   }
00095   if(A == 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 < 0)
00102   throw ReportError("LDA = " + toString(LDA) + ". Should be >= 0", -1);
00103 
00104   if (CV == Copy) {
00105     LDA_ = M_;
00106     const int newsize = LDA_ * N_;
00107     if (newsize > 0) {
00108       A_ = new double[newsize];
00109       CopyMat(A, LDA, 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 Epetra_SerialDenseMatrix& Epetra_SerialDenseMatrix::operator+= ( const Epetra_SerialDenseMatrix & Source) {
00269   if (M() != Source.M()) 
00270     throw ReportError("Row dimension of source = " + toString(Source.M()) +
00271                       " is different than  row dimension of target = " + toString(LDA()), -1);
00272   if (N() != Source.N()) 
00273     throw ReportError("Column dimension of source = " + toString(Source.N()) +
00274                       " is different than column dimension of target = " + toString(N()), -2);
00275 
00276   CopyMat(Source.A(), Source.LDA(), Source.M(), Source.N(), A(), LDA(), true);
00277   return(*this);
00278 }
00279 //=============================================================================
00280 void Epetra_SerialDenseMatrix::CopyMat(double* Source,
00281                                        int Source_LDA,
00282                                        int NumRows,
00283                                        int NumCols, 
00284                                        double* Target,
00285                                        int Target_LDA,
00286                                        bool add)
00287 {
00288   int i;
00289   double* tptr = Target;
00290   double* sptr = Source;
00291   if (add) {
00292     for(int j=0; j<NumCols; ++j) {
00293       for(i=0; i<NumRows; ++i) {
00294         tptr[i] += sptr[i];
00295       }
00296 
00297       tptr += Target_LDA;
00298       sptr += Source_LDA;
00299     }
00300   }
00301   else {
00302     for(int j=0; j<NumCols; ++j) {
00303       for(i=0; i<NumRows; ++i) {
00304         tptr[i] = sptr[i];
00305       }
00306 
00307       tptr += Target_LDA;
00308       sptr += Source_LDA;
00309     }
00310   }
00311 /*
00312   double* targetPtr = Target;
00313   double* sourcePtr = Source;
00314   double* targetEnd = 0;
00315  
00316   for (j = 0; j < NumCols; j++) {
00317     targetPtr = Target + j * Target_LDA;
00318     targetEnd = targetPtr+NumRows;
00319     sourcePtr = Source + j * Source_LDA;
00320     if (add) {
00321       while(targetPtr != targetEnd) 
00322     *targetPtr++ += *sourcePtr++;
00323     }
00324     else {
00325       while(targetPtr != targetEnd) 
00326     *targetPtr++ = *sourcePtr++;
00327     }
00328   }
00329 */
00330   return;
00331 }
00332 //=============================================================================
00333 double Epetra_SerialDenseMatrix::NormOne() const {
00334 
00335   int i, j;
00336 
00337     double anorm = 0.0;
00338     double * ptr;
00339     for (j=0; j<N_; j++) {
00340       double sum=0.0;
00341       ptr = A_ + j*LDA_;
00342       for (i=0; i<M_; i++) sum += std::abs(*ptr++);
00343       anorm = EPETRA_MAX(anorm, sum);
00344     }
00345     UpdateFlops(N_*N_);
00346     return(anorm);
00347 }
00348 //=============================================================================
00349 double Epetra_SerialDenseMatrix::NormInf() const {
00350 
00351   int i, j;
00352 
00353     double anorm = 0.0;
00354     double * ptr;
00355 
00356     // Loop across columns in inner loop.  Most expensive memory access, but 
00357     // requires no extra storage.
00358     for (i=0; i<M_; i++) {
00359       double sum=0.0;
00360       ptr = A_ + i;
00361       for (j=0; j<N_; j++) {
00362   sum += std::abs(*ptr);
00363   ptr += LDA_;
00364       }
00365       anorm = EPETRA_MAX(anorm, sum);
00366     }
00367     UpdateFlops(N_*N_);
00368     return(anorm);
00369 }
00370 //=============================================================================
00371 int Epetra_SerialDenseMatrix::Scale(double ScalarA) {
00372 
00373   int i, j;
00374   
00375   double * ptr;
00376   for (j=0; j<N_; j++) {
00377     ptr = A_ + j*LDA_;
00378     for (i=0; i<M_; i++) { *ptr = ScalarA * (*ptr); ptr++; }
00379   }
00380   UpdateFlops(N_*N_);
00381   return(0);
00382   
00383 }
00384 //=========================================================================
00385 int Epetra_SerialDenseMatrix::Multiply (char TransA, char TransB, double ScalarAB, 
00386           const Epetra_SerialDenseMatrix& A, 
00387           const Epetra_SerialDenseMatrix& B,
00388           double ScalarThis ) {
00389   // Check for compatible dimensions
00390 
00391   if (TransA!='T' && TransA!='N') EPETRA_CHK_ERR(-2); // Return error
00392   if (TransB!='T' && TransB!='N') EPETRA_CHK_ERR(-3);
00393   
00394   int A_nrows = (TransA=='T') ? A.N() : A.M();
00395   int A_ncols = (TransA=='T') ? A.M() : A.N();
00396   int B_nrows = (TransB=='T') ? B.N() : B.M();
00397   int B_ncols = (TransB=='T') ? B.M() : B.N();
00398   
00399   if (M_        != A_nrows     ||
00400       A_ncols   != B_nrows     ||
00401       N_        != B_ncols  ) EPETRA_CHK_ERR(-1); // Return error
00402 
00403     
00404   // Call GEMM function
00405   GEMM(TransA, TransB, M_, N_, A_ncols, ScalarAB, A.A(), A.LDA(), 
00406        B.A(), B.LDA(), ScalarThis, A_, LDA_);
00407   long int nflops = 2*M_;
00408   nflops *= N_;
00409   nflops *= A_ncols;
00410   if (ScalarAB != 1.0) nflops += M_*N_;
00411   if (ScalarThis != 0.0) nflops += M_*N_;
00412   UpdateFlops(nflops);
00413 
00414   return(0);
00415 }
00416 
00417 //=========================================================================
00418 int  Epetra_SerialDenseMatrix::Multiply (bool transA,
00419                                          const Epetra_SerialDenseMatrix& x,
00420                                          Epetra_SerialDenseMatrix& y)
00421 {
00422   int A_nrows = M();
00423   int x_nrows = x.M();
00424   int y_nrows = y.M();
00425   int A_ncols = N();
00426   int x_ncols = x.N();
00427   int y_ncols = y.N();
00428 
00429   if (transA) {
00430     if (x_nrows != A_nrows) EPETRA_CHK_ERR(-1);
00431     if (y_ncols != x_ncols || y_nrows != A_ncols) y.Reshape(A_ncols, x_ncols);
00432   }
00433   else {
00434     if (x_nrows != A_ncols) EPETRA_CHK_ERR(-1);
00435     if (y_ncols != x_ncols || y_nrows != A_nrows) y.Reshape(A_nrows, x_ncols);
00436   }
00437 
00438   double scalar0 = 0.0;
00439   double scalar1 = 1.0;
00440 
00441   int err = 0;
00442   if (transA) {
00443     err = y.Multiply('T', 'N', scalar1, *this, x, scalar0);
00444   }
00445   else {
00446     err = y.Multiply('N', 'N', scalar1, *this, x, scalar0);
00447   }
00448 
00449   return(0);
00450 }
00451 
00452 //=========================================================================
00453 int  Epetra_SerialDenseMatrix::Multiply (char SideA, double ScalarAB, 
00454               const Epetra_SerialSymDenseMatrix& A, 
00455               const Epetra_SerialDenseMatrix& B,
00456               double ScalarThis ) {
00457   // Check for compatible dimensions
00458   
00459   if (SideA=='R') {
00460     if (M_ != B.M() || 
00461   N_ != A.N() ||
00462   B.N() != A.M() ) EPETRA_CHK_ERR(-1); // Return error
00463   }
00464   else if (SideA=='L') {
00465     if (M_ != A.M() || 
00466   N_ != B.N() ||
00467   A.N() != B.M() ) EPETRA_CHK_ERR(-1); // Return error
00468   }
00469   else {
00470     EPETRA_CHK_ERR(-2); // Return error, incorrect value for SideA
00471   }
00472   
00473   // Call SYMM function
00474   SYMM(SideA, A.UPLO(), M_, N_, ScalarAB, A.A(), A.LDA(), 
00475        B.A(), B.LDA(), ScalarThis, A_, LDA_);
00476   long int nflops = 2*M_;
00477   nflops *= N_;
00478   nflops *= A.N();
00479   if (ScalarAB != 1.0) nflops += M_*N_;
00480   if (ScalarThis != 0.0) nflops += M_*N_;
00481   UpdateFlops(nflops);
00482 
00483   return(0);
00484 }
00485 //=========================================================================
00486 void Epetra_SerialDenseMatrix::Print(ostream& os) const {
00487   os << endl;
00488   if(CV_ == Copy)
00489     os << "Data access mode: Copy" << endl;
00490   else
00491     os << "Data access mode: View" << endl;
00492   if(A_Copied_)
00493     os << "A_Copied: yes" << endl;
00494   else
00495     os << "A_Copied: no" << endl;
00496   os << "Rows(M): " << M_ << endl;
00497   os << "Columns(N): " << N_ << endl;
00498   os << "LDA: " << LDA_ << endl;
00499   if(M_ == 0 || N_ == 0)
00500     os << "(matrix is empty, no values to display)" << endl;
00501   else
00502     for(int i = 0; i < M_; i++) {
00503       for(int j = 0; j < N_; j++){
00504         os << (*this)(i,j) << " ";
00505       }
00506       os << endl;
00507     }
00508 }
00509 //=========================================================================
00510 int Epetra_SerialDenseMatrix::Random() {
00511 
00512   Epetra_Util util;
00513 
00514   for(int j = 0; j < N_; j++) {
00515     double* arrayPtr = A_ + (j * LDA_);
00516     for(int i = 0; i < M_; i++) {
00517       *arrayPtr++ = util.RandomDouble();
00518     }
00519   }
00520   
00521   return(0);
00522 }
00523   
00524 //=========================================================================
00525 int Epetra_SerialDenseMatrix::Apply(const Epetra_SerialDenseMatrix& X,
00526             Epetra_SerialDenseMatrix& Y) {
00527   return Multiply(UseTranspose(), X, Y);
00528 }

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