SoftSUSY is hosted by Hepforge, IPPP Durham
SOFTSUSY  4.1
linalg.h
Go to the documentation of this file.
1 
12 #ifndef LINALG_H
13 #define LINALG_H
14 
15 #include <vector>
16 using namespace std;
17 
18 #include "mycomplex.h"
19 #include "xpr-vector.h"
20 #include "xpr-matrix.h"
21 #include "def.h"
22 #include "utils.h"
23 #include <iosfwd>
24 #include <valarray>
25 #include <sstream>
26 using namespace softsusy;
27 
28 /************************************
29  *
30  * *** DOUBLE VECTOR CLASS ***
31  *
32  ************************************/
33 
35 class DoubleVector : public Indexable<double, DoubleVector> {
36 private:
37  int start, end;
38  std::valarray<double> x;
39  DoubleVector(const std::valarray<double> & valarr, int s, int e)
41  : Indexable<double,DoubleVector>(), start(s), end(e), x(valarr) {}
42  friend class DoubleMatrix;
43 
44 public:
45 
46  /*
47  * CONSTRUCTORS
48  */
49 
50 
51  explicit DoubleVector(int e)
52  : Indexable<double,DoubleVector>(), start(1), end(e), x(0.0,e) {}
53  DoubleVector(int s, int e)
54  : Indexable<double,DoubleVector>(), start(s), end(e), x(0.0,e-s+1) {}
55  DoubleVector(double a, double b)
56  : Indexable<double,DoubleVector>(), start(1), end(2), x(0.0,2) {
57  x[0] = a; x[1] = b; }
58  ~DoubleVector() {}
59 
61  template<class E>
63  *this = copy_from(x);
64  return *this;
65  }
66 
67  DoubleVector& operator=(const DoubleVector& other) {
68  x.resize(other.x.size());
69  x = other.x;
70  start = other.start;
71  end = other.end;
72  return *this;
73  }
74 
75  template <typename E>
76  DoubleVector(const Xpr<double,E> & v)
78  start(v.displayStart()), end(v.displayEnd()),
79  x(0.0,v.displayEnd()-v.displayStart()+1)
80  {
81  assign_from(v);
82  }
83 
84  /*
85  * ELEMENT ACCESS
86  */
87 
88  double & operator() (int i);
89  double operator() (int i) const { return x[i-start]; }
90  double display(int i) const { return x[i-start]; }
91  void set(int i, double f);
92 
93  int displayStart() const { return start; }
94  int displayEnd() const { return end; }
95  const DoubleVector & display() const { return *this; }
96  double compare(const DoubleVector & a) const;
97 
99  void append(double);
100 
103  void setEnd(int e);
105  DoubleVector abs() const { return DoubleVector(std::abs(x), start, end); }
106 
107  double norm() const;
109  DoubleVector apply(double (*fn)(double)) const {
110  return DoubleVector(x.apply(fn),start,end);
111  }
112 
113  double nmin(int & p) const;
114  double max() const { return x.max(); }
115 
116  double max(int & p) const {
117  double m = double(0);
118  int i;
119  for (i=displayStart(); i<=displayEnd(); i++)
120  if (display(i) > m) {
121  m = display(i);
122  p = i;
123  }
124  return m;
125  }
126 
127  double min(int & p) const;
128 
130  double absmax(int & p) const { return abs().max(p); }
132  double sumElements() const { return std::abs(x).sum(); }
133  std::size_t size() const { return x.size(); }
134 
135  void swap(int i, int j);
136 
137  double dot(const DoubleVector & v) const {
138  return (x * v.x).sum();
139  }
140 
141  DoubleVector sort() const {
142  DoubleVector temp(*this), saveIt(*this);
143 
144  int j, pos;
145  for (j=start; j<=end; j++) {
146  double d = temp.min(pos);
147  if (d < 1.0e66) saveIt(j) = d;
148 
149  temp(pos) = 6.66e66; // take that element out of the min calculation
150  }
151  return saveIt;
152  }
153 
155  DoubleVector divide(DoubleVector const &v) const;
157  double average() const;
158  int closest(double a) const;
161  void fillArray(double* array, unsigned offset = 0) const;
162 };
163 
164 
165 /*
166  * STANDARD INPUT / OUTPUT
167  */
168 
170 std::ostream & operator<<(std::ostream &left, const DoubleVector &V);
172 std::istream & operator>>(std::istream & left, DoubleVector &V);
173 
174 
175 /*
176  * INLINE FUNCTION DEFINITIONS
177  */
178 
179 inline double & DoubleVector::operator() (int i) {
180 #ifdef ARRAY_BOUNDS_CHECKING
181  if (i>end || i<start) {
182  std::ostringstream ii;
183  ii << "Trying to access " << i << "th element of DoubleVector\n";
184  ii << "start " << start << " end " << end;
185  ii << *this;
186  throw ii.str();
187  }
188 #endif
189  return x[i-start];
190 }
191 
192 inline void DoubleVector::set(int i, double f) {
193 #ifdef ARRAY_BOUNDS_CHECKING
194  if (i < start || i > end)
195  {
196  std::ostringstream ii;
197  ii << "Cannot access " << i << "th element of vector " << *this;
198  throw ii.str();
199  }
200 #endif
201  x[i-start] = f;
202 }
203 
204 
205 /************************************
206  *
207  * *** DOUBLE MATRIX CLASS ***
208  *
209  ************************************/
210 
211 
212 class ComplexVector; class ComplexMatrix;
214 class DoubleMatrix : public MatIndexable<double, DoubleMatrix> {
215 private:
216  int rows, cols;
217  std::valarray<double> x;
218 
221  // for usage see e.g. http://www.roguewave.com/support/docs/sourcepro/stdlibref/2-10.html
222  std::slice_array<double> row(int i) { return x[std::slice((i-1)*cols,cols,1)]; }
223  std::slice_array<double> col(int i) { return x[std::slice(i-1,rows,cols)]; }
224  const std::valarray<double> row(int i) const {
225  return std::valarray<double>(x[std::slice((i-1)*cols,cols,1)]);
226  }
227  const std::valarray<double> col(int i) const {
228  return std::valarray<double>(x[std::slice(i-1,rows,cols)]);
229  }
230 
231  std::slice_array<double> diag() { return x[std::slice(0,rows,cols+1)]; }
232  const std::valarray<double> diag() const {
233  return std::valarray<double>(x[std::slice(0,rows,cols+1)]);
234  }
235 
237  double & elmt(int i, int j) { return x[(i-1) * cols + (j-1)]; }
238  double elmt(int i, int j) const { return x[(i-1) * cols + (j-1)]; }
239 
241  DoubleMatrix(const std::valarray<double> & valarr, int r, int c)
242  : MatIndexable<double,DoubleMatrix>(), rows(r), cols(c), x(valarr) {}
243 
244  friend class ComplexMatrix;
245 
246 public:
247 
248  /*
249  * CONSTRUCTORS
250  */
251 
253  DoubleMatrix(int r, int c)
254  : MatIndexable<double,DoubleMatrix>(), rows(r), cols(c), x(0.0, r*c) {}
256  explicit DoubleMatrix(const DoubleVector &v);
257  ~DoubleMatrix() {}
258 
260  template<class E>
262  *this = copy_from(x);
263  return *this;
264  }
265 
266  DoubleMatrix& operator=(const DoubleMatrix& other) {
267  x.resize(other.x.size());
268  x = other.x;
269  rows = other.rows;
270  cols = other.cols;
271  return *this;
272  }
273 
274  template <typename E>
277  rows(m.displayRows()), cols(m.displayCols()), x(0.0, rows*cols)
278  {
279  assign_from(m);
280  }
281 
282  /*
283  * ELEMENT ACCESS
284  */
285  double operator() (int i, int j) const;
286  double & operator() (int i, int j);
287  double display(int i, int j) const;
288  double sumElements() const { return abs(x).sum(); }
290  int displayRows() const { return rows; };
291  int displayCols() const { return cols; };
292  const DoubleMatrix & display() const { return *this; }
293 
294  DoubleVector displayRow(int i) const {
295  DoubleVector temp(displayCols());
296  int j; for (j=1; j<=displayCols(); j++) temp(j) = display(i, j);
297  return temp;
298  }
299 
300  DoubleVector displayCol(int i) const {
301  DoubleVector temp(displayRows());
302  int j; for (j=1; j<=displayRows(); j++) temp(j) = display(j, i);
303  return temp;
304  }
305 
306 
307  void operator+=(DoubleMatrix & f) {
308  x += f.x;
309  }
310 
311  double nmin(int & k, int & l) const;
312 
314  DoubleMatrix apply(double (*fn)(double)) const {
315  return DoubleMatrix(x.apply(fn),rows,cols);
316  }
317 
319  const DoubleMatrix & operator=(double v);
320 
321  double min(int & k, int & l) const;
322  double max(int & k, int & l) const;
323  void swaprows(int i, int j);
325  void swapcols(int i,int j);
326  std::size_t size() const { return x.size(); }
328  void setCols(int);
330  void setRows(int);
332  void resize(int, int);
333 
334  double trace() const;
335  DoubleMatrix transpose() const;
336  bool testNan() const;
337 
338  /*
339  * NUMERICAL DIAGONALIZATION ROUTINES ETC.
340  */
341 
345  void associateOrderAbs(DoubleVector &v);
346  // Perform on asymmetric matrices M such that
349  void associateOrderAbs(DoubleMatrix & u, DoubleMatrix & v, DoubleVector &w)
350  const;
352  void symmetrise();
355  double compare(const DoubleMatrix & a) const;
356  double compare(const ComplexMatrix & a) const;
358  DoubleVector diagVals() const;
359  DoubleMatrix inverse() const;
360 
364  double diagonalise(DoubleMatrix & u, DoubleMatrix & v, DoubleVector & w) const;
368  double diagonaliseSym(DoubleMatrix & v, DoubleVector & w) const;
369  double diagonaliseSym(ComplexMatrix & v, DoubleVector & w) const;
370 
374  // -[ cos theta sin theta ] A [ cos theta -sin theta ] = diagonal
375  // -[ -sin theta cos theta ] [ sin theta cos theta ]
376  DoubleVector sym2by2(double & theta) const;
378  // -[ cos thetaL sin thetaL ] A [ cos thetaR -sin thetaR ] = diagonal
379  // -[ -sin thetaL cos thetaL ] [ sin thetaR cos thetaR ]
380  DoubleVector asy2by2(double & thetaL, double & thetaR) const;
382  DoubleVector vectorfy() const;
384  double determinant() const;
385 
388  DoubleMatrix ludcmp(double & d) const;
390  void fillArray(double* array, unsigned offset = 0) const;
391 };
392 
393 /*
394  * STANDARD INPUT / OUTPUT
395  */
396 
397 std::istream & operator>>(std::istream & left, DoubleMatrix &M);
398 std::ostream & operator <<(std::ostream &left, const DoubleMatrix &V);
399 
400 /*
401  * NON-MEMBER DIAGONALIZATION ROUTINES ETC.
402  */
403 
405 // [ cos theta sin theta ]
406 // [-sin theta cos theta ]
407 DoubleMatrix rot2d(double theta);
409 // [ -sin(theta) cos(theta) ]
410 // [ cos(theta) sin(theta) ] --
411 DoubleMatrix rot2dTwist(double theta);
412 
415 // [ -cos theta sin theta 0 ]
416 // [ sin theta cos theta 0 ]
417 // [ 0 0 1 ]
418 DoubleMatrix rot3d(double theta);
419 
422 // [ cos thetaL sin thetaL ] A [ cos thetaR -sin thetaR ] = diag
423 // [ -sin thetaL cos thetaL ] [ sin thetaR cos thetaR ]
424 // as given by asy2by2!
426 void positivise(double thetaL, double thetaR, const DoubleVector & diag,
427  ComplexMatrix & u, ComplexMatrix & v);
430 double pythagoras(double a, double b);
431 void diagonaliseJac(DoubleMatrix & a, int n, DoubleVector & d, DoubleMatrix
432  & v, int *nrot);
433 
435 void fillArray(const std::valarray<double>&, double*, unsigned offset = 0);
436 
437 /*
438  * INLINE FUNCTION DEFINITIONS
439  */
440 
441 inline double DoubleMatrix::operator() (int i, int j) const {
442 #ifdef ARRAY_BOUNDS_CHECKING
443  if (i > rows || j > cols || i < 1 || j < 1) {
444  std::ostringstream ii;
445  ii << "Trying to access " << i << "," << j <<
446  "th element of DoubleMatrix\n" << *this;
447  throw ii.str();
448  }
449 #endif
450  return elmt(i,j);
451 }
452 
453 inline double & DoubleMatrix::operator() (int i, int j) {
454 #ifdef ARRAY_BOUNDS_CHECKING
455  if (i > rows || j > cols || i < 1 || j < 1) {
456  std::ostringstream ii;
457  ii << "Trying to access " << i << "," << j <<
458  "th element of DoubleMatrix\n" << *this;
459  throw ii.str();
460  }
461 #endif
462  return elmt(i,j);
463 }
464 
465 inline double DoubleMatrix::display(int i, int j) const {
466 #ifdef ARRAY_BOUNDS_CHECKING
467  if (i < 1 || i > displayRows() || j < 1 || j > displayCols()) {
468  std::ostringstream ii;
469  ii << "Error: Requested display (" << i << "," << j <<
470  ")th element of matrix " << *this;
471  throw ii.str();
472  }
473 #endif
474  return elmt(i,j);
475 }
476 
477 /************************************
478  *
479  * *** COMPLEX VECTOR CLASS ***
480  *
481  ************************************/
482 
484 class ComplexVector : public Indexable<Complex, ComplexVector> {
485 private:
486  int start, end;
487  std::valarray<Complex> x;
488  ComplexVector(const std::valarray<Complex> & valarr, int s, int e)
490  : Indexable<Complex, ComplexVector>(), start(s), end(e), x(valarr) {}
491  friend class ComplexMatrix;
492 
493 
494 public:
495 
496  /*
497  * CONSTRUCTORS
498  */
499 
500  explicit ComplexVector(int e)
501  : Indexable<Complex, ComplexVector>(), start(1), end(e), x(0.0,e) {}
502  ComplexVector(int s, int e)
503  : Indexable<Complex, ComplexVector>(), start(s), end(e), x(0.0,e-s+1) {}
504  ~ComplexVector() {}
505 
507  template<class E>
509  *this = copy_from(x);
510  return *this;
511  }
512 
513  ComplexVector& operator=(const ComplexVector& other) {
514  x.resize(other.x.size());
515  x = other.x;
516  start = other.start;
517  end = other.end;
518  return *this;
519  }
520 
521  template <typename E>
522  ComplexVector(const Xpr<Complex,E> & v)
524  start(v.displayStart()), end(v.displayEnd()),
525  x(0.0,v.displayEnd()-v.displayStart()+1)
526  {
527  assign_from(v);
528  }
529 
530 
531  /*
532  * ELEMENT ACCESS
533  */
534 
535  Complex & operator() (int i);
536  Complex operator() (int i) const { return x[i-start]; }
537  Complex display(int i) const { return x[i-start]; }
538  void set(int i, Complex f);
539 
540  int displayStart() const { return start; }
541  int displayEnd() const { return end; }
542  const ComplexVector & display() const { return *this; }
543 
544 
545 
548  void setEnd(int e);
549 
552  return ComplexVector(x.apply(fn),start,end);
553  }
554 
557  Complex max() const { return x.max(); }
558  Complex min(int & p) const;
559 
560  void swap(int i, int j);
561 
562 };
563 
564 /*
565  * STANDARD INPUT / OUTPUT
566  */
567 
569 std::ostream & operator<<(std::ostream &left, const ComplexVector &V);
571 std::istream & operator>>(std::istream & left, ComplexVector &V);
572 
573 /*
574  * INLINE FUNCTION DEFINITIONS
575  */
576 
578 #ifdef ARRAY_BOUNDS_CHECKING
579  if (i>end || i<start) {
580  std::ostringstream ii;
581  ii << "Trying to access " << i << "th element of DoubleVector\n";
582  ii << "start " << start << " end " << end;
583  ii << *this;
584  throw ii.str();
585  }
586 #endif
587  return x[i-start];
588 }
589 
590 inline void ComplexVector::set(int i, Complex f) {
591 #ifdef ARRAY_BOUNDS_CHECKING
592  if (i < start || i > end) {
593  std::ostringstream ii;
594  ii << "Cannot access " << i << "th element of vector " << *this;
595  throw ii.str();
596  }
597 #endif
598  x[i-start] = f;
599 }
600 
601 /************************************
602  *
603  * *** COMPLEX MATRIX CLASS ***
604  *
605  ************************************/
606 
608 class ComplexMatrix : public MatIndexable<Complex, ComplexMatrix> {
609 private:
610  int rows, cols;
611  std::valarray<Complex> x;
612 
615  // for usage see e.g. http://www.roguewave.com/support/docs/sourcepro/stdlibref/2-10.html
616  std::slice_array<Complex> row(int i) { return x[std::slice((i-1)*cols,cols,1)]; }
617  std::slice_array<Complex> col(int i) { return x[std::slice(i-1,rows,cols)]; }
618  const std::valarray<Complex> row(int i) const {
619  return std::valarray<Complex>(x[std::slice((i-1)*cols,cols,1)]);
620  }
621  const std::valarray<Complex> col(int i) const {
622  return std::valarray<Complex>(x[std::slice(i-1,rows,cols)]);
623  }
624 
625  std::slice_array<Complex> diag() { return x[std::slice(0,rows,cols+1)]; }
626  const std::valarray<Complex> diag() const {
627  return std::valarray<Complex>(x[std::slice(0,rows,cols+1)]);
628  }
629 
631  Complex & elmt(int i, int j) { return x[(i-1)*cols+(j-1)]; }
632  const Complex elmt(int i, int j) const { return x[(i-1)*cols+(j-1)]; }
633 
635  ComplexMatrix(const std::valarray<Complex> & valarr, int r, int c)
636  : MatIndexable<Complex,ComplexMatrix>(), rows(r), cols(c), x(valarr) {}
637 
638  friend class DoubleMatrix;
639 
640 public:
641 
642  /*
643  * CONSTRUCTORS
644  */
645 
647  ComplexMatrix(int r, int c)
648  : MatIndexable<Complex,ComplexMatrix>(), rows(r), cols(c), x(0.0, r*c) {}
650  explicit ComplexMatrix(const DoubleMatrix & m);
652  explicit ComplexMatrix(const ComplexVector &v);
653  ~ComplexMatrix() {}
654 
656  template<class E>
658  *this = copy_from(x);
659  return *this;
660  }
661 
662  ComplexMatrix& operator=(const ComplexMatrix& other) {
663  x.resize(other.x.size());
664  x = other.x;
665  rows = other.rows;
666  cols = other.cols;
667  return *this;
668  }
669 
670  template <typename E>
673  rows(m.displayRows()), cols(m.displayCols()), x(0.0, rows*cols)
674  {
675  assign_from(m);
676  }
677 
678  /*
679  * ELEMENT ACCESS
680  */
681 
682  Complex & operator() (int i, int j);
683  Complex operator() (int i, int j) const;
684  Complex display(int i, int j) const;
685 
686  int displayRows() const { return rows; };
687  int displayCols() const { return cols; };
688  const ComplexMatrix & display() const { return *this; };
689  std::size_t size() const { return x.size(); }
690 
692  const ComplexMatrix & operator=(const Complex &v);
693 
694  Complex min(int & k, int & l) const;
695  void swaprows(int i, int j);
697  void swapcols(int i,int j);
698 
700  void setCols(int);
702  void setRows(int);
704  void resize(int, int);
705 
706  Complex trace() const;
707  ComplexMatrix transpose() const;
708  ComplexMatrix hermitianConjugate() const;
709  ComplexMatrix complexConjugate() const;
711  DoubleMatrix real() const;
713  DoubleMatrix imag() const;
714 
715  /*
716  * NUMERICAL DIAGONALIZATION ROUTINES ETC.
717  */
719  double nonHermiticity() const;
721  void symmetrise();
723  double compare(const ComplexMatrix & a) const;
728  DoubleMatrix makeHermitianRealForDiag() const;
733  double diagonaliseHerm(ComplexMatrix & v, DoubleVector & w) const;
738  double takagi(ComplexMatrix & v, ComplexVector & w) const;
743  double diagonaliseSym2by2(ComplexMatrix & v, ComplexVector & w) const;
748  double diagonalise(ComplexMatrix & u, ComplexMatrix & v, DoubleVector & w)
749  const;
752  ComplexMatrix a(rows, cols);
753  for (int i=1; i<=rows; i++)
754  for (int j=1; j<=rows; j++)
755  a(i, j) = fn(display(i, j));
756 
757  return a;
758  }
759 };
760 
761 /*
762  * STANDARD INPUT / OUTPUT
763  */
764 
766 std::istream & operator>>(std::istream & left, ComplexMatrix &M);
768 std::ostream & operator<<(std::ostream &left, const ComplexMatrix &V);
769 
770 /*
771  * INLINE FUNCTION DEFINITIONS
772  */
773 
774 inline Complex & ComplexMatrix::operator() (int i, int j) {
775 #ifdef ARRAY_BOUNDS_CHECKING
776  if (i > rows || j > cols || i < 0 || j < 0) {
777  std::ostringstream ii;
778  ii << "Trying to access " << i << "," << j <<
779  "th element of ComplexMatrix\n" << *this;
780  throw ii.str();
781  }
782 #endif
783  return elmt(i,j);
784 }
785 
786 inline Complex ComplexMatrix::operator() (int i, int j) const {
787 #ifdef ARRAY_BOUNDS_CHECKING
788  if (i > rows || j > cols || i < 0 || j < 0) {
789  std::ostringstream ii;
790  ii << "Trying to access " << i << "," << j <<
791  "th element of ComplexMatrix\n" << *this;
792  throw ii.str();
793  }
794 #endif
795  return elmt(i,j);
796 }
797 
798 inline Complex ComplexMatrix::display(int i, int j) const {
799 #ifdef ARRAY_BOUNDS_CHECKING
800  if (i < 1 || i > displayRows() || j < 1 || j > displayCols()) {
801  std::ostringstream ii;
802  ii << "Error: Requested display (" << i << "," << j <<
803  ")th element of matrix " << *this;
804  throw ii.str();
805  }
806 #endif
807  return elmt(i,j);
808 }
809 
810 
811 #endif
global variable declaration
Definition: def.cpp:13
ComplexVector(int e)
Definition: linalg.h:500
void setCols(int)
change number of columns (Warning: can be slow because it internally copys a std::valarray<double>) ...
Definition: linalg.cpp:1173
DoubleVector is of variable length, and contains double precision.
Definition: linalg.h:35
matrix form of xpr object
Definition: xpr-base.h:190
DoubleMatrix rot2d(double theta)
Returns a 2x2 orthogonal matrix of rotation by angle theta.
Definition: linalg.cpp:657
int displayEnd() const
returns end of dimension
Definition: linalg.h:94
ComplexVector apply(Complex(*fn)(Complex)) const
Apply fn to every element.
Definition: linalg.h:551
DoubleVector & operator=(const Xpr< double, E > &x)
this is the only required operator to make this class work with ETs
Definition: linalg.h:62
Under-object for speed upgrade of linear algebra.
Definition: xpr-base.h:24
double sumElements() const
Returns sum of absolute values of all elements.
Definition: linalg.h:132
switches (options) and parameters such as default fermion masses, required accuracy etc ...
bool testNan(double f)
Definition: utils.cpp:64
void swaprows(int i, int j)
Obvious elementary row/column operations.
Definition: linalg.cpp:197
void setRows(int)
change number of rows (Warning: can be slow because it internally copys a std::valarray<double>) ...
Definition: linalg.cpp:226
double min(int &k, int &l) const
minimum element
Definition: linalg.cpp:175
void setRows(int)
change number of rows (Warning: can be slow because it internally copys a std::valarray<double>) ...
Definition: linalg.cpp:1190
double dot(const DoubleVector &v) const
Definition: linalg.h:137
std::size_t size() const
returns whole matrix
Definition: linalg.h:689
std::istream & operator>>(std::istream &left, DoubleVector &V)
inputs a vector
int displayRows() const
Routines for outputting size of matrix.
Definition: linalg.h:290
DoubleMatrix(int r, int c)
Constructor for matrix of zeroes 1..r,1..c.
Definition: linalg.h:253
ComplexVector(int s, int e)
Definition: linalg.h:502
provides indexing of xpr type objects
Definition: xpr-base.h:79
double max() const
maximum element in vector
Definition: linalg.h:114
DoubleMatrix & operator=(const MatXpr< double, E > &x)
this is the only required operator to make this class work with ETs
Definition: linalg.h:261
A few handy bits and pieces - little mathematical functions and the like.
void resize(int, int)
resize matrix (Warning: can be slow because it internally copys a std::valarray<double>) ...
Definition: linalg.cpp:239
DoubleVector abs() const
returns absolute of vector
Definition: linalg.h:105
double trace() const
trace must only be performed on a square matrix
Definition: linalg.cpp:246
double display(int i) const
display one element
Definition: linalg.h:90
void fillArray(const std::valarray< double > &, double *, unsigned offset=0)
fill array from valarray, starting at offset
Definition: linalg.cpp:1761
void setCols(int)
change number of columns (Warning: can be slow because it internally copys a std::valarray<double>) ...
Definition: linalg.cpp:209
void diagonaliseSvd(DoubleMatrix &a, DoubleVector &w, DoubleMatrix &v)
Diagonalisation routines.
Definition: linalg.cpp:722
double max(int &k, int &l) const
Definition: linalg.cpp:186
Complex & operator()(int i)
reference one element
Definition: linalg.h:577
matrix of complex double values dimensions (rows x cols)
Definition: linalg.h:608
DoubleVector(int s, int e)
Definition: linalg.h:53
Complex max() const
maximum absolute value
Definition: linalg.h:557
void swapcols(int i, int j)
Swaps column i with column j.
Definition: linalg.cpp:1166
int displayStart() const
start of dimension
Definition: linalg.h:93
Symbolic object for matrices for speed upgrade.
int displayStart() const
displays start of dimension
Definition: linalg.h:540
const ComplexVector & display() const
displays whole vector
Definition: linalg.h:542
Symbolic object for vectors for speed upgrade.
Matrix from 1..rows, 1..cols of double values.
Definition: linalg.h:214
void symmetrise()
fills in bottom-left hand corner of a matrix
Definition: linalg.cpp:327
DoubleVector(int e)
Definition: linalg.h:51
void positivise(double thetaL, double thetaR, const DoubleVector &diag, ComplexMatrix &u, ComplexMatrix &v)
= mdiagpositive
Definition: linalg.cpp:698
double compare(const ComplexMatrix &a) const
Returns the sum of the modulus of the difference of each element.
Definition: linalg.cpp:1485
DoubleMatrix transpose() const
can be any size
Definition: linalg.cpp:260
Complex & operator()(int i, int j)
Returns ijth element.
Definition: linalg.h:774
void set(int i, double f)
sets ith element
Definition: linalg.h:192
int displayEnd() const
displays end of dimension
Definition: linalg.h:541
double absmax(int &p) const
returns absolute maximum
Definition: linalg.h:130
double operator()(int i, int j) const
to reference one element
Definition: linalg.h:441
ComplexMatrix apply(Complex(*fn)(Complex)) const
Applies fn to every element.
Definition: linalg.h:751
Indexable form of xpr matrix.
Definition: xpr-base.h:373
void ludcmp(DoubleMatrix &a, int n, int *indx, double &d)
Get rid of int n.
Definition: numerics.cpp:1889
complex numbers and operators between them
Vector of double complex values.
Definition: linalg.h:484
Complex display(int i) const
display ith element
Definition: linalg.h:537
double diagonalise(ComplexMatrix &u, ComplexMatrix &v, DoubleVector &w) const
Now for any Complex matrix.
Definition: linalg.cpp:1393
std::ostream & operator<<(std::ostream &left, const DoubleVector &V)
prints a vector out formatted, maximum of 5 elements per line
Complex display(int i, int j) const
returns ijth element
Definition: linalg.h:798
void swaprows(int i, int j)
Obvious elementary row/column operations.
Definition: linalg.cpp:1160
Complex min(int &k, int &l) const
Definition: linalg.cpp:1147
const DoubleVector & display() const
returns whole thing
Definition: linalg.h:95
double diagonalise(DoubleMatrix &u, DoubleMatrix &v, DoubleVector &w) const
Definition: linalg.cpp:411
DoubleMatrix rot3d(double theta)
Definition: linalg.cpp:682
const DoubleMatrix & display() const
whole matrix returned
Definition: linalg.h:292
void set(int i, Complex f)
set ith element to f
Definition: linalg.h:590
void resize(int, int)
resize matrix (Warning: can be slow because it internally copys a std::valarray<double>) ...
Definition: linalg.cpp:1204
void symmetrise()
Fills in lower bottom half of a square matrix copying the top right.
Definition: linalg.cpp:1471
ComplexMatrix(int r, int c)
Default constructor: full of zeroes.
Definition: linalg.h:647
DoubleVector apply(double(*fn)(double)) const
applies fn to every element in a vector
Definition: linalg.h:109
ComplexVector & operator=(const Xpr< Complex, E > &x)
this is the only required operator to make this class work with ETs
Definition: linalg.h:508
int theta(double a)
Standard theta function: 1 is a>0, 0 otherwise.
Definition: utils.cpp:25
double compare(const DoubleMatrix &a) const
Definition: linalg.cpp:341
DoubleMatrix rot2dTwist(double theta)
Returns 2 by 2 orthogonal mixing matrix.
Definition: linalg.cpp:668
ComplexMatrix & operator=(const MatXpr< Complex, E > &x)
this is the only required operator to make this class work with ETs
Definition: linalg.h:657
DoubleMatrix apply(double(*fn)(double)) const
Applies fn to every element of a matrix.
Definition: linalg.h:314
double & operator()(int i)
Reference a single element.
Definition: linalg.h:179
drop-in replacement for the original home-grown Complex class
Definition: mycomplex.h:17