linalg.h
Go to the documentation of this file.00001 00013 #ifndef LINALG_H 00014 #define LINALG_H 00015 00016 #include "mycomplex.h" 00017 #include "xpr-vector.h" 00018 #include "xpr-matrix.h" 00019 #include "def.h" 00020 #include <iosfwd> 00021 #include <valarray> 00022 #include <sstream> 00023 00024 using namespace softsusy; 00025 00026 /************************************ 00027 * 00028 * *** DOUBLE VECTOR CLASS *** 00029 * 00030 ************************************/ 00031 00034 class DoubleVector : public Indexable<double, DoubleVector> { 00035 private: 00036 int start, end; 00037 std::valarray<double> x; 00038 00039 DoubleVector(const std::valarray<double> & valarr, int s, int e) 00040 : Indexable<double,DoubleVector>(), start(s), end(e), x(valarr) {} 00041 friend class DoubleMatrix; 00042 00043 00044 public: 00045 00046 /* 00047 * CONSTRUCTORS 00048 */ 00049 00050 explicit DoubleVector(int e) 00051 : Indexable<double,DoubleVector>(), start(1), end(e), x(0.0,e) {} 00052 DoubleVector(int s, int e) 00053 : Indexable<double,DoubleVector>(), start(s), end(e), x(0.0,e-s+1) {} 00054 ~DoubleVector() {} 00055 00057 template<class E> 00058 DoubleVector & operator=(const Xpr<double,E> & x) { 00059 *this = copy_from(x); 00060 return *this; 00061 } 00062 00063 template <typename E> 00064 DoubleVector(const Xpr<double,E> & v) 00065 : Indexable<double,DoubleVector>(), 00066 start(v.displayStart()), end(v.displayEnd()), 00067 x(0.0,v.displayEnd()-v.displayStart()+1) 00068 { 00069 assign_from(v); 00070 } 00071 00072 /* 00073 * ELEMENT ACCESS 00074 */ 00075 00076 double & operator() (int i); 00077 double operator() (int i) const { return x[i-start]; } 00078 double display(int i) const { return x[i-start]; } 00079 void set(int i, double f); 00080 00081 int displayStart() const { return start; } 00082 int displayEnd() const { return end; } 00083 const DoubleVector & display() const { return *this; } 00084 double compare(const DoubleVector & a) const; 00085 00088 void setEnd(int e); 00089 00090 double norm() const; 00092 DoubleVector apply(double (*fn)(double)) const { 00093 return DoubleVector(x.apply(fn),start,end); 00094 } 00095 double nmin(int & p) const; 00096 double max() const { return x.max(); } 00097 00098 double max(int & p) const { 00099 double m = double(0); 00100 int i; 00101 for (i=displayStart(); i<=displayEnd(); i++) 00102 if (display(i) > m) { 00103 m = display(i); 00104 p = i; 00105 } 00106 return m; 00107 } 00108 00109 00110 00111 double min(int & p) const; 00112 00113 double sumElements() const { return abs(x).sum(); } 00114 00115 void swap(int i, int j); 00116 00117 double dot(const DoubleVector & v) const { 00118 return (x * v.x).sum(); 00119 } 00120 00121 DoubleVector sort() const { 00122 DoubleVector temp(*this), saveIt(*this); 00123 00124 int j, pos; 00125 for (j=start; j<=end; j++) { 00126 double d = temp.min(pos); 00127 if (d < 1.0e66) saveIt(j) = d; 00128 00129 temp(pos) = 6.66e66; // take that element out of the min calculation 00130 } 00131 return saveIt; 00132 } 00133 00135 DoubleVector divide(DoubleVector const &v) const; 00137 double average() const; 00138 00139 }; 00140 00141 00142 /* 00143 * STANDARD INPUT / OUTPUT 00144 */ 00145 00147 std::ostream & operator<<(std::ostream &left, const DoubleVector &V); 00149 std::istream & operator>>(std::istream & left, DoubleVector &V); 00150 00151 00152 /* 00153 * INLINE FUNCTION DEFINITIONS 00154 */ 00155 00156 inline double & DoubleVector::operator() (int i) { 00157 #ifdef ARRAY_BOUNDS_CHECKING 00158 if (i>end || i<start) { 00159 std::ostringstream ii; 00160 ii << "Trying to access " << i << "th element of DoubleVector\n"; 00161 ii << "start " << start << " end " << end; 00162 ii << *this; 00163 throw ii.str(); 00164 } 00165 #endif 00166 return x[i-start]; 00167 } 00168 00169 inline void DoubleVector::set(int i, double f) { 00170 #ifdef ARRAY_BOUNDS_CHECKING 00171 if (i < start || i > end) 00172 { 00173 std::ostringstream ii; 00174 ii << "Cannot access " << i << "th element of vector " << *this; 00175 throw ii.str(); 00176 } 00177 #endif 00178 x[i-start] = f; 00179 } 00180 00181 00182 /************************************ 00183 * 00184 * *** DOUBLE MATRIX CLASS *** 00185 * 00186 ************************************/ 00187 00189 class ComplexVector; class ComplexMatrix; 00190 class DoubleMatrix : public MatIndexable<double, DoubleMatrix> { 00191 private: 00192 int rows, cols; 00193 std::valarray<double> x; 00194 00197 // for usage see e.g. http://www.roguewave.com/support/docs/sourcepro/stdlibref/2-10.html 00198 std::slice_array<double> row(int i) { return x[std::slice((i-1)*cols,cols,1)]; } 00199 std::slice_array<double> col(int i) { return x[std::slice(i-1,rows,cols)]; } 00200 const std::valarray<double> row(int i) const { 00201 return std::valarray<double>(x[std::slice((i-1)*cols,cols,1)]); 00202 } 00203 const std::valarray<double> col(int i) const { 00204 return std::valarray<double>(x[std::slice(i-1,rows,cols)]); 00205 } 00206 00207 std::slice_array<double> diag() { return x[std::slice(0,rows,cols+1)]; } 00208 const std::valarray<double> diag() const { 00209 return std::valarray<double>(x[std::slice(0,rows,cols+1)]); 00210 } 00211 00213 double & elmt(int i, int j) { return x[(i-1) * cols + (j-1)]; } 00214 double elmt(int i, int j) const { return x[(i-1) * cols + (j-1)]; } 00215 00217 DoubleMatrix(const std::valarray<double> & valarr, int r, int c) 00218 : MatIndexable<double,DoubleMatrix>(), rows(r), cols(c), x(valarr) {} 00219 00220 friend class ComplexMatrix; 00221 00222 public: 00223 00224 /* 00225 * CONSTRUCTORS 00226 */ 00227 00229 DoubleMatrix(int r, int c) 00230 : MatIndexable<double,DoubleMatrix>(), rows(r), cols(c), x(0.0, r*c) {} 00232 explicit DoubleMatrix(const DoubleVector &v); 00233 ~DoubleMatrix() {} 00234 00235 00236 00238 template<class E> 00239 DoubleMatrix & operator=(const MatXpr<double,E> & x) { 00240 *this = copy_from(x); 00241 return *this; 00242 } 00243 00244 template <typename E> 00245 DoubleMatrix(const MatXpr<double, E> & m) 00246 : MatIndexable<double,DoubleMatrix>(), 00247 rows(m.displayRows()), cols(m.displayCols()), x(0.0, rows*cols) 00248 { 00249 assign_from(m); 00250 } 00251 00252 /* 00253 * ELEMENT ACCESS 00254 */ 00255 double operator() (int i, int j) const; 00256 double & operator() (int i, int j); 00257 double display(int i, int j) const; 00258 double sumElements() const { return abs(x).sum(); } 00260 int displayRows() const { return rows; }; 00261 int displayCols() const { return cols; }; 00262 const DoubleMatrix & display() const { return *this; } 00263 00264 DoubleVector displayRow(int i) const { 00265 DoubleVector temp(displayCols()); 00266 int j; for (j=1; j<=displayCols(); j++) temp(j) = display(i, j); 00267 return temp; 00268 } 00269 00270 DoubleVector displayCol(int i) const { 00271 DoubleVector temp(displayRows()); 00272 int j; for (j=1; j<=displayRows(); j++) temp(j) = display(j, i); 00273 return temp; 00274 } 00275 00276 00277 void operator+=(DoubleMatrix & f) { 00278 x += f.x; 00279 } 00280 00281 double nmin(int & k, int & l) const; 00282 00284 DoubleMatrix apply(double (*fn)(double)) const { 00285 return DoubleMatrix(x.apply(fn),rows,cols); 00286 } 00287 00289 const DoubleMatrix & operator=(double v); 00290 00291 double min(int & k, int & l) const; 00292 double max(int & k, int & l) const; 00293 00294 void swaprows(int i, int j); 00295 void swapcols(int i,int j); 00296 00297 double trace() const; 00298 DoubleMatrix transpose() const; 00299 00300 /* 00301 * NUMERICAL DIAGONALIZATION ROUTINES ETC. 00302 */ 00303 00307 void associateOrderAbs(DoubleVector &v); 00308 // Perform on asymmetric matrices M such that 00311 void associateOrderAbs(DoubleMatrix & u, DoubleMatrix & v, DoubleVector &w) 00312 const; 00314 void symmetrise(); 00317 double compare(const DoubleMatrix & a) const; 00318 double compare(const ComplexMatrix & a) const; 00320 DoubleVector diagVals() const; 00321 DoubleMatrix inverse() const; 00322 00326 double diagonalise(DoubleMatrix & u, DoubleMatrix & v, DoubleVector & w) const; 00330 double diagonaliseSym(DoubleMatrix & v, DoubleVector & w) const; 00331 double diagonaliseSym(ComplexMatrix & v, DoubleVector & w) const; 00332 00338 DoubleVector sym2by2(double & theta) const; 00342 DoubleVector asy2by2(double & thetaL, double & thetaR) const; 00344 DoubleVector vectorfy() const; 00345 }; 00346 00347 /* 00348 * STANDARD INPUT / OUTPUT 00349 */ 00350 00351 std::istream & operator>>(std::istream & left, DoubleMatrix &M); 00352 std::ostream & operator <<(std::ostream &left, const DoubleMatrix &V); 00353 00354 /* 00355 * NON-MEMBER DIAGONALIZATION ROUTINES ETC. 00356 */ 00357 00359 // [ cos theta sin theta ] 00360 // [-sin theta cos theta ] 00361 DoubleMatrix rot2d(double theta); 00363 // [ -sin(theta) cos(theta) ] 00364 // [ cos(theta) sin(theta) ] -- 00365 DoubleMatrix rot2dTwist(double theta); 00368 // [ cos thetaL sin thetaL ] A [ cos thetaR -sin thetaR ] = diag 00369 // [ -sin thetaL cos thetaL ] [ sin thetaR cos thetaR ] 00370 // as given by asy2by2! 00372 void positivise(double thetaL, double thetaR, const DoubleVector & diag, 00373 ComplexMatrix & u, ComplexMatrix & v); 00375 void diagonaliseSvd(DoubleMatrix & a, DoubleVector & w, DoubleMatrix & v); 00376 double pythagoras(double a, double b); 00377 void diagonaliseJac(DoubleMatrix & a, int n, DoubleVector & d, DoubleMatrix 00378 & v, int *nrot); 00379 00380 00381 /* 00382 * INLINE FUNCTION DEFINITIONS 00383 */ 00384 00385 inline double DoubleMatrix::operator() (int i, int j) const { 00386 #ifdef ARRAY_BOUNDS_CHECKING 00387 if (i > rows || j > cols || i < 1 || j < 1) { 00388 std::ostringstream ii; 00389 ii << "Trying to access " << i << "," << j << 00390 "th element of DoubleMatrix\n" << *this; 00391 throw ii.str(); 00392 } 00393 #endif 00394 return elmt(i,j); 00395 } 00396 00397 inline double & DoubleMatrix::operator() (int i, int j) { 00398 #ifdef ARRAY_BOUNDS_CHECKING 00399 if (i > rows || j > cols || i < 1 || j < 1) { 00400 std::ostringstream ii; 00401 ii << "Trying to access " << i << "," << j << 00402 "th element of DoubleMatrix\n" << *this; 00403 throw ii.str(); 00404 } 00405 #endif 00406 return elmt(i,j); 00407 } 00408 00409 inline double DoubleMatrix::display(int i, int j) const { 00410 #ifdef ARRAY_BOUNDS_CHECKING 00411 if (i < 1 || i > displayRows() || j < 1 || j > displayCols()) { 00412 std::ostringstream ii; 00413 ii << "Error: Requested display (" << i << "," << j << 00414 ")th element of matrix " << *this; 00415 throw ii.str(); 00416 } 00417 #endif 00418 return elmt(i,j); 00419 } 00420 00421 /************************************ 00422 * 00423 * *** COMPLEX VECTOR CLASS *** 00424 * 00425 ************************************/ 00426 00428 class ComplexVector : public Indexable<Complex, ComplexVector> { 00429 private: 00430 int start, end; 00431 std::valarray<Complex> x; 00432 00433 ComplexVector(const std::valarray<Complex> & valarr, int s, int e) 00434 : Indexable<Complex, ComplexVector>(), start(s), end(e), x(valarr) {} 00435 friend class ComplexMatrix; 00436 00437 00438 public: 00439 00440 /* 00441 * CONSTRUCTORS 00442 */ 00443 00444 explicit ComplexVector(int e) 00445 : Indexable<Complex, ComplexVector>(), start(1), end(e), x(0.0,e) {} 00446 ComplexVector(int s, int e) 00447 : Indexable<Complex, ComplexVector>(), start(s), end(e), x(0.0,e-s+1) {} 00448 ~ComplexVector() {} 00449 00451 template<class E> 00452 ComplexVector & operator=(const Xpr<Complex,E> & x) { 00453 *this = copy_from(x); 00454 return *this; 00455 } 00456 00457 template <typename E> 00458 ComplexVector(const Xpr<Complex,E> & v) 00459 : Indexable<Complex,ComplexVector>(), 00460 start(v.displayStart()), end(v.displayEnd()), 00461 x(0.0,v.displayEnd()-v.displayStart()+1) 00462 { 00463 assign_from(v); 00464 } 00465 00466 00467 /* 00468 * ELEMENT ACCESS 00469 */ 00470 00471 Complex & operator() (int i); 00472 Complex operator() (int i) const { return x[i-start]; } 00473 Complex display(int i) const { return x[i-start]; } 00474 void set(int i, Complex f); 00475 00476 int displayStart() const { return start; } 00477 int displayEnd() const { return end; } 00478 const ComplexVector & display() const { return *this; } 00479 00482 void setEnd(int e); 00483 00485 ComplexVector apply(Complex (*fn)(Complex)) const { 00486 return ComplexVector(x.apply(fn),start,end); 00487 } 00488 00491 Complex max() const { return x.max(); } 00492 Complex min(int & p) const; 00493 00494 void swap(int i, int j); 00495 00496 }; 00497 00498 /* 00499 * STANDARD INPUT / OUTPUT 00500 */ 00501 00503 std::ostream & operator<<(std::ostream &left, const ComplexVector &V); 00505 std::istream & operator>>(std::istream & left, ComplexVector &V); 00506 00507 /* 00508 * INLINE FUNCTION DEFINITIONS 00509 */ 00510 00511 inline Complex & ComplexVector::operator() (int i) { 00512 #ifdef ARRAY_BOUNDS_CHECKING 00513 if (i>end || i<start) { 00514 std::ostringstream ii; 00515 ii << "Trying to access " << i << "th element of DoubleVector\n"; 00516 ii << "start " << start << " end " << end; 00517 ii << *this; 00518 throw ii.str(); 00519 } 00520 #endif 00521 return x[i-start]; 00522 } 00523 00524 inline void ComplexVector::set(int i, Complex f) { 00525 #ifdef ARRAY_BOUNDS_CHECKING 00526 if (i < start || i > end) { 00527 std::ostringstream ii; 00528 ii << "Cannot access " << i << "th element of vector " << *this; 00529 throw ii.str(); 00530 } 00531 #endif 00532 x[i-start] = f; 00533 } 00534 00535 /************************************ 00536 * 00537 * *** COMPLEX MATRIX CLASS *** 00538 * 00539 ************************************/ 00540 00543 class ComplexMatrix : public MatIndexable<Complex, ComplexMatrix> { 00544 private: 00545 int rows, cols; 00546 std::valarray<Complex> x; 00547 00550 // for usage see e.g. http://www.roguewave.com/support/docs/sourcepro/stdlibref/2-10.html 00551 std::slice_array<Complex> row(int i) { return x[std::slice((i-1)*cols,cols,1)]; } 00552 std::slice_array<Complex> col(int i) { return x[std::slice(i-1,rows,cols)]; } 00553 const std::valarray<Complex> row(int i) const { 00554 return std::valarray<Complex>(x[std::slice((i-1)*cols,cols,1)]); 00555 } 00556 const std::valarray<Complex> col(int i) const { 00557 return std::valarray<Complex>(x[std::slice(i-1,rows,cols)]); 00558 } 00559 00560 std::slice_array<Complex> diag() { return x[std::slice(0,rows,cols+1)]; } 00561 const std::valarray<Complex> diag() const { 00562 return std::valarray<Complex>(x[std::slice(0,rows,cols+1)]); 00563 } 00564 00566 Complex & elmt(int i, int j) { return x[(i-1)*cols+(j-1)]; } 00567 const Complex elmt(int i, int j) const { return x[(i-1)*cols+(j-1)]; } 00568 00570 ComplexMatrix(const std::valarray<Complex> & valarr, int r, int c) 00571 : MatIndexable<Complex,ComplexMatrix>(), rows(r), cols(c), x(valarr) {} 00572 00573 friend class DoubleMatrix; 00574 00575 public: 00576 00577 /* 00578 * CONSTRUCTORS 00579 */ 00580 00582 ComplexMatrix(int r, int c) 00583 : MatIndexable<Complex,ComplexMatrix>(), rows(r), cols(c), x(0.0, r*c) {} 00585 explicit ComplexMatrix(const DoubleMatrix & m); 00587 explicit ComplexMatrix(const ComplexVector &v); 00588 ~ComplexMatrix() {} 00589 00591 template<class E> 00592 ComplexMatrix & operator=(const MatXpr<Complex,E> & x) { 00593 *this = copy_from(x); 00594 return *this; 00595 } 00596 00597 template <typename E> 00598 ComplexMatrix(const MatXpr<Complex, E> & m) 00599 : MatIndexable<Complex,ComplexMatrix>(), 00600 rows(m.displayRows()), cols(m.displayCols()), x(0.0, rows*cols) 00601 { 00602 assign_from(m); 00603 } 00604 00605 /* 00606 * ELEMENT ACCESS 00607 */ 00608 00609 Complex & operator() (int i, int j); 00610 Complex operator() (int i, int j) const; 00611 Complex display(int i, int j) const; 00612 00613 int displayRows() const { return rows; }; 00614 int displayCols() const { return cols; }; 00615 const ComplexMatrix & display() const { return *this; }; 00616 00618 const ComplexMatrix & operator=(const Complex &v); 00619 00620 Complex min(int & k, int & l) const; 00621 00622 void swaprows(int i, int j); 00623 void swapcols(int i,int j); 00624 00625 Complex trace() const; 00626 ComplexMatrix transpose() const; 00627 ComplexMatrix hermitianConjugate() const; 00628 ComplexMatrix complexConjugate() const { 00629 return ComplexMatrix(x.apply(conj),rows,cols); 00630 } 00631 00632 /* 00633 * NUMERICAL DIAGONALIZATION ROUTINES ETC. 00634 */ 00635 00637 void symmetrise(); 00639 double compare(const ComplexMatrix & a) const; 00640 }; 00641 00642 /* 00643 * STANDARD INPUT / OUTPUT 00644 */ 00645 00647 std::istream & operator>>(std::istream & left, ComplexMatrix &M); 00649 std::ostream & operator<<(std::ostream &left, const ComplexMatrix &V); 00650 00651 /* 00652 * INLINE FUNCTION DEFINITIONS 00653 */ 00654 00655 inline Complex & ComplexMatrix::operator() (int i, int j) { 00656 #ifdef ARRAY_BOUNDS_CHECKING 00657 if (i > rows || j > cols || i < 0 || j < 0) { 00658 std::ostringstream ii; 00659 ii << "Trying to access " << i << "," << j << 00660 "th element of ComplexMatrix\n" << *this; 00661 throw ii.str(); 00662 } 00663 #endif 00664 return elmt(i,j); 00665 } 00666 00667 inline Complex ComplexMatrix::operator() (int i, int j) const { 00668 #ifdef ARRAY_BOUNDS_CHECKING 00669 if (i > rows || j > cols || i < 0 || j < 0) { 00670 std::ostringstream ii; 00671 ii << "Trying to access " << i << "," << j << 00672 "th element of ComplexMatrix\n" << *this; 00673 throw ii.str(); 00674 } 00675 #endif 00676 return elmt(i,j); 00677 } 00678 00679 inline Complex ComplexMatrix::display(int i, int j) const { 00680 #ifdef ARRAY_BOUNDS_CHECKING 00681 if (i < 1 || i > displayRows() || j < 1 || j > displayCols()) { 00682 std::ostringstream ii; 00683 ii << "Error: Requested display (" << i << "," << j << 00684 ")th element of matrix " << *this; 00685 throw ii.str(); 00686 } 00687 #endif 00688 return elmt(i,j); 00689 } 00690 00691 #endif
Generated on Mon Jul 19 17:11:53 2010 for SOFTSUSY by
1.5.5