00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #ifndef _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_
00031 #define _TEUCHOS_SERIALSYMDENSEMATRIX_HPP_
00032
00036 #include "Teuchos_CompObject.hpp"
00037 #include "Teuchos_BLAS.hpp"
00038 #include "Teuchos_ScalarTraits.hpp"
00039 #include "Teuchos_DataAccess.hpp"
00040 #include "Teuchos_ConfigDefs.hpp"
00041 #include "Teuchos_TestForException.hpp"
00042
00104 namespace Teuchos {
00105
00106 template<typename OrdinalType, typename ScalarType>
00107 class SerialSymDenseMatrix : public CompObject, public Object, public BLAS<OrdinalType,ScalarType>
00108 {
00109 public:
00110
00112
00113
00114
00122 SerialSymDenseMatrix();
00123
00125
00135 SerialSymDenseMatrix(OrdinalType numRowsCols, bool zeroOut = true);
00136
00138
00149 SerialSymDenseMatrix(DataAccess CV, bool upper, ScalarType* values, OrdinalType stride, OrdinalType numRowsCols);
00150
00152 SerialSymDenseMatrix(const SerialSymDenseMatrix<OrdinalType, ScalarType> &Source);
00153
00155
00164 SerialSymDenseMatrix(DataAccess CV, const SerialSymDenseMatrix<OrdinalType, ScalarType> &Source, OrdinalType numRowCols, OrdinalType startRowCol=0);
00165
00167 virtual ~SerialSymDenseMatrix ();
00169
00171
00172
00174
00183 int shape(OrdinalType numRowsCols);
00184
00186
00195 int shapeUninitialized(OrdinalType numRowsCols);
00196
00198
00208 int reshape(OrdinalType numRowsCols);
00209
00211
00213 void setLower();
00214
00216
00218 void setUpper();
00219
00221
00223
00224
00226
00232 SerialSymDenseMatrix<OrdinalType, ScalarType>& operator= (const SerialSymDenseMatrix<OrdinalType, ScalarType>& Source);
00233
00235
00239 SerialSymDenseMatrix<OrdinalType, ScalarType>& assign (const SerialSymDenseMatrix<OrdinalType, ScalarType>& Source);
00240
00242
00247 int putScalar( const ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero(), bool fullMatrix = false );
00248
00250
00252 int random( const ScalarType bias = 0.1*Teuchos::ScalarTraits<ScalarType>::one() );
00253
00255
00257
00258
00260
00267 ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex);
00268
00270
00277 const ScalarType& operator () (OrdinalType rowIndex, OrdinalType colIndex) const;
00278
00280
00282 ScalarType* values() const { return(values_); }
00283
00285
00287
00288
00290 bool upper() const {return(upper_);};
00291
00293 char UPLO() const {return(UPLO_);};
00295
00297
00298
00300
00306 SerialSymDenseMatrix<OrdinalType, ScalarType>& operator*= (const ScalarType alpha);
00307
00309
00312 SerialSymDenseMatrix<OrdinalType, ScalarType>& operator+= (const SerialSymDenseMatrix<OrdinalType, ScalarType>& Source);
00313
00315
00318 SerialSymDenseMatrix<OrdinalType, ScalarType>& operator-= (const SerialSymDenseMatrix<OrdinalType, ScalarType>& Source);
00319
00321
00323
00324
00326
00329 bool operator== (const SerialSymDenseMatrix<OrdinalType, ScalarType> &Operand);
00330
00332
00335 bool operator!= (const SerialSymDenseMatrix<OrdinalType, ScalarType> &Operand);
00336
00338
00340
00341
00343 OrdinalType numRows() const { return(numRowCols_); }
00344
00346 OrdinalType numCols() const { return(numRowCols_); }
00347
00349 OrdinalType stride() const { return(stride_); }
00350
00352 bool empty() const { return(numRowCols_ == 0); }
00353
00355
00357
00358
00360 typename ScalarTraits<ScalarType>::magnitudeType normOne() const;
00361
00363 typename ScalarTraits<ScalarType>::magnitudeType normInf() const;
00364
00366 typename ScalarTraits<ScalarType>::magnitudeType normFrobenius() const;
00368
00370
00371
00372 virtual void print(std::ostream& os) const;
00373
00375
00376 protected:
00377
00378
00379 void scale( const ScalarType alpha );
00380
00381
00382 void copyMat(bool inputUpper, ScalarType* inputMatrix, OrdinalType inputStride,
00383 OrdinalType numRowCols, bool outputUpper, ScalarType* outputMatrix,
00384 OrdinalType outputStride, OrdinalType startRowCol,
00385 ScalarType alpha = ScalarTraits<ScalarType>::zero() );
00386
00387
00388 void copyUPLOMat(bool inputUpper, ScalarType* inputMatrix,
00389 OrdinalType inputStride, OrdinalType inputRows);
00390
00391 void deleteArrays();
00392 void checkIndex( OrdinalType rowIndex, OrdinalType colIndex = 0 ) const;
00393 OrdinalType numRowCols_;
00394 OrdinalType stride_;
00395 bool valuesCopied_;
00396 ScalarType* values_;
00397 bool upper_;
00398 char UPLO_;
00399
00400
00401 };
00402
00403
00404
00405
00406
00407 template<typename OrdinalType, typename ScalarType>
00408 SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix()
00409 : CompObject(), Object("Teuchos::SerialSymDenseMatrix"), numRowCols_(0), stride_(0), valuesCopied_(false), values_(0), upper_(false), UPLO_('L')
00410 {}
00411
00412 template<typename OrdinalType, typename ScalarType>
00413 SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix(OrdinalType numRowCols_in, bool zeroOut)
00414 : CompObject(), Object("Teuchos::SerialSymDenseMatrix"), numRowCols_(numRowCols_in), stride_(numRowCols_in), valuesCopied_(false), upper_(false), UPLO_('L')
00415 {
00416 values_ = new ScalarType[stride_*numRowCols_];
00417 valuesCopied_ = true;
00418 if (zeroOut == true)
00419 putScalar( Teuchos::ScalarTraits<ScalarType>::zero(), true );
00420 }
00421
00422 template<typename OrdinalType, typename ScalarType>
00423 SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix(
00424 DataAccess CV, bool upper_in, ScalarType* values_in, OrdinalType stride_in, OrdinalType numRowCols_in
00425 )
00426 : CompObject(), Object("Teuchos::SerialSymDenseMatrix"), numRowCols_(numRowCols_in), stride_(stride_in), valuesCopied_(false),
00427 values_(values_in), upper_(upper_in)
00428 {
00429 if (upper_)
00430 UPLO_ = 'U';
00431 else
00432 UPLO_ = 'L';
00433
00434 if(CV == Copy)
00435 {
00436 stride_ = numRowCols_;
00437 values_ = new ScalarType[stride_*numRowCols_];
00438 copyMat(upper_in, values_in, stride_in, numRowCols_, upper_, values_, stride_, 0);
00439 valuesCopied_ = true;
00440 }
00441 }
00442
00443 template<typename OrdinalType, typename ScalarType>
00444 SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix(const SerialSymDenseMatrix<OrdinalType, ScalarType> &Source) : CompObject(), Object("Teuchos::SerialSymDenseMatrix"), numRowCols_(Source.numRowCols_), stride_(Source.numRowCols_), valuesCopied_(true), values_(0), upper_(Source.upper_), UPLO_(Source.UPLO_)
00445 {
00446 values_ = new ScalarType[stride_*numRowCols_];
00447 copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0);
00448 valuesCopied_ = true;
00449 }
00450
00451 template<typename OrdinalType, typename ScalarType>
00452 SerialSymDenseMatrix<OrdinalType, ScalarType>::SerialSymDenseMatrix(
00453 DataAccess CV, const SerialSymDenseMatrix<OrdinalType,
00454 ScalarType> &Source, OrdinalType numRowCols_in, OrdinalType startRowCol )
00455 : CompObject(), Object("Teuchos::SerialSymDenseMatrix"), numRowCols_(numRowCols_in), stride_(Source.stride_), valuesCopied_(false), upper_(Source.upper_), UPLO_(Source.UPLO_)
00456 {
00457 if(CV == Copy)
00458 {
00459 stride_ = numRowCols_in;
00460 values_ = new ScalarType[stride_ * numRowCols_in];
00461 copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_in, upper_, values_, stride_, startRowCol);
00462 valuesCopied_ = true;
00463 }
00464 else
00465 {
00466 values_ = Source.values_ + (stride_ * startRowCol) + startRowCol;
00467 }
00468 }
00469
00470 template<typename OrdinalType, typename ScalarType>
00471 SerialSymDenseMatrix<OrdinalType, ScalarType>::~SerialSymDenseMatrix()
00472 {
00473 deleteArrays();
00474 }
00475
00476
00477
00478
00479
00480 template<typename OrdinalType, typename ScalarType>
00481 int SerialSymDenseMatrix<OrdinalType, ScalarType>::shape( OrdinalType numRowCols_in )
00482 {
00483 deleteArrays();
00484 numRowCols_ = numRowCols_in;
00485 stride_ = numRowCols_;
00486 values_ = new ScalarType[stride_*numRowCols_];
00487 putScalar( Teuchos::ScalarTraits<ScalarType>::zero(), true );
00488 valuesCopied_ = true;
00489 return(0);
00490 }
00491
00492 template<typename OrdinalType, typename ScalarType>
00493 int SerialSymDenseMatrix<OrdinalType, ScalarType>::shapeUninitialized( OrdinalType numRowCols_in )
00494 {
00495 deleteArrays();
00496 numRowCols_ = numRowCols_in;
00497 stride_ = numRowCols_;
00498 values_ = new ScalarType[stride_*numRowCols_];
00499 valuesCopied_ = true;
00500 return(0);
00501 }
00502
00503 template<typename OrdinalType, typename ScalarType>
00504 int SerialSymDenseMatrix<OrdinalType, ScalarType>::reshape( OrdinalType numRowCols_in )
00505 {
00506
00507 ScalarType* values_tmp = new ScalarType[numRowCols_in * numRowCols_in];
00508 ScalarType zero = ScalarTraits<ScalarType>::zero();
00509 for(OrdinalType k = 0; k < numRowCols_in * numRowCols_in; k++)
00510 {
00511 values_tmp[k] = zero;
00512 }
00513 OrdinalType numRowCols_tmp = TEUCHOS_MIN(numRowCols_, numRowCols_in);
00514 if(values_ != 0)
00515 {
00516 copyMat(upper_, values_, stride_, numRowCols_tmp, upper_, values_tmp, numRowCols_in, 0);
00517 }
00518 deleteArrays();
00519 numRowCols_ = numRowCols_in;
00520 stride_ = numRowCols_;
00521 values_ = values_tmp;
00522 valuesCopied_ = true;
00523 return(0);
00524 }
00525
00526
00527
00528
00529
00530 template<typename OrdinalType, typename ScalarType>
00531 void SerialSymDenseMatrix<OrdinalType, ScalarType>::setLower()
00532 {
00533
00534 if (upper_ != false) {
00535 copyUPLOMat( true, values_, stride_, numRowCols_ );
00536 upper_ = false;
00537 UPLO_ = 'L';
00538 }
00539 }
00540
00541 template<typename OrdinalType, typename ScalarType>
00542 void SerialSymDenseMatrix<OrdinalType, ScalarType>::setUpper()
00543 {
00544
00545 if (upper_ == false) {
00546 copyUPLOMat( false, values_, stride_, numRowCols_ );
00547 upper_ = true;
00548 UPLO_ = 'U';
00549 }
00550 }
00551
00552 template<typename OrdinalType, typename ScalarType>
00553 int SerialSymDenseMatrix<OrdinalType, ScalarType>::putScalar( const ScalarType value_in, bool fullMatrix )
00554 {
00555
00556 if (fullMatrix == true) {
00557 for(OrdinalType j = 0; j < numRowCols_; j++)
00558 {
00559 for(OrdinalType i = 0; i < numRowCols_; i++)
00560 {
00561 values_[i + j*stride_] = value_in;
00562 }
00563 }
00564 }
00565
00566 else {
00567 if (upper_) {
00568 for(OrdinalType j = 0; j < numRowCols_; j++) {
00569 for(OrdinalType i = 0; i <= j; i++) {
00570 values_[i + j*stride_] = value_in;
00571 }
00572 }
00573 }
00574 else {
00575 for(OrdinalType j = 0; j < numRowCols_; j++) {
00576 for(OrdinalType i = j; i < numRowCols_; i++) {
00577 values_[i + j*stride_] = value_in;
00578 }
00579 }
00580 }
00581 }
00582 return 0;
00583 }
00584
00585 template<typename OrdinalType, typename ScalarType>
00586 int SerialSymDenseMatrix<OrdinalType, ScalarType>::random( const ScalarType bias )
00587 {
00588 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MT;
00589
00590
00591 std::vector<MT> diagSum( numRowCols_, Teuchos::ScalarTraits<MT>::zero() );
00592 if (upper_) {
00593 for(OrdinalType j = 0; j < numRowCols_; j++) {
00594 for(OrdinalType i = 0; i < j; i++) {
00595 values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
00596 diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
00597 diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
00598 }
00599 }
00600 }
00601 else {
00602 for(OrdinalType j = 0; j < numRowCols_; j++) {
00603 for(OrdinalType i = j+1; i < numRowCols_; i++) {
00604 values_[i + j*stride_] = ScalarTraits<ScalarType>::random();
00605 diagSum[i] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
00606 diagSum[j] += Teuchos::ScalarTraits<ScalarType>::magnitude( values_[i + j*stride_] );
00607 }
00608 }
00609 }
00610
00611
00612 for(OrdinalType i = 0; i < numRowCols_; i++) {
00613 values_[i + i*stride_] = diagSum[i] + bias;
00614 }
00615 return 0;
00616 }
00617
00618 template<typename OrdinalType, typename ScalarType>
00619 SerialSymDenseMatrix<OrdinalType,ScalarType>&
00620 SerialSymDenseMatrix<OrdinalType, ScalarType>::operator=( const SerialSymDenseMatrix<OrdinalType,ScalarType>& Source )
00621 {
00622 if(this == &Source)
00623 return(*this);
00624 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
00625 upper_ = Source.upper_;
00626 return(*this);
00627 }
00628
00629
00630 if (!Source.valuesCopied_) {
00631 if(valuesCopied_) {
00632
00633 deleteArrays();
00634 }
00635 numRowCols_ = Source.numRowCols_;
00636 stride_ = Source.stride_;
00637 values_ = Source.values_;
00638 upper_ = Source.upper_;
00639 UPLO_ = Source.UPLO_;
00640 }
00641 else {
00642
00643 if(!valuesCopied_) {
00644 numRowCols_ = Source.numRowCols_;
00645 stride_ = Source.numRowCols_;
00646 upper_ = Source.upper_;
00647 UPLO_ = Source.UPLO_;
00648 const OrdinalType newsize = stride_ * numRowCols_;
00649 if(newsize > 0) {
00650 values_ = new ScalarType[newsize];
00651 valuesCopied_ = true;
00652 }
00653 else {
00654 values_ = 0;
00655 }
00656 }
00657
00658 else {
00659 if((Source.numRowCols_ <= stride_) && (Source.numRowCols_ == numRowCols_)) {
00660 numRowCols_ = Source.numRowCols_;
00661 upper_ = Source.upper_;
00662 UPLO_ = Source.UPLO_;
00663 }
00664 else {
00665 deleteArrays();
00666 numRowCols_ = Source.numRowCols_;
00667 stride_ = Source.numRowCols_;
00668 upper_ = Source.upper_;
00669 UPLO_ = Source.UPLO_;
00670 const OrdinalType newsize = stride_ * numRowCols_;
00671 if(newsize > 0) {
00672 values_ = new ScalarType[newsize];
00673 valuesCopied_ = true;
00674 }
00675 }
00676 }
00677 copyMat(Source.upper_, Source.values_, Source.stride_, Source.numRowCols_, upper_, values_, stride_, 0);
00678 }
00679 return(*this);
00680 }
00681
00682 template<typename OrdinalType, typename ScalarType>
00683 SerialSymDenseMatrix<OrdinalType, ScalarType>& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator*= (const ScalarType alpha)
00684 {
00685 this->scale(alpha);
00686 return(*this);
00687 }
00688
00689 template<typename OrdinalType, typename ScalarType>
00690 SerialSymDenseMatrix<OrdinalType, ScalarType>& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator+= (const SerialSymDenseMatrix<OrdinalType,ScalarType>& Source )
00691 {
00692
00693 if ((numRowCols_ != Source.numRowCols_))
00694 {
00695 TEUCHOS_CHK_REF(*this);
00696 }
00697 copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0, 1.0);
00698 return(*this);
00699 }
00700
00701 template<typename OrdinalType, typename ScalarType>
00702 SerialSymDenseMatrix<OrdinalType, ScalarType>& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator-= (const SerialSymDenseMatrix<OrdinalType,ScalarType>& Source )
00703 {
00704
00705 if ((numRowCols_ != Source.numRowCols_))
00706 {
00707 TEUCHOS_CHK_REF(*this);
00708 }
00709 copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0, -1.0);
00710 return(*this);
00711 }
00712
00713 template<typename OrdinalType, typename ScalarType>
00714 SerialSymDenseMatrix<OrdinalType,ScalarType>& SerialSymDenseMatrix<OrdinalType, ScalarType>::assign (const SerialSymDenseMatrix<OrdinalType,ScalarType>& Source) {
00715 if(this == &Source)
00716 return(*this);
00717 if((!valuesCopied_) && (!Source.valuesCopied_) && (values_ == Source.values_)) {
00718 upper_ = Source.upper_;
00719 return(*this);
00720 }
00721
00722
00723 if ((numRowCols_ != Source.numRowCols_))
00724 {
00725 TEUCHOS_CHK_REF(*this);
00726 }
00727 copyMat(Source.upper_, Source.values_, Source.stride_, numRowCols_, upper_, values_, stride_, 0 );
00728 return(*this);
00729 }
00730
00731
00732
00733
00734
00735 template<typename OrdinalType, typename ScalarType>
00736 inline ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex)
00737 {
00738 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00739 checkIndex( rowIndex, colIndex );
00740 #endif
00741 if ( rowIndex <= colIndex ) {
00742
00743 if (upper_)
00744 return(values_[colIndex * stride_ + rowIndex]);
00745 else
00746 return(values_[rowIndex * stride_ + colIndex]);
00747 }
00748 else {
00749
00750 if (upper_)
00751 return(values_[rowIndex * stride_ + colIndex]);
00752 else
00753 return(values_[colIndex * stride_ + rowIndex]);
00754 }
00755 }
00756
00757 template<typename OrdinalType, typename ScalarType>
00758 inline const ScalarType& SerialSymDenseMatrix<OrdinalType, ScalarType>::operator () (OrdinalType rowIndex, OrdinalType colIndex) const
00759 {
00760 #ifdef HAVE_TEUCHOS_ARRAY_BOUNDSCHECK
00761 checkIndex( rowIndex, colIndex );
00762 #endif
00763 if ( rowIndex <= colIndex ) {
00764
00765 if (upper_)
00766 return(values_[colIndex * stride_ + rowIndex]);
00767 else
00768 return(values_[rowIndex * stride_ + colIndex]);
00769 }
00770 else {
00771
00772 if (upper_)
00773 return(values_[rowIndex * stride_ + colIndex]);
00774 else
00775 return(values_[colIndex * stride_ + rowIndex]);
00776 }
00777 }
00778
00779
00780
00781
00782
00783 template<typename OrdinalType, typename ScalarType>
00784 typename ScalarTraits<ScalarType>::magnitudeType SerialSymDenseMatrix<OrdinalType, ScalarType>::normOne() const
00785 {
00786 return(normInf());
00787 }
00788
00789 template<typename OrdinalType, typename ScalarType>
00790 typename ScalarTraits<ScalarType>::magnitudeType SerialSymDenseMatrix<OrdinalType, ScalarType>::normInf() const
00791 {
00792 typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
00793
00794 OrdinalType i, j;
00795
00796 MT sum, anorm = ScalarTraits<MT>::zero();
00797 ScalarType* ptr;
00798
00799 if (upper_) {
00800 for (j=0; j<numRowCols_; j++) {
00801 sum = ScalarTraits<MT>::zero();
00802 ptr = values_ + j*stride_;
00803 for (i=0; i<j; i++) {
00804 sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
00805 }
00806 ptr = values_ + j + j*stride_;
00807 for (i=j; i<numRowCols_; i++) {
00808 sum += ScalarTraits<ScalarType>::magnitude( *ptr );
00809 ptr += stride_;
00810 }
00811 anorm = TEUCHOS_MAX( anorm, sum );
00812 }
00813 }
00814 else {
00815 for (j=0; j<numRowCols_; j++) {
00816 sum = ScalarTraits<MT>::zero();
00817 ptr = values_ + j + j*stride_;
00818 for (i=j; i<numRowCols_; i++) {
00819 sum += ScalarTraits<ScalarType>::magnitude( *ptr++ );
00820 }
00821 ptr = values_ + j;
00822 for (i=0; i<j; i++) {
00823 sum += ScalarTraits<ScalarType>::magnitude( *ptr );
00824 ptr += stride_;
00825 }
00826 anorm = TEUCHOS_MAX( anorm, sum );
00827 }
00828 }
00829 return(anorm);
00830 }
00831
00832 template<typename OrdinalType, typename ScalarType>
00833 typename ScalarTraits<ScalarType>::magnitudeType SerialSymDenseMatrix<OrdinalType, ScalarType>::normFrobenius() const
00834 {
00835 typedef typename ScalarTraits<ScalarType>::magnitudeType MT;
00836
00837 OrdinalType i, j;
00838
00839 MT sum = ScalarTraits<MT>::zero(), anorm = ScalarTraits<MT>::zero();
00840
00841 if (upper_) {
00842 for (j = 0; j < numRowCols_; j++) {
00843 for (i = 0; i < j; i++) {
00844 sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
00845 }
00846 sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
00847 }
00848 }
00849 else {
00850 for (j = 0; j < numRowCols_; j++) {
00851 sum += ScalarTraits<ScalarType>::magnitude(values_[j + j*stride_]*values_[j + j*stride_]);
00852 for (i = j+1; i < numRowCols_; i++) {
00853 sum += ScalarTraits<ScalarType>::magnitude(2.0*values_[i+j*stride_]*values_[i+j*stride_]);
00854 }
00855 }
00856 }
00857 anorm = ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::squareroot(sum));
00858 return(anorm);
00859 }
00860
00861
00862
00863
00864
00865 template<typename OrdinalType, typename ScalarType>
00866 bool SerialSymDenseMatrix<OrdinalType, ScalarType>::operator== (const SerialSymDenseMatrix<OrdinalType, ScalarType> &Operand)
00867 {
00868 bool result = 1;
00869 if((numRowCols_ != Operand.numRowCols_)) {
00870 result = 0;
00871 }
00872 else {
00873 OrdinalType i, j;
00874 for(i = 0; i < numRowCols_; i++) {
00875 for(j = 0; j < numRowCols_; j++) {
00876 if((*this)(i, j) != Operand(i, j)) {
00877 return 0;
00878 }
00879 }
00880 }
00881 }
00882 return result;
00883 }
00884
00885 template<typename OrdinalType, typename ScalarType>
00886 bool SerialSymDenseMatrix<OrdinalType, ScalarType>::operator!= (const SerialSymDenseMatrix<OrdinalType, ScalarType> &Operand)
00887 {
00888 return !((*this) == Operand);
00889 }
00890
00891
00892
00893
00894
00895 template<typename OrdinalType, typename ScalarType>
00896 void SerialSymDenseMatrix<OrdinalType, ScalarType>::scale( const ScalarType alpha )
00897 {
00898 OrdinalType i, j;
00899 ScalarType* ptr;
00900
00901 if (upper_) {
00902 for (j=0; j<numRowCols_; j++) {
00903 ptr = values_ + j*stride_;
00904 for (i=0; i<= j; i++) { *ptr = alpha * (*ptr); ptr++; }
00905 }
00906 }
00907 else {
00908 for (j=0; j<numRowCols_; j++) {
00909 ptr = values_ + j*stride_;
00910 for (i=j; i<numRowCols_; i++) { *ptr = alpha * (*ptr); ptr++; }
00911 }
00912 }
00913 }
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944 template<typename OrdinalType, typename ScalarType>
00945 void SerialSymDenseMatrix<OrdinalType, ScalarType>::print(std::ostream& os) const
00946 {
00947 os << std::endl;
00948 if(valuesCopied_)
00949 os << "Values_copied : yes" << std::endl;
00950 else
00951 os << "Values_copied : no" << std::endl;
00952 os << "Rows / Columns : " << numRowCols_ << std::endl;
00953 os << "LDA : " << stride_ << std::endl;
00954 if (upper_)
00955 os << "Storage: Upper" << std::endl;
00956 else
00957 os << "Storage: Lower" << std::endl;
00958 if(numRowCols_ == 0) {
00959 os << "(matrix is empty, no values to display)" << std::endl;
00960 } else {
00961 for(OrdinalType i = 0; i < numRowCols_; i++) {
00962 for(OrdinalType j = 0; j < numRowCols_; j++){
00963 os << (*this)(i,j) << " ";
00964 }
00965 os << std::endl;
00966 }
00967 }
00968 }
00969
00970
00971
00972
00973
00974 template<typename OrdinalType, typename ScalarType>
00975 inline void SerialSymDenseMatrix<OrdinalType, ScalarType>::checkIndex( OrdinalType rowIndex, OrdinalType colIndex ) const {
00976 TEST_FOR_EXCEPTION(rowIndex < 0 || rowIndex >= numRowCols_, std::out_of_range,
00977 "SerialSymDenseMatrix<T>::checkIndex: "
00978 "Row index " << rowIndex << " out of range [0, "<< numRowCols_ << ")");
00979 TEST_FOR_EXCEPTION(colIndex < 0 || colIndex >= numRowCols_, std::out_of_range,
00980 "SerialSymDenseMatrix<T>::checkIndex: "
00981 "Col index " << colIndex << " out of range [0, "<< numRowCols_ << ")");
00982 }
00983
00984 template<typename OrdinalType, typename ScalarType>
00985 void SerialSymDenseMatrix<OrdinalType, ScalarType>::deleteArrays(void)
00986 {
00987 if (valuesCopied_)
00988 {
00989 delete [] values_;
00990 values_ = 0;
00991 valuesCopied_ = false;
00992 }
00993 }
00994
00995 template<typename OrdinalType, typename ScalarType>
00996 void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyMat(
00997 bool inputUpper, ScalarType* inputMatrix,
00998 OrdinalType inputStride, OrdinalType numRowCols_in,
00999 bool outputUpper, ScalarType* outputMatrix,
01000 OrdinalType outputStride, OrdinalType startRowCol,
01001 ScalarType alpha
01002 )
01003 {
01004 OrdinalType i, j;
01005 ScalarType* ptr1 = 0;
01006 ScalarType* ptr2 = 0;
01007
01008 for (j = 0; j < numRowCols_in; j++) {
01009 if (inputUpper == true) {
01010
01011 ptr2 = inputMatrix + (j + startRowCol) * inputStride + startRowCol;
01012 if (outputUpper == true) {
01013
01014 ptr1 = outputMatrix + j*outputStride;
01015 if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
01016 for(i = 0; i <= j; i++) {
01017 *ptr1++ += alpha*(*ptr2++);
01018 }
01019 } else {
01020 for(i = 0; i <= j; i++) {
01021 *ptr1++ = *ptr2++;
01022 }
01023 }
01024 }
01025 else {
01026
01027
01028 ptr1 = outputMatrix + j;
01029 if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
01030 for(i = 0; i <= j; i++) {
01031 *ptr1 += alpha*(*ptr2++);
01032 ptr1 += outputStride;
01033 }
01034 } else {
01035 for(i = 0; i <= j; i++) {
01036 *ptr1 = *ptr2++;
01037 ptr1 += outputStride;
01038 }
01039 }
01040 }
01041 }
01042 else {
01043
01044 ptr2 = inputMatrix + (startRowCol+j) * inputStride + startRowCol + j;
01045 if (outputUpper == true) {
01046
01047
01048 ptr1 = outputMatrix + j*outputStride + j;
01049 if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
01050 for(i = j; i < numRowCols_in; i++) {
01051 *ptr1 += alpha*(*ptr2++);
01052 ptr1 += outputStride;
01053 }
01054 } else {
01055 for(i = j; i < numRowCols_in; i++) {
01056 *ptr1 = *ptr2++;
01057 ptr1 += outputStride;
01058 }
01059 }
01060 }
01061 else {
01062
01063 ptr1 = outputMatrix + j*outputStride + j;
01064 if (alpha != Teuchos::ScalarTraits<ScalarType>::zero() ) {
01065 for(i = j; i < numRowCols_in; i++) {
01066 *ptr1++ += alpha*(*ptr2++);
01067 }
01068 } else {
01069 for(i = j; i < numRowCols_in; i++) {
01070 *ptr1++ = *ptr2++;
01071 }
01072 }
01073 }
01074 }
01075 }
01076 }
01077
01078 template<typename OrdinalType, typename ScalarType>
01079 void SerialSymDenseMatrix<OrdinalType, ScalarType>::copyUPLOMat(
01080 bool inputUpper, ScalarType* inputMatrix,
01081 OrdinalType inputStride, OrdinalType inputRows
01082 )
01083 {
01084 OrdinalType i, j;
01085 ScalarType * ptr1 = 0;
01086 ScalarType * ptr2 = 0;
01087
01088 if (inputUpper) {
01089 for (j=1; j<inputRows; j++) {
01090 ptr1 = inputMatrix + j;
01091 ptr2 = inputMatrix + j*inputStride;
01092 for (i=0; i<j; i++) {
01093 *ptr1 = *ptr2++;
01094 ptr1+=inputStride;
01095 }
01096 }
01097 }
01098 else {
01099 for (i=1; i<inputRows; i++) {
01100 ptr1 = inputMatrix + i;
01101 ptr2 = inputMatrix + i*inputStride;
01102 for (j=0; j<i; j++) {
01103 *ptr2++ = *ptr1;
01104 ptr1+=inputStride;
01105 }
01106 }
01107 }
01108 }
01109
01110 }
01111
01112 #endif