Teuchos Package Browser (Single Doxygen Collection)  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_;
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++) {
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
948  EUplo uplo = (A.upper() ? Teuchos::UPPER_TRI : Teuchos::LOWER_TRI);
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>
996 {
997  if (valuesCopied_)
998  {
999  delete [] values_;
1000  values_ = 0;
1001  valuesCopied_ = false;
1002  }
1003 }
1004 
1005 template<typename OrdinalType, typename ScalarType>
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.
ScalarType * values() const
Data array access method.
int shapeUninitialized(OrdinalType numRows, OrdinalType numCols)
Same as shape() except leaves uninitialized.
Definition: PackageB.cpp:3
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...
void checkIndex(OrdinalType rowIndex, OrdinalType colIndex=0) const
Object for storing data and providing functionality that is common to all computational classes...
Definition: PackageA.cpp:3
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.
#define TEUCHOS_CHK_ERR(a)
SerialDenseMatrix< OrdinalType, ScalarType > & operator=(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
Functionality and data that is common to all computational classes.
static T conjugate(T a)
Returns the conjugate of the scalar type a.
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
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.
#define TEUCHOS_MAX(x, y)
SerialDenseMatrix()
Default Constructor.
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ESideChar[]
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
bool empty() const
Returns whether this matrix is empty.
void copyMat(ScalarType *inputMatrix, OrdinalType strideInput, OrdinalType numRows, OrdinalType numCols, ScalarType *outputMatrix, OrdinalType strideOutput, OrdinalType startRow, OrdinalType startCol, ScalarType alpha=ScalarTraits< ScalarType >::zero())
Templated serial, dense, symmetric matrix class.
OrdinalType ordinalType
Typedef for ordinal type.
ScalarType * operator[](OrdinalType colIndex)
Column access method (non-const).
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
int random()
Set all values in the matrix to be random numbers.
int reshape(OrdinalType numRows, OrdinalType numCols)
Reshaping method for changing the size of a SerialDenseMatrix, keeping the entries.
Defines basic traits for the scalar field type.
static T zero()
Returns representation of zero for this scalar type.
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ETranspChar[]
static T random()
Returns a random number (between -one() and +one()) of this scalar type.
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
int shape(OrdinalType numRows, OrdinalType numCols)
Shape method for changing the size of a SerialDenseMatrix, initializing entries to zero...
#define TEUCHOS_CHK_REF(a)
bool operator!=(const SerialDenseMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
OrdinalType numCols() const
Returns the column dimension of this matrix.
#define TEUCHOS_MIN(x, y)
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
Teuchos::DataAccess Mode enumerable type.
This class creates and provides basic support for dense rectangular matrix of templated type...