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
00036 #ifndef MAT_MatrixSymmetric
00037 #define MAT_MatrixSymmetric
00038
00039 #include <algorithm>
00040
00041 #include "MatrixBase.h"
00042 #include "Interval.h"
00043 #include "LanczosLargestMagnitudeEig.h"
00044 #include "mat_utils.h"
00045 #include "truncation.h"
00046
00047 namespace mat {
00065 template<typename Treal, typename Tmatrix>
00066 class MatrixSymmetric : public MatrixBase<Treal, Tmatrix> {
00067 public:
00068 typedef VectorGeneral<Treal, typename Tmatrix::VectorType> VectorType;
00069 typedef Treal real;
00070
00071 MatrixSymmetric()
00072 :MatrixBase<Treal, Tmatrix>() {}
00073 explicit MatrixSymmetric(const MatrixSymmetric<Treal, Tmatrix>& symm)
00074 :MatrixBase<Treal, Tmatrix>(symm) {}
00075 explicit MatrixSymmetric(const XY<Treal, MatrixSymmetric<Treal, Tmatrix> >& sm)
00076 :MatrixBase<Treal, Tmatrix>() { *this = sm.A * sm.B; }
00077 explicit MatrixSymmetric(const MatrixGeneral<Treal, Tmatrix>& matr)
00078 :MatrixBase<Treal, Tmatrix>(matr) {
00079 this->matrixPtr->nosymToSym();
00080 }
00083 #if 0
00084 template<typename Tfull>
00085 inline void assign_from_full
00086 (Tfull const* const fullmatrix,
00087 int const totnrows, int const totncols) {
00088 assert(totnrows == totncols);
00089 this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
00090 this->matrixPtr->nosym_to_sym();
00091 }
00092 inline void assign_from_full
00093 (Treal const* const fullmatrix,
00094 int const totnrows, int const totncols) {
00095 assert(totnrows == totncols);
00096 this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
00097 this->matrixPtr->nosym_to_sym();
00098 }
00099 #endif
00100
00101 inline void assignFromFull
00102 (std::vector<Treal> const & fullMat) {
00103 assert((int)fullMat.size() == this->get_nrows() * this->get_ncols());
00104 this->matrixPtr->assignFromFull(fullMat);
00105 this->matrixPtr->nosymToSym();
00106 }
00107
00108 inline void assignFromFull
00109 (std::vector<Treal> const & fullMat,
00110 std::vector<int> const & rowPermutation,
00111 std::vector<int> const & colPermutation) {
00112 assert((int)fullMat.size() == this->get_nrows() * this->get_ncols());
00113 std::vector<int> rowind(fullMat.size());
00114 std::vector<int> colind(fullMat.size());
00115 int ind = 0;
00116 for (int col = 0; col < this->get_ncols(); ++col)
00117 for (int row = 0; row < this->get_nrows(); ++row) {
00118 rowind[ind] = row;
00119 colind[ind] = col;
00120 ++ind;
00121 }
00122 this->assign_from_sparse(rowind,
00123 colind,
00124 fullMat,
00125 rowPermutation,
00126 colPermutation);
00127 this->matrixPtr->nosymToSym();
00128 }
00129
00130 inline void fullMatrix(std::vector<Treal> & fullMat) const {
00131 this->matrixPtr->syFullMatrix(fullMat);
00132 }
00139 inline void fullMatrix
00140 (std::vector<Treal> & fullMat,
00141 std::vector<int> const & rowInversePermutation,
00142 std::vector<int> const & colInversePermutation) const {
00143 std::vector<int> rowind;
00144 std::vector<int> colind;
00145 std::vector<Treal> values;
00146 get_all_values(rowind, colind, values,
00147 rowInversePermutation,
00148 colInversePermutation);
00149 fullMat.assign(this->get_nrows()*this->get_ncols(),0);
00150 assert(rowind.size() == values.size());
00151 assert(rowind.size() == colind.size());
00152 for (unsigned int ind = 0; ind < values.size(); ++ind) {
00153 assert(rowind[ind] + this->get_nrows() * colind[ind] <
00154 this->get_nrows()*this->get_ncols());
00155 fullMat[rowind[ind] + this->get_nrows() * colind[ind]] =
00156 values[ind];
00157 fullMat[colind[ind] + this->get_nrows() * rowind[ind]] =
00158 values[ind];
00159 }
00160 }
00167 inline void assign_from_sparse
00168 (std::vector<int> const & rowind,
00169 std::vector<int> const & colind,
00170 std::vector<Treal> const & values) {
00171 this->matrixPtr->syAssignFromSparse(rowind, colind, values);
00172 }
00184 inline void assign_from_sparse
00185 (std::vector<int> const & rowind,
00186 std::vector<int> const & colind,
00187 std::vector<Treal> const & values,
00188 SizesAndBlocks const & newRows,
00189 SizesAndBlocks const & newCols) {
00190 this->resetSizesAndBlocks(newRows, newCols);
00191 this->matrixPtr->syAssignFromSparse(rowind, colind, values);
00192 }
00202 inline void assign_from_sparse
00203 (std::vector<int> const & rowind,
00204 std::vector<int> const & colind,
00205 std::vector<Treal> const & values,
00206 std::vector<int> const & rowPermutation,
00207 std::vector<int> const & colPermutation) {
00208 std::vector<int> newRowind;
00209 std::vector<int> newColind;
00210 this->getPermutedAndSymmetrized(rowind, rowPermutation, newRowind,
00211 colind, colPermutation, newColind);
00212
00213 this->matrixPtr->syAssignFromSparse(newRowind, newColind, values);
00214 }
00220 inline void assign_from_sparse
00221 (std::vector<int> const & rowind,
00222 std::vector<int> const & colind,
00223 std::vector<Treal> const & values,
00224 SizesAndBlocks const & newRows,
00225 SizesAndBlocks const & newCols,
00226 std::vector<int> const & rowPermutation,
00227 std::vector<int> const & colPermutation) {
00228 this->resetSizesAndBlocks(newRows, newCols);
00229 assign_from_sparse(rowind, colind, values,
00230 rowPermutation, colPermutation);
00231 }
00239 inline void add_values
00240 (std::vector<int> const & rowind,
00241 std::vector<int> const & colind,
00242 std::vector<Treal> const & values) {
00243 this->matrixPtr->syAddValues(rowind, colind, values);
00244 }
00245
00249 inline void add_values
00250 (std::vector<int> const & rowind,
00251 std::vector<int> const & colind,
00252 std::vector<Treal> const & values,
00253 std::vector<int> const & rowPermutation,
00254 std::vector<int> const & colPermutation) {
00255 std::vector<int> newRowind;
00256 std::vector<int> newColind;
00257 this->getPermutedAndSymmetrized(rowind, rowPermutation, newRowind,
00258 colind, colPermutation, newColind);
00259 this->matrixPtr->syAddValues(newRowind, newColind, values);
00260 }
00261
00262
00263
00264 inline void get_values
00265 (std::vector<int> const & rowind,
00266 std::vector<int> const & colind,
00267 std::vector<Treal> & values) const {
00268 this->matrixPtr->syGetValues(rowind, colind, values);
00269 }
00277 inline void get_values
00278 (std::vector<int> const & rowind,
00279 std::vector<int> const & colind,
00280 std::vector<Treal> & values,
00281 std::vector<int> const & rowPermutation,
00282 std::vector<int> const & colPermutation) const {
00283 std::vector<int> newRowind;
00284 std::vector<int> newColind;
00285 this->getPermutedAndSymmetrized(rowind, rowPermutation, newRowind,
00286 colind, colPermutation, newColind);
00287 this->matrixPtr->syGetValues(newRowind, newColind, values);
00288 }
00294 inline void get_all_values
00295 (std::vector<int> & rowind,
00296 std::vector<int> & colind,
00297 std::vector<Treal> & values) const {
00298 rowind.resize(0);
00299 colind.resize(0);
00300 values.resize(0);
00301 this->matrixPtr->syGetAllValues(rowind, colind, values);
00302 }
00308 inline void get_all_values
00309 (std::vector<int> & rowind,
00310 std::vector<int> & colind,
00311 std::vector<Treal> & values,
00312 std::vector<int> const & rowInversePermutation,
00313 std::vector<int> const & colInversePermutation) const {
00314 std::vector<int> tmpRowind;
00315 std::vector<int> tmpColind;
00316 tmpRowind.reserve(rowind.capacity());
00317 tmpColind.reserve(colind.capacity());
00318 values.resize(0);
00319 this->matrixPtr->syGetAllValues(tmpRowind, tmpColind, values);
00320 this->getPermutedAndSymmetrized(tmpRowind, rowInversePermutation, rowind,
00321 tmpColind, colInversePermutation, colind);
00322 }
00333 MatrixSymmetric<Treal, Tmatrix>&
00334 operator=(const MatrixSymmetric<Treal, Tmatrix>& symm) {
00335 MatrixBase<Treal, Tmatrix>::operator=(symm);
00336 return *this;
00337 }
00338 MatrixSymmetric<Treal, Tmatrix>&
00339 operator=(const MatrixGeneral<Treal, Tmatrix>& matr) {
00340 MatrixBase<Treal, Tmatrix>::operator=(matr);
00341 this->matrixPtr->nosymToSym();
00342 return *this;
00343 }
00344 inline MatrixSymmetric<Treal, Tmatrix>& operator=(int const k) {
00345 *this->matrixPtr = k;
00346 return *this;
00347 }
00348
00349 inline Treal frob() const {
00350 return this->matrixPtr->syFrob();
00351 }
00352 Treal mixed_norm(Treal const requestedAccuracy,
00353 int maxIter = -1) const;
00354 Treal eucl(Treal const requestedAccuracy,
00355 int maxIter = -1) const;
00356
00357 void quickEuclBounds(Treal & euclLowerBound,
00358 Treal & euclUpperBound) const {
00359 Treal frobTmp = frob();
00360 euclLowerBound = frobTmp / template_blas_sqrt( (Treal)this->get_nrows() );
00361 euclUpperBound = frobTmp;
00362 }
00363
00364
00370 static Interval<Treal>
00371 diff(const MatrixSymmetric<Treal, Tmatrix>& A,
00372 const MatrixSymmetric<Treal, Tmatrix>& B,
00373 normType const norm,
00374 Treal const requestedAccuracy);
00382 static Interval<Treal>
00383 diffIfSmall(const MatrixSymmetric<Treal, Tmatrix>& A,
00384 const MatrixSymmetric<Treal, Tmatrix>& B,
00385 normType const norm,
00386 Treal const requestedAccuracy,
00387 Treal const maxAbsVal);
00391 static inline Treal frob_diff
00392 (const MatrixSymmetric<Treal, Tmatrix>& A,
00393 const MatrixSymmetric<Treal, Tmatrix>& B) {
00394 return Tmatrix::syFrobDiff(*A.matrixPtr, *B.matrixPtr);
00395 }
00396
00400 static Treal eucl_diff
00401 (const MatrixSymmetric<Treal, Tmatrix>& A,
00402 const MatrixSymmetric<Treal, Tmatrix>& B,
00403 Treal const requestedAccuracy);
00404
00408 static Treal mixed_diff
00409 (const MatrixSymmetric<Treal, Tmatrix>& A,
00410 const MatrixSymmetric<Treal, Tmatrix>& B,
00411 Treal const requestedAccuracy);
00412
00419 static Interval<Treal> euclDiffIfSmall
00420 (const MatrixSymmetric<Treal, Tmatrix>& A,
00421 const MatrixSymmetric<Treal, Tmatrix>& B,
00422 Treal const requestedAccuracy,
00423 Treal const maxAbsVal,
00424 VectorType * const eVecPtr = 0);
00425
00426
00437 Treal thresh(Treal const threshold,
00438 normType const norm);
00439
00446 inline Treal frob_thresh(Treal const threshold) {
00447 return 2.0 * this->matrixPtr->frob_thresh(threshold / 2);
00448 }
00449
00450 Treal eucl_thresh(Treal const threshold,
00451 MatrixTriangular<Treal, Tmatrix> const * const Zptr = NULL);
00452
00453 Treal eucl_element_level_thresh(Treal const threshold);
00454
00455 void getSizesAndBlocksForFrobNormMat
00456 ( SizesAndBlocks & rows_new, SizesAndBlocks & cols_new ) const;
00457
00458 Treal mixed_norm_thresh(Treal const threshold);
00459
00460 void simple_blockwise_frob_thresh(Treal const threshold) {
00461 this->matrixPtr->frobThreshLowestLevel(threshold*threshold, 0);
00462 }
00463
00464 inline void gersgorin(Treal& lmin, Treal& lmax) const {
00465 this->matrixPtr->sy_gersgorin(lmin, lmax);
00466 }
00467 static inline Treal trace_ab
00468 (const MatrixSymmetric<Treal, Tmatrix>& A,
00469 const MatrixSymmetric<Treal, Tmatrix>& B) {
00470 return Tmatrix::sy_trace_ab(*A.matrixPtr, *B.matrixPtr);
00471 }
00472 inline size_t nnz() const {
00473 return this->matrixPtr->sy_nnz();
00474 }
00475 inline size_t nvalues() const {
00476 return this->matrixPtr->sy_nvalues();
00477 }
00478 inline void write_to_buffer(void* buffer, const int n_bytes) const {
00479 this->write_to_buffer_base(buffer, n_bytes, matrix_symm);
00480 }
00481 inline void read_from_buffer(void* buffer, const int n_bytes) {
00482 this->read_from_buffer_base(buffer, n_bytes, matrix_symm);
00483 }
00484
00485
00487 MatrixSymmetric<Treal, Tmatrix>& operator=
00488 (XY<Treal, MatrixSymmetric<Treal, Tmatrix> > const & sm);
00490 MatrixSymmetric<Treal, Tmatrix>& operator=
00491 (const XYZpUV<Treal,
00492 MatrixSymmetric<Treal, Tmatrix>,
00493 MatrixSymmetric<Treal, Tmatrix>,
00494 Treal,
00495 MatrixSymmetric<Treal, Tmatrix> >& sm2psm);
00497 MatrixSymmetric<Treal, Tmatrix>& operator=
00498 (const XYZ<Treal,
00499 MatrixSymmetric<Treal, Tmatrix>,
00500 MatrixSymmetric<Treal, Tmatrix> >& sm2);
00502 MatrixSymmetric<Treal, Tmatrix>& operator+=
00503 (const XYZ<Treal,
00504 MatrixSymmetric<Treal, Tmatrix>,
00505 MatrixSymmetric<Treal, Tmatrix> >& sm2);
00507 MatrixSymmetric<Treal, Tmatrix>& operator=
00508 (const XYZpUV<Treal,
00509 MatrixGeneral<Treal, Tmatrix>,
00510 MatrixGeneral<Treal, Tmatrix>,
00511 Treal,
00512 MatrixSymmetric<Treal, Tmatrix> >& smmpsm);
00514 MatrixSymmetric<Treal, Tmatrix>& operator=
00515 (const XYZ<Treal,
00516 MatrixGeneral<Treal, Tmatrix>,
00517 MatrixGeneral<Treal, Tmatrix> >& smm);
00519 MatrixSymmetric<Treal, Tmatrix>& operator+=
00520 (const XYZ<Treal,
00521 MatrixGeneral<Treal, Tmatrix>,
00522 MatrixGeneral<Treal, Tmatrix> >& smm);
00523
00527 MatrixSymmetric<Treal, Tmatrix>& operator=
00528 (const XYZ<MatrixTriangular<Treal, Tmatrix>,
00529 MatrixSymmetric<Treal, Tmatrix>,
00530 MatrixTriangular<Treal, Tmatrix> >& zaz);
00531
00536 static void ssmmUpperTriangleOnly(const Treal alpha,
00537 const MatrixSymmetric<Treal, Tmatrix>& A,
00538 const MatrixSymmetric<Treal, Tmatrix>& B,
00539 const Treal beta,
00540 MatrixSymmetric<Treal, Tmatrix>& C);
00541
00542
00543
00545 MatrixSymmetric<Treal, Tmatrix>& operator=
00546 (XpY<MatrixSymmetric<Treal, Tmatrix>,
00547 MatrixSymmetric<Treal, Tmatrix> > const & mpm);
00549 MatrixSymmetric<Treal, Tmatrix>& operator=
00550 (XmY<MatrixSymmetric<Treal, Tmatrix>,
00551 MatrixSymmetric<Treal, Tmatrix> > const & mm);
00553 MatrixSymmetric<Treal, Tmatrix>& operator+=
00554 (MatrixSymmetric<Treal, Tmatrix> const & A);
00556 MatrixSymmetric<Treal, Tmatrix>& operator-=
00557 (MatrixSymmetric<Treal, Tmatrix> const & A);
00558
00560 MatrixSymmetric<Treal, Tmatrix>& operator+=
00561 (XY<Treal, MatrixSymmetric<Treal, Tmatrix> > const & sm);
00562
00564 MatrixSymmetric<Treal, Tmatrix>& operator-=
00565 (XY<Treal, MatrixSymmetric<Treal, Tmatrix> > const & sm);
00566
00567 template<typename Top>
00568 Treal accumulateWith(Top & op) {
00569 return this->matrixPtr->syAccumulateWith(op);
00570 }
00571
00572 void random() {
00573 this->matrixPtr->syRandom();
00574 }
00575
00576 void randomZeroStructure(Treal probabilityBeingZero) {
00577 this->matrixPtr->syRandomZeroStructure(probabilityBeingZero);
00578 }
00579
00584 template<typename TRule>
00585 void setElementsByRule(TRule & rule) {
00586 this->matrixPtr->sySetElementsByRule(rule);
00587 return;
00588 }
00589
00592 void transfer( MatrixSymmetric<Treal, Tmatrix> & dest ) {
00593 ValidPtr<Tmatrix>::swap( this->matrixPtr, dest.matrixPtr );
00594
00595 this->clear();
00596 }
00597
00598 template<typename Tvector>
00599 void matVecProd(Tvector & y, Tvector const & x) const {
00600 Treal const ONE = 1.0;
00601 y = (ONE * (*this) * x);
00602 }
00603
00604
00605 std::string obj_type_id() const {return "MatrixSymmetric";}
00606 protected:
00607 inline void writeToFileProt(std::ofstream & file) const {
00608 this->writeToFileBase(file, matrix_symm);
00609 }
00610 inline void readFromFileProt(std::ifstream & file) {
00611 this->readFromFileBase(file, matrix_symm);
00612 }
00613
00619 static void getPermutedAndSymmetrized
00620 (std::vector<int> const & rowind,
00621 std::vector<int> const & rowPermutation,
00622 std::vector<int> & newRowind,
00623 std::vector<int> const & colind,
00624 std::vector<int> const & colPermutation,
00625 std::vector<int> & newColind) {
00626 MatrixBase<Treal, Tmatrix>::
00627 getPermutedIndexes(rowind, rowPermutation, newRowind);
00628 MatrixBase<Treal, Tmatrix>::
00629 getPermutedIndexes(colind, colPermutation, newColind);
00630 int tmp;
00631 for (unsigned int i = 0; i < newRowind.size(); ++i) {
00632 if (newRowind[i] > newColind[i]) {
00633 tmp = newRowind[i];
00634 newRowind[i] = newColind[i];
00635 newColind[i] = tmp;
00636 }
00637 }
00638 }
00639 private:
00640 };
00641
00642 template<typename Treal, typename Tmatrix>
00643 Treal MatrixSymmetric<Treal, Tmatrix>::
00644 mixed_norm(Treal const requestedAccuracy,
00645 int maxIter) const {
00646
00647 SizesAndBlocks rows_new;
00648 SizesAndBlocks cols_new;
00649 this->getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
00650
00651
00652 MatrixSymmetric<Treal, typename Tmatrix::ElementType> frobNormMat;
00653 frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
00654 frobNormMat.getMatrix().syAssignFrobNormsLowestLevel(this->getMatrix());
00655 return frobNormMat.eucl(requestedAccuracy, maxIter);
00656 }
00657
00658
00659 template<typename Treal, typename Tmatrix>
00660 Treal MatrixSymmetric<Treal, Tmatrix>::
00661 eucl(Treal const requestedAccuracy,
00662 int maxIter) const {
00663 assert(requestedAccuracy >= 0);
00664
00665 Treal frobTmp = this->frob();
00666 if (frobTmp < requestedAccuracy)
00667 return (Treal)0.0;
00668 if (maxIter < 0)
00669 maxIter = this->get_nrows() * 100;
00670 VectorType guess;
00671 SizesAndBlocks cols;
00672 this->getCols(cols);
00673 guess.resetSizesAndBlocks(cols);
00674 guess.rand();
00675
00676 #if 0 // "new code"
00677 MatrixSymmetric<Treal, Tmatrix> Copy(*this);
00678 Copy.frob_thresh(requestedAccuracy / 2.0);
00679 arn::LanczosLargestMagnitudeEig
00680 <Treal, MatrixSymmetric<Treal, Tmatrix>, VectorType>
00681 lan(Copy, guess, maxIter);
00682 lan.setAbsTol( requestedAccuracy / 2.0 );
00683 #else // "old code"
00684 arn::LanczosLargestMagnitudeEig
00685 <Treal, MatrixSymmetric<Treal, Tmatrix>, VectorType>
00686 lan(*this, guess, maxIter);
00687 lan.setAbsTol( requestedAccuracy );
00688 #endif
00689 lan.run();
00690 Treal eVal = 0;
00691 Treal acc = 0;
00692 lan.getLargestMagnitudeEig(eVal, acc);
00693 return template_blas_fabs(eVal);
00694 }
00695
00696 template<typename Treal, typename Tmatrix>
00697 Interval<Treal> MatrixSymmetric<Treal, Tmatrix>::
00698 diff(const MatrixSymmetric<Treal, Tmatrix>& A,
00699 const MatrixSymmetric<Treal, Tmatrix>& B,
00700 normType const norm, Treal const requestedAccuracy) {
00701 Treal diff;
00702 Treal eNMin;
00703 switch (norm) {
00704 case frobNorm:
00705 diff = frob_diff(A, B);
00706 return Interval<Treal>(diff / template_blas_sqrt((Treal)A.get_nrows()), diff);
00707 break;
00708 case euclNorm:
00709 diff = eucl_diff(A, B, requestedAccuracy);
00710 eNMin = diff - requestedAccuracy;
00711 eNMin = eNMin >= 0 ? eNMin : 0;
00712 return Interval<Treal>(eNMin, diff + requestedAccuracy);
00713 break;
00714 default:
00715 throw Failure("MatrixSymmetric<Treal, Tmatrix>::"
00716 "diff(const MatrixSymmetric<Treal, Tmatrix>&, "
00717 "const MatrixSymmetric<Treal, Tmatrix>&, "
00718 "normType const, Treal): "
00719 "Diff not implemented for selected norm");
00720 }
00721 }
00722
00723 #if 1
00724 template<typename Treal, typename Tmatrix>
00725 Interval<Treal> MatrixSymmetric<Treal, Tmatrix>::
00726 diffIfSmall(const MatrixSymmetric<Treal, Tmatrix>& A,
00727 const MatrixSymmetric<Treal, Tmatrix>& B,
00728 normType const norm,
00729 Treal const requestedAccuracy,
00730 Treal const maxAbsVal) {
00731 Treal diff;
00732 switch (norm) {
00733 case frobNorm:
00734 {
00735 diff = frob_diff(A, B);
00736 return Interval<Treal>(diff / template_blas_sqrt((Treal)A.get_nrows()), diff);
00737 }
00738 break;
00739 case euclNorm:
00740 return euclDiffIfSmall(A, B, requestedAccuracy, maxAbsVal);
00741 break;
00742 default:
00743 throw Failure("MatrixSymmetric<Treal, Tmatrix>::"
00744 "diffIfSmall"
00745 "(const MatrixSymmetric<Treal, Tmatrix>&, "
00746 "const MatrixSymmetric<Treal, Tmatrix>&, "
00747 "normType const, Treal const, Treal const): "
00748 "Diff not implemented for selected norm");
00749 }
00750 }
00751
00752 #endif
00753
00754
00755 template<typename Treal, typename Tmatrix>
00756 Treal MatrixSymmetric<Treal, Tmatrix>::eucl_diff
00757 (const MatrixSymmetric<Treal, Tmatrix>& A,
00758 const MatrixSymmetric<Treal, Tmatrix>& B,
00759 Treal const requestedAccuracy) {
00760
00761 mat::DiffMatrix< MatrixSymmetric<Treal, Tmatrix>, Treal> Diff(A, B);
00762 Treal maxAbsVal = 2 * frob_diff(A,B);
00763
00764
00765 Treal relTol = sqrt(sqrt(std::numeric_limits<Treal>::epsilon()));
00766 Interval<Treal> euclInt =
00767 mat::euclIfSmall(Diff, requestedAccuracy, relTol, maxAbsVal);
00768 return euclInt.midPoint();
00769 }
00770
00771 template<typename Treal, typename Tmatrix>
00772 Treal MatrixSymmetric<Treal, Tmatrix>::mixed_diff
00773 (const MatrixSymmetric<Treal, Tmatrix>& A,
00774 const MatrixSymmetric<Treal, Tmatrix>& B,
00775 Treal const requestedAccuracy) {
00776 MatrixSymmetric<Treal, typename Tmatrix::ElementType> frobNormMat;
00777 {
00778 SizesAndBlocks rows_new;
00779 SizesAndBlocks cols_new;
00780 A.getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
00781 frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
00782 frobNormMat.getMatrix().syAssignDiffFrobNormsLowestLevel(A.getMatrix(),B.getMatrix());
00783 }
00784 return frobNormMat.eucl(requestedAccuracy);
00785 }
00786
00787 #if 1
00788 template<typename Treal, typename Tmatrix>
00789 Interval<Treal> MatrixSymmetric<Treal, Tmatrix>::euclDiffIfSmall
00790 (const MatrixSymmetric<Treal, Tmatrix>& A,
00791 const MatrixSymmetric<Treal, Tmatrix>& B,
00792 Treal const requestedAccuracy,
00793 Treal const maxAbsVal,
00794 VectorType * const eVecPtr) {
00795
00796 mat::DiffMatrix< MatrixSymmetric<Treal, Tmatrix>, Treal> Diff(A, B);
00797
00798 Treal relTol = sqrt(sqrt(std::numeric_limits<Treal>::epsilon()));
00799 Interval<Treal> tmpInterval = mat::euclIfSmall(Diff, requestedAccuracy, relTol, maxAbsVal, eVecPtr);
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809 if ( tmpInterval.length() < 2*requestedAccuracy )
00810 return Interval<Treal>( tmpInterval.midPoint()-requestedAccuracy, tmpInterval.midPoint()+requestedAccuracy );
00811 return tmpInterval;
00812 }
00813
00814 #endif
00815
00816
00817 template<typename Treal, typename Tmatrix>
00818 Treal MatrixSymmetric<Treal, Tmatrix>::
00819 thresh(Treal const threshold,
00820 normType const norm) {
00821 switch (norm) {
00822 case frobNorm:
00823 return this->frob_thresh(threshold);
00824 break;
00825 case euclNorm:
00826 return this->eucl_thresh(threshold);
00827 break;
00828 case mixedNorm:
00829 return this->mixed_norm_thresh(threshold);
00830 break;
00831 default:
00832 throw Failure("MatrixSymmetric<Treal, Tmatrix>::"
00833 "thresh(Treal const, "
00834 "normType const): "
00835 "Thresholding not implemented for selected norm");
00836 }
00837 }
00838
00839 #if 1
00840
00841 template<typename Treal, typename Tmatrix>
00842 Treal MatrixSymmetric<Treal, Tmatrix>::
00843 eucl_thresh(Treal const threshold,
00844 MatrixTriangular<Treal, Tmatrix> const * const Zptr) {
00845 if ( Zptr == NULL ) {
00846 EuclTruncationSymm<MatrixSymmetric<Treal, Tmatrix>, Treal> TruncObj(*this);
00847 return TruncObj.run( threshold );
00848 }
00849 EuclTruncationSymmWithZ<MatrixSymmetric<Treal, Tmatrix>, MatrixTriangular<Treal, Tmatrix>, Treal> TruncObj(*this, *Zptr);
00850 return TruncObj.run( threshold );
00851 }
00852
00853 #endif
00854
00855
00856 template<typename Treal, typename Tmatrix>
00857 void MatrixSymmetric<Treal, Tmatrix>::getSizesAndBlocksForFrobNormMat
00858 ( SizesAndBlocks & rows_new, SizesAndBlocks & cols_new ) const {
00859 std::vector<int> rows_block_sizes;
00860 std::vector<int> cols_block_sizes;
00861 int n_rows;
00862 int n_cols;
00863 {
00864 SizesAndBlocks rows;
00865 SizesAndBlocks cols;
00866 this->getRows(rows);
00867 this->getCols(cols);
00868 rows.getBlockSizeVector( rows_block_sizes );
00869 cols.getBlockSizeVector( cols_block_sizes );
00870 rows_block_sizes.pop_back();
00871 cols_block_sizes.pop_back();
00872 n_rows = rows.getNTotalScalars();
00873 n_cols = cols.getNTotalScalars();
00874 int factor_rows = rows_block_sizes[rows_block_sizes.size()-1];
00875 int factor_cols = cols_block_sizes[cols_block_sizes.size()-1];
00876 for (unsigned int ind = 0; ind < rows_block_sizes.size(); ++ind)
00877 rows_block_sizes[ind] = rows_block_sizes[ind] / factor_rows;
00878 for (unsigned int ind = 0; ind < cols_block_sizes.size(); ++ind)
00879 cols_block_sizes[ind] = cols_block_sizes[ind] / factor_cols;
00880
00881
00882
00883 if (n_rows % factor_rows)
00884 n_rows = n_rows / factor_rows + 1;
00885 else
00886 n_rows = n_rows / factor_rows;
00887 if (n_cols % factor_cols)
00888 n_cols = n_cols / factor_cols + 1;
00889 else
00890 n_cols = n_cols / factor_cols;
00891 }
00892 rows_new = SizesAndBlocks( rows_block_sizes, n_rows );
00893 cols_new = SizesAndBlocks( cols_block_sizes, n_cols );
00894 }
00895
00896
00897 template<typename Treal, typename Tmatrix>
00898 Treal MatrixSymmetric<Treal, Tmatrix>::
00899 mixed_norm_thresh(Treal const threshold) {
00900 assert(threshold >= (Treal)0.0);
00901 if (threshold == (Treal)0.0)
00902 return (Treal)0;
00903
00904 SizesAndBlocks rows_new;
00905 SizesAndBlocks cols_new;
00906 this->getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
00907
00908
00909 MatrixSymmetric<Treal, typename Tmatrix::ElementType> frobNormMat;
00910 frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
00911
00912
00913
00914 frobNormMat.getMatrix().syAssignFrobNormsLowestLevel(this->getMatrix());
00915 EuclTruncationSymmElementLevel<MatrixSymmetric<Treal, typename Tmatrix::ElementType>, Treal> TruncObj( frobNormMat );
00916 Treal mixed_norm_result = TruncObj.run( threshold );
00917 frobNormMat.getMatrix().truncateAccordingToSparsityPattern(this->getMatrix());
00918 return mixed_norm_result;
00919 }
00920
00921
00922
00923 template<typename Treal, typename Tmatrix>
00924 inline MatrixSymmetric<Treal, Tmatrix>&
00925 MatrixSymmetric<Treal, Tmatrix>::operator=
00926 (XY<Treal, MatrixSymmetric<Treal, Tmatrix> > const & sm) {
00927 assert(!sm.tB);
00928
00929 this->matrixPtr.haveDataStructureSet(true);
00930 this->matrixPtr->assign(sm.A, *sm.B.matrixPtr);
00931 return *this;
00932 }
00933
00934 template<typename Treal, typename Tmatrix>
00935 MatrixSymmetric<Treal, Tmatrix>&
00936 MatrixSymmetric<Treal, Tmatrix>::operator=
00937 (const XYZpUV<Treal,
00938 MatrixSymmetric<Treal, Tmatrix>,
00939 MatrixSymmetric<Treal, Tmatrix>,
00940 Treal,
00941 MatrixSymmetric<Treal, Tmatrix> >& sm2psm) {
00942 assert(this != &sm2psm.B);
00943 if (this == &sm2psm.E && &sm2psm.B == &sm2psm.C) {
00944
00945 Tmatrix::sysq('U',
00946 sm2psm.A, *sm2psm.B.matrixPtr,
00947 sm2psm.D, *this->matrixPtr);
00948 return *this;
00949 }
00950 else
00951 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
00952 "(const XYZpUV<Treal, MatrixSymmetric"
00953 "<Treal, Tmatrix> >& sm2psm) : "
00954 "D = alpha * A * B + beta * C not supported for C != D"
00955 " and for symmetric matrices not for A != B since this "
00956 "generally will result in a nonsymmetric matrix");
00957 }
00958
00959
00960 template<typename Treal, typename Tmatrix>
00961 MatrixSymmetric<Treal, Tmatrix>&
00962 MatrixSymmetric<Treal, Tmatrix>::operator=
00963 (const XYZ<Treal,
00964 MatrixSymmetric<Treal, Tmatrix>,
00965 MatrixSymmetric<Treal, Tmatrix> >& sm2) {
00966 assert(this != &sm2.B);
00967 if (&sm2.B == &sm2.C) {
00968 this->matrixPtr.haveDataStructureSet(true);
00969 Tmatrix::sysq('U', sm2.A, *sm2.B.matrixPtr, 0, *this->matrixPtr);
00970 return *this;
00971 }
00972 else
00973 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
00974 "(const XYZ<Treal, MatrixSymmetric<Treal, Tmatrix>,"
00975 " MatrixSymmetric<Treal, Tmatrix> >& sm2) : "
00976 "Operation C = alpha * A * B with only symmetric "
00977 "matrices not supported for A != B");
00978 }
00979
00980
00981 template<typename Treal, typename Tmatrix>
00982 MatrixSymmetric<Treal, Tmatrix>&
00983 MatrixSymmetric<Treal, Tmatrix>::operator+=
00984 (const XYZ<Treal,
00985 MatrixSymmetric<Treal, Tmatrix>,
00986 MatrixSymmetric<Treal, Tmatrix> >& sm2) {
00987 assert(this != &sm2.B);
00988 if (&sm2.B == &sm2.C) {
00989 Tmatrix::sysq('U', sm2.A, *sm2.B.matrixPtr, 1, *this->matrixPtr);
00990 return *this;
00991 }
00992 else
00993 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator+="
00994 "(const XYZ<Treal, MatrixSymmetric<Treal, Tmatrix>,"
00995 " MatrixSymmetric<Treal, Tmatrix> >& sm2) : "
00996 "Operation C += alpha * A * B with only symmetric "
00997 "matrices not supported for A != B");
00998 }
00999
01000
01001 template<typename Treal, typename Tmatrix>
01002 MatrixSymmetric<Treal, Tmatrix>&
01003 MatrixSymmetric<Treal, Tmatrix>::operator=
01004 (const XYZpUV<Treal,
01005 MatrixGeneral<Treal, Tmatrix>,
01006 MatrixGeneral<Treal, Tmatrix>,
01007 Treal,
01008 MatrixSymmetric<Treal, Tmatrix> >& smmpsm) {
01009 if (this == &smmpsm.E)
01010 if (&smmpsm.B == &smmpsm.C)
01011 if (!smmpsm.tB && smmpsm.tC) {
01012 Tmatrix::syrk('U', false,
01013 smmpsm.A, *smmpsm.B.matrixPtr,
01014 smmpsm.D, *this->matrixPtr);
01015 }
01016 else
01017 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01018 "(const XYZpUV<Treal, MatrixGeneral"
01019 "<Treal, Tmatrix>, "
01020 "MatrixGeneral<Treal, Tmatrix>, Treal, "
01021 "MatrixSymmetric<Treal, Tmatrix> >&) : "
01022 "C = alpha * A' * A + beta * C, not implemented"
01023 " only C = alpha * A * A' + beta * C");
01024 else
01025 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01026 "(const XYZpUV<"
01027 "Treal, MatrixGeneral<Treal, Tmatrix>, "
01028 "MatrixGeneral<Treal, Tmatrix>, Treal, "
01029 "MatrixSymmetric<Treal, Tmatrix> >&) : "
01030 "You are trying to call C = alpha * A * A' + beta * C"
01031 " with A and A' being different objects");
01032 else
01033 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01034 "(const XYZpUV<"
01035 "Treal, MatrixGeneral<Treal, Tmatrix>, "
01036 "MatrixGeneral<Treal, Tmatrix>, Treal, "
01037 "MatrixSymmetric<Treal, Tmatrix> >&) : "
01038 "D = alpha * A * A' + beta * C not supported for C != D");
01039 return *this;
01040 }
01041
01042
01043 template<typename Treal, typename Tmatrix>
01044 MatrixSymmetric<Treal, Tmatrix>&
01045 MatrixSymmetric<Treal, Tmatrix>::operator=
01046 (const XYZ<Treal,
01047 MatrixGeneral<Treal, Tmatrix>,
01048 MatrixGeneral<Treal, Tmatrix> >& smm) {
01049 if (&smm.B == &smm.C)
01050 if (!smm.tB && smm.tC) {
01051 Tmatrix::syrk('U', false,
01052 smm.A, *smm.B.matrixPtr,
01053 0, *this->matrixPtr);
01054 }
01055 else
01056 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01057 "(const XYZ<"
01058 "Treal, MatrixGeneral<Treal, Tmatrix>, "
01059 "MatrixGeneral<Treal, Tmatrix> >&) : "
01060 "C = alpha * A' * A, not implemented "
01061 "only C = alpha * A * A'");
01062 else
01063 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01064 "(const XYZ<"
01065 "Treal, MatrixGeneral<Treal, Tmatrix>, "
01066 "MatrixGeneral<Treal, Tmatrix> >&) : "
01067 "You are trying to call C = alpha * A * A' "
01068 "with A and A' being different objects");
01069 return *this;
01070 }
01071
01072
01073 template<typename Treal, typename Tmatrix>
01074 MatrixSymmetric<Treal, Tmatrix>&
01075 MatrixSymmetric<Treal, Tmatrix>::operator+=
01076 (const XYZ<Treal,
01077 MatrixGeneral<Treal, Tmatrix>,
01078 MatrixGeneral<Treal, Tmatrix> >& smm) {
01079 if (&smm.B == &smm.C)
01080 if (!smm.tB && smm.tC) {
01081 Tmatrix::syrk('U', false,
01082 smm.A, *smm.B.matrixPtr,
01083 1, *this->matrixPtr);
01084 }
01085 else
01086 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator+="
01087 "(const XYZ<"
01088 "Treal, MatrixGeneral<Treal, Tmatrix>, "
01089 "MatrixGeneral<Treal, Tmatrix> >&) : "
01090 "C += alpha * A' * A, not implemented "
01091 "only C += alpha * A * A'");
01092 else
01093 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator+="
01094 "(const XYZ<"
01095 "Treal, MatrixGeneral<Treal, Tmatrix>, "
01096 "MatrixGeneral<Treal, Tmatrix> >&) : "
01097 "You are trying to call C += alpha * A * A' "
01098 "with A and A' being different objects");
01099 return *this;
01100 }
01101
01102 #if 1
01103
01104
01105 template<typename Treal, typename Tmatrix>
01106 MatrixSymmetric<Treal, Tmatrix>&
01107 MatrixSymmetric<Treal, Tmatrix>::operator=
01108 (const XYZ<MatrixTriangular<Treal, Tmatrix>,
01109 MatrixSymmetric<Treal, Tmatrix>,
01110 MatrixTriangular<Treal, Tmatrix> >& zaz) {
01111 if (this == &zaz.B) {
01112 if (&zaz.A == &zaz.C) {
01113 if (zaz.tA && !zaz.tC) {
01114 Tmatrix::trsytriplemm('R', *zaz.A.matrixPtr, *this->matrixPtr);
01115 }
01116 else if (!zaz.tA && zaz.tC) {
01117 Tmatrix::trsytriplemm('L', *zaz.A.matrixPtr, *this->matrixPtr);
01118 }
01119 else
01120 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01121 "(const XYZ<MatrixTriangular<Treal, Tmatrix>,"
01122 "MatrixSymmetric<Treal, Tmatrix>,"
01123 "MatrixTriangular<Treal, Tmatrix> >&) : "
01124 "A = op1(Z) * A * op2(Z) : Either op1 xor op2 must be "
01125 "the transpose operation.");
01126 }
01127 else
01128 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01129 "(const XYZ<MatrixTriangular<Treal, Tmatrix>,"
01130 "MatrixSymmetric<Treal, Tmatrix>,"
01131 "MatrixTriangular<Treal, Tmatrix> >&) : "
01132 "A = op1(Z1) * A * op2(Z2) : Z1 and Z2 must be the same "
01133 "object");
01134 }
01135 else
01136 throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
01137 "(const XYZ<MatrixTriangular<Treal, Tmatrix>,"
01138 "MatrixSymmetric<Treal, Tmatrix>,"
01139 "MatrixTriangular<Treal, Tmatrix> >&) : "
01140 "C = op1(Z) * A * op2(Z) : A and C must be the same "
01141 "object");
01142 return *this;
01143 }
01144
01145 #endif
01146
01147
01152 template<typename Treal, typename Tmatrix>
01153 void MatrixSymmetric<Treal, Tmatrix>::
01154 ssmmUpperTriangleOnly(const Treal alpha,
01155 const MatrixSymmetric<Treal, Tmatrix>& A,
01156 const MatrixSymmetric<Treal, Tmatrix>& B,
01157 const Treal beta,
01158 MatrixSymmetric<Treal, Tmatrix>& C) {
01159 Tmatrix::ssmm_upper_tr_only(alpha, *A.matrixPtr, *B.matrixPtr,
01160 beta, *C.matrixPtr);
01161 }
01162
01163
01164
01165
01166
01167 template<typename Treal, typename Tmatrix>
01168 inline MatrixSymmetric<Treal, Tmatrix>&
01169 MatrixSymmetric<Treal, Tmatrix>::operator=
01170 (const XpY<MatrixSymmetric<Treal, Tmatrix>,
01171 MatrixSymmetric<Treal, Tmatrix> >& mpm) {
01172 assert(this != &mpm.A);
01173 (*this) = mpm.B;
01174 Tmatrix::add(1.0, *mpm.A.matrixPtr, *this->matrixPtr);
01175 return *this;
01176 }
01177
01178 template<typename Treal, typename Tmatrix>
01179 inline MatrixSymmetric<Treal, Tmatrix>&
01180 MatrixSymmetric<Treal, Tmatrix>::operator=
01181 (const XmY<MatrixSymmetric<Treal, Tmatrix>,
01182 MatrixSymmetric<Treal, Tmatrix> >& mmm) {
01183 assert(this != &mmm.B);
01184 (*this) = mmm.A;
01185 Tmatrix::add(-1.0, *mmm.B.matrixPtr, *this->matrixPtr);
01186 return *this;
01187 }
01188
01189 template<typename Treal, typename Tmatrix>
01190 inline MatrixSymmetric<Treal, Tmatrix>&
01191 MatrixSymmetric<Treal, Tmatrix>::operator+=
01192 (MatrixSymmetric<Treal, Tmatrix> const & A) {
01193 Tmatrix::add(1.0, *A.matrixPtr, *this->matrixPtr);
01194 return *this;
01195 }
01196
01197 template<typename Treal, typename Tmatrix>
01198 inline MatrixSymmetric<Treal, Tmatrix>&
01199 MatrixSymmetric<Treal, Tmatrix>::operator-=
01200 (MatrixSymmetric<Treal, Tmatrix> const & A) {
01201 Tmatrix::add(-1.0, *A.matrixPtr, *this->matrixPtr);
01202 return *this;
01203 }
01204
01205 template<typename Treal, typename Tmatrix>
01206 inline MatrixSymmetric<Treal, Tmatrix>&
01207 MatrixSymmetric<Treal, Tmatrix>::operator+=
01208 (XY<Treal, MatrixSymmetric<Treal, Tmatrix> > const & sm) {
01209 assert(!sm.tB);
01210 Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
01211 return *this;
01212 }
01213
01214
01215 template<typename Treal, typename Tmatrix>
01216 inline MatrixSymmetric<Treal, Tmatrix>&
01217 MatrixSymmetric<Treal, Tmatrix>::operator-=
01218 (XY<Treal, MatrixSymmetric<Treal, Tmatrix> > const & sm) {
01219 assert(!sm.tB);
01220 Tmatrix::add(-sm.A, *sm.B.matrixPtr, *this->matrixPtr);
01221 return *this;
01222 }
01223
01228 template<typename Treal, typename Tmatrix, typename Top>
01229 Treal accumulate(MatrixSymmetric<Treal, Tmatrix> & A,
01230 Top & op) {
01231 return A.accumulateWith(op);
01232 }
01233
01234
01235
01236 }
01237 #endif
01238