Teuchos - Trilinos Tools Package  Version of the Day
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 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #ifndef _TEUCHOS_SERIALDENSEMATRIX_HPP_
43 #define _TEUCHOS_SERIALDENSEMATRIX_HPP_
44 
48 #include "Teuchos_CompObject.hpp"
49 #include "Teuchos_BLAS.hpp"
50 #include "Teuchos_ScalarTraits.hpp"
51 #include "Teuchos_DataAccess.hpp"
52 #include "Teuchos_ConfigDefs.hpp"
53 #include "Teuchos_Assert.hpp"
55 
64 namespace Teuchos {
65 
66 template<typename OrdinalType, typename ScalarType>
67 class SerialDenseMatrix : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>
68 {
69 public:
70 
72  typedef OrdinalType ordinalType;
74  typedef ScalarType scalarType;
75 
77 
78 
80 
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  int random();
208 
210 
212 
213 
215 
222  ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
223 
225 
232  const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
233 
235 
242  ScalarType* operator [] (OrdinalType colIndex);
243 
245 
252  const ScalarType* operator [] (OrdinalType colIndex) const;
253 
255 
256  ScalarType* values() const { return(values_); }
257 
259 
261 
262 
264 
268 
270 
274 
276 
280 
282 
286  int scale ( const ScalarType alpha );
287 
289 
296 
298 
312  int multiply (ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
313 
315 
326  int multiply (ESide sideA, ScalarType alpha, const SerialSymDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
327 
329 
331 
332 
334 
337  bool operator== (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand) const;
338 
340 
343  bool operator!= (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand) const;
344 
346 
348 
349 
351  OrdinalType numRows() const { return(numRows_); }
352 
354  OrdinalType numCols() const { return(numCols_); }
355 
357  OrdinalType stride() const { return(stride_); }
358 
360  bool empty() const { return(numRows_ == 0 || numCols_ == 0); }
362 
364 
365 
368 
371 
375 
377 
378  virtual void print(std::ostream& os) const;
380 
382 protected:
383  void copyMat(ScalarType* inputMatrix, OrdinalType strideInput,
384  OrdinalType numRows, OrdinalType numCols, ScalarType* outputMatrix,
385  OrdinalType strideOutput, OrdinalType startRow, OrdinalType startCol,
386  ScalarType alpha = ScalarTraits<ScalarType>::zero() );
387  void deleteArrays();
388  void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
389  OrdinalType numRows_;
390  OrdinalType numCols_;
391  OrdinalType stride_;
392  bool valuesCopied_;
393  ScalarType* values_;
394 
395 }; // class Teuchos_SerialDenseMatrix
396 
397 //----------------------------------------------------------------------------------------------------
398 // Constructors and Destructor
399 //----------------------------------------------------------------------------------------------------
400 
401 template<typename OrdinalType, typename ScalarType>
403  : CompObject(), Object(), BLAS<OrdinalType,ScalarType>(), numRows_(0), numCols_(0), stride_(0), valuesCopied_(false), values_(0)
404 {}
405 
406 template<typename OrdinalType, typename ScalarType>
408  OrdinalType numRows_in, OrdinalType numCols_in, bool zeroOut
409  )
410  : CompObject(), Object(), BLAS<OrdinalType,ScalarType>(), numRows_(numRows_in), numCols_(numCols_in), stride_(numRows_in)
411 {
412  values_ = new ScalarType[stride_*numCols_];
413  valuesCopied_ = true;
414  if (zeroOut == true)
415  putScalar();
416 }
417 
418 template<typename OrdinalType, typename ScalarType>
420  DataAccess CV, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRows_in,
421  OrdinalType numCols_in
422  )
423  : CompObject(), Object(), BLAS<OrdinalType,ScalarType>(), numRows_(numRows_in), numCols_(numCols_in), stride_(stride_in),
424  valuesCopied_(false), values_(values_in)
425 {
426  if(CV == Copy)
427  {
428  stride_ = numRows_;
429  values_ = new ScalarType[stride_*numCols_];
430  copyMat(values_in, stride_in, numRows_, numCols_, values_, stride_, 0, 0);
431  valuesCopied_ = true;
432  }
433 }
434 
435 template<typename OrdinalType, typename ScalarType>
437  : CompObject(), Object(), BLAS<OrdinalType,ScalarType>(), numRows_(0), numCols_(0), stride_(0), valuesCopied_(true), values_(0)
438 {
439  if ( trans == Teuchos::NO_TRANS )
440  {
441  numRows_ = Source.numRows_;
442  numCols_ = Source.numCols_;
443 
444  if (!Source.valuesCopied_)
445  {
446  stride_ = Source.stride_;
447  values_ = Source.values_;
448  valuesCopied_ = false;
449  }
450  else
451  {
452  stride_ = numRows_;
453  const OrdinalType newsize = stride_ * numCols_;
454  if(newsize > 0) {
455  values_ = new ScalarType[newsize];
456  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
457  }
458  else {
459  numRows_ = 0; numCols_ = 0; stride_ = 0;
460  valuesCopied_ = false;
461  }
462  }
463  }
465  {
466  numRows_ = Source.numCols_;
467  numCols_ = Source.numRows_;
468  stride_ = numRows_;
469  values_ = new ScalarType[stride_*numCols_];
470  for (OrdinalType j=0; j<numCols_; j++) {
471  for (OrdinalType i=0; i<numRows_; i++) {
472  values_[j*stride_ + i] = Teuchos::ScalarTraits<ScalarType>::conjugate(Source.values_[i*Source.stride_ + j]);
473  }
474  }
475  }
476  else
477  {
478  numRows_ = Source.numCols_;
479  numCols_ = Source.numRows_;
480  stride_ = numRows_;
481  values_ = new ScalarType[stride_*numCols_];
482  for (OrdinalType j=0; j<numCols_; j++) {
483  for (OrdinalType i=0; i<numRows_; i++) {
484  values_[j*stride_ + i] = Source.values_[i*Source.stride_ + j];
485  }
486  }
487  }
488 }
489 
490 
491 template<typename OrdinalType, typename ScalarType>
494  )
495  : CompObject(), numRows_(Source.numRows_), numCols_(Source.numCols_), stride_(Source.stride_),
496  valuesCopied_(false), values_(Source.values_)
497 {
498  if(CV == Copy)
499  {
500  stride_ = numRows_;
501  values_ = new ScalarType[stride_ * numCols_];
502  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
503  valuesCopied_ = true;
504  }
505 }
506 
507 
508 template<typename OrdinalType, typename ScalarType>
511  OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType startRow,
512  OrdinalType startCol
513  )
514  : CompObject(), numRows_(numRows_in), numCols_(numCols_in), stride_(Source.stride_),
515  valuesCopied_(false), values_(Source.values_)
516 {
517  if(CV == Copy)
518  {
519  stride_ = numRows_in;
520  values_ = new ScalarType[stride_ * numCols_in];
521  copyMat(Source.values_, Source.stride_, numRows_in, numCols_in, values_, stride_, startRow, startCol);
522  valuesCopied_ = true;
523  }
524  else // CV == View
525  {
526  values_ = values_ + (stride_ * startCol) + startRow;
527  }
528 }
529 
530 template<typename OrdinalType, typename ScalarType>
532 {
533  deleteArrays();
534 }
535 
536 //----------------------------------------------------------------------------------------------------
537 // Shape methods
538 //----------------------------------------------------------------------------------------------------
539 
540 template<typename OrdinalType, typename ScalarType>
542  OrdinalType numRows_in, OrdinalType numCols_in
543  )
544 {
545  deleteArrays(); // Get rid of anything that might be already allocated
546  numRows_ = numRows_in;
547  numCols_ = numCols_in;
548  stride_ = numRows_;
549  values_ = new ScalarType[stride_*numCols_];
550  putScalar();
551  valuesCopied_ = true;
552  return(0);
553 }
554 
555 template<typename OrdinalType, typename ScalarType>
557  OrdinalType numRows_in, OrdinalType numCols_in
558  )
559 {
560  deleteArrays(); // Get rid of anything that might be already allocated
561  numRows_ = numRows_in;
562  numCols_ = numCols_in;
563  stride_ = numRows_;
564  values_ = new ScalarType[stride_*numCols_];
565  valuesCopied_ = true;
566  return(0);
567 }
568 
569 template<typename OrdinalType, typename ScalarType>
571  OrdinalType numRows_in, OrdinalType numCols_in
572  )
573 {
574  // Allocate space for new matrix
575  ScalarType* values_tmp = new ScalarType[numRows_in * numCols_in];
576  ScalarType zero = ScalarTraits<ScalarType>::zero();
577  for(OrdinalType k = 0; k < numRows_in * numCols_in; k++)
578  {
579  values_tmp[k] = zero;
580  }
581  OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
582  OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
583  if(values_ != 0)
584  {
585  copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
586  numRows_in, 0, 0); // Copy principal submatrix of A to new A
587  }
588  deleteArrays(); // Get rid of anything that might be already allocated
589  numRows_ = numRows_in;
590  numCols_ = numCols_in;
591  stride_ = numRows_;
592  values_ = values_tmp; // Set pointer to new A
593  valuesCopied_ = true;
594  return(0);
595 }
596 
597 //----------------------------------------------------------------------------------------------------
598 // Set methods
599 //----------------------------------------------------------------------------------------------------
600 
601 template<typename OrdinalType, typename ScalarType>
603 {
604  // Set each value of the dense matrix to "value".
605  for(OrdinalType j = 0; j < numCols_; j++)
606  {
607  for(OrdinalType i = 0; i < numRows_; i++)
608  {
609  values_[i + j*stride_] = value_in;
610  }
611  }
612  return 0;
613 }
614 
615 template<typename OrdinalType, typename ScalarType>
617 {
618  // Set each value of the dense matrix to a random value.
619  for(OrdinalType j = 0; j < numCols_; j++)
620  {
621  for(OrdinalType i = 0; i < numRows_; i++)
622  {
623  values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
624  }
625  }
626  return 0;
627 }
628 
629 template<typename OrdinalType, typename ScalarType>
633  )
634 {
635  if(this == &Source)
636  return(*this); // Special case of source same as target
637  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
638  return(*this); // Special case of both are views to same data.
639 
640  // If the source is a view then we will return a view, else we will return a copy.
641  if (!Source.valuesCopied_) {
642  if(valuesCopied_) {
643  // Clean up stored data if this was previously a copy.
644  deleteArrays();
645  }
646  numRows_ = Source.numRows_;
647  numCols_ = Source.numCols_;
648  stride_ = Source.stride_;
649  values_ = Source.values_;
650  }
651  else {
652  // If we were a view, we will now be a copy.
653  if(!valuesCopied_) {
654  numRows_ = Source.numRows_;
655  numCols_ = Source.numCols_;
656  stride_ = Source.numRows_;
657  const OrdinalType newsize = stride_ * numCols_;
658  if(newsize > 0) {
659  values_ = new ScalarType[newsize];
660  valuesCopied_ = true;
661  }
662  else {
663  values_ = 0;
664  }
665  }
666  // If we were a copy, we will stay a copy.
667  else {
668  if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) { // we don't need to reallocate
669  numRows_ = Source.numRows_;
670  numCols_ = Source.numCols_;
671  }
672  else { // we need to allocate more space (or less space)
673  deleteArrays();
674  numRows_ = Source.numRows_;
675  numCols_ = Source.numCols_;
676  stride_ = Source.numRows_;
677  const OrdinalType newsize = stride_ * numCols_;
678  if(newsize > 0) {
679  values_ = new ScalarType[newsize];
680  valuesCopied_ = true;
681  }
682  }
683  }
684  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
685  }
686  return(*this);
687 }
688 
689 template<typename OrdinalType, typename ScalarType>
691 {
692  // Check for compatible dimensions
693  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
694  {
695  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
696  }
697  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, ScalarTraits<ScalarType>::one());
698  return(*this);
699 }
700 
701 template<typename OrdinalType, typename ScalarType>
703 {
704  // Check for compatible dimensions
705  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
706  {
707  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
708  }
709  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, -ScalarTraits<ScalarType>::one());
710  return(*this);
711 }
712 
713 template<typename OrdinalType, typename ScalarType>
715  if(this == &Source)
716  return(*this); // Special case of source same as target
717  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
718  return(*this); // Special case of both are views to same data.
719 
720  // Check for compatible dimensions
721  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
722  {
723  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
724  }
725  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
726  return(*this);
727 }
728 
729 //----------------------------------------------------------------------------------------------------
730 // Accessor methods
731 //----------------------------------------------------------------------------------------------------
732 
733 template<typename OrdinalType, typename ScalarType>
734 inline ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
735 {
736 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
737  checkIndex( rowIndex, colIndex );
738 #endif
739  return(values_[colIndex * stride_ + rowIndex]);
740 }
741 
742 template<typename OrdinalType, typename ScalarType>
743 inline const ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
744 {
745 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
746  checkIndex( rowIndex, colIndex );
747 #endif
748  return(values_[colIndex * stride_ + rowIndex]);
749 }
750 
751 template<typename OrdinalType, typename ScalarType>
752 inline const ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex) const
753 {
754 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
755  checkIndex( 0, colIndex );
756 #endif
757  return(values_ + colIndex * stride_);
758 }
759 
760 template<typename OrdinalType, typename ScalarType>
761 inline ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex)
762 {
763 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
764  checkIndex( 0, colIndex );
765 #endif
766  return(values_ + colIndex * stride_);
767 }
768 
769 //----------------------------------------------------------------------------------------------------
770 // Norm methods
771 //----------------------------------------------------------------------------------------------------
772 
773 template<typename OrdinalType, typename ScalarType>
775 {
776  OrdinalType i, j;
779  ScalarType* ptr;
780  for(j = 0; j < numCols_; j++)
781  {
782  ScalarType sum = 0;
783  ptr = values_ + j * stride_;
784  for(i = 0; i < numRows_; i++)
785  {
787  }
789  if(absSum > anorm)
790  {
791  anorm = absSum;
792  }
793  }
794  updateFlops(numRows_ * numCols_);
795  return(anorm);
796 }
797 
798 template<typename OrdinalType, typename ScalarType>
800 {
801  OrdinalType i, j;
803 
804  for (i = 0; i < numRows_; i++) {
806  for (j=0; j< numCols_; j++) {
807  sum += ScalarTraits<ScalarType>::magnitude(*(values_+i+j*stride_));
808  }
809  anorm = TEUCHOS_MAX( anorm, sum );
810  }
811  updateFlops(numRows_ * numCols_);
812  return(anorm);
813 }
814 
815 template<typename OrdinalType, typename ScalarType>
817 {
818  OrdinalType i, j;
820  for (j = 0; j < numCols_; j++) {
821  for (i = 0; i < numRows_; i++) {
822  anorm += ScalarTraits<ScalarType>::magnitude(values_[i+j*stride_]*values_[i+j*stride_]);
823  }
824  }
826  updateFlops(numRows_ * numCols_);
827  return(anorm);
828 }
829 
830 //----------------------------------------------------------------------------------------------------
831 // Comparison methods
832 //----------------------------------------------------------------------------------------------------
833 
834 template<typename OrdinalType, typename ScalarType>
836 {
837  bool result = 1;
838  if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_))
839  {
840  result = 0;
841  }
842  else
843  {
844  OrdinalType i, j;
845  for(i = 0; i < numRows_; i++)
846  {
847  for(j = 0; j < numCols_; j++)
848  {
849  if((*this)(i, j) != Operand(i, j))
850  {
851  return 0;
852  }
853  }
854  }
855  }
856  return result;
857 }
858 
859 template<typename OrdinalType, typename ScalarType>
861 {
862  return !((*this) == Operand);
863 }
864 
865 //----------------------------------------------------------------------------------------------------
866 // Multiplication method
867 //----------------------------------------------------------------------------------------------------
868 
869 template<typename OrdinalType, typename ScalarType>
871 {
872  this->scale( alpha );
873  return(*this);
874 }
875 
876 template<typename OrdinalType, typename ScalarType>
878 {
879  OrdinalType i, j;
880  ScalarType* ptr;
881 
882  for (j=0; j<numCols_; j++) {
883  ptr = values_ + j*stride_;
884  for (i=0; i<numRows_; i++) { *ptr = alpha * (*ptr); ptr++; }
885  }
886  updateFlops( numRows_*numCols_ );
887  return(0);
888 }
889 
890 template<typename OrdinalType, typename ScalarType>
892 {
893  OrdinalType i, j;
894  ScalarType* ptr;
895 
896  // Check for compatible dimensions
897  if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_))
898  {
899  TEUCHOS_CHK_ERR(-1); // Return error
900  }
901  for (j=0; j<numCols_; j++) {
902  ptr = values_ + j*stride_;
903  for (i=0; i<numRows_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
904  }
905  updateFlops( numRows_*numCols_ );
906  return(0);
907 }
908 
909 template<typename OrdinalType, typename ScalarType>
911 {
912  // Check for compatible dimensions
913  OrdinalType A_nrows = (ETranspChar[transa]!='N') ? A.numCols() : A.numRows();
914  OrdinalType A_ncols = (ETranspChar[transa]!='N') ? A.numRows() : A.numCols();
915  OrdinalType B_nrows = (ETranspChar[transb]!='N') ? B.numCols() : B.numRows();
916  OrdinalType B_ncols = (ETranspChar[transb]!='N') ? B.numRows() : B.numCols();
917  if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols))
918  {
919  TEUCHOS_CHK_ERR(-1); // Return error
920  }
921  // Call GEMM function
922  this->GEMM(transa, transb, numRows_, numCols_, A_ncols, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
923  double nflops = 2 * numRows_;
924  nflops *= numCols_;
925  nflops *= A_ncols;
926  updateFlops(nflops);
927  return(0);
928 }
929 
930 template<typename OrdinalType, typename ScalarType>
932 {
933  // Check for compatible dimensions
934  OrdinalType A_nrows = A.numRows(), A_ncols = A.numCols();
935  OrdinalType B_nrows = B.numRows(), B_ncols = B.numCols();
936 
937  if (ESideChar[sideA]=='L') {
938  if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols)) {
939  TEUCHOS_CHK_ERR(-1); // Return error
940  }
941  } else {
942  if ((numRows_ != B_nrows) || (B_ncols != A_nrows) || (numCols_ != A_ncols)) {
943  TEUCHOS_CHK_ERR(-1); // Return error
944  }
945  }
946 
947  // Call SYMM function
949  this->SYMM(sideA, uplo, numRows_, numCols_, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
950  double nflops = 2 * numRows_;
951  nflops *= numCols_;
952  nflops *= A_ncols;
953  updateFlops(nflops);
954  return(0);
955 }
956 
957 template<typename OrdinalType, typename ScalarType>
959 {
960  os << std::endl;
961  if(valuesCopied_)
962  os << "Values_copied : yes" << std::endl;
963  else
964  os << "Values_copied : no" << std::endl;
965  os << "Rows : " << numRows_ << std::endl;
966  os << "Columns : " << numCols_ << std::endl;
967  os << "LDA : " << stride_ << std::endl;
968  if(numRows_ == 0 || numCols_ == 0) {
969  os << "(matrix is empty, no values to display)" << std::endl;
970  } else {
971  for(OrdinalType i = 0; i < numRows_; i++) {
972  for(OrdinalType j = 0; j < numCols_; j++){
973  os << (*this)(i,j) << " ";
974  }
975  os << std::endl;
976  }
977  }
978 }
979 
980 //----------------------------------------------------------------------------------------------------
981 // Protected methods
982 //----------------------------------------------------------------------------------------------------
983 
984 template<typename OrdinalType, typename ScalarType>
985 inline void SerialDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
986  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRows_, std::out_of_range,
987  "SerialDenseMatrix<T>::checkIndex: "
988  "Row index " << rowIndex << " out of range [0, "<< numRows_ << ")");
989  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numCols_, std::out_of_range,
990  "SerialDenseMatrix<T>::checkIndex: "
991  "Col index " << colIndex << " out of range [0, "<< numCols_ << ")");
992 }
993 
994 template<typename OrdinalType, typename ScalarType>
995 void SerialDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
996 {
997  if (valuesCopied_)
998  {
999  delete [] values_;
1000  values_ = 0;
1001  valuesCopied_ = false;
1002  }
1003 }
1004 
1005 template<typename OrdinalType, typename ScalarType>
1006 void SerialDenseMatrix<OrdinalType, ScalarType>::copyMat(
1007  ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
1008  OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput,
1009  OrdinalType startRow, OrdinalType startCol, ScalarType alpha
1010  )
1011 {
1012  OrdinalType i, j;
1013  ScalarType* ptr1 = 0;
1014  ScalarType* ptr2 = 0;
1015  for(j = 0; j < numCols_in; j++) {
1016  ptr1 = outputMatrix + (j * strideOutput);
1017  ptr2 = inputMatrix + (j + startCol) * strideInput + startRow;
1018  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1019  for(i = 0; i < numRows_in; i++)
1020  {
1021  *ptr1++ += alpha*(*ptr2++);
1022  }
1023  } else {
1024  for(i = 0; i < numRows_in; i++)
1025  {
1026  *ptr1++ = *ptr2++;
1027  }
1028  }
1029  }
1030 }
1031 
1032 } // namespace Teuchos
1033 
1034 
1035 #endif /* _TEUCHOS_SERIALDENSEMATRIX_HPP_ */
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
OrdinalType numCols() const
Returns the column dimension of this matrix.
ScalarType * values() const
Data array access method.
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.
SerialDenseMatrix< OrdinalType, ScalarType > & operator*=(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
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.
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.
virtual void print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream << operator inherited from the Object class...
Object for storing data and providing functionality that is common to all computational classes...
This structure defines some basic traits for a scalar field type.
bool operator==(const SerialDenseMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
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.
OrdinalType numRows() const
Returns the row dimension of this matrix.
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
bool upper() const
Returns true if upper triangular part of this matrix has and will be used.
OrdinalType numRows() const
Returns the row dimension of this matrix.
ScalarType scalarType
Typedef for scalar type.
The base Teuchos class.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Set all values in the matrix to a constant value.
SerialDenseMatrix()
Default Constructor.
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
bool empty() const
Returns whether this matrix is empty.
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.
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
int random()
Set all values in the matrix to be random numbers.
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos, as well as a number of utility routines.
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...
bool operator!=(const SerialDenseMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
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.
ScalarType * values() const
Returns the pointer to the ScalarType data array contained in the object.
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
This class creates and provides basic support for dense rectangular matrix of templated type...