Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Epetra_SerialDenseMatrix.cpp
Go to the documentation of this file.
1
2//@HEADER
3// ************************************************************************
4//
5// Epetra: Linear Algebra Services Package
6// Copyright 2011 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ************************************************************************
41//@HEADER
42
46#include "Epetra_Util.h"
47
48//=============================================================================
51 Epetra_Object(-1, false),
52 M_(0),
53 N_(0),
54 A_Copied_(false),
55 CV_(Copy),
56 LDA_(0),
57 A_(0),
58 UseTranspose_(false)
59{
60 if (set_object_label) {
61 SetLabel("Epetra::SerialDenseMatrix");
62 }
63}
64
65//=============================================================================
67 bool set_object_label)
69 Epetra_Object(-1, false),
70 M_(0),
71 N_(0),
72 A_Copied_(false),
73 CV_(Copy),
74 LDA_(0),
75 A_(0),
76 UseTranspose_(false)
77{
78 if (set_object_label) {
79 SetLabel("Epetra::SerialDenseMatrix");
80 }
81 if(NumRows < 0)
82 throw ReportError("NumRows = " + toString(NumRows) + ". Should be >= 0", -1);
83 if(NumCols < 0)
84 throw ReportError("NumCols = " + toString(NumCols) + ". Should be >= 0", -1);
85
86 int errorcode = Shape(NumRows, NumCols);
87 if(errorcode != 0)
88 throw ReportError("Shape returned non-zero value", errorcode);
89}
90
91//=============================================================================
93 int NumRows, int NumCols,
94 bool set_object_label)
96 Epetra_Object(-1, false),
97 M_(NumRows),
98 N_(NumCols),
99 A_Copied_(false),
100 CV_(CV_in),
101 LDA_(LDA_in),
102 A_(A_in),
103 UseTranspose_(false)
104{
105 if (set_object_label) {
106 SetLabel("Epetra::SerialDenseMatrix");
107 }
108 if(A_in == 0)
109 throw ReportError("Null pointer passed as A parameter.", -3);
110 if(NumRows < 0)
111 throw ReportError("NumRows = " + toString(NumRows) + ". Should be >= 0", -1);
112 if(NumCols < 0)
113 throw ReportError("NumCols = " + toString(NumCols) + ". Should be >= 0", -1);
114 if(LDA_in < 0)
115 throw ReportError("LDA = " + toString(LDA_in) + ". Should be >= 0", -1);
116
117 if (CV_in == Copy) {
118 LDA_ = M_;
119 const int newsize = LDA_ * N_;
120 if (newsize > 0) {
121 A_ = new double[newsize];
122 CopyMat(A_in, LDA_in, M_, N_, A_, LDA_);
123 A_Copied_ = true;
124 }
125 else {
126 A_ = 0;
127 }
128 }
129}
130//=============================================================================
132 : Epetra_CompObject(Source),
133 M_(Source.M_),
134 N_(Source.N_),
135 A_Copied_(false),
136 CV_(Source.CV_),
137 LDA_(Source.LDA_),
138 A_(Source.A_),
139 UseTranspose_(false)
140{
141 SetLabel(Source.Label());
142 if(CV_ == Copy) {
143 LDA_ = M_;
144 const int newsize = LDA_ * N_;
145 if(newsize > 0) {
146 A_ = new double[newsize];
147 CopyMat(Source.A_, Source.LDA_, M_, N_, A_, LDA_);
148 A_Copied_ = true;
149 }
150 else {
151 A_ = 0;
152 }
153 }
154}
155
156//=============================================================================
157int Epetra_SerialDenseMatrix::Reshape(int NumRows, int NumCols) {
158 if(NumRows < 0 || NumCols < 0)
159 return(-1);
160
161 double* A_tmp = 0;
162 const int newsize = NumRows * NumCols;
163
164 if(newsize > 0) {
165 // Allocate space for new matrix
166 A_tmp = new double[newsize];
167 for (int k = 0; k < newsize; k++)
168 A_tmp[k] = 0.0; // Zero out values
169 int M_tmp = EPETRA_MIN(M_, NumRows);
170 int N_tmp = EPETRA_MIN(N_, NumCols);
171 if (A_ != 0)
172 CopyMat(A_, LDA_, M_tmp, N_tmp, A_tmp, NumRows); // Copy principal submatrix of A to new A
173 }
174 CleanupData(); // Get rid of anything that might be already allocated
175 M_ = NumRows;
176 N_ = NumCols;
178 if(newsize > 0) {
179 A_ = A_tmp; // Set pointer to new A
180 A_Copied_ = true;
181 }
182
183 return(0);
184}
185//=============================================================================
186int Epetra_SerialDenseMatrix::Shape(int NumRows, int NumCols) {
187 if(NumRows < 0 || NumCols < 0)
188 return(-1);
189
190 CleanupData(); // Get rid of anything that might be already allocated
191 M_ = NumRows;
192 N_ = NumCols;
194 const int newsize = LDA_ * N_;
195 if(newsize > 0) {
196 A_ = new double[newsize];
197 for (int k = 0; k < newsize; k++)
198 A_[k] = 0.0; // Zero out values
199 A_Copied_ = true;
200 }
201
202 return(0);
203}
204//=============================================================================
209//=============================================================================
211{
212 if (A_Copied_)
213 delete [] A_;
214 A_ = 0;
215 A_Copied_ = false;
216 M_ = 0;
217 N_ = 0;
218 LDA_ = 0;
219}
220//=============================================================================
222 if(this == &Source)
223 return(*this); // Special case of source same as target
224 if((CV_ == View) && (Source.CV_ == View) && (A_ == Source.A_))
225 return(*this); // Special case of both are views to same data.
226
227 if(std::strcmp(Label(), Source.Label()) != 0)
228 throw ReportError("operator= type mismatch (lhs = " + std::string(Label()) +
229 ", rhs = " + std::string(Source.Label()) + ").", -5);
230
231 if(Source.CV_ == View) {
232 if(CV_ == Copy) { // C->V only
233 CleanupData();
234 CV_ = View;
235 }
236 M_ = Source.M_; // C->V and V->V
237 N_ = Source.N_;
238 LDA_ = Source.LDA_;
239 A_ = Source.A_;
240 }
241 else {
242 if(CV_ == View) { // V->C
243 CV_ = Copy;
244 M_ = Source.M_;
245 N_ = Source.N_;
246 LDA_ = Source.M_;
247 const int newsize = LDA_ * N_;
248 if(newsize > 0) {
249 A_ = new double[newsize];
250 A_Copied_ = true;
251 }
252 else {
253 A_ = 0;
254 A_Copied_ = false;
255 }
256 }
257 else { // C->C
258 if((Source.M_ <= LDA_) && (Source.N_ == N_)) { // we don't need to reallocate
259 M_ = Source.M_;
260 N_ = Source.N_;
261 }
262 else { // we need to allocate more space (or less space)
263 CleanupData();
264 M_ = Source.M_;
265 N_ = Source.N_;
266 LDA_ = Source.M_;
267 const int newsize = LDA_ * N_;
268 if(newsize > 0) {
269 A_ = new double[newsize];
270 A_Copied_ = true;
271 }
272 }
273 }
274 CopyMat(Source.A_, Source.LDA_, M_, N_, A_, LDA_); // V->C and C->C
275 }
276
277 return(*this);
278}
279
280
281//=============================================================================
283{
284 if (M_ != rhs.M_ || N_ != rhs.N_) return(false);
285
286 const double* A_tmp = A_;
287 const double* rhsA = rhs.A_;
288
289 for(int j=0; j<N_; ++j) {
290 int offset = j*LDA_;
291 int rhsOffset = j*rhs.LDA_;
292 for(int i=0; i<M_; ++i) {
293 if (std::abs(A_tmp[offset+i] - rhsA[rhsOffset+i]) > Epetra_MinDouble) {
294 return(false);
295 }
296 }
297 }
298
299 return(true);
300}
301
302//=============================================================================
304 if (M() != Source.M())
305 throw ReportError("Row dimension of source = " + toString(Source.M()) +
306 " is different than row dimension of target = " + toString(LDA()), -1);
307 if (N() != Source.N())
308 throw ReportError("Column dimension of source = " + toString(Source.N()) +
309 " is different than column dimension of target = " + toString(N()), -2);
310
311 CopyMat(Source.A(), Source.LDA(), Source.M(), Source.N(), A(), LDA(), true);
312 return(*this);
313}
314//=============================================================================
315void Epetra_SerialDenseMatrix::CopyMat(const double* Source,
316 int Source_LDA,
317 int NumRows,
318 int NumCols,
319 double* Target,
320 int Target_LDA,
321 bool add)
322{
323 int i;
324 double* tptr = Target;
325 const double* sptr = Source;
326 if (add) {
327 for(int j=0; j<NumCols; ++j) {
328 for(i=0; i<NumRows; ++i) {
329 tptr[i] += sptr[i];
330 }
331
332 tptr += Target_LDA;
333 sptr += Source_LDA;
334 }
335 }
336 else {
337 for(int j=0; j<NumCols; ++j) {
338 for(i=0; i<NumRows; ++i) {
339 tptr[i] = sptr[i];
340 }
341
342 tptr += Target_LDA;
343 sptr += Source_LDA;
344 }
345 }
346/*
347 double* targetPtr = Target;
348 double* sourcePtr = Source;
349 double* targetEnd = 0;
350
351 for (j = 0; j < NumCols; j++) {
352 targetPtr = Target + j * Target_LDA;
353 targetEnd = targetPtr+NumRows;
354 sourcePtr = Source + j * Source_LDA;
355 if (add) {
356 while(targetPtr != targetEnd)
357 *targetPtr++ += *sourcePtr++;
358 }
359 else {
360 while(targetPtr != targetEnd)
361 *targetPtr++ = *sourcePtr++;
362 }
363 }
364*/
365 return;
366}
367//=============================================================================
369
370 int i, j;
371
372 double anorm = 0.0;
373 double * ptr;
374 for (j=0; j<N_; j++) {
375 double sum=0.0;
376 ptr = A_ + j*LDA_;
377 for (i=0; i<M_; i++) sum += std::abs(*ptr++);
378 anorm = EPETRA_MAX(anorm, sum);
379 }
380 UpdateFlops((double)N_*(double)N_);
381 return(anorm);
382}
383//=============================================================================
385
386 int i, j;
387
388 double anorm = 0.0;
389 double * ptr;
390
391 // Loop across columns in inner loop. Most expensive memory access, but
392 // requires no extra storage.
393 for (i=0; i<M_; i++) {
394 double sum=0.0;
395 ptr = A_ + i;
396 for (j=0; j<N_; j++) {
397 sum += std::abs(*ptr);
398 ptr += LDA_;
399 }
400 anorm = EPETRA_MAX(anorm, sum);
401 }
402 UpdateFlops((double)N_*(double)N_);
403 return(anorm);
404}
405//=============================================================================
407
408 int i, j;
409
410 double * ptr;
411 for (j=0; j<N_; j++) {
412 ptr = A_ + j*LDA_;
413 for (i=0; i<M_; i++) { *ptr = ScalarA * (*ptr); ptr++; }
414 }
415 UpdateFlops((double)N_*(double)N_);
416 return(0);
417
418}
419//=========================================================================
420int Epetra_SerialDenseMatrix::Multiply (char TransA, char TransB, double ScalarAB,
421 const Epetra_SerialDenseMatrix& A_in,
423 double ScalarThis ) {
424 // Check for compatible dimensions
425
426 if (TransA!='T' && TransA!='N') EPETRA_CHK_ERR(-2); // Return error
427 if (TransB!='T' && TransB!='N') EPETRA_CHK_ERR(-3);
428
429 int A_nrows = (TransA=='T') ? A_in.N() : A_in.M();
430 int A_ncols = (TransA=='T') ? A_in.M() : A_in.N();
431 int B_nrows = (TransB=='T') ? B.N() : B.M();
432 int B_ncols = (TransB=='T') ? B.M() : B.N();
433
434 if (M_ != A_nrows ||
435 A_ncols != B_nrows ||
436 N_ != B_ncols ) EPETRA_CHK_ERR(-1); // Return error
437
438
439 // Call GEMM function
440 GEMM(TransA, TransB, M_, N_, A_ncols, ScalarAB, A_in.A(), A_in.LDA(),
441 B.A(), B.LDA(), ScalarThis, A_, LDA_);
442 long int nflops = 2*M_;
443 nflops *= N_;
444 nflops *= A_ncols;
445 if (ScalarAB != 1.0) nflops += M_*N_;
446 if (ScalarThis != 0.0) nflops += M_*N_;
447 UpdateFlops((double)nflops);
448
449 return(0);
450}
451
452//=========================================================================
456{
457 int A_nrows = M();
458 int x_nrows = x.M();
459 int y_nrows = y.M();
460 int A_ncols = N();
461 int x_ncols = x.N();
462 int y_ncols = y.N();
463
464 if (transA) {
465 if (x_nrows != A_nrows) EPETRA_CHK_ERR(-1);
466 if (y_ncols != x_ncols || y_nrows != A_ncols) y.Reshape(A_ncols, x_ncols);
467 }
468 else {
469 if (x_nrows != A_ncols) EPETRA_CHK_ERR(-1);
470 if (y_ncols != x_ncols || y_nrows != A_nrows) y.Reshape(A_nrows, x_ncols);
471 }
472
473 double scalar0 = 0.0;
474 double scalar1 = 1.0;
475
476 int err = 0;
477 if (transA) {
478 err = y.Multiply('T', 'N', scalar1, *this, x, scalar0);
479 }
480 else {
481 err = y.Multiply('N', 'N', scalar1, *this, x, scalar0);
482 }
483 // FIXME (mfh 06 Mar 2013) Why aren't we returning err instead of 0?
484 // In any case, I'm going to silence the unused value compiler
485 // warning for now. I'm not changing the return value because I
486 // don't want to break any downstream code that depends on this
487 // method always returning 0.
488 (void) err;
489
490 return(0);
491}
492
493//=========================================================================
494int Epetra_SerialDenseMatrix::Multiply (char SideA, double ScalarAB,
495 const Epetra_SerialSymDenseMatrix& A_in,
497 double ScalarThis ) {
498 // Check for compatible dimensions
499
500 if (SideA=='R') {
501 if (M_ != B.M() ||
502 N_ != A_in.N() ||
503 B.N() != A_in.M() ) EPETRA_CHK_ERR(-1); // Return error
504 }
505 else if (SideA=='L') {
506 if (M_ != A_in.M() ||
507 N_ != B.N() ||
508 A_in.N() != B.M() ) EPETRA_CHK_ERR(-1); // Return error
509 }
510 else {
511 EPETRA_CHK_ERR(-2); // Return error, incorrect value for SideA
512 }
513
514 // Call SYMM function
515 SYMM(SideA, A_in.UPLO(), M_, N_, ScalarAB, A_in.A(), A_in.LDA(),
516 B.A(), B.LDA(), ScalarThis, A_, LDA_);
517 long int nflops = 2*M_;
518 nflops *= N_;
519 nflops *= A_in.N();
520 if (ScalarAB != 1.0) nflops += M_*N_;
521 if (ScalarThis != 0.0) nflops += M_*N_;
522 UpdateFlops((double)nflops);
523
524 return(0);
525}
526//=========================================================================
527void Epetra_SerialDenseMatrix::Print(std::ostream& os) const {
528 os << std::endl;
529 if(CV_ == Copy)
530 os << "Data access mode: Copy" << std::endl;
531 else
532 os << "Data access mode: View" << std::endl;
533 if(A_Copied_)
534 os << "A_Copied: yes" << std::endl;
535 else
536 os << "A_Copied: no" << std::endl;
537 os << "Rows(M): " << M_ << std::endl;
538 os << "Columns(N): " << N_ << std::endl;
539 os << "LDA: " << LDA_ << std::endl;
540 if(M_ == 0 || N_ == 0)
541 os << "(matrix is empty, no values to display)" << std::endl;
542 else
543 for(int i = 0; i < M_; i++) {
544 for(int j = 0; j < N_; j++){
545 os << (*this)(i,j) << " ";
546 }
547 os << std::endl;
548 }
549}
550//=========================================================================
552
553 Epetra_Util util;
554
555 for(int j = 0; j < N_; j++) {
556 double* arrayPtr = A_ + (j * LDA_);
557 for(int i = 0; i < M_; i++) {
558 *arrayPtr++ = util.RandomDouble();
559 }
560 }
561
562 return(0);
563}
564
565//=========================================================================
#define EPETRA_MIN(x, y)
#define EPETRA_MAX(x, y)
#define EPETRA_CHK_ERR(a)
const double Epetra_MinDouble
Epetra_DataAccess
void SYMM(const char SIDE, const char UPLO, const int M, const int N, const float ALPHA, const float *A, const int LDA, const float *B, const int LDB, const float BETA, float *C, const int LDC) const
Epetra_BLAS symmetric matrix-matrix multiply function (SSYMM)
void GEMM(const char TRANSA, const char TRANSB, const int M, const int N, const int K, const float ALPHA, const float *A, const int LDA, const float *B, const int LDB, const float BETA, float *C, const int LDC) const
Epetra_BLAS matrix-matrix multiply function (SGEMM)
Epetra_CompObject: Functionality and data that is common to all computational classes.
void UpdateFlops(int Flops_in) const
Increment Flop count for this object.
Epetra_Object: The base Epetra class.
virtual int ReportError(const std::string Message, int ErrorCode) const
Error reporting method.
std::string toString(const int &x) const
virtual void SetLabel(const char *const Label)
Epetra_Object Label definition using char *.
Epetra_SerialDenseMatrix: A class for constructing and using real double precision general dense matr...
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
double * A() const
Returns pointer to the this matrix.
virtual void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
int LDA() const
Returns the leading dimension of the this matrix.
Epetra_SerialDenseMatrix(bool set_object_label=true)
Default constructor; defines a zero size object.
virtual double NormInf() const
Computes the Infinity-Norm of the this matrix.
int Random()
Set matrix values to random numbers.
int Scale(double ScalarA)
Inplace scalar-matrix product A = a A.
void CopyMat(const double *Source, int Source_LDA, int NumRows, int NumCols, double *Target, int Target_LDA, bool add=false)
Epetra_SerialDenseMatrix & operator=(const Epetra_SerialDenseMatrix &Source)
Value copy from one matrix to another.
bool operator==(const Epetra_SerialDenseMatrix &rhs) const
Comparison operator.
Epetra_SerialDenseMatrix & operator+=(const Epetra_SerialDenseMatrix &Source)
Add one matrix to another.
int Shape(int NumRows, int NumCols)
Set dimensions of a Epetra_SerialDenseMatrix object; init values to zero.
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_SerialDenseMatrix &A, const Epetra_SerialDenseMatrix &B, double ScalarThis)
Matrix-Matrix multiplication, this = ScalarThis*this + ScalarAB*A*B.
int M() const
Returns row dimension of system.
int N() const
Returns column dimension of system.
virtual const char * Label() const
Returns a character string describing the operator.
virtual double NormOne() const
Computes the 1-Norm of the this matrix.
virtual ~Epetra_SerialDenseMatrix()
Epetra_SerialDenseMatrix destructor.
virtual int Apply(const Epetra_SerialDenseMatrix &X, Epetra_SerialDenseMatrix &Y)
Returns the result of a Epetra_SerialDenseOperator applied to a Epetra_SerialDenseMatrix X in Y.
int Reshape(int NumRows, int NumCols)
Reshape a Epetra_SerialDenseMatrix object.
Epetra_SerialSymDenseMatrix: A class for constructing and using symmetric positive definite dense mat...
char UPLO() const
Returns character value of UPLO used by LAPACK routines.
Epetra_Util: The Epetra Util Wrapper Class.
Definition Epetra_Util.h:79
double RandomDouble()
Returns a random double on the interval (-1.0,1.0)