Teuchos Package Browser (Single Doxygen Collection)  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 #pragma GCC diagnostic push
368 #pragma GCC diagnostic ignored "-Wvla"
369  return new ScalarType[size];
370 #pragma GCC diagnostic pop
371  }
372 
373  OrdinalType numRows_ = 0;
374  OrdinalType numCols_ = 0;
375  OrdinalType stride_ = 0;
376  bool valuesCopied_ = false;
377  ScalarType* values_ = nullptr;
378 }; // class Teuchos_SerialDenseMatrix
379 
380 //----------------------------------------------------------------------------------------------------
381 // Constructors and Destructor
382 //----------------------------------------------------------------------------------------------------
383 
384 template<typename OrdinalType, typename ScalarType>
386  OrdinalType numRows_in, OrdinalType numCols_in, bool zeroOut
387  )
388  : numRows_(numRows_in),
389  numCols_(numCols_in),
390  stride_(numRows_in),
391  valuesCopied_(true),
392  values_(allocateValues(numRows_in, numCols_in))
393 {
394  if (zeroOut) {
395  putScalar();
396  }
397 }
398 
399 template<typename OrdinalType, typename ScalarType>
401  DataAccess CV, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRows_in,
402  OrdinalType numCols_in
403  )
404  : numRows_(numRows_in),
405  numCols_(numCols_in),
406  stride_(stride_in),
407  valuesCopied_(false),
408  values_(values_in)
409 {
410  if(CV == Copy)
411  {
412  stride_ = numRows_;
414  copyMat(values_in, stride_in, numRows_, numCols_, values_, stride_, 0, 0);
415  valuesCopied_ = true;
416  }
417 }
418 
419 template<typename OrdinalType, typename ScalarType>
421  : valuesCopied_(true)
422 {
423  if ( trans == Teuchos::NO_TRANS )
424  {
425  numRows_ = Source.numRows_;
426  numCols_ = Source.numCols_;
427 
428  if (!Source.valuesCopied_)
429  {
430  stride_ = Source.stride_;
431  values_ = Source.values_;
432  valuesCopied_ = false;
433  }
434  else
435  {
436  stride_ = numRows_;
437  if(stride_ > 0 && numCols_ > 0) {
439  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
440  }
441  else {
442  numRows_ = 0; numCols_ = 0; stride_ = 0;
443  valuesCopied_ = false;
444  }
445  }
446  }
448  {
449  numRows_ = Source.numCols_;
450  numCols_ = Source.numRows_;
451  stride_ = numRows_;
453  for (OrdinalType j=0; j<numCols_; j++) {
454  for (OrdinalType i=0; i<numRows_; i++) {
456  }
457  }
458  }
459  else
460  {
461  numRows_ = Source.numCols_;
462  numCols_ = Source.numRows_;
463  stride_ = numRows_;
465  for (OrdinalType j=0; j<numCols_; j++) {
466  for (OrdinalType i=0; i<numRows_; i++) {
467  values_[j*stride_ + i] = Source.values_[i*Source.stride_ + j];
468  }
469  }
470  }
471 }
472 
473 
474 template<typename OrdinalType, typename ScalarType>
477  )
478  : numRows_(Source.numRows_), numCols_(Source.numCols_), stride_(Source.stride_),
479  valuesCopied_(false), values_(Source.values_)
480 {
481  if(CV == Copy)
482  {
483  stride_ = numRows_;
485  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
486  valuesCopied_ = true;
487  }
488 }
489 
490 
491 template<typename OrdinalType, typename ScalarType>
494  OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType startRow,
495  OrdinalType startCol
496  )
497  : CompObject(), numRows_(numRows_in), numCols_(numCols_in), stride_(Source.stride_),
498  valuesCopied_(false), values_(Source.values_)
499 {
500  if(CV == Copy)
501  {
502  stride_ = numRows_in;
503  values_ = allocateValues(stride_, numCols_in);
504  copyMat(Source.values_, Source.stride_, numRows_in, numCols_in, values_, stride_, startRow, startCol);
505  valuesCopied_ = true;
506  }
507  else // CV == View
508  {
509  values_ = values_ + (stride_ * startCol) + startRow;
510  }
511 }
512 
513 template<typename OrdinalType, typename ScalarType>
515 {
516  deleteArrays();
517 }
518 
519 //----------------------------------------------------------------------------------------------------
520 // Shape methods
521 //----------------------------------------------------------------------------------------------------
522 
523 template<typename OrdinalType, typename ScalarType>
525  OrdinalType numRows_in, OrdinalType numCols_in
526  )
527 {
528  deleteArrays(); // Get rid of anything that might be already allocated
529  numRows_ = numRows_in;
530  numCols_ = numCols_in;
531  stride_ = numRows_;
532  values_ = allocateValues(stride_, numCols_);
533  putScalar();
534  valuesCopied_ = true;
535  return(0);
536 }
537 
538 template<typename OrdinalType, typename ScalarType>
540  OrdinalType numRows_in, OrdinalType numCols_in
541  )
542 {
543  deleteArrays(); // Get rid of anything that might be already allocated
544  numRows_ = numRows_in;
545  numCols_ = numCols_in;
546  stride_ = numRows_;
547  values_ = allocateValues(stride_, numCols_);
548  valuesCopied_ = true;
549  return(0);
550 }
551 
552 template<typename OrdinalType, typename ScalarType>
554  OrdinalType numRows_in, OrdinalType numCols_in
555  )
556 {
557  // Allocate space for new matrix
558  ScalarType* values_tmp = allocateValues(numRows_in, numCols_in);
559  ScalarType zero = ScalarTraits<ScalarType>::zero();
560  for(OrdinalType k = 0; k < numRows_in * numCols_in; k++)
561  {
562  values_tmp[k] = zero;
563  }
564  OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
565  OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
566  if(values_ != 0)
567  {
568  copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
569  numRows_in, 0, 0); // Copy principal submatrix of A to new A
570  }
571  deleteArrays(); // Get rid of anything that might be already allocated
572  numRows_ = numRows_in;
573  numCols_ = numCols_in;
574  stride_ = numRows_;
575  values_ = values_tmp; // Set pointer to new A
576  valuesCopied_ = true;
577  return(0);
578 }
579 
580 //----------------------------------------------------------------------------------------------------
581 // Set methods
582 //----------------------------------------------------------------------------------------------------
583 
584 template<typename OrdinalType, typename ScalarType>
586 {
587  // Set each value of the dense matrix to "value".
588  for(OrdinalType j = 0; j < numCols_; j++)
589  {
590  for(OrdinalType i = 0; i < numRows_; i++)
591  {
592  values_[i + j*stride_] = value_in;
593  }
594  }
595  return 0;
596 }
597 
598 template<typename OrdinalType, typename ScalarType> void
601 {
602  // Notes:
603  // > DefaultBLASImpl::SWAP() uses a deep copy. This fn uses a pointer swap.
604  // > this fn covers both Vector and Matrix, such that some care must be
605  // employed to not swap across types (creating a Vector with non-unitary
606  // numCols_)
607  // > Inherited data that is not currently swapped (since inactive/deprecated):
608  // >> Teuchos::CompObject:
609  // Flops *flopCounter_ [Note: all SerialDenseMatrix ctors initialize a
610  // NULL flop-counter using CompObject(), such that any flop increments
611  // that are computed are not accumulated.]
612  // >> Teuchos::Object: (now removed from inheritance list)
613  // static int tracebackMode (no swap for statics)
614  // std::string label_ (has been reported as a cause of memory overhead)
615 
616  std::swap(values_ , B.values_);
617  std::swap(numRows_, B.numRows_);
618  std::swap(numCols_, B.numCols_);
619  std::swap(stride_, B.stride_);
620  std::swap(valuesCopied_, B.valuesCopied_);
621 }
622 
623 template<typename OrdinalType, typename ScalarType>
625 {
626  // Set each value of the dense matrix to a random value.
627  for(OrdinalType j = 0; j < numCols_; j++)
628  {
629  for(OrdinalType i = 0; i < numRows_; i++)
630  {
631  values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
632  }
633  }
634  return 0;
635 }
636 
637 template<typename OrdinalType, typename ScalarType>
641  )
642 {
643  if(this == &Source)
644  return(*this); // Special case of source same as target
645  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
646  return(*this); // Special case of both are views to same data.
647 
648  // If the source is a view then we will return a view, else we will return a copy.
649  if (!Source.valuesCopied_) {
650  if(valuesCopied_) {
651  // Clean up stored data if this was previously a copy.
652  deleteArrays();
653  }
654  numRows_ = Source.numRows_;
655  numCols_ = Source.numCols_;
656  stride_ = Source.stride_;
657  values_ = Source.values_;
658  }
659  else {
660  // If we were a view, we will now be a copy.
661  if(!valuesCopied_) {
662  numRows_ = Source.numRows_;
663  numCols_ = Source.numCols_;
664  stride_ = Source.numRows_;
665  if(stride_ > 0 && numCols_ > 0) {
666  values_ = allocateValues(stride_, numCols_);
667  valuesCopied_ = true;
668  }
669  else {
670  values_ = 0;
671  }
672  }
673  // If we were a copy, we will stay a copy.
674  else {
675  if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) { // we don't need to reallocate
676  numRows_ = Source.numRows_;
677  numCols_ = Source.numCols_;
678  }
679  else { // we need to allocate more space (or less space)
680  deleteArrays();
681  numRows_ = Source.numRows_;
682  numCols_ = Source.numCols_;
683  stride_ = Source.numRows_;
684  if(stride_ > 0 && numCols_ > 0) {
685  values_ = allocateValues(stride_, numCols_);
686  valuesCopied_ = true;
687  }
688  }
689  }
690  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
691  }
692  return(*this);
693 }
694 
695 template<typename OrdinalType, typename ScalarType>
697 {
698  // Check for compatible dimensions
699  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
700  {
701  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
702  }
703  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, ScalarTraits<ScalarType>::one());
704  return(*this);
705 }
706 
707 template<typename OrdinalType, typename ScalarType>
709 {
710  // Check for compatible dimensions
711  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
712  {
713  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
714  }
715  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, -ScalarTraits<ScalarType>::one());
716  return(*this);
717 }
718 
719 template<typename OrdinalType, typename ScalarType>
721  if(this == &Source)
722  return(*this); // Special case of source same as target
723  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
724  return(*this); // Special case of both are views to same data.
725 
726  // Check for compatible dimensions
727  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
728  {
729  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
730  }
731  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
732  return(*this);
733 }
734 
735 //----------------------------------------------------------------------------------------------------
736 // Accessor methods
737 //----------------------------------------------------------------------------------------------------
738 
739 template<typename OrdinalType, typename ScalarType>
740 inline ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
741 {
742 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
743  checkIndex( rowIndex, colIndex );
744 #endif
745  return(values_[colIndex * stride_ + rowIndex]);
746 }
747 
748 template<typename OrdinalType, typename ScalarType>
749 inline const ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
750 {
751 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
752  checkIndex( rowIndex, colIndex );
753 #endif
754  return(values_[colIndex * stride_ + rowIndex]);
755 }
756 
757 template<typename OrdinalType, typename ScalarType>
758 inline const ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex) const
759 {
760 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
761  checkIndex( 0, colIndex );
762 #endif
763  return(values_ + colIndex * stride_);
764 }
765 
766 template<typename OrdinalType, typename ScalarType>
767 inline ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex)
768 {
769 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
770  checkIndex( 0, colIndex );
771 #endif
772  return(values_ + colIndex * stride_);
773 }
774 
775 //----------------------------------------------------------------------------------------------------
776 // Norm methods
777 //----------------------------------------------------------------------------------------------------
778 
779 template<typename OrdinalType, typename ScalarType>
781 {
782  OrdinalType i, j;
785  ScalarType* ptr;
786  for(j = 0; j < numCols_; j++)
787  {
788  ScalarType sum = 0;
789  ptr = values_ + j * stride_;
790  for(i = 0; i < numRows_; i++)
791  {
793  }
795  if(absSum > anorm)
796  {
797  anorm = absSum;
798  }
799  }
800  // don't compute flop increment unless there is an accumulator
801  if (flopCounter_!=0) updateFlops(numRows_ * numCols_);
802  return(anorm);
803 }
804 
805 template<typename OrdinalType, typename ScalarType>
807 {
808  OrdinalType i, j;
810 
811  for (i = 0; i < numRows_; i++) {
813  for (j=0; j< numCols_; j++) {
814  sum += ScalarTraits<ScalarType>::magnitude(*(values_+i+j*stride_));
815  }
816  anorm = TEUCHOS_MAX( anorm, sum );
817  }
818  // don't compute flop increment unless there is an accumulator
819  if (flopCounter_!=0) updateFlops(numRows_ * numCols_);
820  return(anorm);
821 }
822 
823 template<typename OrdinalType, typename ScalarType>
825 {
826  OrdinalType i, j;
828  for (j = 0; j < numCols_; j++) {
829  for (i = 0; i < numRows_; i++) {
830  anorm += ScalarTraits<ScalarType>::magnitude(values_[i+j*stride_]*values_[i+j*stride_]);
831  }
832  }
834  // don't compute flop increment unless there is an accumulator
835  if (flopCounter_!=0) updateFlops(numRows_ * numCols_);
836  return(anorm);
837 }
838 
839 //----------------------------------------------------------------------------------------------------
840 // Comparison methods
841 //----------------------------------------------------------------------------------------------------
842 
843 template<typename OrdinalType, typename ScalarType>
845 {
846  bool result = 1;
847  if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_))
848  {
849  result = 0;
850  }
851  else
852  {
853  OrdinalType i, j;
854  for(i = 0; i < numRows_; i++)
855  {
856  for(j = 0; j < numCols_; j++)
857  {
858  if((*this)(i, j) != Operand(i, j))
859  {
860  return 0;
861  }
862  }
863  }
864  }
865  return result;
866 }
867 
868 template<typename OrdinalType, typename ScalarType>
870 {
871  return !((*this) == Operand);
872 }
873 
874 //----------------------------------------------------------------------------------------------------
875 // Multiplication method
876 //----------------------------------------------------------------------------------------------------
877 
878 template<typename OrdinalType, typename ScalarType>
880 {
881  this->scale( alpha );
882  return(*this);
883 }
884 
885 template<typename OrdinalType, typename ScalarType>
887 {
888  OrdinalType i, j;
889  ScalarType* ptr;
890 
891  for (j=0; j<numCols_; j++) {
892  ptr = values_ + j*stride_;
893  for (i=0; i<numRows_; i++) { *ptr = alpha * (*ptr); ptr++; }
894  }
895  // don't compute flop increment unless there is an accumulator
896  if (flopCounter_!=0) updateFlops( numRows_*numCols_ );
897  return(0);
898 }
899 
900 template<typename OrdinalType, typename ScalarType>
902 {
903  OrdinalType i, j;
904  ScalarType* ptr;
905 
906  // Check for compatible dimensions
907  if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_))
908  {
909  TEUCHOS_CHK_ERR(-1); // Return error
910  }
911  for (j=0; j<numCols_; j++) {
912  ptr = values_ + j*stride_;
913  for (i=0; i<numRows_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
914  }
915  // don't compute flop increment unless there is an accumulator
916  if (flopCounter_!=0) updateFlops( numRows_*numCols_ );
917  return(0);
918 }
919 
920 template<typename OrdinalType, typename ScalarType>
922 {
923  // Check for compatible dimensions
924  OrdinalType A_nrows = (ETranspChar[transa]!='N') ? A.numCols() : A.numRows();
925  OrdinalType A_ncols = (ETranspChar[transa]!='N') ? A.numRows() : A.numCols();
926  OrdinalType B_nrows = (ETranspChar[transb]!='N') ? B.numCols() : B.numRows();
927  OrdinalType B_ncols = (ETranspChar[transb]!='N') ? B.numRows() : B.numCols();
928  if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols))
929  {
930  TEUCHOS_CHK_ERR(-1); // Return error
931  }
932  // Call GEMM function
933  this->GEMM(transa, transb, numRows_, numCols_, A_ncols, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
934 
935  // don't compute flop increment unless there is an accumulator
936  if (flopCounter_!=0) {
937  double nflops = 2 * numRows_;
938  nflops *= numCols_;
939  nflops *= A_ncols;
940  updateFlops(nflops);
941  }
942  return(0);
943 }
944 
945 template<typename OrdinalType, typename ScalarType>
947 {
948  // Check for compatible dimensions
949  OrdinalType A_nrows = A.numRows(), A_ncols = A.numCols();
950  OrdinalType B_nrows = B.numRows(), B_ncols = B.numCols();
951 
952  if (ESideChar[sideA]=='L') {
953  if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols)) {
954  TEUCHOS_CHK_ERR(-1); // Return error
955  }
956  } else {
957  if ((numRows_ != B_nrows) || (B_ncols != A_nrows) || (numCols_ != A_ncols)) {
958  TEUCHOS_CHK_ERR(-1); // Return error
959  }
960  }
961 
962  // Call SYMM function
964  this->SYMM(sideA, uplo, numRows_, numCols_, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
965 
966  // don't compute flop increment unless there is an accumulator
967  if (flopCounter_!=0) {
968  double nflops = 2 * numRows_;
969  nflops *= numCols_;
970  nflops *= A_ncols;
971  updateFlops(nflops);
972  }
973  return(0);
974 }
975 
976 template<typename OrdinalType, typename ScalarType>
977 std::ostream& SerialDenseMatrix<OrdinalType, ScalarType>::print(std::ostream& os) const
978 {
979  os << std::endl;
980  if(valuesCopied_)
981  os << "Values_copied : yes" << std::endl;
982  else
983  os << "Values_copied : no" << std::endl;
984  os << "Rows : " << numRows_ << std::endl;
985  os << "Columns : " << numCols_ << std::endl;
986  os << "LDA : " << stride_ << std::endl;
987  if(numRows_ == 0 || numCols_ == 0) {
988  os << "(matrix is empty, no values to display)" << std::endl;
989  } else {
990  for(OrdinalType i = 0; i < numRows_; i++) {
991  for(OrdinalType j = 0; j < numCols_; j++){
992  os << (*this)(i,j) << " ";
993  }
994  os << std::endl;
995  }
996  }
997  return os;
998 }
999 
1000 //----------------------------------------------------------------------------------------------------
1001 // Protected methods
1002 //----------------------------------------------------------------------------------------------------
1003 
1004 template<typename OrdinalType, typename ScalarType>
1005 inline void SerialDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
1006  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRows_, std::out_of_range,
1007  "SerialDenseMatrix<T>::checkIndex: "
1008  "Row index " << rowIndex << " out of range [0, "<< numRows_ << ")");
1009  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numCols_, std::out_of_range,
1010  "SerialDenseMatrix<T>::checkIndex: "
1011  "Col index " << colIndex << " out of range [0, "<< numCols_ << ")");
1012 }
1013 
1014 template<typename OrdinalType, typename ScalarType>
1016 {
1017  if (valuesCopied_)
1018  {
1019  delete [] values_;
1020  values_ = 0;
1021  valuesCopied_ = false;
1022  }
1023 }
1024 
1025 template<typename OrdinalType, typename ScalarType>
1027  ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
1028  OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput,
1029  OrdinalType startRow, OrdinalType startCol, ScalarType alpha
1030  )
1031 {
1032  OrdinalType i, j;
1033  ScalarType* ptr1 = 0;
1034  ScalarType* ptr2 = 0;
1035  for(j = 0; j < numCols_in; j++) {
1036  ptr1 = outputMatrix + (j * strideOutput);
1037  ptr2 = inputMatrix + (j + startCol) * strideInput + startRow;
1038  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1039  for(i = 0; i < numRows_in; i++)
1040  {
1041  *ptr1++ += alpha*(*ptr2++);
1042  }
1043  } else {
1044  for(i = 0; i < numRows_in; i++)
1045  {
1046  *ptr1++ = *ptr2++;
1047  }
1048  }
1049  }
1050 }
1051 
1053 template<typename OrdinalType, typename ScalarType>
1055 public:
1059  : obj(obj_in) {}
1060 };
1061 
1063 template<typename OrdinalType, typename ScalarType>
1064 std::ostream&
1065 operator<<(std::ostream &out,
1067 {
1068  printer.obj.print(out);
1069  return out;
1070 }
1071 
1073 template<typename OrdinalType, typename ScalarType>
1074 SerialDenseMatrixPrinter<OrdinalType,ScalarType>
1076 {
1078 }
1079 
1080 
1081 } // namespace Teuchos
1082 
1083 
1084 #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.
void checkIndex(OrdinalType rowIndex, OrdinalType colIndex=0) const
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...
SerialDenseMatrixPrinter(const SerialDenseMatrix< OrdinalType, ScalarType > &obj_in)
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.
#define TEUCHOS_CHK_ERR(a)
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.
const SerialDenseMatrix< OrdinalType, ScalarType > & obj
bool empty() const
Returns whether this matrix is empty.
std::ostream & operator<<(std::ostream &os, BigUInt< n > a)
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.
#define TEUCHOS_MAX(x, y)
bool operator!=(const SerialDenseMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ESideChar[]
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
void copyMat(ScalarType *inputMatrix, OrdinalType strideInput, OrdinalType numRows, OrdinalType numCols, ScalarType *outputMatrix, OrdinalType strideOutput, OrdinalType startRow, OrdinalType startCol, ScalarType alpha=ScalarTraits< ScalarType >::zero())
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.
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ETranspChar[]
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.
#define TEUCHOS_CHK_REF(a)
#define TEUCHOS_MIN(x, y)
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
static ScalarType * allocateValues(const OrdinalType numRows, const OrdinalType numCols)
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...