Teuchos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Teuchos_SerialDenseMatrix.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
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 
56 #include <utility>
57 
66 namespace Teuchos {
67 
68 template<typename OrdinalType, typename ScalarType>
69 class SerialDenseMatrix : public CompObject, public BLAS<OrdinalType, ScalarType>
70 {
71 public:
72 
74  typedef OrdinalType ordinalType;
76  typedef ScalarType scalarType;
77 
79 
80 
82 
86 
88 
96  SerialDenseMatrix(OrdinalType numRows, OrdinalType numCols, bool zeroOut = true);
97 
99 
107  SerialDenseMatrix(DataAccess CV, ScalarType* values, OrdinalType stride, OrdinalType numRows, OrdinalType numCols);
108 
110 
116 
118 
121 
123 
135  SerialDenseMatrix(DataAccess CV, const SerialDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRows, OrdinalType numCols, OrdinalType startRow=0, OrdinalType startCol=0);
136 
138  virtual ~SerialDenseMatrix();
140 
142 
143 
154  int shape(OrdinalType numRows, OrdinalType numCols);
155 
157  int shapeUninitialized(OrdinalType numRows, OrdinalType numCols);
158 
160 
170  int reshape(OrdinalType numRows, OrdinalType numCols);
171 
172 
174 
176 
177 
179 
186 
188 
194 
196 
199  SerialDenseMatrix<OrdinalType, ScalarType>& operator= (const ScalarType value) { putScalar(value); return(*this); }
200 
202 
206  int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero() );
207 
209 
213 
215  int random();
216 
218 
220 
221 
223 
230  ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
231 
233 
240  const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
241 
243 
250  ScalarType* operator [] (OrdinalType colIndex);
251 
253 
260  const ScalarType* operator [] (OrdinalType colIndex) const;
261 
263 
264  ScalarType* values() const { return(values_); }
265 
267 
269 
270 
272 
276 
278 
282 
284 
288 
290 
294  int scale ( const ScalarType alpha );
295 
297 
304 
306 
320  int multiply (ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
321 
323 
334  int multiply (ESide sideA, ScalarType alpha, const SerialSymDenseMatrix<OrdinalType, ScalarType> &A, const SerialDenseMatrix<OrdinalType, ScalarType> &B, ScalarType beta);
335 
337 
339 
340 
342 
345  bool operator== (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand) const;
346 
348 
351  bool operator!= (const SerialDenseMatrix<OrdinalType, ScalarType> &Operand) const;
352 
354 
356 
357 
359  OrdinalType numRows() const { return(numRows_); }
360 
362  OrdinalType numCols() const { return(numCols_); }
363 
365  OrdinalType stride() const { return(stride_); }
366 
368  bool empty() const { return(numRows_ == 0 || numCols_ == 0); }
370 
372 
373 
376 
379 
383 
385 
386  virtual void print(std::ostream& os) const;
388 
390 protected:
391  void copyMat(ScalarType* inputMatrix, OrdinalType strideInput,
392  OrdinalType numRows, OrdinalType numCols, ScalarType* outputMatrix,
393  OrdinalType strideOutput, OrdinalType startRow, OrdinalType startCol,
394  ScalarType alpha = ScalarTraits<ScalarType>::zero() );
395  void deleteArrays();
396  void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
397  OrdinalType numRows_;
398  OrdinalType numCols_;
399  OrdinalType stride_;
401  ScalarType* values_;
402 
403 }; // class Teuchos_SerialDenseMatrix
404 
405 //----------------------------------------------------------------------------------------------------
406 // Constructors and Destructor
407 //----------------------------------------------------------------------------------------------------
408 
409 template<typename OrdinalType, typename ScalarType>
411  : CompObject(), BLAS<OrdinalType,ScalarType>(), numRows_(0), numCols_(0), stride_(0), valuesCopied_(false), values_(0)
412 {}
413 
414 template<typename OrdinalType, typename ScalarType>
416  OrdinalType numRows_in, OrdinalType numCols_in, bool zeroOut
417  )
418  : CompObject(), BLAS<OrdinalType,ScalarType>(), numRows_(numRows_in), numCols_(numCols_in), stride_(numRows_in)
419 {
420  values_ = new ScalarType[stride_*numCols_];
421  valuesCopied_ = true;
422  if (zeroOut == true)
423  putScalar();
424 }
425 
426 template<typename OrdinalType, typename ScalarType>
428  DataAccess CV, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRows_in,
429  OrdinalType numCols_in
430  )
431  : CompObject(), BLAS<OrdinalType,ScalarType>(), numRows_(numRows_in), numCols_(numCols_in), stride_(stride_in),
432  valuesCopied_(false), values_(values_in)
433 {
434  if(CV == Copy)
435  {
436  stride_ = numRows_;
437  values_ = new ScalarType[stride_*numCols_];
438  copyMat(values_in, stride_in, numRows_, numCols_, values_, stride_, 0, 0);
439  valuesCopied_ = true;
440  }
441 }
442 
443 template<typename OrdinalType, typename ScalarType>
445  : CompObject(),BLAS<OrdinalType,ScalarType>(), numRows_(0), numCols_(0), stride_(0), valuesCopied_(true), values_(0)
446 {
447  if ( trans == Teuchos::NO_TRANS )
448  {
449  numRows_ = Source.numRows_;
450  numCols_ = Source.numCols_;
451 
452  if (!Source.valuesCopied_)
453  {
454  stride_ = Source.stride_;
455  values_ = Source.values_;
456  valuesCopied_ = false;
457  }
458  else
459  {
460  stride_ = numRows_;
461  const OrdinalType newsize = stride_ * numCols_;
462  if(newsize > 0) {
463  values_ = new ScalarType[newsize];
464  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
465  }
466  else {
467  numRows_ = 0; numCols_ = 0; stride_ = 0;
468  valuesCopied_ = false;
469  }
470  }
471  }
473  {
474  numRows_ = Source.numCols_;
475  numCols_ = Source.numRows_;
476  stride_ = numRows_;
477  values_ = new ScalarType[stride_*numCols_];
478  for (OrdinalType j=0; j<numCols_; j++) {
479  for (OrdinalType i=0; i<numRows_; i++) {
481  }
482  }
483  }
484  else
485  {
486  numRows_ = Source.numCols_;
487  numCols_ = Source.numRows_;
488  stride_ = numRows_;
489  values_ = new ScalarType[stride_*numCols_];
490  for (OrdinalType j=0; j<numCols_; j++) {
491  for (OrdinalType i=0; i<numRows_; i++) {
492  values_[j*stride_ + i] = Source.values_[i*Source.stride_ + j];
493  }
494  }
495  }
496 }
497 
498 
499 template<typename OrdinalType, typename ScalarType>
502  )
503  : CompObject(), numRows_(Source.numRows_), numCols_(Source.numCols_), stride_(Source.stride_),
504  valuesCopied_(false), values_(Source.values_)
505 {
506  if(CV == Copy)
507  {
508  stride_ = numRows_;
509  values_ = new ScalarType[stride_ * numCols_];
510  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
511  valuesCopied_ = true;
512  }
513 }
514 
515 
516 template<typename OrdinalType, typename ScalarType>
519  OrdinalType numRows_in, OrdinalType numCols_in, OrdinalType startRow,
520  OrdinalType startCol
521  )
522  : CompObject(), numRows_(numRows_in), numCols_(numCols_in), stride_(Source.stride_),
523  valuesCopied_(false), values_(Source.values_)
524 {
525  if(CV == Copy)
526  {
527  stride_ = numRows_in;
528  values_ = new ScalarType[stride_ * numCols_in];
529  copyMat(Source.values_, Source.stride_, numRows_in, numCols_in, values_, stride_, startRow, startCol);
530  valuesCopied_ = true;
531  }
532  else // CV == View
533  {
534  values_ = values_ + (stride_ * startCol) + startRow;
535  }
536 }
537 
538 template<typename OrdinalType, typename ScalarType>
540 {
541  deleteArrays();
542 }
543 
544 //----------------------------------------------------------------------------------------------------
545 // Shape methods
546 //----------------------------------------------------------------------------------------------------
547 
548 template<typename OrdinalType, typename ScalarType>
550  OrdinalType numRows_in, OrdinalType numCols_in
551  )
552 {
553  deleteArrays(); // Get rid of anything that might be already allocated
554  numRows_ = numRows_in;
555  numCols_ = numCols_in;
556  stride_ = numRows_;
557  values_ = new ScalarType[stride_*numCols_];
558  putScalar();
559  valuesCopied_ = true;
560  return(0);
561 }
562 
563 template<typename OrdinalType, typename ScalarType>
565  OrdinalType numRows_in, OrdinalType numCols_in
566  )
567 {
568  deleteArrays(); // Get rid of anything that might be already allocated
569  numRows_ = numRows_in;
570  numCols_ = numCols_in;
571  stride_ = numRows_;
572  values_ = new ScalarType[stride_*numCols_];
573  valuesCopied_ = true;
574  return(0);
575 }
576 
577 template<typename OrdinalType, typename ScalarType>
579  OrdinalType numRows_in, OrdinalType numCols_in
580  )
581 {
582  // Allocate space for new matrix
583  ScalarType* values_tmp = new ScalarType[numRows_in * numCols_in];
584  ScalarType zero = ScalarTraits<ScalarType>::zero();
585  for(OrdinalType k = 0; k < numRows_in * numCols_in; k++)
586  {
587  values_tmp[k] = zero;
588  }
589  OrdinalType numRows_tmp = TEUCHOS_MIN(numRows_, numRows_in);
590  OrdinalType numCols_tmp = TEUCHOS_MIN(numCols_, numCols_in);
591  if(values_ != 0)
592  {
593  copyMat(values_, stride_, numRows_tmp, numCols_tmp, values_tmp,
594  numRows_in, 0, 0); // Copy principal submatrix of A to new A
595  }
596  deleteArrays(); // Get rid of anything that might be already allocated
597  numRows_ = numRows_in;
598  numCols_ = numCols_in;
599  stride_ = numRows_;
600  values_ = values_tmp; // Set pointer to new A
601  valuesCopied_ = true;
602  return(0);
603 }
604 
605 //----------------------------------------------------------------------------------------------------
606 // Set methods
607 //----------------------------------------------------------------------------------------------------
608 
609 template<typename OrdinalType, typename ScalarType>
611 {
612  // Set each value of the dense matrix to "value".
613  for(OrdinalType j = 0; j < numCols_; j++)
614  {
615  for(OrdinalType i = 0; i < numRows_; i++)
616  {
617  values_[i + j*stride_] = value_in;
618  }
619  }
620  return 0;
621 }
622 
623 template<typename OrdinalType, typename ScalarType> void
626 {
627  // Notes:
628  // > DefaultBLASImpl::SWAP() uses a deep copy. This fn uses a pointer swap.
629  // > this fn covers both Vector and Matrix, such that some care must be
630  // employed to not swap across types (creating a Vector with non-unitary
631  // numCols_)
632  // > Inherited data that is not currently swapped (since inactive/deprecated):
633  // >> Teuchos::CompObject:
634  // Flops *flopCounter_ [Note: all SerialDenseMatrix ctors initialize a
635  // NULL flop-counter using CompObject(), such that any flop increments
636  // that are computed are not accumulated.]
637  // >> Teuchos::Object: (now removed from inheritance list)
638  // static int tracebackMode (no swap for statics)
639  // std::string label_ (has been reported as a cause of memory overhead)
640 
641  std::swap(values_ , B.values_);
642  std::swap(numRows_, B.numRows_);
643  std::swap(numCols_, B.numCols_);
644  std::swap(stride_, B.stride_);
645  std::swap(valuesCopied_, B.valuesCopied_);
646 }
647 
648 template<typename OrdinalType, typename ScalarType>
650 {
651  // Set each value of the dense matrix to a random value.
652  for(OrdinalType j = 0; j < numCols_; j++)
653  {
654  for(OrdinalType i = 0; i < numRows_; i++)
655  {
656  values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
657  }
658  }
659  return 0;
660 }
661 
662 template<typename OrdinalType, typename ScalarType>
666  )
667 {
668  if(this == &Source)
669  return(*this); // Special case of source same as target
670  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
671  return(*this); // Special case of both are views to same data.
672 
673  // If the source is a view then we will return a view, else we will return a copy.
674  if (!Source.valuesCopied_) {
675  if(valuesCopied_) {
676  // Clean up stored data if this was previously a copy.
677  deleteArrays();
678  }
679  numRows_ = Source.numRows_;
680  numCols_ = Source.numCols_;
681  stride_ = Source.stride_;
682  values_ = Source.values_;
683  }
684  else {
685  // If we were a view, we will now be a copy.
686  if(!valuesCopied_) {
687  numRows_ = Source.numRows_;
688  numCols_ = Source.numCols_;
689  stride_ = Source.numRows_;
690  const OrdinalType newsize = stride_ * numCols_;
691  if(newsize > 0) {
692  values_ = new ScalarType[newsize];
693  valuesCopied_ = true;
694  }
695  else {
696  values_ = 0;
697  }
698  }
699  // If we were a copy, we will stay a copy.
700  else {
701  if((Source.numRows_ <= stride_) && (Source.numCols_ == numCols_)) { // we don't need to reallocate
702  numRows_ = Source.numRows_;
703  numCols_ = Source.numCols_;
704  }
705  else { // we need to allocate more space (or less space)
706  deleteArrays();
707  numRows_ = Source.numRows_;
708  numCols_ = Source.numCols_;
709  stride_ = Source.numRows_;
710  const OrdinalType newsize = stride_ * numCols_;
711  if(newsize > 0) {
712  values_ = new ScalarType[newsize];
713  valuesCopied_ = true;
714  }
715  }
716  }
717  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
718  }
719  return(*this);
720 }
721 
722 template<typename OrdinalType, typename ScalarType>
724 {
725  // Check for compatible dimensions
726  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
727  {
728  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
729  }
730  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, ScalarTraits<ScalarType>::one());
731  return(*this);
732 }
733 
734 template<typename OrdinalType, typename ScalarType>
736 {
737  // Check for compatible dimensions
738  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
739  {
740  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
741  }
742  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0, -ScalarTraits<ScalarType>::one());
743  return(*this);
744 }
745 
746 template<typename OrdinalType, typename ScalarType>
748  if(this == &Source)
749  return(*this); // Special case of source same as target
750  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_))
751  return(*this); // Special case of both are views to same data.
752 
753  // Check for compatible dimensions
754  if ((numRows_ != Source.numRows_) || (numCols_ != Source.numCols_))
755  {
756  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
757  }
758  copyMat(Source.values_, Source.stride_, numRows_, numCols_, values_, stride_, 0, 0);
759  return(*this);
760 }
761 
762 //----------------------------------------------------------------------------------------------------
763 // Accessor methods
764 //----------------------------------------------------------------------------------------------------
765 
766 template<typename OrdinalType, typename ScalarType>
767 inline ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
768 {
769 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
770  checkIndex( rowIndex, colIndex );
771 #endif
772  return(values_[colIndex * stride_ + rowIndex]);
773 }
774 
775 template<typename OrdinalType, typename ScalarType>
776 inline const ScalarType& SerialDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
777 {
778 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
779  checkIndex( rowIndex, colIndex );
780 #endif
781  return(values_[colIndex * stride_ + rowIndex]);
782 }
783 
784 template<typename OrdinalType, typename ScalarType>
785 inline const ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex) const
786 {
787 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
788  checkIndex( 0, colIndex );
789 #endif
790  return(values_ + colIndex * stride_);
791 }
792 
793 template<typename OrdinalType, typename ScalarType>
794 inline ScalarType* SerialDenseMatrix<OrdinalType, ScalarType>::operator [] (OrdinalType colIndex)
795 {
796 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
797  checkIndex( 0, colIndex );
798 #endif
799  return(values_ + colIndex * stride_);
800 }
801 
802 //----------------------------------------------------------------------------------------------------
803 // Norm methods
804 //----------------------------------------------------------------------------------------------------
805 
806 template<typename OrdinalType, typename ScalarType>
808 {
809  OrdinalType i, j;
812  ScalarType* ptr;
813  for(j = 0; j < numCols_; j++)
814  {
815  ScalarType sum = 0;
816  ptr = values_ + j * stride_;
817  for(i = 0; i < numRows_; i++)
818  {
820  }
822  if(absSum > anorm)
823  {
824  anorm = absSum;
825  }
826  }
827  // don't compute flop increment unless there is an accumulator
828  if (flopCounter_!=0) updateFlops(numRows_ * numCols_);
829  return(anorm);
830 }
831 
832 template<typename OrdinalType, typename ScalarType>
834 {
835  OrdinalType i, j;
837 
838  for (i = 0; i < numRows_; i++) {
840  for (j=0; j< numCols_; j++) {
841  sum += ScalarTraits<ScalarType>::magnitude(*(values_+i+j*stride_));
842  }
843  anorm = TEUCHOS_MAX( anorm, sum );
844  }
845  // don't compute flop increment unless there is an accumulator
846  if (flopCounter_!=0) updateFlops(numRows_ * numCols_);
847  return(anorm);
848 }
849 
850 template<typename OrdinalType, typename ScalarType>
852 {
853  OrdinalType i, j;
855  for (j = 0; j < numCols_; j++) {
856  for (i = 0; i < numRows_; i++) {
857  anorm += ScalarTraits<ScalarType>::magnitude(values_[i+j*stride_]*values_[i+j*stride_]);
858  }
859  }
861  // don't compute flop increment unless there is an accumulator
862  if (flopCounter_!=0) updateFlops(numRows_ * numCols_);
863  return(anorm);
864 }
865 
866 //----------------------------------------------------------------------------------------------------
867 // Comparison methods
868 //----------------------------------------------------------------------------------------------------
869 
870 template<typename OrdinalType, typename ScalarType>
872 {
873  bool result = 1;
874  if((numRows_ != Operand.numRows_) || (numCols_ != Operand.numCols_))
875  {
876  result = 0;
877  }
878  else
879  {
880  OrdinalType i, j;
881  for(i = 0; i < numRows_; i++)
882  {
883  for(j = 0; j < numCols_; j++)
884  {
885  if((*this)(i, j) != Operand(i, j))
886  {
887  return 0;
888  }
889  }
890  }
891  }
892  return result;
893 }
894 
895 template<typename OrdinalType, typename ScalarType>
897 {
898  return !((*this) == Operand);
899 }
900 
901 //----------------------------------------------------------------------------------------------------
902 // Multiplication method
903 //----------------------------------------------------------------------------------------------------
904 
905 template<typename OrdinalType, typename ScalarType>
907 {
908  this->scale( alpha );
909  return(*this);
910 }
911 
912 template<typename OrdinalType, typename ScalarType>
914 {
915  OrdinalType i, j;
916  ScalarType* ptr;
917 
918  for (j=0; j<numCols_; j++) {
919  ptr = values_ + j*stride_;
920  for (i=0; i<numRows_; i++) { *ptr = alpha * (*ptr); ptr++; }
921  }
922  // don't compute flop increment unless there is an accumulator
923  if (flopCounter_!=0) updateFlops( numRows_*numCols_ );
924  return(0);
925 }
926 
927 template<typename OrdinalType, typename ScalarType>
929 {
930  OrdinalType i, j;
931  ScalarType* ptr;
932 
933  // Check for compatible dimensions
934  if ((numRows_ != A.numRows_) || (numCols_ != A.numCols_))
935  {
936  TEUCHOS_CHK_ERR(-1); // Return error
937  }
938  for (j=0; j<numCols_; j++) {
939  ptr = values_ + j*stride_;
940  for (i=0; i<numRows_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
941  }
942  // don't compute flop increment unless there is an accumulator
943  if (flopCounter_!=0) updateFlops( numRows_*numCols_ );
944  return(0);
945 }
946 
947 template<typename OrdinalType, typename ScalarType>
949 {
950  // Check for compatible dimensions
951  OrdinalType A_nrows = (ETranspChar[transa]!='N') ? A.numCols() : A.numRows();
952  OrdinalType A_ncols = (ETranspChar[transa]!='N') ? A.numRows() : A.numCols();
953  OrdinalType B_nrows = (ETranspChar[transb]!='N') ? B.numCols() : B.numRows();
954  OrdinalType B_ncols = (ETranspChar[transb]!='N') ? B.numRows() : B.numCols();
955  if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols))
956  {
957  TEUCHOS_CHK_ERR(-1); // Return error
958  }
959  // Call GEMM function
960  this->GEMM(transa, transb, numRows_, numCols_, A_ncols, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
961 
962  // don't compute flop increment unless there is an accumulator
963  if (flopCounter_!=0) {
964  double nflops = 2 * numRows_;
965  nflops *= numCols_;
966  nflops *= A_ncols;
967  updateFlops(nflops);
968  }
969  return(0);
970 }
971 
972 template<typename OrdinalType, typename ScalarType>
974 {
975  // Check for compatible dimensions
976  OrdinalType A_nrows = A.numRows(), A_ncols = A.numCols();
977  OrdinalType B_nrows = B.numRows(), B_ncols = B.numCols();
978 
979  if (ESideChar[sideA]=='L') {
980  if ((numRows_ != A_nrows) || (A_ncols != B_nrows) || (numCols_ != B_ncols)) {
981  TEUCHOS_CHK_ERR(-1); // Return error
982  }
983  } else {
984  if ((numRows_ != B_nrows) || (B_ncols != A_nrows) || (numCols_ != A_ncols)) {
985  TEUCHOS_CHK_ERR(-1); // Return error
986  }
987  }
988 
989  // Call SYMM function
991  this->SYMM(sideA, uplo, numRows_, numCols_, alpha, A.values(), A.stride(), B.values(), B.stride(), beta, values_, stride_);
992 
993  // don't compute flop increment unless there is an accumulator
994  if (flopCounter_!=0) {
995  double nflops = 2 * numRows_;
996  nflops *= numCols_;
997  nflops *= A_ncols;
998  updateFlops(nflops);
999  }
1000  return(0);
1001 }
1002 
1003 template<typename OrdinalType, typename ScalarType>
1005 {
1006  os << std::endl;
1007  if(valuesCopied_)
1008  os << "Values_copied : yes" << std::endl;
1009  else
1010  os << "Values_copied : no" << std::endl;
1011  os << "Rows : " << numRows_ << std::endl;
1012  os << "Columns : " << numCols_ << std::endl;
1013  os << "LDA : " << stride_ << std::endl;
1014  if(numRows_ == 0 || numCols_ == 0) {
1015  os << "(matrix is empty, no values to display)" << std::endl;
1016  } else {
1017  for(OrdinalType i = 0; i < numRows_; i++) {
1018  for(OrdinalType j = 0; j < numCols_; j++){
1019  os << (*this)(i,j) << " ";
1020  }
1021  os << std::endl;
1022  }
1023  }
1024 }
1025 
1026 //----------------------------------------------------------------------------------------------------
1027 // Protected methods
1028 //----------------------------------------------------------------------------------------------------
1029 
1030 template<typename OrdinalType, typename ScalarType>
1031 inline void SerialDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
1032  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRows_, std::out_of_range,
1033  "SerialDenseMatrix<T>::checkIndex: "
1034  "Row index " << rowIndex << " out of range [0, "<< numRows_ << ")");
1035  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numCols_, std::out_of_range,
1036  "SerialDenseMatrix<T>::checkIndex: "
1037  "Col index " << colIndex << " out of range [0, "<< numCols_ << ")");
1038 }
1039 
1040 template<typename OrdinalType, typename ScalarType>
1042 {
1043  if (valuesCopied_)
1044  {
1045  delete [] values_;
1046  values_ = 0;
1047  valuesCopied_ = false;
1048  }
1049 }
1050 
1051 template<typename OrdinalType, typename ScalarType>
1053  ScalarType* inputMatrix, OrdinalType strideInput, OrdinalType numRows_in,
1054  OrdinalType numCols_in, ScalarType* outputMatrix, OrdinalType strideOutput,
1055  OrdinalType startRow, OrdinalType startCol, ScalarType alpha
1056  )
1057 {
1058  OrdinalType i, j;
1059  ScalarType* ptr1 = 0;
1060  ScalarType* ptr2 = 0;
1061  for(j = 0; j < numCols_in; j++) {
1062  ptr1 = outputMatrix + (j * strideOutput);
1063  ptr2 = inputMatrix + (j + startCol) * strideInput + startRow;
1064  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1065  for(i = 0; i < numRows_in; i++)
1066  {
1067  *ptr1++ += alpha*(*ptr2++);
1068  }
1069  } else {
1070  for(i = 0; i < numRows_in; i++)
1071  {
1072  *ptr1++ = *ptr2++;
1073  }
1074  }
1075  }
1076 }
1077 
1078 #ifndef TEUCHOS_HIDE_DEPRECATED_CODE
1079 template<typename OrdinalType, typename ScalarType>
1081 std::ostream& operator<< (std::ostream& os, const Teuchos::SerialDenseMatrix<OrdinalType, ScalarType>& obj)
1082 {
1083  obj.print (os);
1084  return os;
1085 }
1086 #endif
1087 
1088 } // namespace Teuchos
1089 
1090 
1091 #endif /* _TEUCHOS_SERIALDENSEMATRIX_HPP_ */
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
ScalarType * values() const
Data array access method.
bool upper() const
Returns true if upper triangular part of this matrix has and will be used.
void checkIndex(OrdinalType rowIndex, OrdinalType colIndex=0) const
bool operator==(const SerialDenseMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
int shapeUninitialized(OrdinalType numRows, OrdinalType numCols)
Same as shape() except leaves uninitialized.
virtual void print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream &lt;&lt; operator.
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.
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.
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.
#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.
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.
#define TEUCHOS_MAX(x, y)
SerialDenseMatrix()
Default Constructor.
bool operator!=(const SerialDenseMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ESideChar[]
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
void copyMat(ScalarType *inputMatrix, OrdinalType strideInput, OrdinalType numRows, OrdinalType numCols, ScalarType *outputMatrix, OrdinalType strideOutput, OrdinalType startRow, OrdinalType startCol, ScalarType alpha=ScalarTraits< ScalarType >::zero())
Templated serial, dense, symmetric matrix class.
OrdinalType ordinalType
Typedef for ordinal type.
ScalarType * operator[](OrdinalType colIndex)
Column access method (non-const).
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
OrdinalType numCols() const
Returns the column dimension of this matrix.
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
int random()
Set all values in the matrix to be random numbers.
int reshape(OrdinalType numRows, OrdinalType numCols)
Reshaping method for changing the size of a SerialDenseMatrix, keeping the entries.
Defines basic traits for the scalar field type.
static T zero()
Returns representation of zero for this scalar type.
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ETranspChar[]
static T random()
Returns a random number (between -one() and +one()) of this scalar type.
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
int shape(OrdinalType numRows, OrdinalType numCols)
Shape method for changing the size of a SerialDenseMatrix, initializing entries to zero...
OrdinalType numCols() const
Returns the column dimension of this matrix.
#define TEUCHOS_CHK_REF(a)
#define TEUCHOS_MIN(x, y)
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
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...