SoftSUSY is hosted by Hepforge, IPPP Durham
SOFTSUSY  4.0
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 
100  void setEnd(int e);
102  DoubleVector abs() const { return DoubleVector(std::abs(x), start, end); }
103 
104  double norm() const;
106  DoubleVector apply(double (*fn)(double)) const {
107  return DoubleVector(x.apply(fn),start,end);
108  }
109 
110  double nmin(int & p) const;
111  double max() const { return x.max(); }
112 
113  double max(int & p) const {
114  double m = double(0);
115  int i;
116  for (i=displayStart(); i<=displayEnd(); i++)
117  if (display(i) > m) {
118  m = display(i);
119  p = i;
120  }
121  return m;
122  }
123 
124  double min(int & p) const;
125 
127  double absmax(int & p) const { return abs().max(p); }
129  double sumElements() const { return std::abs(x).sum(); }
130  std::size_t size() const { return x.size(); }
131 
132  void swap(int i, int j);
133 
134  double dot(const DoubleVector & v) const {
135  return (x * v.x).sum();
136  }
137 
138  DoubleVector sort() const {
139  DoubleVector temp(*this), saveIt(*this);
140 
141  int j, pos;
142  for (j=start; j<=end; j++) {
143  double d = temp.min(pos);
144  if (d < 1.0e66) saveIt(j) = d;
145 
146  temp(pos) = 6.66e66; // take that element out of the min calculation
147  }
148  return saveIt;
149  }
150 
152  DoubleVector divide(DoubleVector const &v) const;
154  double average() const;
155  int closest(double a) const;
158  void fillArray(double* array, unsigned offset = 0) const;
159 };
160 
161 
162 /*
163  * STANDARD INPUT / OUTPUT
164  */
165 
167 std::ostream & operator<<(std::ostream &left, const DoubleVector &V);
169 std::istream & operator>>(std::istream & left, DoubleVector &V);
170 
171 
172 /*
173  * INLINE FUNCTION DEFINITIONS
174  */
175 
176 inline double & DoubleVector::operator() (int i) {
177 #ifdef ARRAY_BOUNDS_CHECKING
178  if (i>end || i<start) {
179  std::ostringstream ii;
180  ii << "Trying to access " << i << "th element of DoubleVector\n";
181  ii << "start " << start << " end " << end;
182  ii << *this;
183  throw ii.str();
184  }
185 #endif
186  return x[i-start];
187 }
188 
189 inline void DoubleVector::set(int i, double f) {
190 #ifdef ARRAY_BOUNDS_CHECKING
191  if (i < start || i > end)
192  {
193  std::ostringstream ii;
194  ii << "Cannot access " << i << "th element of vector " << *this;
195  throw ii.str();
196  }
197 #endif
198  x[i-start] = f;
199 }
200 
201 
202 /************************************
203  *
204  * *** DOUBLE MATRIX CLASS ***
205  *
206  ************************************/
207 
208 
209 class ComplexVector; class ComplexMatrix;
211 class DoubleMatrix : public MatIndexable<double, DoubleMatrix> {
212 private:
213  int rows, cols;
214  std::valarray<double> x;
215 
218  // for usage see e.g. http://www.roguewave.com/support/docs/sourcepro/stdlibref/2-10.html
219  std::slice_array<double> row(int i) { return x[std::slice((i-1)*cols,cols,1)]; }
220  std::slice_array<double> col(int i) { return x[std::slice(i-1,rows,cols)]; }
221  const std::valarray<double> row(int i) const {
222  return std::valarray<double>(x[std::slice((i-1)*cols,cols,1)]);
223  }
224  const std::valarray<double> col(int i) const {
225  return std::valarray<double>(x[std::slice(i-1,rows,cols)]);
226  }
227 
228  std::slice_array<double> diag() { return x[std::slice(0,rows,cols+1)]; }
229  const std::valarray<double> diag() const {
230  return std::valarray<double>(x[std::slice(0,rows,cols+1)]);
231  }
232 
234  double & elmt(int i, int j) { return x[(i-1) * cols + (j-1)]; }
235  double elmt(int i, int j) const { return x[(i-1) * cols + (j-1)]; }
236 
238  DoubleMatrix(const std::valarray<double> & valarr, int r, int c)
239  : MatIndexable<double,DoubleMatrix>(), rows(r), cols(c), x(valarr) {}
240 
241  friend class ComplexMatrix;
242 
243 public:
244 
245  /*
246  * CONSTRUCTORS
247  */
248 
250  DoubleMatrix(int r, int c)
251  : MatIndexable<double,DoubleMatrix>(), rows(r), cols(c), x(0.0, r*c) {}
253  explicit DoubleMatrix(const DoubleVector &v);
254  ~DoubleMatrix() {}
255 
257  template<class E>
259  *this = copy_from(x);
260  return *this;
261  }
262 
263  DoubleMatrix& operator=(const DoubleMatrix& other) {
264  x.resize(other.x.size());
265  x = other.x;
266  rows = other.rows;
267  cols = other.cols;
268  return *this;
269  }
270 
271  template <typename E>
274  rows(m.displayRows()), cols(m.displayCols()), x(0.0, rows*cols)
275  {
276  assign_from(m);
277  }
278 
279  /*
280  * ELEMENT ACCESS
281  */
282  double operator() (int i, int j) const;
283  double & operator() (int i, int j);
284  double display(int i, int j) const;
285  double sumElements() const { return abs(x).sum(); }
287  int displayRows() const { return rows; };
288  int displayCols() const { return cols; };
289  const DoubleMatrix & display() const { return *this; }
290 
291  DoubleVector displayRow(int i) const {
292  DoubleVector temp(displayCols());
293  int j; for (j=1; j<=displayCols(); j++) temp(j) = display(i, j);
294  return temp;
295  }
296 
297  DoubleVector displayCol(int i) const {
298  DoubleVector temp(displayRows());
299  int j; for (j=1; j<=displayRows(); j++) temp(j) = display(j, i);
300  return temp;
301  }
302 
303 
304  void operator+=(DoubleMatrix & f) {
305  x += f.x;
306  }
307 
308  double nmin(int & k, int & l) const;
309 
311  DoubleMatrix apply(double (*fn)(double)) const {
312  return DoubleMatrix(x.apply(fn),rows,cols);
313  }
314 
316  const DoubleMatrix & operator=(double v);
317 
318  double min(int & k, int & l) const;
319  double max(int & k, int & l) const;
320  void swaprows(int i, int j);
322  void swapcols(int i,int j);
323  std::size_t size() const { return x.size(); }
325  void setCols(int);
327  void setRows(int);
329  void resize(int, int);
330 
331  double trace() const;
332  DoubleMatrix transpose() const;
333  bool testNan() const;
334 
335  /*
336  * NUMERICAL DIAGONALIZATION ROUTINES ETC.
337  */
338 
342  void associateOrderAbs(DoubleVector &v);
343  // Perform on asymmetric matrices M such that
346  void associateOrderAbs(DoubleMatrix & u, DoubleMatrix & v, DoubleVector &w)
347  const;
349  void symmetrise();
352  double compare(const DoubleMatrix & a) const;
353  double compare(const ComplexMatrix & a) const;
355  DoubleVector diagVals() const;
356  DoubleMatrix inverse() const;
357 
361  double diagonalise(DoubleMatrix & u, DoubleMatrix & v, DoubleVector & w) const;
365  double diagonaliseSym(DoubleMatrix & v, DoubleVector & w) const;
366  double diagonaliseSym(ComplexMatrix & v, DoubleVector & w) const;
367 
371  // -[ cos theta sin theta ] A [ cos theta -sin theta ] = diagonal
372  // -[ -sin theta cos theta ] [ sin theta cos theta ]
373  DoubleVector sym2by2(double & theta) const;
375  // -[ cos thetaL sin thetaL ] A [ cos thetaR -sin thetaR ] = diagonal
376  // -[ -sin thetaL cos thetaL ] [ sin thetaR cos thetaR ]
377  DoubleVector asy2by2(double & thetaL, double & thetaR) const;
379  DoubleVector vectorfy() const;
381  double determinant() const;
382 
385  DoubleMatrix ludcmp(double & d) const;
387  void fillArray(double* array, unsigned offset = 0) const;
388 };
389 
390 /*
391  * STANDARD INPUT / OUTPUT
392  */
393 
394 std::istream & operator>>(std::istream & left, DoubleMatrix &M);
395 std::ostream & operator <<(std::ostream &left, const DoubleMatrix &V);
396 
397 /*
398  * NON-MEMBER DIAGONALIZATION ROUTINES ETC.
399  */
400 
402 // [ cos theta sin theta ]
403 // [-sin theta cos theta ]
404 DoubleMatrix rot2d(double theta);
406 // [ -sin(theta) cos(theta) ]
407 // [ cos(theta) sin(theta) ] --
408 DoubleMatrix rot2dTwist(double theta);
409 
412 // [ -cos theta sin theta 0 ]
413 // [ sin theta cos theta 0 ]
414 // [ 0 0 1 ]
415 DoubleMatrix rot3d(double theta);
416 
419 // [ cos thetaL sin thetaL ] A [ cos thetaR -sin thetaR ] = diag
420 // [ -sin thetaL cos thetaL ] [ sin thetaR cos thetaR ]
421 // as given by asy2by2!
423 void positivise(double thetaL, double thetaR, const DoubleVector & diag,
424  ComplexMatrix & u, ComplexMatrix & v);
427 double pythagoras(double a, double b);
428 void diagonaliseJac(DoubleMatrix & a, int n, DoubleVector & d, DoubleMatrix
429  & v, int *nrot);
430 
432 void fillArray(const std::valarray<double>&, double*, unsigned offset = 0);
433 
434 /*
435  * INLINE FUNCTION DEFINITIONS
436  */
437 
438 inline double DoubleMatrix::operator() (int i, int j) const {
439 #ifdef ARRAY_BOUNDS_CHECKING
440  if (i > rows || j > cols || i < 1 || j < 1) {
441  std::ostringstream ii;
442  ii << "Trying to access " << i << "," << j <<
443  "th element of DoubleMatrix\n" << *this;
444  throw ii.str();
445  }
446 #endif
447  return elmt(i,j);
448 }
449 
450 inline double & DoubleMatrix::operator() (int i, int j) {
451 #ifdef ARRAY_BOUNDS_CHECKING
452  if (i > rows || j > cols || i < 1 || j < 1) {
453  std::ostringstream ii;
454  ii << "Trying to access " << i << "," << j <<
455  "th element of DoubleMatrix\n" << *this;
456  throw ii.str();
457  }
458 #endif
459  return elmt(i,j);
460 }
461 
462 inline double DoubleMatrix::display(int i, int j) const {
463 #ifdef ARRAY_BOUNDS_CHECKING
464  if (i < 1 || i > displayRows() || j < 1 || j > displayCols()) {
465  std::ostringstream ii;
466  ii << "Error: Requested display (" << i << "," << j <<
467  ")th element of matrix " << *this;
468  throw ii.str();
469  }
470 #endif
471  return elmt(i,j);
472 }
473 
474 /************************************
475  *
476  * *** COMPLEX VECTOR CLASS ***
477  *
478  ************************************/
479 
481 class ComplexVector : public Indexable<Complex, ComplexVector> {
482 private:
483  int start, end;
484  std::valarray<Complex> x;
485  ComplexVector(const std::valarray<Complex> & valarr, int s, int e)
487  : Indexable<Complex, ComplexVector>(), start(s), end(e), x(valarr) {}
488  friend class ComplexMatrix;
489 
490 
491 public:
492 
493  /*
494  * CONSTRUCTORS
495  */
496 
497  explicit ComplexVector(int e)
498  : Indexable<Complex, ComplexVector>(), start(1), end(e), x(0.0,e) {}
499  ComplexVector(int s, int e)
500  : Indexable<Complex, ComplexVector>(), start(s), end(e), x(0.0,e-s+1) {}
501  ~ComplexVector() {}
502 
504  template<class E>
506  *this = copy_from(x);
507  return *this;
508  }
509 
510  ComplexVector& operator=(const ComplexVector& other) {
511  x.resize(other.x.size());
512  x = other.x;
513  start = other.start;
514  end = other.end;
515  return *this;
516  }
517 
518  template <typename E>
519  ComplexVector(const Xpr<Complex,E> & v)
521  start(v.displayStart()), end(v.displayEnd()),
522  x(0.0,v.displayEnd()-v.displayStart()+1)
523  {
524  assign_from(v);
525  }
526 
527 
528  /*
529  * ELEMENT ACCESS
530  */
531 
532  Complex & operator() (int i);
533  Complex operator() (int i) const { return x[i-start]; }
534  Complex display(int i) const { return x[i-start]; }
535  void set(int i, Complex f);
536 
537  int displayStart() const { return start; }
538  int displayEnd() const { return end; }
539  const ComplexVector & display() const { return *this; }
540 
541 
542 
545  void setEnd(int e);
546 
549  return ComplexVector(x.apply(fn),start,end);
550  }
551 
554  Complex max() const { return x.max(); }
555  Complex min(int & p) const;
556 
557  void swap(int i, int j);
558 
559 };
560 
561 /*
562  * STANDARD INPUT / OUTPUT
563  */
564 
566 std::ostream & operator<<(std::ostream &left, const ComplexVector &V);
568 std::istream & operator>>(std::istream & left, ComplexVector &V);
569 
570 /*
571  * INLINE FUNCTION DEFINITIONS
572  */
573 
575 #ifdef ARRAY_BOUNDS_CHECKING
576  if (i>end || i<start) {
577  std::ostringstream ii;
578  ii << "Trying to access " << i << "th element of DoubleVector\n";
579  ii << "start " << start << " end " << end;
580  ii << *this;
581  throw ii.str();
582  }
583 #endif
584  return x[i-start];
585 }
586 
587 inline void ComplexVector::set(int i, Complex f) {
588 #ifdef ARRAY_BOUNDS_CHECKING
589  if (i < start || i > end) {
590  std::ostringstream ii;
591  ii << "Cannot access " << i << "th element of vector " << *this;
592  throw ii.str();
593  }
594 #endif
595  x[i-start] = f;
596 }
597 
598 /************************************
599  *
600  * *** COMPLEX MATRIX CLASS ***
601  *
602  ************************************/
603 
605 class ComplexMatrix : public MatIndexable<Complex, ComplexMatrix> {
606 private:
607  int rows, cols;
608  std::valarray<Complex> x;
609 
612  // for usage see e.g. http://www.roguewave.com/support/docs/sourcepro/stdlibref/2-10.html
613  std::slice_array<Complex> row(int i) { return x[std::slice((i-1)*cols,cols,1)]; }
614  std::slice_array<Complex> col(int i) { return x[std::slice(i-1,rows,cols)]; }
615  const std::valarray<Complex> row(int i) const {
616  return std::valarray<Complex>(x[std::slice((i-1)*cols,cols,1)]);
617  }
618  const std::valarray<Complex> col(int i) const {
619  return std::valarray<Complex>(x[std::slice(i-1,rows,cols)]);
620  }
621 
622  std::slice_array<Complex> diag() { return x[std::slice(0,rows,cols+1)]; }
623  const std::valarray<Complex> diag() const {
624  return std::valarray<Complex>(x[std::slice(0,rows,cols+1)]);
625  }
626 
628  Complex & elmt(int i, int j) { return x[(i-1)*cols+(j-1)]; }
629  const Complex elmt(int i, int j) const { return x[(i-1)*cols+(j-1)]; }
630 
632  ComplexMatrix(const std::valarray<Complex> & valarr, int r, int c)
633  : MatIndexable<Complex,ComplexMatrix>(), rows(r), cols(c), x(valarr) {}
634 
635  friend class DoubleMatrix;
636 
637 public:
638 
639  /*
640  * CONSTRUCTORS
641  */
642 
644  ComplexMatrix(int r, int c)
645  : MatIndexable<Complex,ComplexMatrix>(), rows(r), cols(c), x(0.0, r*c) {}
647  explicit ComplexMatrix(const DoubleMatrix & m);
649  explicit ComplexMatrix(const ComplexVector &v);
650  ~ComplexMatrix() {}
651 
653  template<class E>
655  *this = copy_from(x);
656  return *this;
657  }
658 
659  ComplexMatrix& operator=(const ComplexMatrix& other) {
660  x.resize(other.x.size());
661  x = other.x;
662  rows = other.rows;
663  cols = other.cols;
664  return *this;
665  }
666 
667  template <typename E>
670  rows(m.displayRows()), cols(m.displayCols()), x(0.0, rows*cols)
671  {
672  assign_from(m);
673  }
674 
675  /*
676  * ELEMENT ACCESS
677  */
678 
679  Complex & operator() (int i, int j);
680  Complex operator() (int i, int j) const;
681  Complex display(int i, int j) const;
682 
683  int displayRows() const { return rows; };
684  int displayCols() const { return cols; };
685  const ComplexMatrix & display() const { return *this; };
686  std::size_t size() const { return x.size(); }
687 
689  const ComplexMatrix & operator=(const Complex &v);
690 
691  Complex min(int & k, int & l) const;
692  void swaprows(int i, int j);
694  void swapcols(int i,int j);
695 
697  void setCols(int);
699  void setRows(int);
701  void resize(int, int);
702 
703  Complex trace() const;
704  ComplexMatrix transpose() const;
705  ComplexMatrix hermitianConjugate() const;
706  ComplexMatrix complexConjugate() const;
708  DoubleMatrix real() const;
710  DoubleMatrix imag() const;
711 
712  /*
713  * NUMERICAL DIAGONALIZATION ROUTINES ETC.
714  */
716  double nonHermiticity() const;
718  void symmetrise();
720  double compare(const ComplexMatrix & a) const;
725  DoubleMatrix makeHermitianRealForDiag() const;
730  double diagonaliseHerm(ComplexMatrix & v, DoubleVector & w) const;
735  double takagi(ComplexMatrix & v, ComplexVector & w) const;
740  double diagonaliseSym2by2(ComplexMatrix & v, ComplexVector & w) const;
745  double diagonalise(ComplexMatrix & u, ComplexMatrix & v, DoubleVector & w)
746  const;
749  ComplexMatrix a(rows, cols);
750  for (int i=1; i<=rows; i++)
751  for (int j=1; j<=rows; j++)
752  a(i, j) = fn(display(i, j));
753 
754  return a;
755  }
756 };
757 
758 /*
759  * STANDARD INPUT / OUTPUT
760  */
761 
763 std::istream & operator>>(std::istream & left, ComplexMatrix &M);
765 std::ostream & operator<<(std::ostream &left, const ComplexMatrix &V);
766 
767 /*
768  * INLINE FUNCTION DEFINITIONS
769  */
770 
771 inline Complex & ComplexMatrix::operator() (int i, int j) {
772 #ifdef ARRAY_BOUNDS_CHECKING
773  if (i > rows || j > cols || i < 0 || j < 0) {
774  std::ostringstream ii;
775  ii << "Trying to access " << i << "," << j <<
776  "th element of ComplexMatrix\n" << *this;
777  throw ii.str();
778  }
779 #endif
780  return elmt(i,j);
781 }
782 
783 inline Complex ComplexMatrix::operator() (int i, int j) const {
784 #ifdef ARRAY_BOUNDS_CHECKING
785  if (i > rows || j > cols || i < 0 || j < 0) {
786  std::ostringstream ii;
787  ii << "Trying to access " << i << "," << j <<
788  "th element of ComplexMatrix\n" << *this;
789  throw ii.str();
790  }
791 #endif
792  return elmt(i,j);
793 }
794 
795 inline Complex ComplexMatrix::display(int i, int j) const {
796 #ifdef ARRAY_BOUNDS_CHECKING
797  if (i < 1 || i > displayRows() || j < 1 || j > displayCols()) {
798  std::ostringstream ii;
799  ii << "Error: Requested display (" << i << "," << j <<
800  ")th element of matrix " << *this;
801  throw ii.str();
802  }
803 #endif
804  return elmt(i,j);
805 }
806 
807 
808 #endif
global variable declaration
Definition: def.cpp:13
ComplexVector(int e)
Definition: linalg.h:497
void setCols(int)
change number of columns (Warning: can be slow because it internally copys a std::valarray<double>) ...
Definition: linalg.cpp:1168
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:652
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:548
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:129
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:192
void setRows(int)
change number of rows (Warning: can be slow because it internally copys a std::valarray<double>) ...
Definition: linalg.cpp:221
double min(int &k, int &l) const
minimum element
Definition: linalg.cpp:170
void setRows(int)
change number of rows (Warning: can be slow because it internally copys a std::valarray<double>) ...
Definition: linalg.cpp:1185
double dot(const DoubleVector &v) const
Definition: linalg.h:134
std::size_t size() const
returns whole matrix
Definition: linalg.h:686
std::istream & operator>>(std::istream &left, DoubleVector &V)
inputs a vector
int displayRows() const
Routines for outputting size of matrix.
Definition: linalg.h:287
DoubleMatrix(int r, int c)
Constructor for matrix of zeroes 1..r,1..c.
Definition: linalg.h:250
ComplexVector(int s, int e)
Definition: linalg.h:499
provides indexing of xpr type objects
Definition: xpr-base.h:79
double max() const
maximum element in vector
Definition: linalg.h:111
DoubleMatrix & operator=(const MatXpr< double, E > &x)
this is the only required operator to make this class work with ETs
Definition: linalg.h:258
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:234
DoubleVector abs() const
returns absolute of vector
Definition: linalg.h:102
double trace() const
trace must only be performed on a square matrix
Definition: linalg.cpp:241
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:1733
void setCols(int)
change number of columns (Warning: can be slow because it internally copys a std::valarray<double>) ...
Definition: linalg.cpp:204
void diagonaliseSvd(DoubleMatrix &a, DoubleVector &w, DoubleMatrix &v)
Diagonalisation routines.
Definition: linalg.cpp:717
double max(int &k, int &l) const
Definition: linalg.cpp:181
Complex & operator()(int i)
reference one element
Definition: linalg.h:574
matrix of complex double values dimensions (rows x cols)
Definition: linalg.h:605
DoubleVector(int s, int e)
Definition: linalg.h:53
Complex max() const
maximum absolute value
Definition: linalg.h:554
void swapcols(int i, int j)
Swaps column i with column j.
Definition: linalg.cpp:1161
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:537
const ComplexVector & display() const
displays whole vector
Definition: linalg.h:539
Symbolic object for vectors for speed upgrade.
Matrix from 1..rows, 1..cols of double values.
Definition: linalg.h:211
void symmetrise()
fills in bottom-left hand corner of a matrix
Definition: linalg.cpp:322
DoubleVector(int e)
Definition: linalg.h:51
void positivise(double thetaL, double thetaR, const DoubleVector &diag, ComplexMatrix &u, ComplexMatrix &v)
= mdiagpositive
Definition: linalg.cpp:693
double compare(const ComplexMatrix &a) const
Returns the sum of the modulus of the difference of each element.
Definition: linalg.cpp:1480
DoubleMatrix transpose() const
can be any size
Definition: linalg.cpp:255
Complex & operator()(int i, int j)
Returns ijth element.
Definition: linalg.h:771
void set(int i, double f)
sets ith element
Definition: linalg.h:189
int displayEnd() const
displays end of dimension
Definition: linalg.h:538
double absmax(int &p) const
returns absolute maximum
Definition: linalg.h:127
double operator()(int i, int j) const
to reference one element
Definition: linalg.h:438
ComplexMatrix apply(Complex(*fn)(Complex)) const
Applies fn to every element.
Definition: linalg.h:748
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:1682
complex numbers and operators between them
Vector of double complex values.
Definition: linalg.h:481
Complex display(int i) const
display ith element
Definition: linalg.h:534
double diagonalise(ComplexMatrix &u, ComplexMatrix &v, DoubleVector &w) const
Now for any Complex matrix.
Definition: linalg.cpp:1388
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:795
void swaprows(int i, int j)
Obvious elementary row/column operations.
Definition: linalg.cpp:1155
Complex min(int &k, int &l) const
Definition: linalg.cpp:1142
const DoubleVector & display() const
returns whole thing
Definition: linalg.h:95
double diagonalise(DoubleMatrix &u, DoubleMatrix &v, DoubleVector &w) const
Definition: linalg.cpp:406
DoubleMatrix rot3d(double theta)
Definition: linalg.cpp:677
const DoubleMatrix & display() const
whole matrix returned
Definition: linalg.h:289
void set(int i, Complex f)
set ith element to f
Definition: linalg.h:587
void resize(int, int)
resize matrix (Warning: can be slow because it internally copys a std::valarray<double>) ...
Definition: linalg.cpp:1199
void symmetrise()
Fills in lower bottom half of a square matrix copying the top right.
Definition: linalg.cpp:1466
ComplexMatrix(int r, int c)
Default constructor: full of zeroes.
Definition: linalg.h:644
DoubleVector apply(double(*fn)(double)) const
applies fn to every element in a vector
Definition: linalg.h:106
ComplexVector & operator=(const Xpr< Complex, E > &x)
this is the only required operator to make this class work with ETs
Definition: linalg.h:505
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:336
DoubleMatrix rot2dTwist(double theta)
Returns 2 by 2 orthogonal mixing matrix.
Definition: linalg.cpp:663
ComplexMatrix & operator=(const MatXpr< Complex, E > &x)
this is the only required operator to make this class work with ETs
Definition: linalg.h:654
DoubleMatrix apply(double(*fn)(double)) const
Applies fn to every element of a matrix.
Definition: linalg.h:311
double & operator()(int i)
Reference a single element.
Definition: linalg.h:176
drop-in replacement for the original home-grown Complex class
Definition: mycomplex.h:17