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