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 //
4 // Teuchos: Common Tools Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // ***********************************************************************
38 // @HEADER
39 
40 #ifndef _TEUCHOS_SERIALDENSEMATRIX_HPP_
41 #define _TEUCHOS_SERIALDENSEMATRIX_HPP_
42 
46 #include "Teuchos_CompObject.hpp"
47 #include "Teuchos_BLAS.hpp"
48 #include "Teuchos_ScalarTraits.hpp"
49 #include "Teuchos_DataAccess.hpp"
50 #include "Teuchos_ConfigDefs.hpp"
51 #include "Teuchos_Assert.hpp"
53 #include <cstddef>
54 #include <utility>
55 
64 namespace Teuchos {
65 
66 template<typename OrdinalType, typename ScalarType>
67 class SerialDenseMatrix : public CompObject, public BLAS<OrdinalType, ScalarType>
68 {
69 public:
70 
72  typedef OrdinalType ordinalType;
74  typedef ScalarType scalarType;
75 
77 
78 
80 
83  SerialDenseMatrix() = default;
84 
86 
94  SerialDenseMatrix(OrdinalType numRows, OrdinalType numCols, bool zeroOut = true);
95 
97 
105  SerialDenseMatrix(DataAccess CV, ScalarType* values, OrdinalType stride, OrdinalType numRows, OrdinalType numCols);
106 
108 
114 
116 
119 
121 
133  SerialDenseMatrix(DataAccess CV, const SerialDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRows, OrdinalType numCols, OrdinalType startRow=0, OrdinalType startCol=0);
134 
136  virtual ~SerialDenseMatrix();
138 
140 
141 
152  int shape(OrdinalType numRows, OrdinalType numCols);
153 
155  int shapeUninitialized(OrdinalType numRows, OrdinalType numCols);
156 
158 
168  int reshape(OrdinalType numRows, OrdinalType numCols);
169 
170 
172 
174 
175 
177 
184 
186 
192 
194 
197  SerialDenseMatrix<OrdinalType, ScalarType>& operator= (const ScalarType value) { putScalar(value); return(*this); }
198 
200 
204  int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero() );
205 
207 
211 
213  int random();
214 
216 
218 
219 
221 
228  ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
229 
231 
238  const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
239 
241 
246  ScalarType* operator [] (OrdinalType colIndex);
247 
249 
254  const ScalarType* operator [] (OrdinalType colIndex) const;
255 
257 
258  ScalarType* values() const { return values_; }
259 
261 
263 
264 
266 
270 
272 
276 
278 
282 
284 
288  int scale ( const ScalarType alpha );
289 
291 
298 
300 
314  int multiply (ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
315 
317 
328  int multiply (ESide sideA, ScalarType alpha, const SerialSymDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
329 
331 
333 
334 
336 
339  bool operator== (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand) const;
340 
342 
345  bool operator!= (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand) const;
346 
348 
350 
351 
353  OrdinalType numRows() const { return(numRows_); }
354 
356  OrdinalType numCols() const { return(numCols_); }
357 
359  OrdinalType stride() const { return(stride_); }
360 
362  bool empty() const { return(numRows_ == 0 || numCols_ == 0); }
364 
366 
367 
370 
373 
377 
379 
380  virtual std::ostream& print(std::ostream& os) const;
382 
384 protected:
385  void copyMat(ScalarType* inputMatrix, OrdinalType strideInput,
386  OrdinalType numRows, OrdinalType numCols, ScalarType* outputMatrix,
387  OrdinalType strideOutput, OrdinalType startRow, OrdinalType startCol,
388  ScalarType alpha = ScalarTraits<ScalarType>::zero() );
389  void deleteArrays();
390  void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
391 
392  static ScalarType*
393  allocateValues(const OrdinalType numRows,
394  const OrdinalType numCols)
395  {
396  const size_t size = size_t(numRows) * size_t(numCols);
397 #pragma GCC diagnostic push
398 #pragma GCC diagnostic ignored "-Wvla"
399  return new ScalarType[size];
400 #pragma GCC diagnostic pop
401  }
402 
403  OrdinalType numRows_ = 0;
404  OrdinalType numCols_ = 0;
405  OrdinalType stride_ = 0;
406  bool valuesCopied_ = false;
407  ScalarType* values_ = nullptr;
408 }; // class Teuchos_SerialDenseMatrix
409 
410 //----------------------------------------------------------------------------------------------------
411 // Constructors and Destructor
412 //----------------------------------------------------------------------------------------------------
413 
414 template<typename OrdinalType, typename ScalarType>
416  OrdinalType numRows_in, OrdinalType numCols_in, bool zeroOut
417  )
418  : numRows_(numRows_in),
419  numCols_(numCols_in),
420  stride_(numRows_in),
421  valuesCopied_(true),
422  values_(allocateValues(numRows_in, numCols_in))
423 {
424  if (zeroOut) {
425  putScalar();
426  }
427 }
428 
429 template<typename OrdinalType, typename ScalarType>
431  DataAccess CV, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRows_in,
432  OrdinalType numCols_in
433  )
434  : numRows_(numRows_in),
435  numCols_(numCols_in),
436  stride_(stride_in),
437  valuesCopied_(false),
438  values_(values_in)
439 {
440  if(CV == Copy)
441  {
442  stride_ = numRows_;
443  values_ = allocateValues(stride_, numCols_);
444  copyMat(values_in, stride_in, numRows_, numCols_, values_, stride_, 0, 0);
445  valuesCopied_ = true;
446  }
447 }
448 
449 template<typename OrdinalType, typename ScalarType>
451  : valuesCopied_(true)
452 {
453  if ( trans == Teuchos::NO_TRANS )
454  {
455  numRows_ = Source.numRows_;
456  numCols_ = Source.numCols_;
457 
458  if (!Source.valuesCopied_)
459  {
460  stride_ = Source.stride_;
461  values_ = Source.values_;
462  valuesCopied_ = false;
463  }
464  else
465  {
466  stride_ = numRows_;
467  if(stride_ > 0 && numCols_ > 0) {
468  values_ = allocateValues(stride_, numCols_);
469  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
470  }
471  else {
472  numRows_ = 0; numCols_ = 0; stride_ = 0;
473  valuesCopied_ = false;
474  }
475  }
476  }
478  {
479  numRows_ = Source.numCols_;
480  numCols_ = Source.numRows_;
481  stride_ = numRows_;
482  values_ = allocateValues(stride_, numCols_);
483  for (OrdinalType j=0; j<numCols_; j++) {
484  for (OrdinalType i=0; i<numRows_; i++) {
485  values_[j*stride_ + i] = Teuchos::ScalarTraits<ScalarType>::conjugate(Source.values_[i*Source.stride_ + j]);
486  }
487  }
488  }
489  else
490  {
491  numRows_ = Source.numCols_;
492  numCols_ = Source.numRows_;
493  stride_ = numRows_;
494  values_ = allocateValues(stride_, numCols_);
495  for (OrdinalType j=0; j<numCols_; j++) {
496  for (OrdinalType i=0; i<numRows_; i++) {
497  values_[j*stride_ + i] = Source.values_[i*Source.stride_ + j];
498  }
499  }
500  }
501 }
502 
503 
504 template<typename OrdinalType, typename ScalarType>
507  )
508  : numRows_(Source.numRows_), numCols_(Source.numCols_), stride_(Source.stride_),
509  valuesCopied_(false), values_(Source.values_)
510 {
511  if(CV == Copy)
512  {
513  stride_ = numRows_;
514  values_ = allocateValues(stride_, numCols_);
515  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
516  valuesCopied_ = true;
517  }
518 }
519 
520 
521 template<typename OrdinalType, typename ScalarType>
524  OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType startRow,
525  OrdinalType startCol
526  )
527  : CompObject(), numRows_(numRows_in), numCols_(numCols_in), stride_(Source.stride_),
528  valuesCopied_(false), values_(Source.values_)
529 {
530  if(CV == Copy)
531  {
532  stride_ = numRows_in;
533  values_ = allocateValues(stride_, numCols_in);
534  copyMat(Source.values_, Source.stride_, numRows_in, numCols_in, values_, stride_, startRow, startCol);
535  valuesCopied_ = true;
536  }
537  else // CV == View
538  {
539  values_ = values_ + (stride_ * startCol) + startRow;
540  }
541 }
542 
543 template<typename OrdinalType, typename ScalarType>
545 {
546  deleteArrays();
547 }
548 
549 //----------------------------------------------------------------------------------------------------
550 // Shape methods
551 //----------------------------------------------------------------------------------------------------
552 
553 template<typename OrdinalType, typename ScalarType>
555  OrdinalType numRows_in, OrdinalType numCols_in
556  )
557 {
558  deleteArrays(); // Get rid of anything that might be already allocated
559  numRows_ = numRows_in;
560  numCols_ = numCols_in;
561  stride_ = numRows_;
562  values_ = allocateValues(stride_, numCols_);
563  putScalar();
564  valuesCopied_ = true;
565  return(0);
566 }
567 
568 template<typename OrdinalType, typename ScalarType>
570  OrdinalType numRows_in, OrdinalType numCols_in
571  )
572 {
573  deleteArrays(); // Get rid of anything that might be already allocated
574  numRows_ = numRows_in;
575  numCols_ = numCols_in;
576  stride_ = numRows_;
577  values_ = allocateValues(stride_, numCols_);
578  valuesCopied_ = true;
579  return(0);
580 }
581 
582 template<typename OrdinalType, typename ScalarType>
584  OrdinalType numRows_in, OrdinalType numCols_in
585  )
586 {
587  // Allocate space for new matrix
588  ScalarType* values_tmp = allocateValues(numRows_in, numCols_in);
589  ScalarType zero = ScalarTraits<ScalarType>::zero();
590  for(OrdinalType k = 0; k < numRows_in * numCols_in; k++)
591  {
592  values_tmp[k] = zero;
593  }
594  OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
595  OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
596  if(values_ != 0)
597  {
598  copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
599  numRows_in, 0, 0); // Copy principal submatrix of A to new A
600  }
601  deleteArrays(); // Get rid of anything that might be already allocated
602  numRows_ = numRows_in;
603  numCols_ = numCols_in;
604  stride_ = numRows_;
605  values_ = values_tmp; // Set pointer to new A
606  valuesCopied_ = true;
607  return(0);
608 }
609 
610 //----------------------------------------------------------------------------------------------------
611 // Set methods
612 //----------------------------------------------------------------------------------------------------
613 
614 template<typename OrdinalType, typename ScalarType>
616 {
617  // Set each value of the dense matrix to "value".
618  for(OrdinalType j = 0; j < numCols_; j++)
619  {
620  for(OrdinalType i = 0; i < numRows_; i++)
621  {
622  values_[i + j*stride_] = value_in;
623  }
624  }
625  return 0;
626 }
627 
628 template<typename OrdinalType, typename ScalarType> void
631 {
632  // Notes:
633  // > DefaultBLASImpl::SWAP() uses a deep copy. This fn uses a pointer swap.
634  // > this fn covers both Vector and Matrix, such that some care must be
635  // employed to not swap across types (creating a Vector with non-unitary
636  // numCols_)
637  // > Inherited data that is not currently swapped (since inactive/deprecated):
638  // >> Teuchos::CompObject:
639  // Flops *flopCounter_ [Note: all SerialDenseMatrix ctors initialize a
640  // NULL flop-counter using CompObject(), such that any flop increments
641  // that are computed are not accumulated.]
642  // >> Teuchos::Object: (now removed from inheritance list)
643  // static int tracebackMode (no swap for statics)
644  // std::string label_ (has been reported as a cause of memory overhead)
645 
646  std::swap(values_ , B.values_);
647  std::swap(numRows_, B.numRows_);
648  std::swap(numCols_, B.numCols_);
649  std::swap(stride_, B.stride_);
650  std::swap(valuesCopied_, B.valuesCopied_);
651 }
652 
653 template<typename OrdinalType, typename ScalarType>
655 {
656  // Set each value of the dense matrix to a random value.
657  for(OrdinalType j = 0; j < numCols_; j++)
658  {
659  for(OrdinalType i = 0; i < numRows_; i++)
660  {
661  values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
662  }
663  }
664  return 0;
665 }
666 
667 template<typename OrdinalType, typename ScalarType>
671  )
672 {
673  if(this == &Source)
674  return(*this); // Special case of source same as target
675  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
676  return(*this); // Special case of both are views to same data.
677 
678  // If the source is a view then we will return a view, else we will return a copy.
679  if (!Source.valuesCopied_) {
680  if(valuesCopied_) {
681  // Clean up stored data if this was previously a copy.
682  deleteArrays();
683  }
684  numRows_ = Source.numRows_;
685  numCols_ = Source.numCols_;
686  stride_ = Source.stride_;
687  values_ = Source.values_;
688  }
689  else {
690  // If we were a view, we will now be a copy.
691  if(!valuesCopied_) {
692  numRows_ = Source.numRows_;
693  numCols_ = Source.numCols_;
694  stride_ = Source.numRows_;
695  if(stride_ > 0 && numCols_ > 0) {
696  values_ = allocateValues(stride_, numCols_);
697  valuesCopied_ = true;
698  }
699  else {
700  values_ = 0;
701  }
702  }
703  // If we were a copy, we will stay a copy.
704  else {
705  if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) { // we don't need to reallocate
706  numRows_ = Source.numRows_;
707  numCols_ = Source.numCols_;
708  }
709  else { // we need to allocate more space (or less space)
710  deleteArrays();
711  numRows_ = Source.numRows_;
712  numCols_ = Source.numCols_;
713  stride_ = Source.numRows_;
714  if(stride_ > 0 && numCols_ > 0) {
715  values_ = allocateValues(stride_, numCols_);
716  valuesCopied_ = true;
717  }
718  }
719  }
720  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
721  }
722  return(*this);
723 }
724 
725 template<typename OrdinalType, typename ScalarType>
727 {
728  // Check for compatible dimensions
729  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
730  {
731  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
732  }
733  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, ScalarTraits<ScalarType>::one());
734  return(*this);
735 }
736 
737 template<typename OrdinalType, typename ScalarType>
739 {
740  // Check for compatible dimensions
741  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
742  {
743  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
744  }
745  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, -ScalarTraits<ScalarType>::one());
746  return(*this);
747 }
748 
749 template<typename OrdinalType, typename ScalarType>
751  if(this == &Source)
752  return(*this); // Special case of source same as target
753  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
754  return(*this); // Special case of both are views to same data.
755 
756  // Check for compatible dimensions
757  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
758  {
759  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
760  }
761  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
762  return(*this);
763 }
764 
765 //----------------------------------------------------------------------------------------------------
766 // Accessor methods
767 //----------------------------------------------------------------------------------------------------
768 
769 template<typename OrdinalType, typename ScalarType>
770 inline ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
771 {
772 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
773  checkIndex( rowIndex, colIndex );
774 #endif
775  return(values_[colIndex * stride_ + rowIndex]);
776 }
777 
778 template<typename OrdinalType, typename ScalarType>
779 inline const ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
780 {
781 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
782  checkIndex( rowIndex, colIndex );
783 #endif
784  return(values_[colIndex * stride_ + rowIndex]);
785 }
786 
787 template<typename OrdinalType, typename ScalarType>
788 inline const ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex) const
789 {
790 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
791  checkIndex( 0, colIndex );
792 #endif
793  return(values_ + colIndex * stride_);
794 }
795 
796 template<typename OrdinalType, typename ScalarType>
797 inline ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex)
798 {
799 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
800  checkIndex( 0, colIndex );
801 #endif
802  return(values_ + colIndex * stride_);
803 }
804 
805 //----------------------------------------------------------------------------------------------------
806 // Norm methods
807 //----------------------------------------------------------------------------------------------------
808 
809 template<typename OrdinalType, typename ScalarType>
811 {
812  OrdinalType i, j;
815  ScalarType* ptr;
816  for(j = 0; j < numCols_; j++)
817  {
818  ScalarType sum = 0;
819  ptr = values_ + j * stride_;
820  for(i = 0; i < numRows_; i++)
821  {
823  }
825  if(absSum > anorm)
826  {
827  anorm = absSum;
828  }
829  }
830  // don't compute flop increment unless there is an accumulator
831  if (flopCounter_!=0) updateFlops(numRows_ * numCols_);
832  return(anorm);
833 }
834 
835 template<typename OrdinalType, typename ScalarType>
837 {
838  OrdinalType i, j;
840 
841  for (i = 0; i < numRows_; i++) {
843  for (j=0; j< numCols_; j++) {
844  sum += ScalarTraits<ScalarType>::magnitude(*(values_+i+j*stride_));
845  }
846  anorm = TEUCHOS_MAX( anorm, sum );
847  }
848  // don't compute flop increment unless there is an accumulator
849  if (flopCounter_!=0) updateFlops(numRows_ * numCols_);
850  return(anorm);
851 }
852 
853 template<typename OrdinalType, typename ScalarType>
855 {
856  OrdinalType i, j;
858  for (j = 0; j < numCols_; j++) {
859  for (i = 0; i < numRows_; i++) {
860  anorm += ScalarTraits<ScalarType>::magnitude(values_[i+j*stride_]*values_[i+j*stride_]);
861  }
862  }
864  // don't compute flop increment unless there is an accumulator
865  if (flopCounter_!=0) updateFlops(numRows_ * numCols_);
866  return(anorm);
867 }
868 
869 //----------------------------------------------------------------------------------------------------
870 // Comparison methods
871 //----------------------------------------------------------------------------------------------------
872 
873 template<typename OrdinalType, typename ScalarType>
875 {
876  bool result = 1;
877  if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_))
878  {
879  result = 0;
880  }
881  else
882  {
883  OrdinalType i, j;
884  for(i = 0; i < numRows_; i++)
885  {
886  for(j = 0; j < numCols_; j++)
887  {
888  if((*this)(i, j) != Operand(i, j))
889  {
890  return 0;
891  }
892  }
893  }
894  }
895  return result;
896 }
897 
898 template<typename OrdinalType, typename ScalarType>
900 {
901  return !((*this) == Operand);
902 }
903 
904 //----------------------------------------------------------------------------------------------------
905 // Multiplication method
906 //----------------------------------------------------------------------------------------------------
907 
908 template<typename OrdinalType, typename ScalarType>
910 {
911  this->scale( alpha );
912  return(*this);
913 }
914 
915 template<typename OrdinalType, typename ScalarType>
917 {
918  OrdinalType i, j;
919  ScalarType* ptr;
920 
921  for (j=0; j<numCols_; j++) {
922  ptr = values_ + j*stride_;
923  for (i=0; i<numRows_; i++) { *ptr = alpha * (*ptr); ptr++; }
924  }
925  // don't compute flop increment unless there is an accumulator
926  if (flopCounter_!=0) updateFlops( numRows_*numCols_ );
927  return(0);
928 }
929 
930 template<typename OrdinalType, typename ScalarType>
932 {
933  OrdinalType i, j;
934  ScalarType* ptr;
935 
936  // Check for compatible dimensions
937  if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_))
938  {
939  TEUCHOS_CHK_ERR(-1); // Return error
940  }
941  for (j=0; j<numCols_; j++) {
942  ptr = values_ + j*stride_;
943  for (i=0; i<numRows_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
944  }
945  // don't compute flop increment unless there is an accumulator
946  if (flopCounter_!=0) updateFlops( numRows_*numCols_ );
947  return(0);
948 }
949 
950 template<typename OrdinalType, typename ScalarType>
952 {
953  // Check for compatible dimensions
954  OrdinalType A_nrows = (ETranspChar[transa]!='N') ? A.numCols() : A.numRows();
955  OrdinalType A_ncols = (ETranspChar[transa]!='N') ? A.numRows() : A.numCols();
956  OrdinalType B_nrows = (ETranspChar[transb]!='N') ? B.numCols() : B.numRows();
957  OrdinalType B_ncols = (ETranspChar[transb]!='N') ? B.numRows() : B.numCols();
958  if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols))
959  {
960  TEUCHOS_CHK_ERR(-1); // Return error
961  }
962  // Call GEMM function
963  this->GEMM(transa, transb, numRows_, numCols_, A_ncols, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
964 
965  // don't compute flop increment unless there is an accumulator
966  if (flopCounter_!=0) {
967  double nflops = 2 * numRows_;
968  nflops *= numCols_;
969  nflops *= A_ncols;
970  updateFlops(nflops);
971  }
972  return(0);
973 }
974 
975 template<typename OrdinalType, typename ScalarType>
977 {
978  // Check for compatible dimensions
979  OrdinalType A_nrows = A.numRows(), A_ncols = A.numCols();
980  OrdinalType B_nrows = B.numRows(), B_ncols = B.numCols();
981 
982  if (ESideChar[sideA]=='L') {
983  if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols)) {
984  TEUCHOS_CHK_ERR(-1); // Return error
985  }
986  } else {
987  if ((numRows_ != B_nrows) || (B_ncols != A_nrows) || (numCols_ != A_ncols)) {
988  TEUCHOS_CHK_ERR(-1); // Return error
989  }
990  }
991 
992  // Call SYMM function
994  this->SYMM(sideA, uplo, numRows_, numCols_, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
995 
996  // don't compute flop increment unless there is an accumulator
997  if (flopCounter_!=0) {
998  double nflops = 2 * numRows_;
999  nflops *= numCols_;
1000  nflops *= A_ncols;
1001  updateFlops(nflops);
1002  }
1003  return(0);
1004 }
1005 
1006 template<typename OrdinalType, typename ScalarType>
1007 std::ostream& SerialDenseMatrix<OrdinalType, ScalarType>::print(std::ostream& os) const
1008 {
1009  os << std::endl;
1010  if(valuesCopied_)
1011  os << "Values_copied : yes" << std::endl;
1012  else
1013  os << "Values_copied : no" << std::endl;
1014  os << "Rows : " << numRows_ << std::endl;
1015  os << "Columns : " << numCols_ << std::endl;
1016  os << "LDA : " << stride_ << std::endl;
1017  if(numRows_ == 0 || numCols_ == 0) {
1018  os << "(matrix is empty, no values to display)" << std::endl;
1019  } else {
1020  for(OrdinalType i = 0; i < numRows_; i++) {
1021  for(OrdinalType j = 0; j < numCols_; j++){
1022  os << (*this)(i,j) << " ";
1023  }
1024  os << std::endl;
1025  }
1026  }
1027  return os;
1028 }
1029 
1030 //----------------------------------------------------------------------------------------------------
1031 // Protected methods
1032 //----------------------------------------------------------------------------------------------------
1033 
1034 template<typename OrdinalType, typename ScalarType>
1035 inline void SerialDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
1036  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRows_, std::out_of_range,
1037  "SerialDenseMatrix<T>::checkIndex: "
1038  "Row index " << rowIndex << " out of range [0, "<< numRows_ << ")");
1039  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numCols_, std::out_of_range,
1040  "SerialDenseMatrix<T>::checkIndex: "
1041  "Col index " << colIndex << " out of range [0, "<< numCols_ << ")");
1042 }
1043 
1044 template<typename OrdinalType, typename ScalarType>
1045 void SerialDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
1046 {
1047  if (valuesCopied_)
1048  {
1049  delete [] values_;
1050  values_ = 0;
1051  valuesCopied_ = false;
1052  }
1053 }
1054 
1055 template<typename OrdinalType, typename ScalarType>
1056 void SerialDenseMatrix<OrdinalType, ScalarType>::copyMat(
1057  ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
1058  OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput,
1059  OrdinalType startRow, OrdinalType startCol, ScalarType alpha
1060  )
1061 {
1062  OrdinalType i, j;
1063  ScalarType* ptr1 = 0;
1064  ScalarType* ptr2 = 0;
1065  for(j = 0; j < numCols_in; j++) {
1066  ptr1 = outputMatrix + (j * strideOutput);
1067  ptr2 = inputMatrix + (j + startCol) * strideInput + startRow;
1068  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1069  for(i = 0; i < numRows_in; i++)
1070  {
1071  *ptr1++ += alpha*(*ptr2++);
1072  }
1073  } else {
1074  for(i = 0; i < numRows_in; i++)
1075  {
1076  *ptr1++ = *ptr2++;
1077  }
1078  }
1079  }
1080 }
1081 
1083 template<typename OrdinalType, typename ScalarType>
1085 public:
1089  : obj(obj_in) {}
1090 };
1091 
1093 template<typename OrdinalType, typename ScalarType>
1094 std::ostream&
1095 operator<<(std::ostream &out,
1097 {
1098  printer.obj.print(out);
1099  return out;
1100 }
1101 
1103 template<typename OrdinalType, typename ScalarType>
1104 SerialDenseMatrixPrinter<OrdinalType,ScalarType>
1106 {
1108 }
1109 
1110 
1111 } // namespace Teuchos
1112 
1113 
1114 #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...