Teuchos - Trilinos Tools Package  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Teuchos_SerialDenseMatrix.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Teuchos: Common Tools Package
4 //
5 // Copyright 2004 NTESS and the Teuchos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef _TEUCHOS_SERIALDENSEMATRIX_HPP_
11 #define _TEUCHOS_SERIALDENSEMATRIX_HPP_
12 
16 #include "Teuchos_CompObject.hpp"
17 #include "Teuchos_BLAS.hpp"
18 #include "Teuchos_ScalarTraits.hpp"
19 #include "Teuchos_DataAccess.hpp"
20 #include "Teuchos_ConfigDefs.hpp"
21 #include "Teuchos_Assert.hpp"
23 #include <cstddef>
24 #include <utility>
25 
34 namespace Teuchos {
35 
36 template<typename OrdinalType, typename ScalarType>
37 class SerialDenseMatrix : public CompObject, public BLAS<OrdinalType, ScalarType>
38 {
39 public:
40 
42  typedef OrdinalType ordinalType;
44  typedef ScalarType scalarType;
45 
47 
48 
50 
53  SerialDenseMatrix() = default;
54 
56 
64  SerialDenseMatrix(OrdinalType numRows, OrdinalType numCols, bool zeroOut = true);
65 
67 
75  SerialDenseMatrix(DataAccess CV, ScalarType* values, OrdinalType stride, OrdinalType numRows, OrdinalType numCols);
76 
78 
84 
86 
89 
91 
103  SerialDenseMatrix(DataAccess CV, const SerialDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRows, OrdinalType numCols, OrdinalType startRow=0, OrdinalType startCol=0);
104 
106  virtual ~SerialDenseMatrix();
108 
110 
111 
122  int shape(OrdinalType numRows, OrdinalType numCols);
123 
125  int shapeUninitialized(OrdinalType numRows, OrdinalType numCols);
126 
128 
138  int reshape(OrdinalType numRows, OrdinalType numCols);
139 
140 
142 
144 
145 
147 
154 
156 
162 
164 
167  SerialDenseMatrix<OrdinalType, ScalarType>& operator= (const ScalarType value) { putScalar(value); return(*this); }
168 
170 
174  int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero() );
175 
177 
181 
183  int random();
184 
186 
188 
189 
191 
198  ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
199 
201 
208  const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
209 
211 
216  ScalarType* operator [] (OrdinalType colIndex);
217 
219 
224  const ScalarType* operator [] (OrdinalType colIndex) const;
225 
227 
228  ScalarType* values() const { return values_; }
229 
231 
233 
234 
236 
240 
242 
246 
248 
252 
254 
258  int scale ( const ScalarType alpha );
259 
261 
268 
270 
284  int multiply (ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
285 
287 
298  int multiply (ESide sideA, ScalarType alpha, const SerialSymDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
299 
301 
303 
304 
306 
309  bool operator== (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand) const;
310 
312 
315  bool operator!= (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand) const;
316 
318 
320 
321 
323  OrdinalType numRows() const { return(numRows_); }
324 
326  OrdinalType numCols() const { return(numCols_); }
327 
329  OrdinalType stride() const { return(stride_); }
330 
332  bool empty() const { return(numRows_ == 0 || numCols_ == 0); }
334 
336 
337 
340 
343 
347 
349 
350  virtual std::ostream& print(std::ostream& os) const;
352 
354 protected:
355  void copyMat(ScalarType* inputMatrix, OrdinalType strideInput,
356  OrdinalType numRows, OrdinalType numCols, ScalarType* outputMatrix,
357  OrdinalType strideOutput, OrdinalType startRow, OrdinalType startCol,
358  ScalarType alpha = ScalarTraits<ScalarType>::zero() );
359  void deleteArrays();
360  void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
361 
362  static ScalarType*
363  allocateValues(const OrdinalType numRows,
364  const OrdinalType numCols)
365  {
366  const size_t size = size_t(numRows) * size_t(numCols);
367  return new ScalarType[size];
368  }
369 
370  OrdinalType numRows_ = 0;
371  OrdinalType numCols_ = 0;
372  OrdinalType stride_ = 0;
373  bool valuesCopied_ = false;
374  ScalarType* values_ = nullptr;
375 }; // class Teuchos_SerialDenseMatrix
376 
377 //----------------------------------------------------------------------------------------------------
378 // Constructors and Destructor
379 //----------------------------------------------------------------------------------------------------
380 
381 template<typename OrdinalType, typename ScalarType>
383  OrdinalType numRows_in, OrdinalType numCols_in, bool zeroOut
384  )
385  : numRows_(numRows_in),
386  numCols_(numCols_in),
387  stride_(numRows_in),
388  valuesCopied_(true),
389  values_(allocateValues(numRows_in, numCols_in))
390 {
391  if (zeroOut) {
392  putScalar();
393  }
394 }
395 
396 template<typename OrdinalType, typename ScalarType>
398  DataAccess CV, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRows_in,
399  OrdinalType numCols_in
400  )
401  : numRows_(numRows_in),
402  numCols_(numCols_in),
403  stride_(stride_in),
404  valuesCopied_(false),
405  values_(values_in)
406 {
407  if(CV == Copy)
408  {
409  stride_ = numRows_;
410  values_ = allocateValues(stride_, numCols_);
411  copyMat(values_in, stride_in, numRows_, numCols_, values_, stride_, 0, 0);
412  valuesCopied_ = true;
413  }
414 }
415 
416 template<typename OrdinalType, typename ScalarType>
418  : valuesCopied_(true)
419 {
420  if ( trans == Teuchos::NO_TRANS )
421  {
422  numRows_ = Source.numRows_;
423  numCols_ = Source.numCols_;
424 
425  if (!Source.valuesCopied_)
426  {
427  stride_ = Source.stride_;
428  values_ = Source.values_;
429  valuesCopied_ = false;
430  }
431  else
432  {
433  stride_ = numRows_;
434  if(stride_ > 0 && numCols_ > 0) {
435  values_ = allocateValues(stride_, numCols_);
436  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
437  }
438  else {
439  numRows_ = 0; numCols_ = 0; stride_ = 0;
440  valuesCopied_ = false;
441  }
442  }
443  }
445  {
446  numRows_ = Source.numCols_;
447  numCols_ = Source.numRows_;
448  stride_ = numRows_;
449  values_ = allocateValues(stride_, numCols_);
450  for (OrdinalType j=0; j<numCols_; j++) {
451  for (OrdinalType i=0; i<numRows_; i++) {
452  values_[j*stride_ + i] = Teuchos::ScalarTraits<ScalarType>::conjugate(Source.values_[i*Source.stride_ + j]);
453  }
454  }
455  }
456  else
457  {
458  numRows_ = Source.numCols_;
459  numCols_ = Source.numRows_;
460  stride_ = numRows_;
461  values_ = allocateValues(stride_, numCols_);
462  for (OrdinalType j=0; j<numCols_; j++) {
463  for (OrdinalType i=0; i<numRows_; i++) {
464  values_[j*stride_ + i] = Source.values_[i*Source.stride_ + j];
465  }
466  }
467  }
468 }
469 
470 
471 template<typename OrdinalType, typename ScalarType>
474  )
475  : numRows_(Source.numRows_), numCols_(Source.numCols_), stride_(Source.stride_),
476  valuesCopied_(false), values_(Source.values_)
477 {
478  if(CV == Copy)
479  {
480  stride_ = numRows_;
481  values_ = allocateValues(stride_, numCols_);
482  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
483  valuesCopied_ = true;
484  }
485 }
486 
487 
488 template<typename OrdinalType, typename ScalarType>
491  OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType startRow,
492  OrdinalType startCol
493  )
494  : CompObject(), numRows_(numRows_in), numCols_(numCols_in), stride_(Source.stride_),
495  valuesCopied_(false), values_(Source.values_)
496 {
497  if(CV == Copy)
498  {
499  stride_ = numRows_in;
500  values_ = allocateValues(stride_, numCols_in);
501  copyMat(Source.values_, Source.stride_, numRows_in, numCols_in, values_, stride_, startRow, startCol);
502  valuesCopied_ = true;
503  }
504  else // CV == View
505  {
506  values_ = values_ + (stride_ * startCol) + startRow;
507  }
508 }
509 
510 template<typename OrdinalType, typename ScalarType>
512 {
513  deleteArrays();
514 }
515 
516 //----------------------------------------------------------------------------------------------------
517 // Shape methods
518 //----------------------------------------------------------------------------------------------------
519 
520 template<typename OrdinalType, typename ScalarType>
522  OrdinalType numRows_in, OrdinalType numCols_in
523  )
524 {
525  deleteArrays(); // Get rid of anything that might be already allocated
526  numRows_ = numRows_in;
527  numCols_ = numCols_in;
528  stride_ = numRows_;
529  values_ = allocateValues(stride_, numCols_);
530  putScalar();
531  valuesCopied_ = true;
532  return(0);
533 }
534 
535 template<typename OrdinalType, typename ScalarType>
537  OrdinalType numRows_in, OrdinalType numCols_in
538  )
539 {
540  deleteArrays(); // Get rid of anything that might be already allocated
541  numRows_ = numRows_in;
542  numCols_ = numCols_in;
543  stride_ = numRows_;
544  values_ = allocateValues(stride_, numCols_);
545  valuesCopied_ = true;
546  return(0);
547 }
548 
549 template<typename OrdinalType, typename ScalarType>
551  OrdinalType numRows_in, OrdinalType numCols_in
552  )
553 {
554  // Allocate space for new matrix
555  ScalarType* values_tmp = allocateValues(numRows_in, numCols_in);
556  ScalarType zero = ScalarTraits<ScalarType>::zero();
557  for(OrdinalType k = 0; k < numRows_in * numCols_in; k++)
558  {
559  values_tmp[k] = zero;
560  }
561  OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
562  OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
563  if(values_ != 0)
564  {
565  copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
566  numRows_in, 0, 0); // Copy principal submatrix of A to new A
567  }
568  deleteArrays(); // Get rid of anything that might be already allocated
569  numRows_ = numRows_in;
570  numCols_ = numCols_in;
571  stride_ = numRows_;
572  values_ = values_tmp; // Set pointer to new A
573  valuesCopied_ = true;
574  return(0);
575 }
576 
577 //----------------------------------------------------------------------------------------------------
578 // Set methods
579 //----------------------------------------------------------------------------------------------------
580 
581 template<typename OrdinalType, typename ScalarType>
583 {
584  // Set each value of the dense matrix to "value".
585  for(OrdinalType j = 0; j < numCols_; j++)
586  {
587  for(OrdinalType i = 0; i < numRows_; i++)
588  {
589  values_[i + j*stride_] = value_in;
590  }
591  }
592  return 0;
593 }
594 
595 template<typename OrdinalType, typename ScalarType> void
598 {
599  // Notes:
600  // > DefaultBLASImpl::SWAP() uses a deep copy. This fn uses a pointer swap.
601  // > this fn covers both Vector and Matrix, such that some care must be
602  // employed to not swap across types (creating a Vector with non-unitary
603  // numCols_)
604  // > Inherited data that is not currently swapped (since inactive/deprecated):
605  // >> Teuchos::CompObject:
606  // Flops *flopCounter_ [Note: all SerialDenseMatrix ctors initialize a
607  // NULL flop-counter using CompObject(), such that any flop increments
608  // that are computed are not accumulated.]
609  // >> Teuchos::Object: (now removed from inheritance list)
610  // static int tracebackMode (no swap for statics)
611  // std::string label_ (has been reported as a cause of memory overhead)
612 
613  std::swap(values_ , B.values_);
614  std::swap(numRows_, B.numRows_);
615  std::swap(numCols_, B.numCols_);
616  std::swap(stride_, B.stride_);
617  std::swap(valuesCopied_, B.valuesCopied_);
618 }
619 
620 template<typename OrdinalType, typename ScalarType>
622 {
623  // Set each value of the dense matrix to a random value.
624  for(OrdinalType j = 0; j < numCols_; j++)
625  {
626  for(OrdinalType i = 0; i < numRows_; i++)
627  {
628  values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
629  }
630  }
631  return 0;
632 }
633 
634 template<typename OrdinalType, typename ScalarType>
638  )
639 {
640  if(this == &Source)
641  return(*this); // Special case of source same as target
642  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
643  return(*this); // Special case of both are views to same data.
644 
645  // If the source is a view then we will return a view, else we will return a copy.
646  if (!Source.valuesCopied_) {
647  if(valuesCopied_) {
648  // Clean up stored data if this was previously a copy.
649  deleteArrays();
650  }
651  numRows_ = Source.numRows_;
652  numCols_ = Source.numCols_;
653  stride_ = Source.stride_;
654  values_ = Source.values_;
655  }
656  else {
657  // If we were a view, we will now be a copy.
658  if(!valuesCopied_) {
659  numRows_ = Source.numRows_;
660  numCols_ = Source.numCols_;
661  stride_ = Source.numRows_;
662  if(stride_ > 0 && numCols_ > 0) {
663  values_ = allocateValues(stride_, numCols_);
664  valuesCopied_ = true;
665  }
666  else {
667  values_ = 0;
668  }
669  }
670  // If we were a copy, we will stay a copy.
671  else {
672  if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) { // we don't need to reallocate
673  numRows_ = Source.numRows_;
674  numCols_ = Source.numCols_;
675  }
676  else { // we need to allocate more space (or less space)
677  deleteArrays();
678  numRows_ = Source.numRows_;
679  numCols_ = Source.numCols_;
680  stride_ = Source.numRows_;
681  if(stride_ > 0 && numCols_ > 0) {
682  values_ = allocateValues(stride_, numCols_);
683  valuesCopied_ = true;
684  }
685  }
686  }
687  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
688  }
689  return(*this);
690 }
691 
692 template<typename OrdinalType, typename ScalarType>
694 {
695  // Check for compatible dimensions
696  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
697  {
698  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
699  }
700  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, ScalarTraits<ScalarType>::one());
701  return(*this);
702 }
703 
704 template<typename OrdinalType, typename ScalarType>
706 {
707  // Check for compatible dimensions
708  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
709  {
710  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
711  }
712  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, -ScalarTraits<ScalarType>::one());
713  return(*this);
714 }
715 
716 template<typename OrdinalType, typename ScalarType>
718  if(this == &Source)
719  return(*this); // Special case of source same as target
720  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
721  return(*this); // Special case of both are views to same data.
722 
723  // Check for compatible dimensions
724  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
725  {
726  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
727  }
728  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
729  return(*this);
730 }
731 
732 //----------------------------------------------------------------------------------------------------
733 // Accessor methods
734 //----------------------------------------------------------------------------------------------------
735 
736 template<typename OrdinalType, typename ScalarType>
737 inline ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
738 {
739 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
740  checkIndex( rowIndex, colIndex );
741 #endif
742  return(values_[colIndex * stride_ + rowIndex]);
743 }
744 
745 template<typename OrdinalType, typename ScalarType>
746 inline const ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
747 {
748 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
749  checkIndex( rowIndex, colIndex );
750 #endif
751  return(values_[colIndex * stride_ + rowIndex]);
752 }
753 
754 template<typename OrdinalType, typename ScalarType>
755 inline const ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex) const
756 {
757 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
758  checkIndex( 0, colIndex );
759 #endif
760  return(values_ + colIndex * stride_);
761 }
762 
763 template<typename OrdinalType, typename ScalarType>
764 inline ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex)
765 {
766 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
767  checkIndex( 0, colIndex );
768 #endif
769  return(values_ + colIndex * stride_);
770 }
771 
772 //----------------------------------------------------------------------------------------------------
773 // Norm methods
774 //----------------------------------------------------------------------------------------------------
775 
776 template<typename OrdinalType, typename ScalarType>
778 {
779  OrdinalType i, j;
782  ScalarType* ptr;
783  for(j = 0; j < numCols_; j++)
784  {
785  ScalarType sum = 0;
786  ptr = values_ + j * stride_;
787  for(i = 0; i < numRows_; i++)
788  {
790  }
792  if(absSum > anorm)
793  {
794  anorm = absSum;
795  }
796  }
797  // don't compute flop increment unless there is an accumulator
798  if (flopCounter_!=0) updateFlops(numRows_ * numCols_);
799  return(anorm);
800 }
801 
802 template<typename OrdinalType, typename ScalarType>
804 {
805  OrdinalType i, j;
807 
808  for (i = 0; i < numRows_; i++) {
810  for (j=0; j< numCols_; j++) {
811  sum += ScalarTraits<ScalarType>::magnitude(*(values_+i+j*stride_));
812  }
813  anorm = TEUCHOS_MAX( anorm, sum );
814  }
815  // don't compute flop increment unless there is an accumulator
816  if (flopCounter_!=0) updateFlops(numRows_ * numCols_);
817  return(anorm);
818 }
819 
820 template<typename OrdinalType, typename ScalarType>
822 {
823  OrdinalType i, j;
825  for (j = 0; j < numCols_; j++) {
826  for (i = 0; i < numRows_; i++) {
827  anorm += ScalarTraits<ScalarType>::magnitude(values_[i+j*stride_]*values_[i+j*stride_]);
828  }
829  }
831  // don't compute flop increment unless there is an accumulator
832  if (flopCounter_!=0) updateFlops(numRows_ * numCols_);
833  return(anorm);
834 }
835 
836 //----------------------------------------------------------------------------------------------------
837 // Comparison methods
838 //----------------------------------------------------------------------------------------------------
839 
840 template<typename OrdinalType, typename ScalarType>
842 {
843  bool result = 1;
844  if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_))
845  {
846  result = 0;
847  }
848  else
849  {
850  OrdinalType i, j;
851  for(i = 0; i < numRows_; i++)
852  {
853  for(j = 0; j < numCols_; j++)
854  {
855  if((*this)(i, j) != Operand(i, j))
856  {
857  return 0;
858  }
859  }
860  }
861  }
862  return result;
863 }
864 
865 template<typename OrdinalType, typename ScalarType>
867 {
868  return !((*this) == Operand);
869 }
870 
871 //----------------------------------------------------------------------------------------------------
872 // Multiplication method
873 //----------------------------------------------------------------------------------------------------
874 
875 template<typename OrdinalType, typename ScalarType>
877 {
878  this->scale( alpha );
879  return(*this);
880 }
881 
882 template<typename OrdinalType, typename ScalarType>
884 {
885  OrdinalType i, j;
886  ScalarType* ptr;
887 
888  for (j=0; j<numCols_; j++) {
889  ptr = values_ + j*stride_;
890  for (i=0; i<numRows_; i++) { *ptr = alpha * (*ptr); ptr++; }
891  }
892  // don't compute flop increment unless there is an accumulator
893  if (flopCounter_!=0) updateFlops( numRows_*numCols_ );
894  return(0);
895 }
896 
897 template<typename OrdinalType, typename ScalarType>
899 {
900  OrdinalType i, j;
901  ScalarType* ptr;
902 
903  // Check for compatible dimensions
904  if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_))
905  {
906  TEUCHOS_CHK_ERR(-1); // Return error
907  }
908  for (j=0; j<numCols_; j++) {
909  ptr = values_ + j*stride_;
910  for (i=0; i<numRows_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
911  }
912  // don't compute flop increment unless there is an accumulator
913  if (flopCounter_!=0) updateFlops( numRows_*numCols_ );
914  return(0);
915 }
916 
917 template<typename OrdinalType, typename ScalarType>
919 {
920  // Check for compatible dimensions
921  OrdinalType A_nrows = (ETranspChar[transa]!='N') ? A.numCols() : A.numRows();
922  OrdinalType A_ncols = (ETranspChar[transa]!='N') ? A.numRows() : A.numCols();
923  OrdinalType B_nrows = (ETranspChar[transb]!='N') ? B.numCols() : B.numRows();
924  OrdinalType B_ncols = (ETranspChar[transb]!='N') ? B.numRows() : B.numCols();
925  if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols))
926  {
927  TEUCHOS_CHK_ERR(-1); // Return error
928  }
929  // Call GEMM function
930  this->GEMM(transa, transb, numRows_, numCols_, A_ncols, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
931 
932  // don't compute flop increment unless there is an accumulator
933  if (flopCounter_!=0) {
934  double nflops = 2 * numRows_;
935  nflops *= numCols_;
936  nflops *= A_ncols;
937  updateFlops(nflops);
938  }
939  return(0);
940 }
941 
942 template<typename OrdinalType, typename ScalarType>
944 {
945  // Check for compatible dimensions
946  OrdinalType A_nrows = A.numRows(), A_ncols = A.numCols();
947  OrdinalType B_nrows = B.numRows(), B_ncols = B.numCols();
948 
949  if (ESideChar[sideA]=='L') {
950  if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols)) {
951  TEUCHOS_CHK_ERR(-1); // Return error
952  }
953  } else {
954  if ((numRows_ != B_nrows) || (B_ncols != A_nrows) || (numCols_ != A_ncols)) {
955  TEUCHOS_CHK_ERR(-1); // Return error
956  }
957  }
958 
959  // Call SYMM function
961  this->SYMM(sideA, uplo, numRows_, numCols_, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
962 
963  // don't compute flop increment unless there is an accumulator
964  if (flopCounter_!=0) {
965  double nflops = 2 * numRows_;
966  nflops *= numCols_;
967  nflops *= A_ncols;
968  updateFlops(nflops);
969  }
970  return(0);
971 }
972 
973 template<typename OrdinalType, typename ScalarType>
974 std::ostream& SerialDenseMatrix<OrdinalType, ScalarType>::print(std::ostream& os) const
975 {
976  os << std::endl;
977  if(valuesCopied_)
978  os << "Values_copied : yes" << std::endl;
979  else
980  os << "Values_copied : no" << std::endl;
981  os << "Rows : " << numRows_ << std::endl;
982  os << "Columns : " << numCols_ << std::endl;
983  os << "LDA : " << stride_ << std::endl;
984  if(numRows_ == 0 || numCols_ == 0) {
985  os << "(matrix is empty, no values to display)" << std::endl;
986  } else {
987  for(OrdinalType i = 0; i < numRows_; i++) {
988  for(OrdinalType j = 0; j < numCols_; j++){
989  os << (*this)(i,j) << " ";
990  }
991  os << std::endl;
992  }
993  }
994  return os;
995 }
996 
997 //----------------------------------------------------------------------------------------------------
998 // Protected methods
999 //----------------------------------------------------------------------------------------------------
1000 
1001 template<typename OrdinalType, typename ScalarType>
1002 inline void SerialDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
1003  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRows_, std::out_of_range,
1004  "SerialDenseMatrix<T>::checkIndex: "
1005  "Row index " << rowIndex << " out of range [0, "<< numRows_ << ")");
1006  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numCols_, std::out_of_range,
1007  "SerialDenseMatrix<T>::checkIndex: "
1008  "Col index " << colIndex << " out of range [0, "<< numCols_ << ")");
1009 }
1010 
1011 template<typename OrdinalType, typename ScalarType>
1012 void SerialDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
1013 {
1014  if (valuesCopied_)
1015  {
1016  delete [] values_;
1017  values_ = 0;
1018  valuesCopied_ = false;
1019  }
1020 }
1021 
1022 template<typename OrdinalType, typename ScalarType>
1023 void SerialDenseMatrix<OrdinalType, ScalarType>::copyMat(
1024  ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
1025  OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput,
1026  OrdinalType startRow, OrdinalType startCol, ScalarType alpha
1027  )
1028 {
1029  OrdinalType i, j;
1030  ScalarType* ptr1 = 0;
1031  ScalarType* ptr2 = 0;
1032  for(j = 0; j < numCols_in; j++) {
1033  ptr1 = outputMatrix + (j * strideOutput);
1034  ptr2 = inputMatrix + (j + startCol) * strideInput + startRow;
1035  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1036  for(i = 0; i < numRows_in; i++)
1037  {
1038  *ptr1++ += alpha*(*ptr2++);
1039  }
1040  } else {
1041  for(i = 0; i < numRows_in; i++)
1042  {
1043  *ptr1++ = *ptr2++;
1044  }
1045  }
1046  }
1047 }
1048 
1050 template<typename OrdinalType, typename ScalarType>
1052 public:
1056  : obj(obj_in) {}
1057 };
1058 
1060 template<typename OrdinalType, typename ScalarType>
1061 std::ostream&
1062 operator<<(std::ostream &out,
1064 {
1065  printer.obj.print(out);
1066  return out;
1067 }
1068 
1070 template<typename OrdinalType, typename ScalarType>
1071 SerialDenseMatrixPrinter<OrdinalType,ScalarType>
1073 {
1075 }
1076 
1077 
1078 } // namespace Teuchos
1079 
1080 
1081 #endif /* _TEUCHOS_SERIALDENSEMATRIX_HPP_ */
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
ScalarType * values() const
Data array access method.
bool upper() const
Returns true if upper triangular part of this matrix has and will be used.
bool operator==(const SerialDenseMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
int shapeUninitialized(OrdinalType numRows, OrdinalType numCols)
Same as shape() except leaves uninitialized.
Templated interface class to BLAS routines.
SerialDenseMatrix< OrdinalType, ScalarType > & operator-=(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Subtract another matrix from this matrix.
SerialBandDenseMatrixPrinter< OrdinalType, ScalarType > printMat(const SerialBandDenseMatrix< OrdinalType, ScalarType > &obj)
Return SerialBandDenseMatrix ostream manipulator Use as:
SerialDenseMatrix< OrdinalType, ScalarType > & operator*=(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
void swap(SerialDenseMatrix< OrdinalType, ScalarType > &B)
Swap values between this matrix and incoming matrix.
ScalarType * values() const
Returns the pointer to the ScalarType data array contained in the object.
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
Multiply A * B and add them to this; this = beta * this + alpha*A*B.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
virtual std::ostream & print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream &lt;&lt; operator.
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
Templated BLAS wrapper.
This class creates and provides basic support for symmetric, positive-definite dense matrices of temp...
int scale(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
OrdinalType numRows() const
Returns the row dimension of this matrix.
Object for storing data and providing functionality that is common to all computational classes...
SerialDenseMatrix()=default
Default Constructor.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
This structure defines some basic traits for a scalar field type.
SerialDenseMatrix< OrdinalType, ScalarType > & operator+=(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Add another matrix to this matrix.
SerialDenseMatrix< OrdinalType, ScalarType > & operator=(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
Functionality and data that is common to all computational classes.
static T conjugate(T a)
Returns the conjugate of the scalar type a.
bool empty() const
Returns whether this matrix is empty.
ScalarType scalarType
Typedef for scalar type.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Set all values in the matrix to a constant value.
bool operator!=(const SerialDenseMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
Templated serial, dense, symmetric matrix class.
OrdinalType ordinalType
Typedef for ordinal type.
ScalarType * operator[](OrdinalType colIndex)
Column access method (non-const).
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
OrdinalType numCols() const
Returns the column dimension of this matrix.
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
int random()
Set all values in the matrix to be random numbers.
Ostream manipulator for SerialDenseMatrix.
int reshape(OrdinalType numRows, OrdinalType numCols)
Reshaping method for changing the size of a SerialDenseMatrix, keeping the entries.
Defines basic traits for the scalar field type.
static T zero()
Returns representation of zero for this scalar type.
static T random()
Returns a random number (between -one() and +one()) of this scalar type.
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
int shape(OrdinalType numRows, OrdinalType numCols)
Shape method for changing the size of a SerialDenseMatrix, initializing entries to zero...
OrdinalType numCols() const
Returns the column dimension of this matrix.
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
Teuchos::DataAccess Mode enumerable type.
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
OrdinalType numRows() const
Returns the row dimension of this matrix.
This class creates and provides basic support for dense rectangular matrix of templated type...