Teuchos - Trilinos Tools Package Version of the Day
Loading...
Searching...
No Matches
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_
48#include "Teuchos_BLAS.hpp"
52#include "Teuchos_Assert.hpp"
53#include <utility>
54#include <vector>
55
118namespace Teuchos {
119
120template<typename OrdinalType, typename ScalarType>
121class SerialSymDenseMatrix : public CompObject, public BLAS<OrdinalType,ScalarType>
122{
123 public:
124
129
131
132
142
144
154 SerialSymDenseMatrix(OrdinalType numRowsCols, bool zeroOut = true);
155
157
169
172
174
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
265
267
272 int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero(), bool fullMatrix = false );
273
275
279
281
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
403 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.
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,
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
448template<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) {
456 putScalar(ZERO, true);
457 }
458}
459
460template<typename OrdinalType, typename ScalarType>
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
481template<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
507template<typename OrdinalType, typename 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
527template<typename OrdinalType, typename ScalarType>
532
533//----------------------------------------------------------------------------------------------------
534// Shape methods
535//----------------------------------------------------------------------------------------------------
536
537template<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
549template<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
560template<typename OrdinalType, typename ScalarType>
562{
563 // Allocate space for new matrix
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
587template<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
598template<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
609template<typename OrdinalType, typename ScalarType>
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
642template<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
654template<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
687template<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
749template<typename OrdinalType, typename ScalarType>
755
756template<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
768template<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
780template<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
802template<typename OrdinalType, typename ScalarType>
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
824template<typename OrdinalType, typename ScalarType>
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
850template<typename OrdinalType, typename ScalarType>
855
856template<typename OrdinalType, typename ScalarType>
858{
860
861 OrdinalType i, j;
862
865
866 if (upper_) {
867 for (j=0; j<numRowCols_; j++) {
869 ptr = values_ + j*stride_;
870 for (i=0; i<j; i++) {
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++) {
884 ptr = values_ + j + j*stride_;
885 for (i=j; i<numRowCols_; i++) {
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
899template<typename OrdinalType, typename ScalarType>
901{
903
904 OrdinalType i, j;
905
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
932template<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
952template<typename OrdinalType, typename ScalarType>
957
958//----------------------------------------------------------------------------------------------------
959// Multiplication method
960//----------------------------------------------------------------------------------------------------
961
962template<typename OrdinalType, typename ScalarType>
964{
965 OrdinalType i, j;
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/*
983template<typename OrdinalType, typename ScalarType>
984int 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
1011template<typename OrdinalType, typename ScalarType>
1012std::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
1042template<typename OrdinalType, typename ScalarType>
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
1052template<typename OrdinalType, typename ScalarType>
1053void SerialSymDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
1054{
1055 if (valuesCopied_)
1056 {
1057 delete [] values_;
1058 values_ = 0;
1059 valuesCopied_ = false;
1060 }
1061}
1062
1063template<typename OrdinalType, typename ScalarType>
1064void 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.
1080 if (outputUpper == true) {
1081 // The output matrix matches the same pattern as the input matrix.
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;
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.
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.
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.
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
1146template<typename OrdinalType, typename ScalarType>
1147void 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;
1160 for (i=0; i<j; i++) {
1161 *ptr1 = *ptr2++;
1163 }
1164 }
1165 }
1166 else {
1167 for (i=1; i<inputRows; i++) {
1168 ptr1 = inputMatrix + i;
1170 for (j=0; j<i; j++) {
1171 *ptr2++ = *ptr1;
1173 }
1174 }
1175 }
1176}
1177
1179template<typename OrdinalType, typename ScalarType>
1187
1189template<typename OrdinalType, typename ScalarType>
1190std::ostream&
1191operator<<(std::ostream &out,
1193{
1194 printer.obj.print(out);
1195 return out;
1196}
1197
1199template<typename OrdinalType, typename ScalarType>
1200SerialSymDenseMatrixPrinter<OrdinalType,ScalarType>
1205
1206
1207} // namespace Teuchos
1208
1209#endif /* _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_ */
Templated interface class to BLAS routines.
Object for storing data and providing functionality that is common to all computational classes.
Teuchos header file which uses auto-configuration information to include necessary C++ headers.
Teuchos::DataAccess Mode enumerable type.
Defines basic traits for the scalar field type.
Templated BLAS wrapper.
Functionality and data that is common to all computational classes.
Smart reference counting pointer class for automatic garbage collection.
RCP(ENull null_arg=null)
Initialize RCP<T> to NULL.
Ptr< T > ptr() const
Get a safer wrapper raw C++ pointer to the underlying object.
This class creates and provides basic support for symmetric, positive-definite dense matrices of temp...
SerialSymDenseMatrix()=default
Default constructor; defines a zero size object.
bool operator==(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
int shapeUninitialized(OrdinalType numRowsCols)
Set dimensions of a Teuchos::SerialSymDenseMatrix object; don't initialize values.
int shape(OrdinalType numRowsCols)
Set dimensions of a Teuchos::SerialSymDenseMatrix object; init values to zero.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator-=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Subtract another matrix from this matrix.
virtual std::ostream & print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream << operator.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator+=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Add another matrix to this matrix.
bool operator!=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.
void swap(SerialSymDenseMatrix< OrdinalType, ScalarType > &B)
Swap values between this matrix and incoming matrix.
OrdinalType numCols() const
Returns the column dimension of this matrix.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero(), bool fullMatrix=false)
Set all values in the matrix to a constant value.
virtual ~SerialSymDenseMatrix()
Teuchos::SerialSymDenseMatrix destructor.
char UPLO() const
Returns character value of UPLO used by LAPACK routines.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator*=(const ScalarType alpha)
Inplace scalar-matrix product A = alpha*A.
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.
ScalarType scalarType
Typedef for scalar type.
bool upper() const
Returns true if upper triangular part of this matrix has and will be used.
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
ScalarType * values() const
Returns the pointer to the ScalarType data array contained in the object.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
bool empty() const
Returns whether this matrix is empty.
SerialSymDenseMatrix< OrdinalType, ScalarType > & assign(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
void setLower()
Specify that the lower triangle of the this matrix should be used.
int reshape(OrdinalType numRowsCols)
Reshape a Teuchos::SerialSymDenseMatrix object.
OrdinalType numRows() const
Returns the row dimension of this matrix.
void setUpper()
Specify that the upper triangle of the this matrix should be used.
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
OrdinalType ordinalType
Typedef for ordinal type.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos,...
SerialBandDenseMatrixPrinter< OrdinalType, ScalarType > printMat(const SerialBandDenseMatrix< OrdinalType, ScalarType > &obj)
Return SerialBandDenseMatrix ostream manipulator Use as:
This structure defines some basic traits for a scalar field type.
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
T magnitudeType
Mandatory typedef for result of magnitude.
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.
Ostream manipulator for SerialSymDenseMatrix.