Teuchos - Trilinos Tools Package  Version of the Day
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 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ***********************************************************************
41 // @HEADER
42 
43 #ifndef _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_
44 #define _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_
45 
49 #include "Teuchos_CompObject.hpp"
50 #include "Teuchos_BLAS.hpp"
51 #include "Teuchos_ScalarTraits.hpp"
52 #include "Teuchos_DataAccess.hpp"
53 #include "Teuchos_ConfigDefs.hpp"
54 #include "Teuchos_Assert.hpp"
55 
118 namespace Teuchos {
119 
120 template<typename OrdinalType, typename ScalarType>
121 class SerialSymDenseMatrix : public CompObject, public Object, public BLAS<OrdinalType,ScalarType>
122 {
123  public:
124 
126  typedef OrdinalType ordinalType;
128  typedef ScalarType scalarType;
129 
131 
132 
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 
277  int random( const ScalarType bias = 0.1*Teuchos::ScalarTraits<ScalarType>::one() );
278 
280 
282 
283 
285 
292  ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
293 
295 
302  const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
303 
305 
307  ScalarType* values() const { return(values_); }
308 
310 
312 
313 
315  bool upper() const {return(upper_);};
316 
318  char UPLO() const {return(UPLO_);};
320 
322 
323 
325 
332 
334 
338 
340 
344 
346 
348 
349 
351 
355 
357 
361 
363 
365 
366 
368  OrdinalType numRows() const { return(numRowCols_); }
369 
371  OrdinalType numCols() const { return(numRowCols_); }
372 
374  OrdinalType stride() const { return(stride_); }
375 
377  bool empty() const { return(numRowCols_ == 0); }
378 
380 
382 
383 
386 
389 
393 
395 
396  virtual void print(std::ostream& os) const;
398 
400 
401  protected:
402 
403  // In-place scaling of matrix.
404  void scale( const ScalarType alpha );
405 
406  // Copy the values from one matrix to the other.
407  void copyMat(bool inputUpper, ScalarType* inputMatrix, OrdinalType inputStride,
408  OrdinalType numRowCols, bool outputUpper, ScalarType* outputMatrix,
409  OrdinalType outputStride, OrdinalType startRowCol,
410  ScalarType alpha = ScalarTraits<ScalarType>::zero() );
411 
412  // Copy the values from the active triangle of the matrix to the other to make the matrix full symmetric.
413  void copyUPLOMat(bool inputUpper, ScalarType* inputMatrix,
414  OrdinalType inputStride, OrdinalType inputRows);
415 
416  void deleteArrays();
417  void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
418  OrdinalType numRowCols_;
419  OrdinalType stride_;
420  bool valuesCopied_;
421  ScalarType* values_;
422  bool upper_;
423  char UPLO_;
424 
425 
426 };
427 
428 //----------------------------------------------------------------------------------------------------
429 // Constructors and Destructor
430 //----------------------------------------------------------------------------------------------------
431 
432 template<typename OrdinalType, typename ScalarType>
434  : CompObject(), Object(), BLAS<OrdinalType,ScalarType>(),
435  numRowCols_(0), stride_(0), valuesCopied_(false), values_(0), upper_(false), UPLO_('L')
436 {}
437 
438 template<typename OrdinalType, typename ScalarType>
440  : CompObject(), Object(), BLAS<OrdinalType,ScalarType>(),
441  numRowCols_(numRowCols_in), stride_(numRowCols_in), valuesCopied_(false), upper_(false), UPLO_('L')
442 {
443  values_ = new ScalarType[stride_*numRowCols_];
444  valuesCopied_ = true;
445  if (zeroOut == true)
447 }
448 
449 template<typename OrdinalType, typename ScalarType>
451  DataAccess CV, bool upper_in, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRowCols_in
452  )
453  : CompObject(), Object(), BLAS<OrdinalType,ScalarType>(),
454  numRowCols_(numRowCols_in), stride_(stride_in), valuesCopied_(false),
455  values_(values_in), upper_(upper_in)
456 {
457  if (upper_)
458  UPLO_ = 'U';
459  else
460  UPLO_ = 'L';
461 
462  if(CV == Copy)
463  {
464  stride_ = numRowCols_;
465  values_ = new ScalarType[stride_*numRowCols_];
466  copyMat(upper_in, values_in, stride_in, numRowCols_, upper_, values_, stride_, 0);
467  valuesCopied_ = true;
468  }
469 }
470 
471 template<typename OrdinalType, typename ScalarType>
473  : CompObject(), Object(), BLAS<OrdinalType,ScalarType>(),
474  numRowCols_(Source.numRowCols_), stride_(0), valuesCopied_(true),
475  values_(0), upper_(Source.upper_), UPLO_(Source.UPLO_)
476 {
477  if (!Source.valuesCopied_)
478  {
479  stride_ = Source.stride_;
480  values_ = Source.values_;
481  valuesCopied_ = false;
482  }
483  else
484  {
485  stride_ = numRowCols_;
486  const OrdinalType newsize = stride_ * numRowCols_;
487  if(newsize > 0) {
488  values_ = new ScalarType[newsize];
489  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0);
490  }
491  else {
492  numRowCols_ = 0; stride_ = 0;
493  valuesCopied_ = false;
494  }
495  }
496 }
497 
498 template<typename OrdinalType, typename ScalarType>
500  DataAccess CV, const SerialSymDenseMatrix<OrdinalType,
501  ScalarType> &Source, OrdinalType numRowCols_in, OrdinalType startRowCol )
502  : CompObject(), Object(), BLAS<OrdinalType,ScalarType>(),
503  numRowCols_(numRowCols_in), stride_(Source.stride_), valuesCopied_(false), upper_(Source.upper_), UPLO_(Source.UPLO_)
504 {
505  if(CV == Copy)
506  {
507  stride_ = numRowCols_in;
508  values_ = new ScalarType[stride_ * numRowCols_in];
509  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_in, upper_, values_, stride_, startRowCol);
510  valuesCopied_ = true;
511  }
512  else // CV == View
513  {
514  values_ = Source.values_ + (stride_ * startRowCol) + startRowCol;
515  }
516 }
517 
518 template<typename OrdinalType, typename ScalarType>
520 {
521  deleteArrays();
522 }
523 
524 //----------------------------------------------------------------------------------------------------
525 // Shape methods
526 //----------------------------------------------------------------------------------------------------
527 
528 template<typename OrdinalType, typename ScalarType>
530 {
531  deleteArrays(); // Get rid of anything that might be already allocated
532  numRowCols_ = numRowCols_in;
533  stride_ = numRowCols_;
534  values_ = new ScalarType[stride_*numRowCols_];
535  putScalar( Teuchos::ScalarTraits<ScalarType>::zero(), true );
536  valuesCopied_ = true;
537  return(0);
538 }
539 
540 template<typename OrdinalType, typename ScalarType>
542 {
543  deleteArrays(); // Get rid of anything that might be already allocated
544  numRowCols_ = numRowCols_in;
545  stride_ = numRowCols_;
546  values_ = new ScalarType[stride_*numRowCols_];
547  valuesCopied_ = true;
548  return(0);
549 }
550 
551 template<typename OrdinalType, typename ScalarType>
553 {
554  // Allocate space for new matrix
555  ScalarType* values_tmp = new ScalarType[numRowCols_in * numRowCols_in];
556  ScalarType zero = ScalarTraits<ScalarType>::zero();
557  for(OrdinalType k = 0; k < numRowCols_in * numRowCols_in; k++)
558  {
559  values_tmp[k] = zero;
560  }
561  OrdinalType numRowCols_tmp = TEUCHOS_MIN(numRowCols_, numRowCols_in);
562  if(values_ != 0)
563  {
564  copyMat(upper_, values_, stride_, numRowCols_tmp, upper_, values_tmp, numRowCols_in, 0); // Copy principal submatrix of A to new A
565  }
566  deleteArrays(); // Get rid of anything that might be already allocated
567  numRowCols_ = numRowCols_in;
568  stride_ = numRowCols_;
569  values_ = values_tmp; // Set pointer to new A
570  valuesCopied_ = true;
571  return(0);
572 }
573 
574 //----------------------------------------------------------------------------------------------------
575 // Set methods
576 //----------------------------------------------------------------------------------------------------
577 
578 template<typename OrdinalType, typename ScalarType>
580 {
581  // Do nothing if the matrix is already a lower triangular matrix
582  if (upper_ != false) {
583  copyUPLOMat( true, values_, stride_, numRowCols_ );
584  upper_ = false;
585  UPLO_ = 'L';
586  }
587 }
588 
589 template<typename OrdinalType, typename ScalarType>
591 {
592  // Do nothing if the matrix is already an upper triangular matrix
593  if (upper_ == false) {
594  copyUPLOMat( false, values_, stride_, numRowCols_ );
595  upper_ = true;
596  UPLO_ = 'U';
597  }
598 }
599 
600 template<typename OrdinalType, typename ScalarType>
601 int SerialSymDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in, bool fullMatrix )
602 {
603  // Set each value of the dense matrix to "value".
604  if (fullMatrix == true) {
605  for(OrdinalType j = 0; j < numRowCols_; j++)
606  {
607  for(OrdinalType i = 0; i < numRowCols_; i++)
608  {
609  values_[i + j*stride_] = value_in;
610  }
611  }
612  }
613  // Set the active upper or lower triangular part of the matrix to "value"
614  else {
615  if (upper_) {
616  for(OrdinalType j = 0; j < numRowCols_; j++) {
617  for(OrdinalType i = 0; i <= j; i++) {
618  values_[i + j*stride_] = value_in;
619  }
620  }
621  }
622  else {
623  for(OrdinalType j = 0; j < numRowCols_; j++) {
624  for(OrdinalType i = j; i < numRowCols_; i++) {
625  values_[i + j*stride_] = value_in;
626  }
627  }
628  }
629  }
630  return 0;
631 }
632 
633 template<typename OrdinalType, typename ScalarType>
635 {
637 
638  // Set each value of the dense matrix to a random value.
639  std::vector<MT> diagSum( numRowCols_, Teuchos::ScalarTraits<MT>::zero() );
640  if (upper_) {
641  for(OrdinalType j = 0; j < numRowCols_; j++) {
642  for(OrdinalType i = 0; i < j; 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  else {
650  for(OrdinalType j = 0; j < numRowCols_; j++) {
651  for(OrdinalType i = j+1; i < numRowCols_; i++) {
652  values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
653  diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
654  diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
655  }
656  }
657  }
658 
659  // Set the diagonal.
660  for(OrdinalType i = 0; i < numRowCols_; i++) {
661  values_[i + i*stride_] = diagSum[i] + bias;
662  }
663  return 0;
664 }
665 
666 template<typename OrdinalType, typename ScalarType>
669 {
670  if(this == &Source)
671  return(*this); // Special case of source same as target
672  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
673  upper_ = Source.upper_; // Might have to change the active part of the matrix.
674  return(*this); // Special case of both are views to same data.
675  }
676 
677  // If the source is a view then we will return a view, else we will return a copy.
678  if (!Source.valuesCopied_) {
679  if(valuesCopied_) {
680  // Clean up stored data if this was previously a copy.
681  deleteArrays();
682  }
683  numRowCols_ = Source.numRowCols_;
684  stride_ = Source.stride_;
685  values_ = Source.values_;
686  upper_ = Source.upper_;
687  UPLO_ = Source.UPLO_;
688  }
689  else {
690  // If we were a view, we will now be a copy.
691  if(!valuesCopied_) {
692  numRowCols_ = Source.numRowCols_;
693  stride_ = Source.numRowCols_;
694  upper_ = Source.upper_;
695  UPLO_ = Source.UPLO_;
696  const OrdinalType newsize = stride_ * numRowCols_;
697  if(newsize > 0) {
698  values_ = new ScalarType[newsize];
699  valuesCopied_ = true;
700  }
701  else {
702  values_ = 0;
703  }
704  }
705  // If we were a copy, we will stay a copy.
706  else {
707  if((Source.numRowCols_ <= stride_) && (Source.numRowCols_ == numRowCols_)) { // we don't need to reallocate
708  numRowCols_ = Source.numRowCols_;
709  upper_ = Source.upper_;
710  UPLO_ = Source.UPLO_;
711  }
712  else { // we need to allocate more space (or less space)
713  deleteArrays();
714  numRowCols_ = Source.numRowCols_;
715  stride_ = Source.numRowCols_;
716  upper_ = Source.upper_;
717  UPLO_ = Source.UPLO_;
718  const OrdinalType newsize = stride_ * numRowCols_;
719  if(newsize > 0) {
720  values_ = new ScalarType[newsize];
721  valuesCopied_ = true;
722  }
723  }
724  }
725  copyMat(Source.upper_, Source.values_, Source.stride_, Source.numRowCols_, upper_, values_, stride_, 0);
726  }
727  return(*this);
728 }
729 
730 template<typename OrdinalType, typename ScalarType>
732 {
733  this->scale(alpha);
734  return(*this);
735 }
736 
737 template<typename OrdinalType, typename ScalarType>
739 {
740  // Check for compatible dimensions
741  if ((numRowCols_ != Source.numRowCols_))
742  {
743  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
744  }
745  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0, ScalarTraits<ScalarType>::one());
746  return(*this);
747 }
748 
749 template<typename OrdinalType, typename ScalarType>
751 {
752  // Check for compatible dimensions
753  if ((numRowCols_ != Source.numRowCols_))
754  {
755  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
756  }
757  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0, -ScalarTraits<ScalarType>::one());
758  return(*this);
759 }
760 
761 template<typename OrdinalType, typename ScalarType>
763  if(this == &Source)
764  return(*this); // Special case of source same as target
765  if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
766  upper_ = Source.upper_; // We may have to change the active part of the matrix.
767  return(*this); // Special case of both are views to same data.
768  }
769 
770  // Check for compatible dimensions
771  if ((numRowCols_ != Source.numRowCols_))
772  {
773  TEUCHOS_CHK_REF(*this); // Return *this without altering it.
774  }
775  copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0 );
776  return(*this);
777 }
778 
779 //----------------------------------------------------------------------------------------------------
780 // Accessor methods
781 //----------------------------------------------------------------------------------------------------
782 
783 template<typename OrdinalType, typename ScalarType>
784 inline ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
785 {
786 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
787  checkIndex( rowIndex, colIndex );
788 #endif
789  if ( rowIndex <= colIndex ) {
790  // Accessing upper triangular part of matrix
791  if (upper_)
792  return(values_[colIndex * stride_ + rowIndex]);
793  else
794  return(values_[rowIndex * stride_ + colIndex]);
795  }
796  else {
797  // Accessing lower triangular part of matrix
798  if (upper_)
799  return(values_[rowIndex * stride_ + colIndex]);
800  else
801  return(values_[colIndex * stride_ + rowIndex]);
802  }
803 }
804 
805 template<typename OrdinalType, typename ScalarType>
806 inline const ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
807 {
808 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
809  checkIndex( rowIndex, colIndex );
810 #endif
811  if ( rowIndex <= colIndex ) {
812  // Accessing upper triangular part of matrix
813  if (upper_)
814  return(values_[colIndex * stride_ + rowIndex]);
815  else
816  return(values_[rowIndex * stride_ + colIndex]);
817  }
818  else {
819  // Accessing lower triangular part of matrix
820  if (upper_)
821  return(values_[rowIndex * stride_ + colIndex]);
822  else
823  return(values_[colIndex * stride_ + rowIndex]);
824  }
825 }
826 
827 //----------------------------------------------------------------------------------------------------
828 // Norm methods
829 //----------------------------------------------------------------------------------------------------
830 
831 template<typename OrdinalType, typename ScalarType>
833 {
834  return(normInf());
835 }
836 
837 template<typename OrdinalType, typename ScalarType>
839 {
840  typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
841 
842  OrdinalType i, j;
843 
844  MT sum, anorm = ScalarTraits<MT>::zero();
845  ScalarType* ptr;
846 
847  if (upper_) {
848  for (j=0; j<numRowCols_; j++) {
849  sum = ScalarTraits<MT>::zero();
850  ptr = values_ + j*stride_;
851  for (i=0; i<j; i++) {
852  sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
853  }
854  ptr = values_ + j + j*stride_;
855  for (i=j; i<numRowCols_; i++) {
857  ptr += stride_;
858  }
859  anorm = TEUCHOS_MAX( anorm, sum );
860  }
861  }
862  else {
863  for (j=0; j<numRowCols_; j++) {
864  sum = ScalarTraits<MT>::zero();
865  ptr = values_ + j + j*stride_;
866  for (i=j; i<numRowCols_; i++) {
867  sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
868  }
869  ptr = values_ + j;
870  for (i=0; i<j; i++) {
872  ptr += stride_;
873  }
874  anorm = TEUCHOS_MAX( anorm, sum );
875  }
876  }
877  return(anorm);
878 }
879 
880 template<typename OrdinalType, typename ScalarType>
882 {
883  typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
884 
885  OrdinalType i, j;
886 
887  MT sum = ScalarTraits<MT>::zero(), anorm = ScalarTraits<MT>::zero();
888 
889  if (upper_) {
890  for (j = 0; j < numRowCols_; j++) {
891  for (i = 0; i < j; i++) {
892  sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
893  }
894  sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
895  }
896  }
897  else {
898  for (j = 0; j < numRowCols_; j++) {
899  sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
900  for (i = j+1; i < numRowCols_; i++) {
901  sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
902  }
903  }
904  }
906  return(anorm);
907 }
908 
909 //----------------------------------------------------------------------------------------------------
910 // Comparison methods
911 //----------------------------------------------------------------------------------------------------
912 
913 template<typename OrdinalType, typename ScalarType>
915 {
916  bool result = 1;
917  if((numRowCols_ != Operand.numRowCols_)) {
918  result = 0;
919  }
920  else {
921  OrdinalType i, j;
922  for(i = 0; i < numRowCols_; i++) {
923  for(j = 0; j < numRowCols_; j++) {
924  if((*this)(i, j) != Operand(i, j)) {
925  return 0;
926  }
927  }
928  }
929  }
930  return result;
931 }
932 
933 template<typename OrdinalType, typename ScalarType>
935 {
936  return !((*this) == Operand);
937 }
938 
939 //----------------------------------------------------------------------------------------------------
940 // Multiplication method
941 //----------------------------------------------------------------------------------------------------
942 
943 template<typename OrdinalType, typename ScalarType>
944 void SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
945 {
946  OrdinalType i, j;
947  ScalarType* ptr;
948 
949  if (upper_) {
950  for (j=0; j<numRowCols_; j++) {
951  ptr = values_ + j*stride_;
952  for (i=0; i<= j; i++) { *ptr = alpha * (*ptr); ptr++; }
953  }
954  }
955  else {
956  for (j=0; j<numRowCols_; j++) {
957  ptr = values_ + j*stride_ + j;
958  for (i=j; i<numRowCols_; i++) { *ptr = alpha * (*ptr); ptr++; }
959  }
960  }
961 }
962 
963 /*
964 template<typename OrdinalType, typename ScalarType>
965 int SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const SerialSymDenseMatrix<OrdinalType,ScalarType>& A )
966 {
967  OrdinalType i, j;
968  ScalarType* ptr;
969 
970  // Check for compatible dimensions
971  if ((numRowCols_ != A.numRowCols_)) {
972  TEUCHOS_CHK_ERR(-1); // Return error
973  }
974 
975  if (upper_) {
976  for (j=0; j<numRowCols_; j++) {
977  ptr = values_ + j*stride_;
978  for (i=0; i<= j; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
979  }
980  }
981  else {
982  for (j=0; j<numRowCols_; j++) {
983  ptr = values_ + j*stride_;
984  for (i=j; i<numRowCols_; i++) { *ptr = A(i,j) * (*ptr); ptr++; }
985  }
986  }
987 
988  return(0);
989 }
990 */
991 
992 template<typename OrdinalType, typename ScalarType>
994 {
995  os << std::endl;
996  if(valuesCopied_)
997  os << "Values_copied : yes" << std::endl;
998  else
999  os << "Values_copied : no" << std::endl;
1000  os << "Rows / Columns : " << numRowCols_ << std::endl;
1001  os << "LDA : " << stride_ << std::endl;
1002  if (upper_)
1003  os << "Storage: Upper" << std::endl;
1004  else
1005  os << "Storage: Lower" << std::endl;
1006  if(numRowCols_ == 0) {
1007  os << "(matrix is empty, no values to display)" << std::endl;
1008  } else {
1009  for(OrdinalType i = 0; i < numRowCols_; i++) {
1010  for(OrdinalType j = 0; j < numRowCols_; j++){
1011  os << (*this)(i,j) << " ";
1012  }
1013  os << std::endl;
1014  }
1015  }
1016 }
1017 
1018 //----------------------------------------------------------------------------------------------------
1019 // Protected methods
1020 //----------------------------------------------------------------------------------------------------
1021 
1022 template<typename OrdinalType, typename ScalarType>
1023 inline void SerialSymDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
1024  TEUCHOS_TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRowCols_, std::out_of_range,
1025  "SerialSymDenseMatrix<T>::checkIndex: "
1026  "Row index " << rowIndex << " out of range [0, "<< numRowCols_ << ")");
1027  TEUCHOS_TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numRowCols_, std::out_of_range,
1028  "SerialSymDenseMatrix<T>::checkIndex: "
1029  "Col index " << colIndex << " out of range [0, "<< numRowCols_ << ")");
1030 }
1031 
1032 template<typename OrdinalType, typename ScalarType>
1033 void SerialSymDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
1034 {
1035  if (valuesCopied_)
1036  {
1037  delete [] values_;
1038  values_ = 0;
1039  valuesCopied_ = false;
1040  }
1041 }
1042 
1043 template<typename OrdinalType, typename ScalarType>
1044 void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyMat(
1045  bool inputUpper, ScalarType* inputMatrix,
1046  OrdinalType inputStride, OrdinalType numRowCols_in,
1047  bool outputUpper, ScalarType* outputMatrix,
1048  OrdinalType outputStride, OrdinalType startRowCol,
1049  ScalarType alpha
1050  )
1051 {
1052  OrdinalType i, j;
1053  ScalarType* ptr1 = 0;
1054  ScalarType* ptr2 = 0;
1055 
1056  for (j = 0; j < numRowCols_in; j++) {
1057  if (inputUpper == true) {
1058  // The input matrix is upper triangular, start at the beginning of each column.
1059  ptr2 = inputMatrix + (j + startRowCol) * inputStride + startRowCol;
1060  if (outputUpper == true) {
1061  // The output matrix matches the same pattern as the input matrix.
1062  ptr1 = outputMatrix + j*outputStride;
1063  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1064  for(i = 0; i <= j; i++) {
1065  *ptr1++ += alpha*(*ptr2++);
1066  }
1067  } else {
1068  for(i = 0; i <= j; i++) {
1069  *ptr1++ = *ptr2++;
1070  }
1071  }
1072  }
1073  else {
1074  // The output matrix has the opposite pattern as the input matrix.
1075  // Copy down across rows of the output matrix, but down columns of the input matrix.
1076  ptr1 = outputMatrix + j;
1077  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1078  for(i = 0; i <= j; i++) {
1079  *ptr1 += alpha*(*ptr2++);
1080  ptr1 += outputStride;
1081  }
1082  } else {
1083  for(i = 0; i <= j; i++) {
1084  *ptr1 = *ptr2++;
1085  ptr1 += outputStride;
1086  }
1087  }
1088  }
1089  }
1090  else {
1091  // The input matrix is lower triangular, start at the diagonal of each row.
1092  ptr2 = inputMatrix + (startRowCol+j) * inputStride + startRowCol + j;
1093  if (outputUpper == true) {
1094  // The output matrix has the opposite pattern as the input matrix.
1095  // Copy across rows of the output matrix, but down columns of the input matrix.
1096  ptr1 = outputMatrix + j*outputStride + j;
1097  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1098  for(i = j; i < numRowCols_in; i++) {
1099  *ptr1 += alpha*(*ptr2++);
1100  ptr1 += outputStride;
1101  }
1102  } else {
1103  for(i = j; i < numRowCols_in; i++) {
1104  *ptr1 = *ptr2++;
1105  ptr1 += outputStride;
1106  }
1107  }
1108  }
1109  else {
1110  // The output matrix matches the same pattern as the input matrix.
1111  ptr1 = outputMatrix + j*outputStride + j;
1112  if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
1113  for(i = j; i < numRowCols_in; i++) {
1114  *ptr1++ += alpha*(*ptr2++);
1115  }
1116  } else {
1117  for(i = j; i < numRowCols_in; i++) {
1118  *ptr1++ = *ptr2++;
1119  }
1120  }
1121  }
1122  }
1123  }
1124 }
1125 
1126 template<typename OrdinalType, typename ScalarType>
1127 void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyUPLOMat(
1128  bool inputUpper, ScalarType* inputMatrix,
1129  OrdinalType inputStride, OrdinalType inputRows
1130  )
1131 {
1132  OrdinalType i, j;
1133  ScalarType * ptr1 = 0;
1134  ScalarType * ptr2 = 0;
1135 
1136  if (inputUpper) {
1137  for (j=1; j<inputRows; j++) {
1138  ptr1 = inputMatrix + j;
1139  ptr2 = inputMatrix + j*inputStride;
1140  for (i=0; i<j; i++) {
1141  *ptr1 = *ptr2++;
1142  ptr1+=inputStride;
1143  }
1144  }
1145  }
1146  else {
1147  for (i=1; i<inputRows; i++) {
1148  ptr1 = inputMatrix + i;
1149  ptr2 = inputMatrix + i*inputStride;
1150  for (j=0; j<i; j++) {
1151  *ptr2++ = *ptr1;
1152  ptr1+=inputStride;
1153  }
1154  }
1155  }
1156 }
1157 
1158 } // namespace Teuchos
1159 
1160 #endif /* _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_ */
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
virtual void print(std::ostream &os) const
Print method. Defines the behavior of the std::ostream << operator inherited from the Object class...
OrdinalType numCols() const
Returns the column dimension of this matrix.
T magnitudeType
Mandatory typedef for result of magnitude.
SerialSymDenseMatrix()
Default constructor; defines a zero size object.
Templated interface class to BLAS routines.
bool empty() const
Returns whether this matrix is empty.
ScalarType & operator()(OrdinalType rowIndex, OrdinalType colIndex)
Element access method (non-const).
#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...
bool operator==(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Operand) const
Equality of two matrices.
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.
OrdinalType numRows() const
Returns the row dimension of this matrix.
bool upper() const
Returns true if upper triangular part of this matrix has and will be used.
SerialSymDenseMatrix< OrdinalType, ScalarType > & operator*=(const ScalarType alpha)
Inplace scalar-matrix product A = alpha*A.
The base Teuchos class.
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.
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...
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.
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos, as well as a number of utility routines.
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.
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
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.
char UPLO() const
Returns character value of UPLO used by LAPACK routines.
SerialSymDenseMatrix< OrdinalType, ScalarType > & assign(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Source)
Copies values from one matrix to another.
Teuchos::DataAccess Mode enumerable type.
ScalarType * values() const
Returns the pointer to the ScalarType data array contained in the object.
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
ScalarType scalarType
Typedef for scalar type.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
bool operator!=(const SerialSymDenseMatrix< OrdinalType, ScalarType > &Operand) const
Inequality of two matrices.