Teuchos - Trilinos Tools Package  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Teuchos_SerialSymDenseMatrix.hpp
Go to the documentation of this file.
1 
2 // @HEADER
3 // *****************************************************************************
4 // Teuchos: Common Tools Package
5 //
6 // Copyright 2004 NTESS and the Teuchos contributors.
7 // SPDX-License-Identifier: BSD-3-Clause
8 // *****************************************************************************
9 // @HEADER
10 
11 #ifndef _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_
12 #define _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_
13 
17 #include "Teuchos_CompObject.hpp"
18 #include "Teuchos_BLAS.hpp"
19 #include "Teuchos_ScalarTraits.hpp"
20 #include "Teuchos_DataAccess.hpp"
21 #include "Teuchos_ConfigDefs.hpp"
22 #include "Teuchos_Assert.hpp"
23 #include <utility>
24 #include <vector>
25 
88 namespace Teuchos {
89 
90 template<typename OrdinalType, typename ScalarType>
91 class SerialSymDenseMatrix : public CompObject, public BLAS<OrdinalType,ScalarType>
92 {
93  public:
94 
96  typedef OrdinalType ordinalType;
98  typedef ScalarType scalarType;
99 
101 
102 
111  SerialSymDenseMatrix() = default;
112 
114 
124  SerialSymDenseMatrix(OrdinalType numRowsCols, bool zeroOut = true);
125 
127 
138  SerialSymDenseMatrix(DataAccess CV, bool upper, ScalarType* values, OrdinalType stride, OrdinalType numRowsCols);
139 
142 
144 
153  SerialSymDenseMatrix(DataAccess CV, const SerialSymDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRowCols, OrdinalType startRowCol=0);
154 
156  virtual ~SerialSymDenseMatrix ();
158 
160 
161 
163 
172  int shape(OrdinalType numRowsCols);
173 
175 
184  int shapeUninitialized(OrdinalType numRowsCols);
185 
187 
197  int reshape(OrdinalType numRowsCols);
198 
200 
202  void setLower();
203 
205 
207  void setUpper();
208 
210 
212 
213 
215 
222 
224 
229 
231 
234  SerialSymDenseMatrix<OrdinalType, ScalarType>& operator= (const ScalarType value) { putScalar(value); return(*this); }
235 
237 
242  int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero(), bool fullMatrix = false );
243 
245 
249 
251 
253  int random( const ScalarType bias = 0.1*Teuchos::ScalarTraits<ScalarType>::one() );
254 
256 
258 
259 
261 
268  ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
269 
271 
278  const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
279 
281 
283  ScalarType* values() const { return(values_); }
284 
286 
288 
289 
291  bool upper() const {return(upper_);};
292 
294  char UPLO() const {return(UPLO_);};
296 
298 
299 
301 
308 
310 
314 
316 
320 
322 
324 
325 
327 
331 
333 
337 
339 
341 
342 
344  OrdinalType numRows() const { return(numRowCols_); }
345 
347  OrdinalType numCols() const { return(numRowCols_); }
348 
350  OrdinalType stride() const { return(stride_); }
351 
353  bool empty() const { return(numRowCols_ == 0); }
354 
356 
358 
359 
362 
365 
369 
371 
372  virtual std::ostream& print(std::ostream& os) const;
374 
376 
377  protected:
378 
379  // In-place scaling of matrix.
380  void scale( const ScalarType alpha );
381 
382  // Copy the values from one matrix to the other.
383  void copyMat(bool inputUpper, ScalarType* inputMatrix, OrdinalType inputStride,
384  OrdinalType numRowCols, bool outputUpper, ScalarType* outputMatrix,
385  OrdinalType outputStride, OrdinalType startRowCol,
386  ScalarType alpha = ScalarTraits<ScalarType>::zero() );
387 
388  // Copy the values from the active triangle of the matrix to the other to make the matrix full symmetric.
389  void copyUPLOMat(bool inputUpper, ScalarType* inputMatrix,
390  OrdinalType inputStride, OrdinalType inputRows);
391 
392  void deleteArrays();
393  void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
394 
395  static ScalarType*
396  allocateValues(const OrdinalType numRows,
397  const OrdinalType numCols)
398  {
399  const size_t size = size_t(numRows) * size_t(numCols);
400 #pragma GCC diagnostic push
401 #pragma GCC diagnostic ignored "-Wvla"
402  return new ScalarType[size];
403 #pragma GCC diagnostic pop
404  }
405 
406  OrdinalType numRowCols_ = 0;
407  OrdinalType stride_ = 0;
408  bool valuesCopied_ = false;
409  ScalarType* values_ = nullptr;
410  bool upper_ = false;
411  char UPLO_ {'L'};
412 };
413 
414 //----------------------------------------------------------------------------------------------------
415 // Constructors and Destructor
416 //----------------------------------------------------------------------------------------------------
417 
418 template<typename OrdinalType, typename ScalarType>
420  : numRowCols_(numRowCols_in), stride_(numRowCols_in), valuesCopied_(false), upper_(false), UPLO_('L')
421 {
422  values_ = allocateValues(stride_, numRowCols_);
423  valuesCopied_ = true;
424  if (zeroOut) {
425  const ScalarType ZERO = Teuchos::ScalarTraits<ScalarType>::zero();
426  putScalar(ZERO, true);
427  }
428 }
429 
430 template<typename OrdinalType, typename ScalarType>
432  DataAccess CV, bool upper_in, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRowCols_in
433  )
434  : numRowCols_(numRowCols_in), stride_(stride_in), valuesCopied_(false),
435  values_(values_in), upper_(upper_in)
436 {
437  if (upper_)
438  UPLO_ = 'U';
439  else
440  UPLO_ = 'L';
441 
442  if(CV == Copy)
443  {
444  stride_ = numRowCols_;
445  values_ = allocateValues(stride_, numRowCols_);
446  copyMat(upper_in, values_in, stride_in, numRowCols_, upper_, values_, stride_, 0);
447  valuesCopied_ = true;
448  }
449 }
450 
451 template<typename OrdinalType, typename ScalarType>
453  : numRowCols_(Source.numRowCols_), stride_(0), valuesCopied_(true),
454  values_(0), upper_(Source.upper_), UPLO_(Source.UPLO_)
455 {
456  if (!Source.valuesCopied_)
457  {
458  stride_ = Source.stride_;
459  values_ = Source.values_;
460  valuesCopied_ = false;
461  }
462  else
463  {
464  stride_ = numRowCols_;
465  if(stride_ > 0 && numRowCols_ > 0) {
466  values_ = allocateValues(stride_, numRowCols_);
467  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0);
468  }
469  else {
470  numRowCols_ = 0;
471  stride_ = 0;
472  valuesCopied_ = false;
473  }
474  }
475 }
476 
477 template<typename OrdinalType, typename ScalarType>
479  DataAccess CV, const SerialSymDenseMatrix<OrdinalType,
480  ScalarType> &Source, OrdinalType numRowCols_in, OrdinalType startRowCol )
481  : CompObject(), BLAS<OrdinalType,ScalarType>(),
482  numRowCols_(numRowCols_in), stride_(Source.stride_), valuesCopied_(false), upper_(Source.upper_), UPLO_(Source.UPLO_)
483 {
484  if(CV == Copy)
485  {
486  stride_ = numRowCols_in;
487  values_ = allocateValues(stride_, numRowCols_in);
488  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_in, upper_, values_, stride_, startRowCol);
489  valuesCopied_ = true;
490  }
491  else // CV == View
492  {
493  values_ = Source.values_ + (stride_ * startRowCol) + startRowCol;
494  }
495 }
496 
497 template<typename OrdinalType, typename ScalarType>
499 {
500  deleteArrays();
501 }
502 
503 //----------------------------------------------------------------------------------------------------
504 // Shape methods
505 //----------------------------------------------------------------------------------------------------
506 
507 template<typename OrdinalType, typename ScalarType>
509 {
510  deleteArrays(); // Get rid of anything that might be already allocated
511  numRowCols_ = numRowCols_in;
512  stride_ = numRowCols_;
513  values_ = allocateValues(stride_, numRowCols_);
514  putScalar( Teuchos::ScalarTraits<ScalarType>::zero(), true );
515  valuesCopied_ = true;
516  return(0);
517 }
518 
519 template<typename OrdinalType, typename ScalarType>
521 {
522  deleteArrays(); // Get rid of anything that might be already allocated
523  numRowCols_ = numRowCols_in;
524  stride_ = numRowCols_;
525  values_ = allocateValues(stride_, numRowCols_);
526  valuesCopied_ = true;
527  return(0);
528 }
529 
530 template<typename OrdinalType, typename ScalarType>
532 {
533  // Allocate space for new matrix
534  ScalarType* values_tmp = allocateValues(numRowCols_in, numRowCols_in);
535  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
536  for(OrdinalType k = 0; k < numRowCols_in * numRowCols_in; k++)
537  {
538  values_tmp[k] = zero;
539  }
540  OrdinalType numRowCols_tmp = TEUCHOS_MIN(numRowCols_, numRowCols_in);
541  if(values_ != 0)
542  {
543  copyMat(upper_, values_, stride_, numRowCols_tmp, upper_, values_tmp, numRowCols_in, 0); // Copy principal submatrix of A to new A
544  }
545  deleteArrays(); // Get rid of anything that might be already allocated
546  numRowCols_ = numRowCols_in;
547  stride_ = numRowCols_;
548  values_ = values_tmp; // Set pointer to new A
549  valuesCopied_ = true;
550  return(0);
551 }
552 
553 //----------------------------------------------------------------------------------------------------
554 // Set methods
555 //----------------------------------------------------------------------------------------------------
556 
557 template<typename OrdinalType, typename ScalarType>
559 {
560  // Do nothing if the matrix is already a lower triangular matrix
561  if (upper_ != false) {
562  copyUPLOMat( true, values_, stride_, numRowCols_ );
563  upper_ = false;
564  UPLO_ = 'L';
565  }
566 }
567 
568 template<typename OrdinalType, typename ScalarType>
570 {
571  // Do nothing if the matrix is already an upper triangular matrix
572  if (upper_ == false) {
573  copyUPLOMat( false, values_, stride_, numRowCols_ );
574  upper_ = true;
575  UPLO_ = 'U';
576  }
577 }
578 
579 template<typename OrdinalType, typename ScalarType>
580 int SerialSymDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in, bool fullMatrix )
581 {
582  // Set each value of the dense matrix to "value".
583  if (fullMatrix == true) {
584  for(OrdinalType j = 0; j < numRowCols_; j++)
585  {
586  for(OrdinalType i = 0; i < numRowCols_; i++)
587  {
588  values_[i + j*stride_] = value_in;
589  }
590  }
591  }
592  // Set the active upper or lower triangular part of the matrix to "value"
593  else {
594  if (upper_) {
595  for(OrdinalType j = 0; j < numRowCols_; j++) {
596  for(OrdinalType i = 0; i <= j; i++) {
597  values_[i + j*stride_] = value_in;
598  }
599  }
600  }
601  else {
602  for(OrdinalType j = 0; j < numRowCols_; j++) {
603  for(OrdinalType i = j; i < numRowCols_; i++) {
604  values_[i + j*stride_] = value_in;
605  }
606  }
607  }
608  }
609  return 0;
610 }
611 
612 template<typename OrdinalType, typename ScalarType> void
615 {
616  std::swap(values_ , B.values_);
617  std::swap(numRowCols_, B.numRowCols_);
618  std::swap(stride_, B.stride_);
619  std::swap(valuesCopied_, B.valuesCopied_);
620  std::swap(upper_, B.upper_);
621  std::swap(UPLO_, B.UPLO_);
622 }
623 
624 template<typename OrdinalType, typename ScalarType>
626 {
628 
629  // Set each value of the dense matrix to a random value.
630  std::vector<MT> diagSum( numRowCols_, Teuchos::ScalarTraits<MT>::zero() );
631  if (upper_) {
632  for(OrdinalType j = 0; j < numRowCols_; j++) {
633  for(OrdinalType i = 0; i < j; i++) {
634  values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
635  diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
636  diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
637  }
638  }
639  }
640  else {
641  for(OrdinalType j = 0; j < numRowCols_; j++) {
642  for(OrdinalType i = j+1; i < numRowCols_; i++) {
643  values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
644  diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
645  diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
646  }
647  }
648  }
649 
650  // Set the diagonal.
651  for(OrdinalType i = 0; i < numRowCols_; i++) {
652  values_[i + i*stride_] = diagSum[i] + bias;
653  }
654  return 0;
655 }
656 
657 template<typename OrdinalType, typename ScalarType>
660 {
661  if(this == &Source)
662  return(*this); // Special case of source same as target
663  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
664  upper_ = Source.upper_; // Might have to change the active part of the matrix.
665  return(*this); // Special case of both are views to same data.
666  }
667 
668  // If the source is a view then we will return a view, else we will return a copy.
669  if (!Source.valuesCopied_) {
670  if(valuesCopied_) {
671  // Clean up stored data if this was previously a copy.
672  deleteArrays();
673  }
674  numRowCols_ = Source.numRowCols_;
675  stride_ = Source.stride_;
676  values_ = Source.values_;
677  upper_ = Source.upper_;
678  UPLO_ = Source.UPLO_;
679  }
680  else {
681  // If we were a view, we will now be a copy.
682  if(!valuesCopied_) {
683  numRowCols_ = Source.numRowCols_;
684  stride_ = Source.numRowCols_;
685  upper_ = Source.upper_;
686  UPLO_ = Source.UPLO_;
687  if(stride_ > 0 && numRowCols_ > 0) {
688  values_ = allocateValues(stride_, numRowCols_);
689  valuesCopied_ = true;
690  }
691  else {
692  values_ = 0;
693  }
694  }
695  // If we were a copy, we will stay a copy.
696  else {
697  if((Source.numRowCols_ <= stride_) && (Source.numRowCols_ == numRowCols_)) { // we don't need to reallocate
698  numRowCols_ = Source.numRowCols_;
699  upper_ = Source.upper_;
700  UPLO_ = Source.UPLO_;
701  }
702  else { // we need to allocate more space (or less space)
703  deleteArrays();
704  numRowCols_ = Source.numRowCols_;
705  stride_ = Source.numRowCols_;
706  upper_ = Source.upper_;
707  UPLO_ = Source.UPLO_;
708  if(stride_ > 0 && numRowCols_ > 0) {
709  values_ = allocateValues(stride_, numRowCols_);
710  valuesCopied_ = true;
711  }
712  }
713  }
714  copyMat(Source.upper_, Source.values_, Source.stride_, Source.numRowCols_, upper_, values_, stride_, 0);
715  }
716  return(*this);
717 }
718 
719 template<typename OrdinalType, typename ScalarType>
721 {
722  this->scale(alpha);
723  return(*this);
724 }
725 
726 template<typename OrdinalType, typename ScalarType>
728 {
729  // Check for compatible dimensions
730  if ((numRowCols_ != Source.numRowCols_))
731  {
732  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
733  }
734  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0, ScalarTraits<ScalarType>::one());
735  return(*this);
736 }
737 
738 template<typename OrdinalType, typename ScalarType>
740 {
741  // Check for compatible dimensions
742  if ((numRowCols_ != Source.numRowCols_))
743  {
744  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
745  }
746  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0, -ScalarTraits<ScalarType>::one());
747  return(*this);
748 }
749 
750 template<typename OrdinalType, typename ScalarType>
752  if(this == &Source)
753  return(*this); // Special case of source same as target
754  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
755  upper_ = Source.upper_; // We may have to change the active part of the matrix.
756  return(*this); // Special case of both are views to same data.
757  }
758 
759  // Check for compatible dimensions
760  if ((numRowCols_ != Source.numRowCols_))
761  {
762  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
763  }
764  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0 );
765  return(*this);
766 }
767 
768 //----------------------------------------------------------------------------------------------------
769 // Accessor methods
770 //----------------------------------------------------------------------------------------------------
771 
772 template<typename OrdinalType, typename ScalarType>
773 inline ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
774 {
775 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
776  checkIndex( rowIndex, colIndex );
777 #endif
778  if ( rowIndex <= colIndex ) {
779  // Accessing upper triangular part of matrix
780  if (upper_)
781  return(values_[colIndex * stride_ + rowIndex]);
782  else
783  return(values_[rowIndex * stride_ + colIndex]);
784  }
785  else {
786  // Accessing lower triangular part of matrix
787  if (upper_)
788  return(values_[rowIndex * stride_ + colIndex]);
789  else
790  return(values_[colIndex * stride_ + rowIndex]);
791  }
792 }
793 
794 template<typename OrdinalType, typename ScalarType>
795 inline const ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
796 {
797 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
798  checkIndex( rowIndex, colIndex );
799 #endif
800  if ( rowIndex <= colIndex ) {
801  // Accessing upper triangular part of matrix
802  if (upper_)
803  return(values_[colIndex * stride_ + rowIndex]);
804  else
805  return(values_[rowIndex * stride_ + colIndex]);
806  }
807  else {
808  // Accessing lower triangular part of matrix
809  if (upper_)
810  return(values_[rowIndex * stride_ + colIndex]);
811  else
812  return(values_[colIndex * stride_ + rowIndex]);
813  }
814 }
815 
816 //----------------------------------------------------------------------------------------------------
817 // Norm methods
818 //----------------------------------------------------------------------------------------------------
819 
820 template<typename OrdinalType, typename ScalarType>
822 {
823  return(normInf());
824 }
825 
826 template<typename OrdinalType, typename ScalarType>
828 {
829  typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
830 
831  OrdinalType i, j;
832 
833  MT sum, anorm = ScalarTraits<MT>::zero();
834  ScalarType* ptr;
835 
836  if (upper_) {
837  for (j=0; j<numRowCols_; j++) {
838  sum = ScalarTraits<MT>::zero();
839  ptr = values_ + j*stride_;
840  for (i=0; i<j; i++) {
841  sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
842  }
843  ptr = values_ + j + j*stride_;
844  for (i=j; i<numRowCols_; i++) {
846  ptr += stride_;
847  }
848  anorm = TEUCHOS_MAX( anorm, sum );
849  }
850  }
851  else {
852  for (j=0; j<numRowCols_; j++) {
853  sum = ScalarTraits<MT>::zero();
854  ptr = values_ + j + j*stride_;
855  for (i=j; i<numRowCols_; i++) {
856  sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
857  }
858  ptr = values_ + j;
859  for (i=0; i<j; i++) {
861  ptr += stride_;
862  }
863  anorm = TEUCHOS_MAX( anorm, sum );
864  }
865  }
866  return(anorm);
867 }
868 
869 template<typename OrdinalType, typename ScalarType>
871 {
872  typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
873 
874  OrdinalType i, j;
875 
876  MT sum = ScalarTraits<MT>::zero(), anorm = ScalarTraits<MT>::zero();
877 
878  if (upper_) {
879  for (j = 0; j < numRowCols_; j++) {
880  for (i = 0; i < j; i++) {
881  sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
882  }
883  sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
884  }
885  }
886  else {
887  for (j = 0; j < numRowCols_; j++) {
888  sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
889  for (i = j+1; i < numRowCols_; i++) {
890  sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
891  }
892  }
893  }
895  return(anorm);
896 }
897 
898 //----------------------------------------------------------------------------------------------------
899 // Comparison methods
900 //----------------------------------------------------------------------------------------------------
901 
902 template<typename OrdinalType, typename ScalarType>
904 {
905  bool result = 1;
906  if((numRowCols_ != Operand.numRowCols_)) {
907  result = 0;
908  }
909  else {
910  OrdinalType i, j;
911  for(i = 0; i < numRowCols_; i++) {
912  for(j = 0; j < numRowCols_; j++) {
913  if((*this)(i, j) != Operand(i, j)) {
914  return 0;
915  }
916  }
917  }
918  }
919  return result;
920 }
921 
922 template<typename OrdinalType, typename ScalarType>
924 {
925  return !((*this) == Operand);
926 }
927 
928 //----------------------------------------------------------------------------------------------------
929 // Multiplication method
930 //----------------------------------------------------------------------------------------------------
931 
932 template<typename OrdinalType, typename ScalarType>
933 void SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
934 {
935  OrdinalType i, j;
936  ScalarType* ptr;
937 
938  if (upper_) {
939  for (j=0; j<numRowCols_; j++) {
940  ptr = values_ + j*stride_;
941  for (i=0; i<= j; i++) { *ptr = alpha * (*ptr); ptr++; }
942  }
943  }
944  else {
945  for (j=0; j<numRowCols_; j++) {
946  ptr = values_ + j*stride_ + j;
947  for (i=j; i<numRowCols_; i++) { *ptr = alpha * (*ptr); ptr++; }
948  }
949  }
950 }
951 
952 /*
953 template<typename OrdinalType, typename ScalarType>
954 int SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const SerialSymDenseMatrix<OrdinalType,ScalarType>& A )
955 {
956  OrdinalType i, j;
957  ScalarType* ptr;
958 
959  // Check for compatible dimensions
960  if ((numRowCols_ != A.numRowCols_)) {
961  TEUCHOS_CHK_ERR(-1); // Return error
962  }
963 
964  if (upper_) {
965  for (j=0; j<numRowCols_; j++) {
966  ptr = values_ + j*stride_;
967  for (i=0; i<= j; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
968  }
969  }
970  else {
971  for (j=0; j<numRowCols_; j++) {
972  ptr = values_ + j*stride_;
973  for (i=j; i<numRowCols_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
974  }
975  }
976 
977  return(0);
978 }
979 */
980 
981 template<typename OrdinalType, typename ScalarType>
982 std::ostream& SerialSymDenseMatrix<OrdinalType, ScalarType>::print(std::ostream& os) const
983 {
984  os << std::endl;
985  if(valuesCopied_)
986  os << "Values_copied : yes" << std::endl;
987  else
988  os << "Values_copied : no" << std::endl;
989  os << "Rows / Columns : " << numRowCols_ << std::endl;
990  os << "LDA : " << stride_ << std::endl;
991  if (upper_)
992  os << "Storage: Upper" << std::endl;
993  else
994  os << "Storage: Lower" << std::endl;
995  if(numRowCols_ == 0) {
996  os << "(matrix is empty, no values to display)" << std::endl;
997  } else {
998  for(OrdinalType i = 0; i < numRowCols_; i++) {
999  for(OrdinalType j = 0; j < numRowCols_; j++){
1000  os << (*this)(i,j) << " ";
1001  }
1002  os << std::endl;
1003  }
1004  }
1005  return os;
1006 }
1007 
1008 //----------------------------------------------------------------------------------------------------
1009 // Protected methods
1010 //----------------------------------------------------------------------------------------------------
1011 
1012 template<typename OrdinalType, typename ScalarType>
1013 inline void SerialSymDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
1014  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRowCols_, std::out_of_range,
1015  "SerialSymDenseMatrix<T>::checkIndex: "
1016  "Row index " << rowIndex << " out of range [0, "<< numRowCols_ << ")");
1017  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numRowCols_, std::out_of_range,
1018  "SerialSymDenseMatrix<T>::checkIndex: "
1019  "Col index " << colIndex << " out of range [0, "<< numRowCols_ << ")");
1020 }
1021 
1022 template<typename OrdinalType, typename ScalarType>
1023 void SerialSymDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
1024 {
1025  if (valuesCopied_)
1026  {
1027  delete [] values_;
1028  values_ = 0;
1029  valuesCopied_ = false;
1030  }
1031 }
1032 
1033 template<typename OrdinalType, typename ScalarType>
1034 void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyMat(
1035  bool inputUpper, ScalarType* inputMatrix,
1036  OrdinalType inputStride, OrdinalType numRowCols_in,
1037  bool outputUpper, ScalarType* outputMatrix,
1038  OrdinalType outputStride, OrdinalType startRowCol,
1039  ScalarType alpha
1040  )
1041 {
1042  OrdinalType i, j;
1043  ScalarType* ptr1 = 0;
1044  ScalarType* ptr2 = 0;
1045 
1046  for (j = 0; j < numRowCols_in; j++) {
1047  if (inputUpper == true) {
1048  // The input matrix is upper triangular, start at the beginning of each column.
1049  ptr2 = inputMatrix + (j + startRowCol) * inputStride + startRowCol;
1050  if (outputUpper == true) {
1051  // The output matrix matches the same pattern as the input matrix.
1052  ptr1 = outputMatrix + j*outputStride;
1053  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1054  for(i = 0; i <= j; i++) {
1055  *ptr1++ += alpha*(*ptr2++);
1056  }
1057  } else {
1058  for(i = 0; i <= j; i++) {
1059  *ptr1++ = *ptr2++;
1060  }
1061  }
1062  }
1063  else {
1064  // The output matrix has the opposite pattern as the input matrix.
1065  // Copy down across rows of the output matrix, but down columns of the input matrix.
1066  ptr1 = outputMatrix + j;
1067  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1068  for(i = 0; i <= j; i++) {
1069  *ptr1 += alpha*(*ptr2++);
1070  ptr1 += outputStride;
1071  }
1072  } else {
1073  for(i = 0; i <= j; i++) {
1074  *ptr1 = *ptr2++;
1075  ptr1 += outputStride;
1076  }
1077  }
1078  }
1079  }
1080  else {
1081  // The input matrix is lower triangular, start at the diagonal of each row.
1082  ptr2 = inputMatrix + (startRowCol+j) * inputStride + startRowCol + j;
1083  if (outputUpper == true) {
1084  // The output matrix has the opposite pattern as the input matrix.
1085  // Copy across rows of the output matrix, but down columns of the input matrix.
1086  ptr1 = outputMatrix + j*outputStride + j;
1087  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1088  for(i = j; i < numRowCols_in; i++) {
1089  *ptr1 += alpha*(*ptr2++);
1090  ptr1 += outputStride;
1091  }
1092  } else {
1093  for(i = j; i < numRowCols_in; i++) {
1094  *ptr1 = *ptr2++;
1095  ptr1 += outputStride;
1096  }
1097  }
1098  }
1099  else {
1100  // The output matrix matches the same pattern as the input matrix.
1101  ptr1 = outputMatrix + j*outputStride + j;
1102  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1103  for(i = j; i < numRowCols_in; i++) {
1104  *ptr1++ += alpha*(*ptr2++);
1105  }
1106  } else {
1107  for(i = j; i < numRowCols_in; i++) {
1108  *ptr1++ = *ptr2++;
1109  }
1110  }
1111  }
1112  }
1113  }
1114 }
1115 
1116 template<typename OrdinalType, typename ScalarType>
1117 void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyUPLOMat(
1118  bool inputUpper, ScalarType* inputMatrix,
1119  OrdinalType inputStride, OrdinalType inputRows
1120  )
1121 {
1122  OrdinalType i, j;
1123  ScalarType * ptr1 = 0;
1124  ScalarType * ptr2 = 0;
1125 
1126  if (inputUpper) {
1127  for (j=1; j<inputRows; j++) {
1128  ptr1 = inputMatrix + j;
1129  ptr2 = inputMatrix + j*inputStride;
1130  for (i=0; i<j; i++) {
1131  *ptr1 = *ptr2++;
1132  ptr1+=inputStride;
1133  }
1134  }
1135  }
1136  else {
1137  for (i=1; i<inputRows; i++) {
1138  ptr1 = inputMatrix + i;
1139  ptr2 = inputMatrix + i*inputStride;
1140  for (j=0; j<i; j++) {
1141  *ptr2++ = *ptr1;
1142  ptr1+=inputStride;
1143  }
1144  }
1145  }
1146 }
1147 
1149 template<typename OrdinalType, typename ScalarType>
1151 public:
1155  : obj(obj_in) {}
1156 };
1157 
1159 template<typename OrdinalType, typename ScalarType>
1160 std::ostream&
1161 operator<<(std::ostream &out,
1163 {
1164  printer.obj.print(out);
1165  return out;
1166 }
1167 
1169 template<typename OrdinalType, typename ScalarType>
1170 SerialSymDenseMatrixPrinter<OrdinalType,ScalarType>
1172 {
1174 }
1175 
1176 
1177 } // namespace Teuchos
1178 
1179 #endif /* _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_ */
bool upper() const
Returns true if upper triangular part of this matrix has and will be used.
Ostream manipulator for SerialSymDenseMatrix.
SerialSymDenseMatrix()=default
Default constructor; defines a zero size object.
T magnitudeType
Mandatory typedef for result of magnitude.
Templated interface class to BLAS routines.
virtual std::ostream & print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream &lt;&lt; operator.
SerialBandDenseMatrixPrinter< OrdinalType, ScalarType > printMat(const SerialBandDenseMatrix< OrdinalType, ScalarType > &obj)
Return SerialBandDenseMatrix ostream manipulator Use as:
ScalarType * values() const
Returns the pointer to the ScalarType data array contained in the object.
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
void swap(SerialSymDenseMatrix< OrdinalType, ScalarType > &B)
Swap values between this matrix and incoming matrix.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
int shape(OrdinalType numRowsCols)
Set dimensions of a Teuchos::SerialSymDenseMatrix object; init values to zero.
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...
OrdinalType numRows() const
Returns the row dimension of this matrix.
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
Object for storing data and providing functionality that is common to all computational classes...
int reshape(OrdinalType numRowsCols)
Reshape a Teuchos::SerialSymDenseMatrix object.
This structure defines some basic traits for a scalar field type.
void setLower()
Specify that the lower triangle of the this matrix should be used.
void setUpper()
Specify that the upper triangle of the this matrix should be used.
Functionality and data that is common to all computational classes.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator*=(const ScalarType alpha)
Inplace scalar-matrix product A = alpha*A.
bool operator==(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero(), bool fullMatrix=false)
Set all values in the matrix to a constant value.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator+=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Add another matrix to this matrix.
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
int random(const ScalarType bias=0.1 *Teuchos::ScalarTraits< ScalarType >::one())
Set all values in the active area (upper/lower triangle) of this matrix to be random numbers...
char UPLO() const
Returns character value of UPLO used by LAPACK routines.
virtual ~SerialSymDenseMatrix()
Teuchos::SerialSymDenseMatrix destructor.
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator-=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Subtract another matrix from this matrix.
OrdinalType ordinalType
Typedef for ordinal type.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
int shapeUninitialized(OrdinalType numRowsCols)
Set dimensions of a Teuchos::SerialSymDenseMatrix object; don&#39;t initialize values.
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.
OrdinalType numCols() const
Returns the column dimension of this matrix.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
SerialSymDenseMatrix< OrdinalType, ScalarType > & assign(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
bool operator!=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
Teuchos::DataAccess Mode enumerable type.
bool empty() const
Returns whether this matrix is empty.
ScalarType scalarType
Typedef for scalar type.