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(), 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(), 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(), 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>
436 SerialDenseMatrix<OrdinalType, ScalarType>::SerialDenseMatrix(const SerialDenseMatrix<OrdinalType, ScalarType> &Source, ETransp trans) : CompObject(), BLAS<OrdinalType,ScalarType>(), numRows_(0), numCols_(0), stride_(0), valuesCopied_(true), values_(0)
437 {
438  if ( trans == Teuchos::NO_TRANS )
439  {
440  numRows_ = Source.numRows_;
441  numCols_ = Source.numCols_;
442 
443  if (!Source.valuesCopied_)
444  {
445  stride_ = Source.stride_;
446  values_ = Source.values_;
447  valuesCopied_ = false;
448  }
449  else
450  {
451  stride_ = numRows_;
452  const OrdinalType newsize = stride_ * numCols_;
453  if(newsize > 0) {
454  values_ = new ScalarType[newsize];
455  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
456  }
457  else {
458  numRows_ = 0; numCols_ = 0; stride_ = 0;
459  valuesCopied_ = false;
460  }
461  }
462  }
464  {
465  numRows_ = Source.numCols_;
466  numCols_ = Source.numRows_;
467  stride_ = numRows_;
468  values_ = new ScalarType[stride_*numCols_];
469  for (OrdinalType j=0; j<numCols_; j++) {
470  for (OrdinalType i=0; i<numRows_; i++) {
471  values_[j*stride_ + i] = Teuchos::ScalarTraits<ScalarType>::conjugate(Source.values_[i*Source.stride_ + j]);
472  }
473  }
474  }
475  else
476  {
477  numRows_ = Source.numCols_;
478  numCols_ = Source.numRows_;
479  stride_ = numRows_;
480  values_ = new ScalarType[stride_*numCols_];
481  for (OrdinalType j=0; j<numCols_; j++) {
482  for (OrdinalType i=0; i<numRows_; i++) {
483  values_[j*stride_ + i] = Source.values_[i*Source.stride_ + j];
484  }
485  }
486  }
487 }
488 
489 
490 template<typename OrdinalType, typename ScalarType>
493  )
494  : CompObject(), numRows_(Source.numRows_), numCols_(Source.numCols_), stride_(Source.stride_),
495  valuesCopied_(false), values_(Source.values_)
496 {
497  if(CV == Copy)
498  {
499  stride_ = numRows_;
500  values_ = new ScalarType[stride_ * numCols_];
501  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
502  valuesCopied_ = true;
503  }
504 }
505 
506 
507 template<typename OrdinalType, typename ScalarType>
510  OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType startRow,
511  OrdinalType startCol
512  )
513  : CompObject(), numRows_(numRows_in), numCols_(numCols_in), stride_(Source.stride_),
514  valuesCopied_(false), values_(Source.values_)
515 {
516  if(CV == Copy)
517  {
518  stride_ = numRows_in;
519  values_ = new ScalarType[stride_ * numCols_in];
520  copyMat(Source.values_, Source.stride_, numRows_in, numCols_in, values_, stride_, startRow, startCol);
521  valuesCopied_ = true;
522  }
523  else // CV == View
524  {
525  values_ = values_ + (stride_ * startCol) + startRow;
526  }
527 }
528 
529 template<typename OrdinalType, typename ScalarType>
531 {
532  deleteArrays();
533 }
534 
535 //----------------------------------------------------------------------------------------------------
536 // Shape methods
537 //----------------------------------------------------------------------------------------------------
538 
539 template<typename OrdinalType, typename ScalarType>
541  OrdinalType numRows_in, OrdinalType numCols_in
542  )
543 {
544  deleteArrays(); // Get rid of anything that might be already allocated
545  numRows_ = numRows_in;
546  numCols_ = numCols_in;
547  stride_ = numRows_;
548  values_ = new ScalarType[stride_*numCols_];
549  putScalar();
550  valuesCopied_ = true;
551  return(0);
552 }
553 
554 template<typename OrdinalType, typename ScalarType>
556  OrdinalType numRows_in, OrdinalType numCols_in
557  )
558 {
559  deleteArrays(); // Get rid of anything that might be already allocated
560  numRows_ = numRows_in;
561  numCols_ = numCols_in;
562  stride_ = numRows_;
563  values_ = new ScalarType[stride_*numCols_];
564  valuesCopied_ = true;
565  return(0);
566 }
567 
568 template<typename OrdinalType, typename ScalarType>
570  OrdinalType numRows_in, OrdinalType numCols_in
571  )
572 {
573  // Allocate space for new matrix
574  ScalarType* values_tmp = new ScalarType[numRows_in * numCols_in];
575  ScalarType zero = ScalarTraits<ScalarType>::zero();
576  for(OrdinalType k = 0; k < numRows_in * numCols_in; k++)
577  {
578  values_tmp[k] = zero;
579  }
580  OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
581  OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
582  if(values_ != 0)
583  {
584  copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
585  numRows_in, 0, 0); // Copy principal submatrix of A to new A
586  }
587  deleteArrays(); // Get rid of anything that might be already allocated
588  numRows_ = numRows_in;
589  numCols_ = numCols_in;
590  stride_ = numRows_;
591  values_ = values_tmp; // Set pointer to new A
592  valuesCopied_ = true;
593  return(0);
594 }
595 
596 //----------------------------------------------------------------------------------------------------
597 // Set methods
598 //----------------------------------------------------------------------------------------------------
599 
600 template<typename OrdinalType, typename ScalarType>
602 {
603  // Set each value of the dense matrix to "value".
604  for(OrdinalType j = 0; j < numCols_; j++)
605  {
606  for(OrdinalType i = 0; i < numRows_; i++)
607  {
608  values_[i + j*stride_] = value_in;
609  }
610  }
611  return 0;
612 }
613 
614 template<typename OrdinalType, typename ScalarType>
616 {
617  // Set each value of the dense matrix to a random value.
618  for(OrdinalType j = 0; j < numCols_; j++)
619  {
620  for(OrdinalType i = 0; i < numRows_; i++)
621  {
622  values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
623  }
624  }
625  return 0;
626 }
627 
628 template<typename OrdinalType, typename ScalarType>
632  )
633 {
634  if(this == &Source)
635  return(*this); // Special case of source same as target
636  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
637  return(*this); // Special case of both are views to same data.
638 
639  // If the source is a view then we will return a view, else we will return a copy.
640  if (!Source.valuesCopied_) {
641  if(valuesCopied_) {
642  // Clean up stored data if this was previously a copy.
643  deleteArrays();
644  }
645  numRows_ = Source.numRows_;
646  numCols_ = Source.numCols_;
647  stride_ = Source.stride_;
648  values_ = Source.values_;
649  }
650  else {
651  // If we were a view, we will now be a copy.
652  if(!valuesCopied_) {
653  numRows_ = Source.numRows_;
654  numCols_ = Source.numCols_;
655  stride_ = Source.numRows_;
656  const OrdinalType newsize = stride_ * numCols_;
657  if(newsize > 0) {
658  values_ = new ScalarType[newsize];
659  valuesCopied_ = true;
660  }
661  else {
662  values_ = 0;
663  }
664  }
665  // If we were a copy, we will stay a copy.
666  else {
667  if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) { // we don't need to reallocate
668  numRows_ = Source.numRows_;
669  numCols_ = Source.numCols_;
670  }
671  else { // we need to allocate more space (or less space)
672  deleteArrays();
673  numRows_ = Source.numRows_;
674  numCols_ = Source.numCols_;
675  stride_ = Source.numRows_;
676  const OrdinalType newsize = stride_ * numCols_;
677  if(newsize > 0) {
678  values_ = new ScalarType[newsize];
679  valuesCopied_ = true;
680  }
681  }
682  }
683  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
684  }
685  return(*this);
686 }
687 
688 template<typename OrdinalType, typename ScalarType>
690 {
691  // Check for compatible dimensions
692  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
693  {
694  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
695  }
696  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, ScalarTraits<ScalarType>::one());
697  return(*this);
698 }
699 
700 template<typename OrdinalType, typename ScalarType>
702 {
703  // Check for compatible dimensions
704  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
705  {
706  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
707  }
708  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, -ScalarTraits<ScalarType>::one());
709  return(*this);
710 }
711 
712 template<typename OrdinalType, typename ScalarType>
714  if(this == &Source)
715  return(*this); // Special case of source same as target
716  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
717  return(*this); // Special case of both are views to same data.
718 
719  // Check for compatible dimensions
720  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
721  {
722  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
723  }
724  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
725  return(*this);
726 }
727 
728 //----------------------------------------------------------------------------------------------------
729 // Accessor methods
730 //----------------------------------------------------------------------------------------------------
731 
732 template<typename OrdinalType, typename ScalarType>
733 inline ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
734 {
735 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
736  checkIndex( rowIndex, colIndex );
737 #endif
738  return(values_[colIndex * stride_ + rowIndex]);
739 }
740 
741 template<typename OrdinalType, typename ScalarType>
742 inline const ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
743 {
744 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
745  checkIndex( rowIndex, colIndex );
746 #endif
747  return(values_[colIndex * stride_ + rowIndex]);
748 }
749 
750 template<typename OrdinalType, typename ScalarType>
751 inline const ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex) const
752 {
753 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
754  checkIndex( 0, colIndex );
755 #endif
756  return(values_ + colIndex * stride_);
757 }
758 
759 template<typename OrdinalType, typename ScalarType>
760 inline ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex)
761 {
762 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
763  checkIndex( 0, colIndex );
764 #endif
765  return(values_ + colIndex * stride_);
766 }
767 
768 //----------------------------------------------------------------------------------------------------
769 // Norm methods
770 //----------------------------------------------------------------------------------------------------
771 
772 template<typename OrdinalType, typename ScalarType>
774 {
775  OrdinalType i, j;
778  ScalarType* ptr;
779  for(j = 0; j < numCols_; j++)
780  {
781  ScalarType sum = 0;
782  ptr = values_ + j * stride_;
783  for(i = 0; i < numRows_; i++)
784  {
786  }
788  if(absSum > anorm)
789  {
790  anorm = absSum;
791  }
792  }
793  updateFlops(numRows_ * numCols_);
794  return(anorm);
795 }
796 
797 template<typename OrdinalType, typename ScalarType>
799 {
800  OrdinalType i, j;
802 
803  for (i = 0; i < numRows_; i++) {
805  for (j=0; j< numCols_; j++) {
806  sum += ScalarTraits<ScalarType>::magnitude(*(values_+i+j*stride_));
807  }
808  anorm = TEUCHOS_MAX( anorm, sum );
809  }
810  updateFlops(numRows_ * numCols_);
811  return(anorm);
812 }
813 
814 template<typename OrdinalType, typename ScalarType>
816 {
817  OrdinalType i, j;
819  for (j = 0; j < numCols_; j++) {
820  for (i = 0; i < numRows_; i++) {
821  anorm += ScalarTraits<ScalarType>::magnitude(values_[i+j*stride_]*values_[i+j*stride_]);
822  }
823  }
825  updateFlops(numRows_ * numCols_);
826  return(anorm);
827 }
828 
829 //----------------------------------------------------------------------------------------------------
830 // Comparison methods
831 //----------------------------------------------------------------------------------------------------
832 
833 template<typename OrdinalType, typename ScalarType>
835 {
836  bool result = 1;
837  if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_))
838  {
839  result = 0;
840  }
841  else
842  {
843  OrdinalType i, j;
844  for(i = 0; i < numRows_; i++)
845  {
846  for(j = 0; j < numCols_; j++)
847  {
848  if((*this)(i, j) != Operand(i, j))
849  {
850  return 0;
851  }
852  }
853  }
854  }
855  return result;
856 }
857 
858 template<typename OrdinalType, typename ScalarType>
860 {
861  return !((*this) == Operand);
862 }
863 
864 //----------------------------------------------------------------------------------------------------
865 // Multiplication method
866 //----------------------------------------------------------------------------------------------------
867 
868 template<typename OrdinalType, typename ScalarType>
870 {
871  this->scale( alpha );
872  return(*this);
873 }
874 
875 template<typename OrdinalType, typename ScalarType>
877 {
878  OrdinalType i, j;
879  ScalarType* ptr;
880 
881  for (j=0; j<numCols_; j++) {
882  ptr = values_ + j*stride_;
883  for (i=0; i<numRows_; i++) { *ptr = alpha * (*ptr); ptr++; }
884  }
885  updateFlops( numRows_*numCols_ );
886  return(0);
887 }
888 
889 template<typename OrdinalType, typename ScalarType>
891 {
892  OrdinalType i, j;
893  ScalarType* ptr;
894 
895  // Check for compatible dimensions
896  if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_))
897  {
898  TEUCHOS_CHK_ERR(-1); // Return error
899  }
900  for (j=0; j<numCols_; j++) {
901  ptr = values_ + j*stride_;
902  for (i=0; i<numRows_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
903  }
904  updateFlops( numRows_*numCols_ );
905  return(0);
906 }
907 
908 template<typename OrdinalType, typename ScalarType>
910 {
911  // Check for compatible dimensions
912  OrdinalType A_nrows = (ETranspChar[transa]!='N') ? A.numCols() : A.numRows();
913  OrdinalType A_ncols = (ETranspChar[transa]!='N') ? A.numRows() : A.numCols();
914  OrdinalType B_nrows = (ETranspChar[transb]!='N') ? B.numCols() : B.numRows();
915  OrdinalType B_ncols = (ETranspChar[transb]!='N') ? B.numRows() : B.numCols();
916  if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols))
917  {
918  TEUCHOS_CHK_ERR(-1); // Return error
919  }
920  // Call GEMM function
921  this->GEMM(transa, transb, numRows_, numCols_, A_ncols, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
922  double nflops = 2 * numRows_;
923  nflops *= numCols_;
924  nflops *= A_ncols;
925  updateFlops(nflops);
926  return(0);
927 }
928 
929 template<typename OrdinalType, typename ScalarType>
931 {
932  // Check for compatible dimensions
933  OrdinalType A_nrows = A.numRows(), A_ncols = A.numCols();
934  OrdinalType B_nrows = B.numRows(), B_ncols = B.numCols();
935 
936  if (ESideChar[sideA]=='L') {
937  if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols)) {
938  TEUCHOS_CHK_ERR(-1); // Return error
939  }
940  } else {
941  if ((numRows_ != B_nrows) || (B_ncols != A_nrows) || (numCols_ != A_ncols)) {
942  TEUCHOS_CHK_ERR(-1); // Return error
943  }
944  }
945 
946  // Call SYMM function
948  this->SYMM(sideA, uplo, numRows_, numCols_, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
949  double nflops = 2 * numRows_;
950  nflops *= numCols_;
951  nflops *= A_ncols;
952  updateFlops(nflops);
953  return(0);
954 }
955 
956 template<typename OrdinalType, typename ScalarType>
958 {
959  os << std::endl;
960  if(valuesCopied_)
961  os << "Values_copied : yes" << std::endl;
962  else
963  os << "Values_copied : no" << std::endl;
964  os << "Rows : " << numRows_ << std::endl;
965  os << "Columns : " << numCols_ << std::endl;
966  os << "LDA : " << stride_ << std::endl;
967  if(numRows_ == 0 || numCols_ == 0) {
968  os << "(matrix is empty, no values to display)" << std::endl;
969  } else {
970  for(OrdinalType i = 0; i < numRows_; i++) {
971  for(OrdinalType j = 0; j < numCols_; j++){
972  os << (*this)(i,j) << " ";
973  }
974  os << std::endl;
975  }
976  }
977 }
978 
979 //----------------------------------------------------------------------------------------------------
980 // Protected methods
981 //----------------------------------------------------------------------------------------------------
982 
983 template<typename OrdinalType, typename ScalarType>
984 inline void SerialDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
985  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRows_, std::out_of_range,
986  "SerialDenseMatrix<T>::checkIndex: "
987  "Row index " << rowIndex << " out of range [0, "<< numRows_ << ")");
988  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numCols_, std::out_of_range,
989  "SerialDenseMatrix<T>::checkIndex: "
990  "Col index " << colIndex << " out of range [0, "<< numCols_ << ")");
991 }
992 
993 template<typename OrdinalType, typename ScalarType>
995 {
996  if (valuesCopied_)
997  {
998  delete [] values_;
999  values_ = 0;
1000  valuesCopied_ = false;
1001  }
1002 }
1003 
1004 template<typename OrdinalType, typename ScalarType>
1006  ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
1007  OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput,
1008  OrdinalType startRow, OrdinalType startCol, ScalarType alpha
1009  )
1010 {
1011  OrdinalType i, j;
1012  ScalarType* ptr1 = 0;
1013  ScalarType* ptr2 = 0;
1014  for(j = 0; j < numCols_in; j++) {
1015  ptr1 = outputMatrix + (j * strideOutput);
1016  ptr2 = inputMatrix + (j + startCol) * strideInput + startRow;
1017  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1018  for(i = 0; i < numRows_in; i++)
1019  {
1020  *ptr1++ += alpha*(*ptr2++);
1021  }
1022  } else {
1023  for(i = 0; i < numRows_in; i++)
1024  {
1025  *ptr1++ = *ptr2++;
1026  }
1027  }
1028  }
1029 }
1030 
1031 } // namespace Teuchos
1032 
1033 
1034 #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.
virtual void print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream << operator inherited from the Object class...
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.
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.
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...
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.
void SYMM(ESide side, EUplo uplo, const OrdinalType m, const OrdinalType n, const alpha_type alpha, const A_type *A, const OrdinalType lda, const B_type *B, const OrdinalType ldb, const beta_type beta, ScalarType *C, const OrdinalType ldc) const
Performs the matrix-matrix operation: C <- alpha*A*B+beta*C or C <- alpha*B*A+beta*C where A is an m ...
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.
The base Teuchos class.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Set all values in the matrix to a constant value.
void GEMM(ETransp transa, ETransp transb, const OrdinalType m, const OrdinalType n, const OrdinalType k, const alpha_type alpha, const A_type *A, const OrdinalType lda, const B_type *B, const OrdinalType ldb, const beta_type beta, ScalarType *C, const OrdinalType ldc) const
General matrix-matrix multiply.
SerialDenseMatrix()
Default Constructor.
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.
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.
void updateFlops(int addflops) const
Increment Flop count for this object.
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...