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_MATRIXGENERAL
00037 #define MAT_MATRIXGENERAL
00038 #include "MatrixBase.h"
00039 namespace mat {
00056 template<typename Treal, typename Tmatrix>
00057 class MatrixGeneral : public MatrixBase<Treal, Tmatrix> {
00058 public:
00059 typedef VectorGeneral<Treal, typename Tmatrix::VectorType> VectorType;
00060
00061 MatrixGeneral()
00062 :MatrixBase<Treal, Tmatrix>() {}
00063 explicit MatrixGeneral(const MatrixGeneral<Treal, Tmatrix>& matr)
00064 :MatrixBase<Treal, Tmatrix>(matr) {}
00065 explicit MatrixGeneral(const MatrixSymmetric<Treal, Tmatrix>& symm)
00066 :MatrixBase<Treal, Tmatrix>(symm) {
00067 this->matrixPtr->symToNosym();
00068 }
00069 explicit MatrixGeneral(const MatrixTriangular<Treal, Tmatrix>& triang)
00070 :MatrixBase<Treal, Tmatrix>(triang) {}
00073 #if 0
00074 template<typename Tfull>
00075 inline void assign_from_full
00076 (Tfull const* const fullmatrix,
00077 const int totnrows, const int totncols) {
00078 this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
00079 }
00080 inline void assign_from_full
00081 (Treal const* const fullmatrix,
00082 const int totnrows, const int totncols) {
00083 this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
00084 }
00085 #endif
00086
00087 inline void assignFromFull
00088 (std::vector<Treal> const & fullMat) {
00089 assert((int)fullMat.size() == this->get_nrows() * this->get_ncols());
00090 this->matrixPtr->assignFromFull(fullMat);
00091 }
00092
00093 inline void fullMatrix(std::vector<Treal> & fullMat) const {
00094 this->matrixPtr->fullMatrix(fullMat);
00095 }
00096
00097 inline void fullMatrix
00098 (std::vector<Treal> & fullMat,
00099 std::vector<int> const & rowInversePermutation,
00100 std::vector<int> const & colInversePermutation) const {
00101 std::vector<int> rowind;
00102 std::vector<int> colind;
00103 std::vector<Treal> values;
00104 get_all_values(rowind, colind, values,
00105 rowInversePermutation,
00106 colInversePermutation);
00107 fullMat.assign(this->get_nrows()*this->get_ncols(),0);
00108 for (unsigned int ind = 0; ind < values.size(); ++ind)
00109 fullMat[rowind[ind] + this->get_nrows() * colind[ind]] =
00110 values[ind];
00111 }
00119 inline void assign_from_sparse
00120 (std::vector<int> const & rowind,
00121 std::vector<int> const & colind,
00122 std::vector<Treal> const & values,
00123 SizesAndBlocks const & newRows,
00124 SizesAndBlocks const & newCols) {
00125 this->resetSizesAndBlocks(newRows, newCols);
00126 this->matrixPtr->assignFromSparse(rowind, colind, values);
00127 }
00135 inline void assign_from_sparse
00136 (std::vector<int> const & rowind,
00137 std::vector<int> const & colind,
00138 std::vector<Treal> const & values,
00139 std::vector<int> const & rowPermutation,
00140 std::vector<int> const & colPermutation) {
00141 std::vector<int> newRowind;
00142 std::vector<int> newColind;
00143 MatrixBase<Treal, Tmatrix>::
00144 getPermutedIndexes(rowind, rowPermutation, newRowind);
00145 MatrixBase<Treal, Tmatrix>::
00146 getPermutedIndexes(colind, colPermutation, newColind);
00147 this->matrixPtr->assignFromSparse(newRowind, newColind, values);
00148 }
00154 inline void assign_from_sparse
00155 (std::vector<int> const & rowind,
00156 std::vector<int> const & colind,
00157 std::vector<Treal> const & values,
00158 SizesAndBlocks const & newRows,
00159 SizesAndBlocks const & newCols,
00160 std::vector<int> const & rowPermutation,
00161 std::vector<int> const & colPermutation) {
00162 this->resetSizesAndBlocks(newRows, newCols);
00163 assign_from_sparse(rowind, colind, values,
00164 rowPermutation, colPermutation);
00165 }
00170 inline void get_values
00171 (std::vector<int> const & rowind,
00172 std::vector<int> const & colind,
00173 std::vector<Treal> & values) const {
00174 this->matrixPtr->getValues(rowind, colind, values);
00175 }
00183 inline void get_values
00184 (std::vector<int> const & rowind,
00185 std::vector<int> const & colind,
00186 std::vector<Treal> & values,
00187 std::vector<int> const & rowPermutation,
00188 std::vector<int> const & colPermutation) const {
00189 std::vector<int> newRowind;
00190 std::vector<int> newColind;
00191 MatrixBase<Treal, Tmatrix>::
00192 getPermutedIndexes(rowind, rowPermutation, newRowind);
00193 MatrixBase<Treal, Tmatrix>::
00194 getPermutedIndexes(colind, colPermutation, newColind);
00195 this->matrixPtr->getValues(newRowind, newColind, values);
00196 }
00201 inline void get_all_values
00202 (std::vector<int> & rowind,
00203 std::vector<int> & colind,
00204 std::vector<Treal> & values) const {
00205 rowind.resize(0);
00206 colind.resize(0);
00207 values.resize(0);
00208 this->matrixPtr->getAllValues(rowind, colind, values);
00209 }
00219 inline void get_all_values
00220 (std::vector<int> & rowind,
00221 std::vector<int> & colind,
00222 std::vector<Treal> & values,
00223 std::vector<int> const & rowInversePermutation,
00224 std::vector<int> const & colInversePermutation) const {
00225 std::vector<int> tmpRowind;
00226 std::vector<int> tmpColind;
00227 tmpRowind.reserve(rowind.capacity());
00228 tmpColind.reserve(colind.capacity());
00229 values.resize(0);
00230 this->matrixPtr->getAllValues(tmpRowind, tmpColind, values);
00231 MatrixBase<Treal, Tmatrix>::
00232 getPermutedIndexes(tmpRowind, rowInversePermutation, rowind);
00233 MatrixBase<Treal, Tmatrix>::
00234 getPermutedIndexes(tmpColind, colInversePermutation, colind);
00235
00236 }
00246 #if 0
00247 inline void fullmatrix(Treal* const full,
00248 const int totnrows, const int totncols) const {
00249 this->matrixPtr->fullmatrix(full, totnrows, totncols);
00250 }
00251 #endif
00252 MatrixGeneral<Treal, Tmatrix>&
00253 operator=(const MatrixGeneral<Treal, Tmatrix>& mat) {
00254 MatrixBase<Treal, Tmatrix>::operator=(mat);
00255 return *this;
00256 }
00257 inline MatrixGeneral<Treal, Tmatrix>&
00258 operator=(const Xtrans<MatrixGeneral<Treal, Tmatrix> >& mt) {
00259 if (mt.tA)
00260 MatrixBase<Treal, Tmatrix>::operator=(transpose(mt.A));
00261 else
00262 MatrixBase<Treal, Tmatrix>::operator=(mt.A);
00263 return *this;
00264 }
00265
00266 MatrixGeneral<Treal, Tmatrix>&
00267 operator=(const MatrixSymmetric<Treal, Tmatrix>& symm) {
00268 MatrixBase<Treal, Tmatrix>::operator=(symm);
00269 this->matrixPtr->symToNosym();
00270 return *this;
00271 }
00272 MatrixGeneral<Treal, Tmatrix>&
00273 operator=(const MatrixTriangular<Treal, Tmatrix>& triang) {
00274 MatrixBase<Treal, Tmatrix>::operator=(triang);
00275 return *this;
00276 }
00277
00278 inline MatrixGeneral<Treal, Tmatrix>& operator=(int const k) {
00279 *this->matrixPtr = k;
00280 return *this;
00281 }
00282 inline Treal frob() const {
00283 return this->matrixPtr->frob();
00284 }
00285 static inline Treal frob_diff
00286 (const MatrixGeneral<Treal, Tmatrix>& A,
00287 const MatrixGeneral<Treal, Tmatrix>& B) {
00288 return Tmatrix::frobDiff(*A.matrixPtr, *B.matrixPtr);
00289 }
00290 Treal eucl(Treal const requestedAccuracy,
00291 int maxIter = -1) const;
00292
00293
00294 void thresh(Treal const threshold,
00295 normType const norm);
00296
00297 inline void frob_thresh(Treal threshold) {
00298 this->matrixPtr->frob_thresh(threshold);
00299 }
00300 Treal eucl_thresh(Treal const threshold);
00301
00302 inline void gersgorin(Treal& lmin, Treal& lmax) {
00303 this->matrixPtr->gersgorin(lmin, lmax);
00304 }
00305 static inline Treal trace_ab
00306 (const MatrixGeneral<Treal, Tmatrix>& A,
00307 const MatrixGeneral<Treal, Tmatrix>& B) {
00308 return Tmatrix::trace_ab(*A.matrixPtr, *B.matrixPtr);
00309 }
00310 static inline Treal trace_aTb
00311 (const MatrixGeneral<Treal, Tmatrix>& A,
00312 const MatrixGeneral<Treal, Tmatrix>& B) {
00313 return Tmatrix::trace_aTb(*A.matrixPtr, *B.matrixPtr);
00314 }
00315 inline size_t nnz() const {
00316 return this->matrixPtr->nnz();
00317 }
00318 inline size_t nvalues() const {
00319 return this->matrixPtr->nvalues();
00320 }
00321
00322 inline void write_to_buffer(void* buffer, const int n_bytes) const {
00323 this->write_to_buffer_base(buffer, n_bytes, matrix_matr);
00324 }
00325 inline void read_from_buffer(void* buffer, const int n_bytes) {
00326 this->read_from_buffer_base(buffer, n_bytes, matrix_matr);
00327 }
00328
00329
00331 MatrixGeneral<Treal, Tmatrix>& operator=
00332 (const XYZ<Treal,
00333 MatrixGeneral<Treal, Tmatrix>,
00334 MatrixGeneral<Treal, Tmatrix> >& smm);
00335
00337 MatrixGeneral<Treal, Tmatrix>& operator=
00338 (const XY<MatrixGeneral<Treal, Tmatrix>,
00339 MatrixGeneral<Treal, Tmatrix> >& mm);
00340
00342 MatrixGeneral<Treal, Tmatrix>& operator+=
00343 (const XYZ<Treal,
00344 MatrixGeneral<Treal, Tmatrix>,
00345 MatrixGeneral<Treal, Tmatrix> >& smm);
00346
00348 MatrixGeneral<Treal, Tmatrix>& operator=
00349 (const XYZpUV<Treal,
00350 MatrixGeneral<Treal, Tmatrix>,
00351 MatrixGeneral<Treal, Tmatrix>,
00352 Treal,
00353 MatrixGeneral<Treal, Tmatrix> >& smmpsm);
00354
00356 MatrixGeneral<Treal, Tmatrix>& operator=
00357 (XpY<MatrixGeneral<Treal, Tmatrix>,
00358 MatrixGeneral<Treal, Tmatrix> > const & mpm);
00360 MatrixGeneral<Treal, Tmatrix>& operator+=
00361 (MatrixGeneral<Treal, Tmatrix> const & A);
00362 MatrixGeneral<Treal, Tmatrix>& operator-=
00363 (MatrixGeneral<Treal, Tmatrix> const & A);
00365 MatrixGeneral<Treal, Tmatrix>& operator+=
00366 (XY<Treal, MatrixGeneral<Treal, Tmatrix> > const & sm);
00367
00368
00369
00371 MatrixGeneral<Treal, Tmatrix>& operator=
00372 (const XYZ<Treal,
00373 MatrixSymmetric<Treal, Tmatrix>,
00374 MatrixGeneral<Treal, Tmatrix> >& smm);
00376 MatrixGeneral<Treal, Tmatrix>& operator=
00377 (const XY<MatrixSymmetric<Treal, Tmatrix>,
00378 MatrixGeneral<Treal, Tmatrix> >& mm);
00380 MatrixGeneral<Treal, Tmatrix>& operator+=
00381 (const XYZ<Treal,
00382 MatrixSymmetric<Treal, Tmatrix>,
00383 MatrixGeneral<Treal, Tmatrix> >& smm);
00385 MatrixGeneral<Treal, Tmatrix>& operator=
00386 (const XYZpUV<Treal,
00387 MatrixSymmetric<Treal, Tmatrix>,
00388 MatrixGeneral<Treal, Tmatrix>,
00389 Treal,
00390 MatrixGeneral<Treal, Tmatrix> >& smmpsm);
00392 MatrixGeneral<Treal, Tmatrix>& operator=
00393 (const XYZ<Treal,
00394 MatrixGeneral<Treal, Tmatrix>,
00395 MatrixSymmetric<Treal, Tmatrix> >& smm);
00397 MatrixGeneral<Treal, Tmatrix>& operator=
00398 (const XY<MatrixGeneral<Treal, Tmatrix>,
00399 MatrixSymmetric<Treal, Tmatrix> >& mm);
00401 MatrixGeneral<Treal, Tmatrix>& operator+=
00402 (const XYZ<Treal,
00403 MatrixGeneral<Treal, Tmatrix>,
00404 MatrixSymmetric<Treal, Tmatrix> >& smm);
00406 MatrixGeneral<Treal, Tmatrix>& operator=
00407 (const XYZpUV<Treal,
00408 MatrixGeneral<Treal, Tmatrix>,
00409 MatrixSymmetric<Treal, Tmatrix>,
00410 Treal,
00411 MatrixGeneral<Treal, Tmatrix> >& smmpsm);
00413 MatrixGeneral<Treal, Tmatrix>& operator=
00414 (const XYZ<Treal,
00415 MatrixSymmetric<Treal, Tmatrix>,
00416 MatrixSymmetric<Treal, Tmatrix> >& smm);
00418 MatrixGeneral<Treal, Tmatrix>& operator=
00419 (const XY<MatrixSymmetric<Treal, Tmatrix>,
00420 MatrixSymmetric<Treal, Tmatrix> >& mm);
00422 MatrixGeneral<Treal, Tmatrix>& operator+=
00423 (const XYZ<Treal,
00424 MatrixSymmetric<Treal, Tmatrix>,
00425 MatrixSymmetric<Treal, Tmatrix> >& smm);
00427 MatrixGeneral<Treal, Tmatrix>& operator=
00428 (const XYZpUV<Treal,
00429 MatrixSymmetric<Treal, Tmatrix>,
00430 MatrixSymmetric<Treal, Tmatrix>,
00431 Treal,
00432 MatrixGeneral<Treal, Tmatrix> >& smmpsm);
00433
00434
00436 MatrixGeneral<Treal, Tmatrix>& operator=
00437 (const XYZ<Treal,
00438 MatrixTriangular<Treal, Tmatrix>,
00439 MatrixGeneral<Treal, Tmatrix> >& smm);
00441 MatrixGeneral<Treal, Tmatrix>& operator=
00442 (const XYZ<Treal,
00443 MatrixGeneral<Treal, Tmatrix>,
00444 MatrixTriangular<Treal, Tmatrix> >& smm);
00445
00446 void random() {
00447 this->matrixPtr->random();
00448 }
00449
00450 void randomZeroStructure(Treal probabilityBeingZero) {
00451 this->matrixPtr->randomZeroStructure(probabilityBeingZero);
00452 }
00453
00454 template<typename TRule>
00455 void setElementsByRule(TRule & rule) {
00456 this->matrixPtr->setElementsByRule(rule);
00457 return;
00458 }
00459
00460 std::string obj_type_id() const {return "MatrixGeneral";}
00461 protected:
00462 inline void writeToFileProt(std::ofstream & file) const {
00463 this->writeToFileBase(file, matrix_matr);
00464 }
00465 inline void readFromFileProt(std::ifstream & file) {
00466 this->readFromFileBase(file, matrix_matr);
00467 }
00468 private:
00469
00470 };
00471
00472 template<typename Treal, typename Tmatrix>
00473 Treal MatrixGeneral<Treal, Tmatrix>::
00474 eucl(Treal const requestedAccuracy,
00475 int maxIter) const {
00476 VectorType guess;
00477 SizesAndBlocks cols;
00478 this->getCols(cols);
00479 guess.resetSizesAndBlocks(cols);
00480 guess.rand();
00481 mat::ATAMatrix<MatrixGeneral<Treal, Tmatrix>, Treal> ata(*this);
00482 if (maxIter < 0)
00483 maxIter = this->get_nrows() * 100;
00484 arn::LanczosLargestMagnitudeEig
00485 <Treal, ATAMatrix<MatrixGeneral<Treal, Tmatrix>, Treal>, VectorType>
00486 lan(ata, guess, maxIter);
00487 lan.setRelTol( 100 * std::numeric_limits<Treal>::epsilon() );
00488 lan.run();
00489 Treal eVal = 0;
00490 Treal acc = 0;
00491 lan.getLargestMagnitudeEig(eVal, acc);
00492 Interval<Treal> euclInt( sqrt(eVal-acc),
00493 sqrt(eVal+acc) );
00494 if ( euclInt.low() < 0 )
00495 euclInt = Interval<Treal>( 0, sqrt(eVal+acc) );
00496 if ( euclInt.length() / 2.0 > requestedAccuracy ) {
00497 std::cout << "req: " << requestedAccuracy
00498 << " obt: " << euclInt.length() / 2.0 << std::endl;
00499 throw std::runtime_error("Desired accuracy not obtained in Lanczos.");
00500 }
00501 return euclInt.midPoint();
00502 }
00503
00504
00505 template<typename Treal, typename Tmatrix>
00506 void MatrixGeneral<Treal, Tmatrix>::
00507 thresh(Treal const threshold,
00508 normType const norm) {
00509 switch (norm) {
00510 case frobNorm:
00511 this->frob_thresh(threshold);
00512 break;
00513 default:
00514 throw Failure("MatrixGeneral<Treal, Tmatrix>::"
00515 "thresh(Treal const, "
00516 "normType const): "
00517 "Thresholding not imlpemented for selected norm");
00518 }
00519 }
00520
00521 template<typename Treal, typename Tmatrix>
00522 Treal MatrixGeneral<Treal, Tmatrix>::
00523 eucl_thresh(Treal const threshold) {
00524 EuclTruncationGeneral<MatrixGeneral<Treal, Tmatrix>, Treal> TruncObj( *this );
00525 return TruncObj.run( threshold );
00526 }
00527
00528
00529
00530
00531 template<typename Treal, typename Tmatrix>
00532 MatrixGeneral<Treal, Tmatrix>&
00533 MatrixGeneral<Treal, Tmatrix>::operator=
00534 (const XYZ<Treal,
00535 MatrixGeneral<Treal, Tmatrix>,
00536 MatrixGeneral<Treal, Tmatrix> >& smm) {
00537 assert(this != &smm.B && this != &smm.C);
00538 this->matrixPtr.haveDataStructureSet(true);
00539 Tmatrix::gemm(smm.tB, smm.tC, smm.A,
00540 *smm.B.matrixPtr, *smm.C.matrixPtr,
00541 0, *this->matrixPtr);
00542 return *this;
00543 }
00544
00545
00546 template<typename Treal, typename Tmatrix>
00547 MatrixGeneral<Treal, Tmatrix>&
00548 MatrixGeneral<Treal, Tmatrix>::operator=
00549 (const XY<MatrixGeneral<Treal, Tmatrix>,
00550 MatrixGeneral<Treal, Tmatrix> >& mm) {
00551 assert(this != &mm.A && this != &mm.B);
00552 Tmatrix::gemm(mm.tA, mm.tB, 1.0,
00553 *mm.A.matrixPtr, *mm.B.matrixPtr,
00554 0, *this->matrixPtr);
00555 return *this;
00556 }
00557
00558
00559 template<typename Treal, typename Tmatrix>
00560 MatrixGeneral<Treal, Tmatrix>&
00561 MatrixGeneral<Treal, Tmatrix>::operator+=
00562 (const XYZ<Treal,
00563 MatrixGeneral<Treal, Tmatrix>,
00564 MatrixGeneral<Treal, Tmatrix> >& smm) {
00565 assert(this != &smm.B && this != &smm.C);
00566 Tmatrix::gemm(smm.tB, smm.tC, smm.A,
00567 *smm.B.matrixPtr, *smm.C.matrixPtr,
00568 1, *this->matrixPtr);
00569 return *this;
00570 }
00571
00572
00573 template<typename Treal, typename Tmatrix>
00574 MatrixGeneral<Treal, Tmatrix>&
00575 MatrixGeneral<Treal, Tmatrix>::operator=
00576 (const XYZpUV<Treal,
00577 MatrixGeneral<Treal, Tmatrix>,
00578 MatrixGeneral<Treal, Tmatrix>,
00579 Treal,
00580 MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
00581 assert(this != &smmpsm.B && this != &smmpsm.C);
00582 assert(!smmpsm.tE);
00583 if (this == &smmpsm.E)
00584 Tmatrix::gemm(smmpsm.tB, smmpsm.tC, smmpsm.A,
00585 *smmpsm.B.matrixPtr, *smmpsm.C.matrixPtr,
00586 smmpsm.D, *this->matrixPtr);
00587 else
00588 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
00589 "(const XYZpUV<Treal, MatrixGeneral"
00590 "<Treal, Tmatrix> >&) : D = alpha "
00591 "* op(A) * op(B) + beta * C not supported for C != D");
00592 return *this;
00593 }
00594
00595
00596
00597 template<typename Treal, typename Tmatrix>
00598 inline MatrixGeneral<Treal, Tmatrix>&
00599 MatrixGeneral<Treal, Tmatrix>::operator=
00600 (const XpY<MatrixGeneral<Treal, Tmatrix>,
00601 MatrixGeneral<Treal, Tmatrix> >& mpm) {
00602 assert(this != &mpm.A);
00603 (*this) = mpm.B;
00604 Tmatrix::add(1.0, *mpm.A.matrixPtr, *this->matrixPtr);
00605 return *this;
00606 }
00607
00608 template<typename Treal, typename Tmatrix>
00609 inline MatrixGeneral<Treal, Tmatrix>&
00610 MatrixGeneral<Treal, Tmatrix>::operator+=
00611 (MatrixGeneral<Treal, Tmatrix> const & A) {
00612 Tmatrix::add(1.0, *A.matrixPtr, *this->matrixPtr);
00613 return *this;
00614 }
00615
00616 template<typename Treal, typename Tmatrix>
00617 inline MatrixGeneral<Treal, Tmatrix>&
00618 MatrixGeneral<Treal, Tmatrix>::operator-=
00619 (MatrixGeneral<Treal, Tmatrix> const & A) {
00620 Tmatrix::add(-1.0, *A.matrixPtr, *this->matrixPtr);
00621 return *this;
00622 }
00623
00624 template<typename Treal, typename Tmatrix>
00625 inline MatrixGeneral<Treal, Tmatrix>&
00626 MatrixGeneral<Treal, Tmatrix>::operator+=
00627 (XY<Treal, MatrixGeneral<Treal, Tmatrix> > const & sm) {
00628 assert(!sm.tB);
00629 Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
00630 return *this;
00631 }
00632
00633
00634
00635
00636 template<typename Treal, typename Tmatrix>
00637 MatrixGeneral<Treal, Tmatrix>&
00638 MatrixGeneral<Treal, Tmatrix>::operator=
00639 (const XYZ<Treal,
00640 MatrixSymmetric<Treal, Tmatrix>,
00641 MatrixGeneral<Treal, Tmatrix> >& smm) {
00642 assert(this != &smm.C);
00643 assert(!smm.tB && !smm.tC);
00644 this->matrixPtr.haveDataStructureSet(true);
00645 Tmatrix::symm('L', 'U', smm.A,
00646 *smm.B.matrixPtr, *smm.C.matrixPtr,
00647 0, *this->matrixPtr);
00648 return *this;
00649 }
00650
00651
00652 template<typename Treal, typename Tmatrix>
00653 MatrixGeneral<Treal, Tmatrix>&
00654 MatrixGeneral<Treal, Tmatrix>::operator=
00655 (const XY<MatrixSymmetric<Treal, Tmatrix>,
00656 MatrixGeneral<Treal, Tmatrix> >& mm) {
00657 assert(this != &mm.B);
00658 assert(!mm.tB);
00659 Tmatrix::symm('L', 'U', 1.0,
00660 *mm.A.matrixPtr, *mm.B.matrixPtr,
00661 0, *this->matrixPtr);
00662 return *this;
00663 }
00664
00665
00666 template<typename Treal, typename Tmatrix>
00667 MatrixGeneral<Treal, Tmatrix>&
00668 MatrixGeneral<Treal, Tmatrix>::operator+=
00669 (const XYZ<Treal,
00670 MatrixSymmetric<Treal, Tmatrix>,
00671 MatrixGeneral<Treal, Tmatrix> >& smm) {
00672 assert(this != &smm.C);
00673 assert(!smm.tB && !smm.tC);
00674 Tmatrix::symm('L', 'U', smm.A,
00675 *smm.B.matrixPtr, *smm.C.matrixPtr,
00676 1, *this->matrixPtr);
00677 return *this;
00678 }
00679
00680
00681 template<typename Treal, typename Tmatrix>
00682 MatrixGeneral<Treal, Tmatrix>&
00683 MatrixGeneral<Treal, Tmatrix>::operator=
00684 (const XYZpUV<Treal,
00685 MatrixSymmetric<Treal, Tmatrix>,
00686 MatrixGeneral<Treal, Tmatrix>,
00687 Treal,
00688 MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
00689 assert(this != &smmpsm.C);
00690 assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
00691 if (this == &smmpsm.E)
00692 Tmatrix::symm('L', 'U', smmpsm.A,
00693 *smmpsm.B.matrixPtr, *smmpsm.C.matrixPtr,
00694 smmpsm.D, *this->matrixPtr);
00695 else
00696 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
00697 "(const XYZpUV<Treal, MatrixGeneral"
00698 "<Treal, Tmatrix>, MatrixSymmetric<Treal, "
00699 "Tmatrix>, Treal, MatrixGeneral"
00700 "<Treal, Tmatrix> >&) "
00701 ": D = alpha * A * B + beta * C (with A symmetric)"
00702 " not supported for C != D");
00703 return *this;
00704 }
00705
00706
00707 template<typename Treal, typename Tmatrix>
00708 MatrixGeneral<Treal, Tmatrix>&
00709 MatrixGeneral<Treal, Tmatrix>::operator=
00710 (const XYZ<Treal,
00711 MatrixGeneral<Treal, Tmatrix>,
00712 MatrixSymmetric<Treal, Tmatrix> >& smm) {
00713 assert(this != &smm.B);
00714 assert(!smm.tB && !smm.tC);
00715 this->matrixPtr.haveDataStructureSet(true);
00716 Tmatrix::symm('R', 'U', smm.A,
00717 *smm.C.matrixPtr, *smm.B.matrixPtr,
00718 0, *this->matrixPtr);
00719 return *this;
00720 }
00721
00722
00723 template<typename Treal, typename Tmatrix>
00724 MatrixGeneral<Treal, Tmatrix>&
00725 MatrixGeneral<Treal, Tmatrix>::operator=
00726 (const XY<MatrixGeneral<Treal, Tmatrix>,
00727 MatrixSymmetric<Treal, Tmatrix> >& mm) {
00728 assert(this != &mm.A);
00729 assert(!mm.tA && !mm.tB);
00730 Tmatrix::symm('R', 'U', 1.0,
00731 *mm.B.matrixPtr, *mm.A.matrixPtr,
00732 0, *this->matrixPtr);
00733 return *this;
00734 }
00735
00736
00737 template<typename Treal, typename Tmatrix>
00738 MatrixGeneral<Treal, Tmatrix>&
00739 MatrixGeneral<Treal, Tmatrix>::operator+=
00740 (const XYZ<Treal,
00741 MatrixGeneral<Treal, Tmatrix>,
00742 MatrixSymmetric<Treal, Tmatrix> >& smm) {
00743 assert(this != &smm.B);
00744 assert(!smm.tB && !smm.tC);
00745 Tmatrix::symm('R', 'U', smm.A,
00746 *smm.C.matrixPtr, *smm.B.matrixPtr,
00747 1, *this->matrixPtr);
00748 return *this;
00749 }
00750
00751
00752 template<typename Treal, typename Tmatrix>
00753 MatrixGeneral<Treal, Tmatrix>&
00754 MatrixGeneral<Treal, Tmatrix>::operator=
00755 (const XYZpUV<Treal,
00756 MatrixGeneral<Treal, Tmatrix>,
00757 MatrixSymmetric<Treal, Tmatrix>,
00758 Treal,
00759 MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
00760 assert(this != &smmpsm.B);
00761 assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
00762 if (this == &smmpsm.E)
00763 Tmatrix::symm('R', 'U', smmpsm.A,
00764 *smmpsm.C.matrixPtr, *smmpsm.B.matrixPtr,
00765 smmpsm.D, *this->matrixPtr);
00766 else
00767 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
00768 "(const XYZpUV<Treal, MatrixSymmetric"
00769 "<Treal, Tmatrix>, MatrixGeneral<Treal, "
00770 "Tmatrix>, Treal, MatrixGeneral"
00771 "<Treal, Tmatrix> >&) "
00772 ": D = alpha * B * A + beta * C (with A symmetric)"
00773 " not supported for C != D");
00774 return *this;
00775 }
00776
00777
00779 template<typename Treal, typename Tmatrix>
00780 MatrixGeneral<Treal, Tmatrix>&
00781 MatrixGeneral<Treal, Tmatrix>::operator=
00782 (const XYZ<Treal,
00783 MatrixSymmetric<Treal, Tmatrix>,
00784 MatrixSymmetric<Treal, Tmatrix> >& smm) {
00785 assert(!smm.tB && !smm.tC);
00786 this->matrixPtr.haveDataStructureSet(true);
00787 Tmatrix::ssmm(smm.A,
00788 *smm.B.matrixPtr,
00789 *smm.C.matrixPtr,
00790 0,
00791 *this->matrixPtr);
00792 return *this;
00793 }
00794
00796 template<typename Treal, typename Tmatrix>
00797 MatrixGeneral<Treal, Tmatrix>&
00798 MatrixGeneral<Treal, Tmatrix>::operator=
00799 (const XY<MatrixSymmetric<Treal, Tmatrix>,
00800 MatrixSymmetric<Treal, Tmatrix> >& mm) {
00801 assert(!mm.tB);
00802 Tmatrix::ssmm(1.0,
00803 *mm.A.matrixPtr,
00804 *mm.B.matrixPtr,
00805 0,
00806 *this->matrixPtr);
00807 return *this;
00808 }
00809
00811 template<typename Treal, typename Tmatrix>
00812 MatrixGeneral<Treal, Tmatrix>&
00813 MatrixGeneral<Treal, Tmatrix>::operator+=
00814 (const XYZ<Treal,
00815 MatrixSymmetric<Treal, Tmatrix>,
00816 MatrixSymmetric<Treal, Tmatrix> >& smm) {
00817 assert(!smm.tB && !smm.tC);
00818 Tmatrix::ssmm(smm.A,
00819 *smm.B.matrixPtr,
00820 *smm.C.matrixPtr,
00821 1,
00822 *this->matrixPtr);
00823 return *this;
00824 }
00825
00826
00828 template<typename Treal, typename Tmatrix>
00829 MatrixGeneral<Treal, Tmatrix>&
00830 MatrixGeneral<Treal, Tmatrix>::operator=
00831 (const XYZpUV<Treal,
00832 MatrixSymmetric<Treal, Tmatrix>,
00833 MatrixSymmetric<Treal, Tmatrix>,
00834 Treal,
00835 MatrixGeneral<Treal, Tmatrix> >& smmpsm) {
00836 assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
00837 if (this == &smmpsm.E)
00838 Tmatrix::ssmm(smmpsm.A,
00839 *smmpsm.B.matrixPtr,
00840 *smmpsm.C.matrixPtr,
00841 smmpsm.D,
00842 *this->matrixPtr);
00843 else
00844 throw Failure("MatrixGeneral<Treal, Tmatrix>::"
00845 "operator=(const XYZpUV<"
00846 "Treal, MatrixSymmetric<Treal, Tmatrix>, "
00847 "MatrixSymmetric<Treal, Tmatrix>, Treal, "
00848 "MatrixGeneral<Treal, Tmatrix> >&) "
00849 ": D = alpha * A * B + beta * C (with A and B symmetric)"
00850 " not supported for C != D");
00851 return *this;
00852 }
00853
00854
00855
00856
00857
00858
00859 template<typename Treal, typename Tmatrix>
00860 MatrixGeneral<Treal, Tmatrix>&
00861 MatrixGeneral<Treal, Tmatrix>::operator=
00862 (const XYZ<Treal,
00863 MatrixTriangular<Treal, Tmatrix>,
00864 MatrixGeneral<Treal, Tmatrix> >& smm) {
00865 assert(!smm.tC);
00866 if (this == &smm.C)
00867 Tmatrix::trmm('L', 'U', smm.tB, smm.A,
00868 *smm.B.matrixPtr, *this->matrixPtr);
00869 else
00870 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
00871 "(const XYZ<Treal, MatrixTriangular"
00872 "<Treal, Tmatrix>, MatrixGeneral<Treal,"
00873 " Tmatrix> >& : D = alpha * op(A) * B (with"
00874 " A upper triangular) not supported for B != D");
00875 return *this;
00876 }
00877
00878
00879
00880 template<typename Treal, typename Tmatrix>
00881 MatrixGeneral<Treal, Tmatrix>&
00882 MatrixGeneral<Treal, Tmatrix>::operator=
00883 (const XYZ<Treal,
00884 MatrixGeneral<Treal, Tmatrix>,
00885 MatrixTriangular<Treal, Tmatrix> >& smm) {
00886 assert(!smm.tB);
00887 if (this == &smm.B)
00888 Tmatrix::trmm('R', 'U', smm.tC, smm.A,
00889 *smm.C.matrixPtr, *this->matrixPtr);
00890 else
00891 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
00892 "(const XYZ<Treal, MatrixGeneral"
00893 "<Treal, Tmatrix>, MatrixTriangular<Treal,"
00894 " Tmatrix> >& : D = alpha * A * op(B) (with"
00895 " B upper triangular) not supported for A != D");
00896 return *this;
00897 }
00898
00899
00900
00901 template<typename Treal, typename Tmatrix>
00902 Treal trace(const XYZ<Treal,
00903 MatrixGeneral<Treal, Tmatrix>,
00904 MatrixGeneral<Treal, Tmatrix> >& smm) {
00905 if ((!smm.tB && !smm.tC) || (smm.tB && smm.tC))
00906 return smm.A * MatrixGeneral<Treal, Tmatrix>::
00907 trace_ab(smm.B, smm.C);
00908 else if (smm.tB)
00909 return smm.A * MatrixGeneral<Treal, Tmatrix>::
00910 trace_aTb(smm.B, smm.C);
00911 else
00912 return smm.A * MatrixGeneral<Treal, Tmatrix>::
00913 trace_aTb(smm.C, smm.B);
00914 }
00915
00916
00917 }
00918 #endif
00919
00920