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 //
5 // Teuchos: Common Tools Package
6 // Copyright (2004) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // ***********************************************************************
39 // @HEADER
40 
41 #ifndef _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_
42 #define _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_
43 
47 #include "Teuchos_CompObject.hpp"
48 #include "Teuchos_BLAS.hpp"
49 #include "Teuchos_ScalarTraits.hpp"
50 #include "Teuchos_DataAccess.hpp"
51 #include "Teuchos_ConfigDefs.hpp"
52 #include "Teuchos_Assert.hpp"
53 #include <utility>
54 #include <vector>
55 
118 namespace Teuchos {
119 
120 template<typename OrdinalType, typename ScalarType>
121 class SerialSymDenseMatrix : public CompObject, public BLAS<OrdinalType,ScalarType>
122 {
123  public:
124 
126  typedef OrdinalType ordinalType;
128  typedef ScalarType scalarType;
129 
131 
132 
141  SerialSymDenseMatrix() = default;
142 
144 
154  SerialSymDenseMatrix(OrdinalType numRowsCols, bool zeroOut = true);
155 
157 
168  SerialSymDenseMatrix(DataAccess CV, bool upper, ScalarType* values, OrdinalType stride, OrdinalType numRowsCols);
169 
172 
174 
183  SerialSymDenseMatrix(DataAccess CV, const SerialSymDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRowCols, OrdinalType startRowCol=0);
184 
186  virtual ~SerialSymDenseMatrix ();
188 
190 
191 
193 
202  int shape(OrdinalType numRowsCols);
203 
205 
214  int shapeUninitialized(OrdinalType numRowsCols);
215 
217 
227  int reshape(OrdinalType numRowsCols);
228 
230 
232  void setLower();
233 
235 
237  void setUpper();
238 
240 
242 
243 
245 
252 
254 
259 
261 
264  SerialSymDenseMatrix<OrdinalType, ScalarType>& operator= (const ScalarType value) { putScalar(value); return(*this); }
265 
267 
272  int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero(), bool fullMatrix = false );
273 
275 
279 
281 
283  int random( const ScalarType bias = 0.1*Teuchos::ScalarTraits<ScalarType>::one() );
284 
286 
288 
289 
291 
298  ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
299 
301 
308  const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
309 
311 
313  ScalarType* values() const { return(values_); }
314 
316 
318 
319 
321  bool upper() const {return(upper_);};
322 
324  char UPLO() const {return(UPLO_);};
326 
328 
329 
331 
338 
340 
344 
346 
350 
352 
354 
355 
357 
361 
363 
367 
369 
371 
372 
374  OrdinalType numRows() const { return(numRowCols_); }
375 
377  OrdinalType numCols() const { return(numRowCols_); }
378 
380  OrdinalType stride() const { return(stride_); }
381 
383  bool empty() const { return(numRowCols_ == 0); }
384 
386 
388 
389 
392 
395 
399 
401 
402  virtual std::ostream& print(std::ostream& os) const;
404 
406 
407  protected:
408 
409  // In-place scaling of matrix.
410  void scale( const ScalarType alpha );
411 
412  // Copy the values from one matrix to the other.
413  void copyMat(bool inputUpper, ScalarType* inputMatrix, OrdinalType inputStride,
414  OrdinalType numRowCols, bool outputUpper, ScalarType* outputMatrix,
415  OrdinalType outputStride, OrdinalType startRowCol,
416  ScalarType alpha = ScalarTraits<ScalarType>::zero() );
417 
418  // Copy the values from the active triangle of the matrix to the other to make the matrix full symmetric.
419  void copyUPLOMat(bool inputUpper, ScalarType* inputMatrix,
420  OrdinalType inputStride, OrdinalType inputRows);
421 
422  void deleteArrays();
423  void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
424 
425  static ScalarType*
426  allocateValues(const OrdinalType numRows,
427  const OrdinalType numCols)
428  {
429  const size_t size = size_t(numRows) * size_t(numCols);
430 #pragma GCC diagnostic push
431 #pragma GCC diagnostic ignored "-Wvla"
432  return new ScalarType[size];
433 #pragma GCC diagnostic pop
434  }
435 
436  OrdinalType numRowCols_ = 0;
437  OrdinalType stride_ = 0;
438  bool valuesCopied_ = false;
439  ScalarType* values_ = nullptr;
440  bool upper_ = false;
441  char UPLO_ {'L'};
442 };
443 
444 //----------------------------------------------------------------------------------------------------
445 // Constructors and Destructor
446 //----------------------------------------------------------------------------------------------------
447 
448 template<typename OrdinalType, typename ScalarType>
450  : numRowCols_(numRowCols_in), stride_(numRowCols_in), valuesCopied_(false), upper_(false), UPLO_('L')
451 {
452  values_ = allocateValues(stride_, numRowCols_);
453  valuesCopied_ = true;
454  if (zeroOut) {
455  const ScalarType ZERO = Teuchos::ScalarTraits<ScalarType>::zero();
456  putScalar(ZERO, true);
457  }
458 }
459 
460 template<typename OrdinalType, typename ScalarType>
462  DataAccess CV, bool upper_in, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRowCols_in
463  )
464  : numRowCols_(numRowCols_in), stride_(stride_in), valuesCopied_(false),
465  values_(values_in), upper_(upper_in)
466 {
467  if (upper_)
468  UPLO_ = 'U';
469  else
470  UPLO_ = 'L';
471 
472  if(CV == Copy)
473  {
474  stride_ = numRowCols_;
475  values_ = allocateValues(stride_, numRowCols_);
476  copyMat(upper_in, values_in, stride_in, numRowCols_, upper_, values_, stride_, 0);
477  valuesCopied_ = true;
478  }
479 }
480 
481 template<typename OrdinalType, typename ScalarType>
483  : numRowCols_(Source.numRowCols_), stride_(0), valuesCopied_(true),
484  values_(0), upper_(Source.upper_), UPLO_(Source.UPLO_)
485 {
486  if (!Source.valuesCopied_)
487  {
488  stride_ = Source.stride_;
489  values_ = Source.values_;
490  valuesCopied_ = false;
491  }
492  else
493  {
494  stride_ = numRowCols_;
495  if(stride_ > 0 && numRowCols_ > 0) {
496  values_ = allocateValues(stride_, numRowCols_);
497  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0);
498  }
499  else {
500  numRowCols_ = 0;
501  stride_ = 0;
502  valuesCopied_ = false;
503  }
504  }
505 }
506 
507 template<typename OrdinalType, typename ScalarType>
509  DataAccess CV, const SerialSymDenseMatrix<OrdinalType,
510  ScalarType> &Source, OrdinalType numRowCols_in, OrdinalType startRowCol )
511  : CompObject(), BLAS<OrdinalType,ScalarType>(),
512  numRowCols_(numRowCols_in), stride_(Source.stride_), valuesCopied_(false), upper_(Source.upper_), UPLO_(Source.UPLO_)
513 {
514  if(CV == Copy)
515  {
516  stride_ = numRowCols_in;
517  values_ = allocateValues(stride_, numRowCols_in);
518  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_in, upper_, values_, stride_, startRowCol);
519  valuesCopied_ = true;
520  }
521  else // CV == View
522  {
523  values_ = Source.values_ + (stride_ * startRowCol) + startRowCol;
524  }
525 }
526 
527 template<typename OrdinalType, typename ScalarType>
529 {
530  deleteArrays();
531 }
532 
533 //----------------------------------------------------------------------------------------------------
534 // Shape methods
535 //----------------------------------------------------------------------------------------------------
536 
537 template<typename OrdinalType, typename ScalarType>
539 {
540  deleteArrays(); // Get rid of anything that might be already allocated
541  numRowCols_ = numRowCols_in;
542  stride_ = numRowCols_;
543  values_ = allocateValues(stride_, numRowCols_);
544  putScalar( Teuchos::ScalarTraits<ScalarType>::zero(), true );
545  valuesCopied_ = true;
546  return(0);
547 }
548 
549 template<typename OrdinalType, typename ScalarType>
551 {
552  deleteArrays(); // Get rid of anything that might be already allocated
553  numRowCols_ = numRowCols_in;
554  stride_ = numRowCols_;
555  values_ = allocateValues(stride_, numRowCols_);
556  valuesCopied_ = true;
557  return(0);
558 }
559 
560 template<typename OrdinalType, typename ScalarType>
562 {
563  // Allocate space for new matrix
564  ScalarType* values_tmp = allocateValues(numRowCols_in, numRowCols_in);
565  const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
566  for(OrdinalType k = 0; k < numRowCols_in * numRowCols_in; k++)
567  {
568  values_tmp[k] = zero;
569  }
570  OrdinalType numRowCols_tmp = TEUCHOS_MIN(numRowCols_, numRowCols_in);
571  if(values_ != 0)
572  {
573  copyMat(upper_, values_, stride_, numRowCols_tmp, upper_, values_tmp, numRowCols_in, 0); // Copy principal submatrix of A to new A
574  }
575  deleteArrays(); // Get rid of anything that might be already allocated
576  numRowCols_ = numRowCols_in;
577  stride_ = numRowCols_;
578  values_ = values_tmp; // Set pointer to new A
579  valuesCopied_ = true;
580  return(0);
581 }
582 
583 //----------------------------------------------------------------------------------------------------
584 // Set methods
585 //----------------------------------------------------------------------------------------------------
586 
587 template<typename OrdinalType, typename ScalarType>
589 {
590  // Do nothing if the matrix is already a lower triangular matrix
591  if (upper_ != false) {
592  copyUPLOMat( true, values_, stride_, numRowCols_ );
593  upper_ = false;
594  UPLO_ = 'L';
595  }
596 }
597 
598 template<typename OrdinalType, typename ScalarType>
600 {
601  // Do nothing if the matrix is already an upper triangular matrix
602  if (upper_ == false) {
603  copyUPLOMat( false, values_, stride_, numRowCols_ );
604  upper_ = true;
605  UPLO_ = 'U';
606  }
607 }
608 
609 template<typename OrdinalType, typename ScalarType>
610 int SerialSymDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in, bool fullMatrix )
611 {
612  // Set each value of the dense matrix to "value".
613  if (fullMatrix == true) {
614  for(OrdinalType j = 0; j < numRowCols_; j++)
615  {
616  for(OrdinalType i = 0; i < numRowCols_; i++)
617  {
618  values_[i + j*stride_] = value_in;
619  }
620  }
621  }
622  // Set the active upper or lower triangular part of the matrix to "value"
623  else {
624  if (upper_) {
625  for(OrdinalType j = 0; j < numRowCols_; j++) {
626  for(OrdinalType i = 0; i <= j; i++) {
627  values_[i + j*stride_] = value_in;
628  }
629  }
630  }
631  else {
632  for(OrdinalType j = 0; j < numRowCols_; j++) {
633  for(OrdinalType i = j; i < numRowCols_; i++) {
634  values_[i + j*stride_] = value_in;
635  }
636  }
637  }
638  }
639  return 0;
640 }
641 
642 template<typename OrdinalType, typename ScalarType> void
645 {
646  std::swap(values_ , B.values_);
647  std::swap(numRowCols_, B.numRowCols_);
648  std::swap(stride_, B.stride_);
649  std::swap(valuesCopied_, B.valuesCopied_);
650  std::swap(upper_, B.upper_);
651  std::swap(UPLO_, B.UPLO_);
652 }
653 
654 template<typename OrdinalType, typename ScalarType>
656 {
658 
659  // Set each value of the dense matrix to a random value.
660  std::vector<MT> diagSum( numRowCols_, Teuchos::ScalarTraits<MT>::zero() );
661  if (upper_) {
662  for(OrdinalType j = 0; j < numRowCols_; j++) {
663  for(OrdinalType i = 0; i < j; i++) {
664  values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
665  diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
666  diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
667  }
668  }
669  }
670  else {
671  for(OrdinalType j = 0; j < numRowCols_; j++) {
672  for(OrdinalType i = j+1; i < numRowCols_; i++) {
673  values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
674  diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
675  diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
676  }
677  }
678  }
679 
680  // Set the diagonal.
681  for(OrdinalType i = 0; i < numRowCols_; i++) {
682  values_[i + i*stride_] = diagSum[i] + bias;
683  }
684  return 0;
685 }
686 
687 template<typename OrdinalType, typename ScalarType>
690 {
691  if(this == &Source)
692  return(*this); // Special case of source same as target
693  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
694  upper_ = Source.upper_; // Might have to change the active part of the matrix.
695  return(*this); // Special case of both are views to same data.
696  }
697 
698  // If the source is a view then we will return a view, else we will return a copy.
699  if (!Source.valuesCopied_) {
700  if(valuesCopied_) {
701  // Clean up stored data if this was previously a copy.
702  deleteArrays();
703  }
704  numRowCols_ = Source.numRowCols_;
705  stride_ = Source.stride_;
706  values_ = Source.values_;
707  upper_ = Source.upper_;
708  UPLO_ = Source.UPLO_;
709  }
710  else {
711  // If we were a view, we will now be a copy.
712  if(!valuesCopied_) {
713  numRowCols_ = Source.numRowCols_;
714  stride_ = Source.numRowCols_;
715  upper_ = Source.upper_;
716  UPLO_ = Source.UPLO_;
717  if(stride_ > 0 && numRowCols_ > 0) {
718  values_ = allocateValues(stride_, numRowCols_);
719  valuesCopied_ = true;
720  }
721  else {
722  values_ = 0;
723  }
724  }
725  // If we were a copy, we will stay a copy.
726  else {
727  if((Source.numRowCols_ <= stride_) && (Source.numRowCols_ == numRowCols_)) { // we don't need to reallocate
728  numRowCols_ = Source.numRowCols_;
729  upper_ = Source.upper_;
730  UPLO_ = Source.UPLO_;
731  }
732  else { // we need to allocate more space (or less space)
733  deleteArrays();
734  numRowCols_ = Source.numRowCols_;
735  stride_ = Source.numRowCols_;
736  upper_ = Source.upper_;
737  UPLO_ = Source.UPLO_;
738  if(stride_ > 0 && numRowCols_ > 0) {
739  values_ = allocateValues(stride_, numRowCols_);
740  valuesCopied_ = true;
741  }
742  }
743  }
744  copyMat(Source.upper_, Source.values_, Source.stride_, Source.numRowCols_, upper_, values_, stride_, 0);
745  }
746  return(*this);
747 }
748 
749 template<typename OrdinalType, typename ScalarType>
751 {
752  this->scale(alpha);
753  return(*this);
754 }
755 
756 template<typename OrdinalType, typename ScalarType>
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, ScalarTraits<ScalarType>::one());
765  return(*this);
766 }
767 
768 template<typename OrdinalType, typename ScalarType>
770 {
771  // Check for compatible dimensions
772  if ((numRowCols_ != Source.numRowCols_))
773  {
774  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
775  }
776  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0, -ScalarTraits<ScalarType>::one());
777  return(*this);
778 }
779 
780 template<typename OrdinalType, typename ScalarType>
782  if(this == &Source)
783  return(*this); // Special case of source same as target
784  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
785  upper_ = Source.upper_; // We may have to change the active part of the matrix.
786  return(*this); // Special case of both are views to same data.
787  }
788 
789  // Check for compatible dimensions
790  if ((numRowCols_ != Source.numRowCols_))
791  {
792  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
793  }
794  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0 );
795  return(*this);
796 }
797 
798 //----------------------------------------------------------------------------------------------------
799 // Accessor methods
800 //----------------------------------------------------------------------------------------------------
801 
802 template<typename OrdinalType, typename ScalarType>
803 inline ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
804 {
805 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
806  checkIndex( rowIndex, colIndex );
807 #endif
808  if ( rowIndex <= colIndex ) {
809  // Accessing upper triangular part of matrix
810  if (upper_)
811  return(values_[colIndex * stride_ + rowIndex]);
812  else
813  return(values_[rowIndex * stride_ + colIndex]);
814  }
815  else {
816  // Accessing lower triangular part of matrix
817  if (upper_)
818  return(values_[rowIndex * stride_ + colIndex]);
819  else
820  return(values_[colIndex * stride_ + rowIndex]);
821  }
822 }
823 
824 template<typename OrdinalType, typename ScalarType>
825 inline const ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
826 {
827 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
828  checkIndex( rowIndex, colIndex );
829 #endif
830  if ( rowIndex <= colIndex ) {
831  // Accessing upper triangular part of matrix
832  if (upper_)
833  return(values_[colIndex * stride_ + rowIndex]);
834  else
835  return(values_[rowIndex * stride_ + colIndex]);
836  }
837  else {
838  // Accessing lower triangular part of matrix
839  if (upper_)
840  return(values_[rowIndex * stride_ + colIndex]);
841  else
842  return(values_[colIndex * stride_ + rowIndex]);
843  }
844 }
845 
846 //----------------------------------------------------------------------------------------------------
847 // Norm methods
848 //----------------------------------------------------------------------------------------------------
849 
850 template<typename OrdinalType, typename ScalarType>
852 {
853  return(normInf());
854 }
855 
856 template<typename OrdinalType, typename ScalarType>
858 {
859  typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
860 
861  OrdinalType i, j;
862 
863  MT sum, anorm = ScalarTraits<MT>::zero();
864  ScalarType* ptr;
865 
866  if (upper_) {
867  for (j=0; j<numRowCols_; j++) {
868  sum = ScalarTraits<MT>::zero();
869  ptr = values_ + j*stride_;
870  for (i=0; i<j; i++) {
871  sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
872  }
873  ptr = values_ + j + j*stride_;
874  for (i=j; i<numRowCols_; i++) {
876  ptr += stride_;
877  }
878  anorm = TEUCHOS_MAX( anorm, sum );
879  }
880  }
881  else {
882  for (j=0; j<numRowCols_; j++) {
883  sum = ScalarTraits<MT>::zero();
884  ptr = values_ + j + j*stride_;
885  for (i=j; i<numRowCols_; i++) {
886  sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
887  }
888  ptr = values_ + j;
889  for (i=0; i<j; i++) {
891  ptr += stride_;
892  }
893  anorm = TEUCHOS_MAX( anorm, sum );
894  }
895  }
896  return(anorm);
897 }
898 
899 template<typename OrdinalType, typename ScalarType>
901 {
902  typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
903 
904  OrdinalType i, j;
905 
906  MT sum = ScalarTraits<MT>::zero(), anorm = ScalarTraits<MT>::zero();
907 
908  if (upper_) {
909  for (j = 0; j < numRowCols_; j++) {
910  for (i = 0; i < j; i++) {
911  sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
912  }
913  sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
914  }
915  }
916  else {
917  for (j = 0; j < numRowCols_; j++) {
918  sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
919  for (i = j+1; i < numRowCols_; i++) {
920  sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
921  }
922  }
923  }
925  return(anorm);
926 }
927 
928 //----------------------------------------------------------------------------------------------------
929 // Comparison methods
930 //----------------------------------------------------------------------------------------------------
931 
932 template<typename OrdinalType, typename ScalarType>
934 {
935  bool result = 1;
936  if((numRowCols_ != Operand.numRowCols_)) {
937  result = 0;
938  }
939  else {
940  OrdinalType i, j;
941  for(i = 0; i < numRowCols_; i++) {
942  for(j = 0; j < numRowCols_; j++) {
943  if((*this)(i, j) != Operand(i, j)) {
944  return 0;
945  }
946  }
947  }
948  }
949  return result;
950 }
951 
952 template<typename OrdinalType, typename ScalarType>
954 {
955  return !((*this) == Operand);
956 }
957 
958 //----------------------------------------------------------------------------------------------------
959 // Multiplication method
960 //----------------------------------------------------------------------------------------------------
961 
962 template<typename OrdinalType, typename ScalarType>
963 void SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
964 {
965  OrdinalType i, j;
966  ScalarType* ptr;
967 
968  if (upper_) {
969  for (j=0; j<numRowCols_; j++) {
970  ptr = values_ + j*stride_;
971  for (i=0; i<= j; i++) { *ptr = alpha * (*ptr); ptr++; }
972  }
973  }
974  else {
975  for (j=0; j<numRowCols_; j++) {
976  ptr = values_ + j*stride_ + j;
977  for (i=j; i<numRowCols_; i++) { *ptr = alpha * (*ptr); ptr++; }
978  }
979  }
980 }
981 
982 /*
983 template<typename OrdinalType, typename ScalarType>
984 int SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const SerialSymDenseMatrix<OrdinalType,ScalarType>& A )
985 {
986  OrdinalType i, j;
987  ScalarType* ptr;
988 
989  // Check for compatible dimensions
990  if ((numRowCols_ != A.numRowCols_)) {
991  TEUCHOS_CHK_ERR(-1); // Return error
992  }
993 
994  if (upper_) {
995  for (j=0; j<numRowCols_; j++) {
996  ptr = values_ + j*stride_;
997  for (i=0; i<= j; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
998  }
999  }
1000  else {
1001  for (j=0; j<numRowCols_; j++) {
1002  ptr = values_ + j*stride_;
1003  for (i=j; i<numRowCols_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
1004  }
1005  }
1006 
1007  return(0);
1008 }
1009 */
1010 
1011 template<typename OrdinalType, typename ScalarType>
1012 std::ostream& SerialSymDenseMatrix<OrdinalType, ScalarType>::print(std::ostream& os) const
1013 {
1014  os << std::endl;
1015  if(valuesCopied_)
1016  os << "Values_copied : yes" << std::endl;
1017  else
1018  os << "Values_copied : no" << std::endl;
1019  os << "Rows / Columns : " << numRowCols_ << std::endl;
1020  os << "LDA : " << stride_ << std::endl;
1021  if (upper_)
1022  os << "Storage: Upper" << std::endl;
1023  else
1024  os << "Storage: Lower" << std::endl;
1025  if(numRowCols_ == 0) {
1026  os << "(matrix is empty, no values to display)" << std::endl;
1027  } else {
1028  for(OrdinalType i = 0; i < numRowCols_; i++) {
1029  for(OrdinalType j = 0; j < numRowCols_; j++){
1030  os << (*this)(i,j) << " ";
1031  }
1032  os << std::endl;
1033  }
1034  }
1035  return os;
1036 }
1037 
1038 //----------------------------------------------------------------------------------------------------
1039 // Protected methods
1040 //----------------------------------------------------------------------------------------------------
1041 
1042 template<typename OrdinalType, typename ScalarType>
1043 inline void SerialSymDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
1044  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRowCols_, std::out_of_range,
1045  "SerialSymDenseMatrix<T>::checkIndex: "
1046  "Row index " << rowIndex << " out of range [0, "<< numRowCols_ << ")");
1047  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numRowCols_, std::out_of_range,
1048  "SerialSymDenseMatrix<T>::checkIndex: "
1049  "Col index " << colIndex << " out of range [0, "<< numRowCols_ << ")");
1050 }
1051 
1052 template<typename OrdinalType, typename ScalarType>
1053 void SerialSymDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
1054 {
1055  if (valuesCopied_)
1056  {
1057  delete [] values_;
1058  values_ = 0;
1059  valuesCopied_ = false;
1060  }
1061 }
1062 
1063 template<typename OrdinalType, typename ScalarType>
1064 void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyMat(
1065  bool inputUpper, ScalarType* inputMatrix,
1066  OrdinalType inputStride, OrdinalType numRowCols_in,
1067  bool outputUpper, ScalarType* outputMatrix,
1068  OrdinalType outputStride, OrdinalType startRowCol,
1069  ScalarType alpha
1070  )
1071 {
1072  OrdinalType i, j;
1073  ScalarType* ptr1 = 0;
1074  ScalarType* ptr2 = 0;
1075 
1076  for (j = 0; j < numRowCols_in; j++) {
1077  if (inputUpper == true) {
1078  // The input matrix is upper triangular, start at the beginning of each column.
1079  ptr2 = inputMatrix + (j + startRowCol) * inputStride + startRowCol;
1080  if (outputUpper == true) {
1081  // The output matrix matches the same pattern as the input matrix.
1082  ptr1 = outputMatrix + j*outputStride;
1083  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1084  for(i = 0; i <= j; i++) {
1085  *ptr1++ += alpha*(*ptr2++);
1086  }
1087  } else {
1088  for(i = 0; i <= j; i++) {
1089  *ptr1++ = *ptr2++;
1090  }
1091  }
1092  }
1093  else {
1094  // The output matrix has the opposite pattern as the input matrix.
1095  // Copy down across rows of the output matrix, but down columns of the input matrix.
1096  ptr1 = outputMatrix + j;
1097  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1098  for(i = 0; i <= j; i++) {
1099  *ptr1 += alpha*(*ptr2++);
1100  ptr1 += outputStride;
1101  }
1102  } else {
1103  for(i = 0; i <= j; i++) {
1104  *ptr1 = *ptr2++;
1105  ptr1 += outputStride;
1106  }
1107  }
1108  }
1109  }
1110  else {
1111  // The input matrix is lower triangular, start at the diagonal of each row.
1112  ptr2 = inputMatrix + (startRowCol+j) * inputStride + startRowCol + j;
1113  if (outputUpper == true) {
1114  // The output matrix has the opposite pattern as the input matrix.
1115  // Copy across rows of the output matrix, but down columns of the input matrix.
1116  ptr1 = outputMatrix + j*outputStride + j;
1117  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1118  for(i = j; i < numRowCols_in; i++) {
1119  *ptr1 += alpha*(*ptr2++);
1120  ptr1 += outputStride;
1121  }
1122  } else {
1123  for(i = j; i < numRowCols_in; i++) {
1124  *ptr1 = *ptr2++;
1125  ptr1 += outputStride;
1126  }
1127  }
1128  }
1129  else {
1130  // The output matrix matches the same pattern as the input matrix.
1131  ptr1 = outputMatrix + j*outputStride + j;
1132  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1133  for(i = j; i < numRowCols_in; i++) {
1134  *ptr1++ += alpha*(*ptr2++);
1135  }
1136  } else {
1137  for(i = j; i < numRowCols_in; i++) {
1138  *ptr1++ = *ptr2++;
1139  }
1140  }
1141  }
1142  }
1143  }
1144 }
1145 
1146 template<typename OrdinalType, typename ScalarType>
1147 void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyUPLOMat(
1148  bool inputUpper, ScalarType* inputMatrix,
1149  OrdinalType inputStride, OrdinalType inputRows
1150  )
1151 {
1152  OrdinalType i, j;
1153  ScalarType * ptr1 = 0;
1154  ScalarType * ptr2 = 0;
1155 
1156  if (inputUpper) {
1157  for (j=1; j<inputRows; j++) {
1158  ptr1 = inputMatrix + j;
1159  ptr2 = inputMatrix + j*inputStride;
1160  for (i=0; i<j; i++) {
1161  *ptr1 = *ptr2++;
1162  ptr1+=inputStride;
1163  }
1164  }
1165  }
1166  else {
1167  for (i=1; i<inputRows; i++) {
1168  ptr1 = inputMatrix + i;
1169  ptr2 = inputMatrix + i*inputStride;
1170  for (j=0; j<i; j++) {
1171  *ptr2++ = *ptr1;
1172  ptr1+=inputStride;
1173  }
1174  }
1175  }
1176 }
1177 
1179 template<typename OrdinalType, typename ScalarType>
1181 public:
1185  : obj(obj_in) {}
1186 };
1187 
1189 template<typename OrdinalType, typename ScalarType>
1190 std::ostream&
1191 operator<<(std::ostream &out,
1193 {
1194  printer.obj.print(out);
1195  return out;
1196 }
1197 
1199 template<typename OrdinalType, typename ScalarType>
1200 SerialSymDenseMatrixPrinter<OrdinalType,ScalarType>
1202 {
1204 }
1205 
1206 
1207 } // namespace Teuchos
1208 
1209 #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.