ampsci
c++ program for high-precision atomic structure calculations of single-valence systems
Loading...
Searching...
No Matches
RDMatrix.hpp
1#pragma once
2#include "LinAlg/Matrix.hpp"
3#include "Maths/Grid.hpp"
4#include "Maths/Interpolator.hpp"
5#include "Wavefunction/DiracSpinor.hpp"
6#include <cassert>
7#include <iostream>
8#include <type_traits>
9
10// XXX NOTE: with multiplication defined as int... identity is NOT 1!
11
12namespace MBPT {
13
43template <typename T>
44class RDMatrix {
45
46 std::size_t m_i0, m_stride;
47 std::size_t m_size;
48 std::size_t m_g_size;
49 LinAlg::Matrix<T> m_ff, m_fg, m_gf, m_gg;
50 bool m_incl_g;
51 std::shared_ptr<const Grid> m_rgrid; // "full" grid
52 std::vector<double> sub_r{}; // sub grid
53
54public:
55 //============================================================================
56
57 RDMatrix(std::size_t i0, std::size_t stride, std::size_t size, bool incl_g,
58 std::shared_ptr<const Grid> rgrid)
59 : m_i0(i0),
60 m_stride(stride),
61 m_size(size),
62 m_g_size(incl_g ? size : 0),
63 m_ff(m_size),
64 m_fg(m_g_size),
65 m_gf(m_g_size),
66 m_gg(m_g_size),
67 m_incl_g(incl_g),
68 m_rgrid(rgrid) {
69 //------------------
70 // create vector of r on sub-grid, used to interpolate values onto full
71 const auto &r = m_rgrid->r();
72 sub_r.reserve(m_size);
73 assert(m_i0 + m_stride * m_size <= r.size());
74 for (std::size_t i = 0; i < m_size; ++i) {
75 sub_r.push_back(r[index_to_fullgrid(i)]);
76 }
77 assert(m_i0 + m_stride * m_size <= r.size());
78 assert(sub_r[1] == r[index_to_fullgrid(1)]);
79 assert(sub_r[m_size - 1] == r[index_to_fullgrid(m_size - 1)]);
80 //------------------
81 }
82
83 //============================================================================
85 T &ff(std::size_t i, std::size_t j) { return m_ff(i, j); }
86 T &fg(std::size_t i, std::size_t j) { return m_fg(i, j); }
87 T &gf(std::size_t i, std::size_t j) { return m_gf(i, j); }
88 T &gg(std::size_t i, std::size_t j) { return m_gg(i, j); }
89 const T ff(std::size_t i, std::size_t j) const { return m_ff(i, j); }
90 const T fg(std::size_t i, std::size_t j) const { return m_fg(i, j); }
91 const T gf(std::size_t i, std::size_t j) const { return m_gf(i, j); }
92 const T gg(std::size_t i, std::size_t j) const { return m_gg(i, j); }
93
95 const LinAlg::Matrix<T> &ff() const { return m_ff; }
96 const LinAlg::Matrix<T> &fg() const { return m_fg; }
97 const LinAlg::Matrix<T> &gf() const { return m_gf; }
98 const LinAlg::Matrix<T> &gg() const { return m_gg; }
99 LinAlg::Matrix<T> &ff() { return m_ff; }
100 LinAlg::Matrix<T> &fg() { return m_fg; }
101 LinAlg::Matrix<T> &gf() { return m_gf; }
102 LinAlg::Matrix<T> &gg() { return m_gg; }
103
104 std::size_t size() const { return m_size; }
105 std::size_t g_size() const { return m_g_size; }
106 bool includes_g() const { return m_g_size == m_size; };
107 std::size_t i0() const { return m_i0; }
108 std::size_t stride() const { return m_stride; }
109
110 //============================================================================
112 void zero() {
113 m_ff.zero();
114 m_fg.zero();
115 m_gf.zero();
116 m_gg.zero();
117 }
118
119 // Makes 1
120 [[deprecated]] void make_identity() {
121 m_ff.make_identity();
122 m_fg.make_identity();
123 m_gf.make_identity();
124 m_gg.make_identity();
125 }
126
127 // only for transition, kill!
128 [[deprecated]] RDMatrix<T> &plusIdent(T a = T{1.0}) {
129 (*this) += a;
130 return *this;
131 }
132
133 //============================================================================
136 m_ff += rhs.m_ff;
137 m_fg += rhs.m_fg;
138 m_gf += rhs.m_gf;
139 m_gg += rhs.m_gg;
140 return *this;
141 }
144 m_ff -= rhs.m_ff;
145 m_fg -= rhs.m_fg;
146 m_gf -= rhs.m_gf;
147 m_gg -= rhs.m_gg;
148 return *this;
149 }
152 m_ff *= x;
153 m_fg *= x;
154 m_gf *= x;
155 m_gg *= x;
156 return *this;
157 }
158
160 [[nodiscard]] friend RDMatrix<T> operator+(RDMatrix<T> lhs,
161 const RDMatrix<T> &rhs) {
162 return (lhs += rhs);
163 }
165 [[nodiscard]] friend RDMatrix<T> operator-(RDMatrix<T> lhs,
166 const RDMatrix<T> &rhs) {
167 return (lhs -= rhs);
168 }
170 [[nodiscard]] friend RDMatrix<T> operator*(const T x, RDMatrix<T> rhs) {
171 return (rhs *= x);
172 }
173
176 m_ff += aI;
177 m_gg += aI;
178 return *this;
179 }
182 m_ff -= aI;
183 m_gg -= aI;
184 return *this;
185 }
186
188 [[nodiscard]] friend RDMatrix<T> operator+(RDMatrix<T> M, T aI) {
189 return (M += aI);
190 }
192 [[nodiscard]] friend RDMatrix<T> operator-(RDMatrix<T> M, T aI) {
193 return (M -= aI);
194 }
195
196 //============================================================================
197
200 [[nodiscard]] friend RDMatrix<T> operator*(const RDMatrix<T> &a,
201 const RDMatrix<T> &b) {
202
203 RDMatrix<T> out(a.m_i0, a.m_stride, a.m_size, a.m_incl_g, a.m_rgrid);
204
205 // FF = FF*FF + FG*GF
206 // FG = FF*FG + FG*GG
207 // GF = GF*FF + GG*GF
208 // GG = GF*FG + GG*GG
209 out.ff() = a.ff() * b.ff();
210 if (a.m_incl_g) {
211 out.ff() += a.fg() * b.gf();
212 out.fg() = a.ff() * b.fg() + a.fg() * b.gg();
213 out.gf() = a.gf() * b.ff() + a.gg() * b.gf();
214 out.gg() = a.gf() * b.fg() + a.gg() * b.gg();
215 }
216 return out;
217 }
218
219 //============================================================================
222 m_ff.mult_elements_by(rhs.ff());
223 m_fg.mult_elements_by(rhs.fg());
224 m_gf.mult_elements_by(rhs.gf());
225 m_gg.mult_elements_by(rhs.gg());
226 return *this;
227 }
229 [[nodiscard]] friend RDMatrix<T> mult_elements(RDMatrix<T> lhs,
230 const RDMatrix<T> &rhs) {
231 lhs.mult_elements_by(rhs);
232 return lhs;
233 }
234
235 //============================================================================
236
238 [[nodiscard]] RDMatrix<T> conj() const {
239 auto out = *this;
240 out.ff().conj_in_place();
241 out.fg().conj_in_place();
242 out.gf().conj_in_place();
243 out.gg().conj_in_place();
244 return out;
245 }
248 [[nodiscard]] RDMatrix<double> real() const {
249 RDMatrix<double> out(m_i0, m_stride, m_size, m_incl_g, m_rgrid);
250 out.ff() = m_ff.real();
251 out.fg() = m_fg.real();
252 out.gf() = m_gf.real();
253 out.gg() = m_gg.real();
254 return out;
255 }
257 [[nodiscard]] RDMatrix<double> imag() const {
258 RDMatrix<double> out(m_i0, m_stride, m_size, m_incl_g, m_rgrid);
259 out.ff() = m_ff.imag();
260 out.fg() = m_fg.imag();
261 out.gf() = m_gf.imag();
262 out.gg() = m_gg.imag();
263 return out;
264 }
267 [[nodiscard]] RDMatrix<std::complex<double>> complex() const {
268 RDMatrix<std::complex<double>> out(m_i0, m_stride, m_size, m_incl_g,
269 m_rgrid);
270 out.ff() = m_ff.complex();
271 out.fg() = m_fg.complex();
272 out.gf() = m_gf.complex();
273 out.gg() = m_gg.complex();
274 return out;
275 }
276
277 //============================================================================
280 m_ff.invert_in_place();
281 if (m_incl_g) {
282 const auto &ai = m_ff; // already inverted
283 const auto &b = m_fg;
284 const auto &c = m_gf;
285 const auto &d = m_gg;
286 const auto cai = c * ai;
287 const auto dmcaib = (d - cai * b).invert_in_place();
288 const auto aib_dmcaib = ai * b * dmcaib;
289 m_ff += aib_dmcaib * cai;
290 m_fg = -1.0 * aib_dmcaib;
291 m_gf = -1.0 * dmcaib * cai;
292 m_gg = dmcaib;
293 }
294 return *this;
295 }
297 [[nodiscard]] RDMatrix<T> inverse() const {
298 auto out = *this; //
299 return out.invert_in_place();
300 }
301
302 //============================================================================
305 const auto dus = m_rgrid->du() * double(m_stride);
306 for (auto i = 0ul; i < m_size; ++i) {
307 for (auto j = 0ul; j < m_size; ++j) {
308 const auto sj = index_to_fullgrid(j);
309 const auto dr = m_rgrid->drdu(sj) * dus;
310 m_ff[i][j] *= dr;
311 }
312 }
313 if (m_incl_g) {
314 for (auto i = 0ul; i < m_size; ++i) {
315 for (auto j = 0ul; j < m_size; ++j) {
316 const auto sj = index_to_fullgrid(j);
317 const auto dr = m_rgrid->drdu(sj) * dus;
318 m_fg[i][j] *= dr;
319 m_gf[i][j] *= dr;
320 m_gg[i][j] *= dr;
321 }
322 }
323 }
324 return *this;
325 }
328 const auto dus = m_rgrid->du() * double(m_stride);
329 for (auto i = 0ul; i < m_size; ++i) {
330 const auto si = index_to_fullgrid(i);
331 const auto dr = m_rgrid->drdu(si) * dus;
332 for (auto j = 0ul; j < m_size; ++j) {
333 m_ff[i][j] *= dr;
334 }
335 }
336 if (m_incl_g) {
337 for (auto i = 0ul; i < m_size; ++i) {
338 const auto si = index_to_fullgrid(i);
339 const auto dr = m_rgrid->drdu(si) * dus;
340 for (auto j = 0ul; j < m_size; ++j) {
341 m_fg[i][j] *= dr;
342 m_gf[i][j] *= dr;
343 m_gg[i][j] *= dr;
344 }
345 }
346 }
347 return *this;
348 }
350 RDMatrix<T> drj() const {
351 auto out = *this;
352 return out.drj_in_place();
353 }
355 RDMatrix<T> dri() const {
356 auto out = *this;
357 return out.dri_in_place();
358 }
359
361 double dr(std::size_t sub_index) const {
362 const auto full_index = index_to_fullgrid(sub_index);
363 return m_rgrid->drdu(full_index) * m_rgrid->du() * double(m_stride);
364 }
365
366 //============================================================================
368 std::size_t index_to_fullgrid(std::size_t i) const {
369 return m_i0 + i * m_stride;
370 }
371
372 //============================================================================
374 void add(const DiracSpinor &ket, const DiracSpinor &bra, T k = T(1.0)) {
375 // Adds (k)*|ket><bra| to G matrix
376 // G_ij = f * Q_i * W_j
377 // Q = Q(1) = ket, W = W(2) = bra
378 // Takes sub-grid into account; ket,bra are on full grid, G on sub-grid
379 for (auto i = 0ul; i < m_size; ++i) {
380 const auto si = index_to_fullgrid(i);
381 for (auto j = 0ul; j < m_size; ++j) {
382 const auto sj = index_to_fullgrid(j);
383 m_ff[i][j] += k * ket.f(si) * bra.f(sj);
384 } // j
385 } // i
386
387 if (m_incl_g) {
388 for (auto i = 0ul; i < m_size; ++i) {
389 const auto si = index_to_fullgrid(i);
390 for (auto j = 0ul; j < m_size; ++j) {
391 const auto sj = index_to_fullgrid(j);
392 // XXX Double check fg/gf right way!
393 m_fg[i][j] += k * ket.f(si) * bra.g(sj);
394 m_gf[i][j] += k * ket.g(si) * bra.f(sj); // symmetric, transpose?
395 m_gg[i][j] += k * ket.g(si) * bra.g(sj);
396 } // j
397 } // i
398 }
399 }
400
401 //============================================================================
405
406 const auto &r = Fn.grid().r();
407 // const auto &drdu = Fn.grid().drdu();
408 // const double s_du = double(m_stride) * Fn.grid().du();
409
410 // include dr?? No, not by default?
411 std::vector<double> f(m_size), g;
412 for (auto i = 0ul; i < m_size; ++i) {
413 for (auto j = 0ul; j < m_size; ++j) {
414 const auto j_f = index_to_fullgrid(j);
415 f[i] += m_ff(i, j) * Fn.f(j_f); // * drdu[j_f] * s_du;
416 }
417 }
418 if (m_incl_g) {
419 g.resize(m_size);
420 for (auto i = 0ul; i < m_size; ++i) {
421 for (auto j = 0ul; j < m_size; ++j) {
422 const auto j_f = index_to_fullgrid(j);
423 f[i] += m_fg(i, j) * Fn.g(j_f); // * drdu[j_f] * s_du;
424 g[i] += (m_gf(i, j) * Fn.f(j_f) + m_gg(i, j) * Fn.g(j_f)); // *
425 // drdu[j_f] * s_du;
426 }
427 }
428 }
429
430 DiracSpinor out = Fn * 0.0;
431 // Interpolate from sub-grid to full grid
432 out.f() = Interpolator::interpolate(sub_r, f, r);
433 if (m_incl_g) {
434 out.g() = Interpolator::interpolate(sub_r, g, r);
435 }
436 return out;
437 }
438
439 //============================================================================
440 // For testing only:
441 friend std::ostream &operator<<(std::ostream &os, const RDMatrix<T> &a) {
442 os << "FF:\n";
443 os << a.m_ff;
444 if (a.m_incl_g) {
445 os << "FG:\n";
446 os << a.m_fg;
447 os << "GF:\n";
448 os << a.m_gf;
449 os << "GG:\n";
450 os << a.m_gg;
451 }
452 return os;
453 }
454};
455
457template <typename T>
458bool equal(const RDMatrix<T> &lhs, const RDMatrix<T> &rhs) {
459 return equal(lhs.ff(), rhs.ff()) && equal(lhs.fg(), rhs.fg()) &&
460 equal(lhs.gf(), rhs.gf()) && equal(lhs.gg(), rhs.gg());
461}
462
464template <typename T>
465double max_element(const RDMatrix<T> &a) {
466 double xff = 0.0, xfg = 0.0, xgf = 0.0, xgg = 0.0;
467 for (auto i = 0ul; i < a.size(); ++i) {
468 for (auto j = 0ul; j < a.size(); ++j) {
469 if (std::abs(a.ff(i, j)) > xff)
470 xff = std::abs(a.ff(i, j));
471 if (a.g_size() != 0) {
472 if (std::abs(a.fg(i, j)) > xfg)
473 xfg = std::abs(a.fg(i, j));
474 if (std::abs(a.gf(i, j)) > xgf)
475 xgf = std::abs(a.gf(i, j));
476 if (std::abs(a.gg(i, j)) > xgg)
477 xgg = std::abs(a.gg(i, j));
478 }
479 }
480 }
481 return std::max({xff, xfg, xgf, xgg});
482}
483
485template <typename T>
486double max_delta(const RDMatrix<T> &a, const RDMatrix<T> &b) {
487 double xff = 0.0, xfg = 0.0, xgf = 0.0, xgg = 0.0;
488 for (auto i = 0ul; i < a.size(); ++i) {
489 for (auto j = 0ul; j < a.size(); ++j) {
490 if (std::abs(a.ff(i, j) - b.ff(i, j)) > xff)
491 xff = std::abs(a.ff(i, j) - b.ff(i, j));
492 if (a.g_size() != 0) {
493 if (std::abs(a.fg(i, j) - b.fg(i, j)) > xfg)
494 xfg = std::abs(a.fg(i, j) - b.fg(i, j));
495 if (std::abs(a.gf(i, j) - b.gf(i, j)) > xgf)
496 xgf = std::abs(a.gf(i, j) - b.gf(i, j));
497 if (std::abs(a.gg(i, j) - b.gg(i, j)) > xgg)
498 xgg = std::abs(a.gg(i, j) - b.gg(i, j));
499 }
500 }
501 }
502 return std::max({xff, xfg, xgf, xgg});
503}
504
507template <typename T>
508double max_epsilon(const RDMatrix<T> &a, const RDMatrix<T> &b) {
509 double xff = 0.0, xfg = 0.0, xgf = 0.0, xgg = 0.0;
510 for (auto i = 0ul; i < a.size(); ++i) {
511 for (auto j = 0ul; j < a.size(); ++j) {
512 const auto eps_ff =
513 std::abs((a.ff(i, j) - b.ff(i, j)) / (a.ff(i, j) + b.ff(i, j)));
514 if (eps_ff > xff)
515 xff = eps_ff;
516 if (a.g_size() != 0) {
517 const auto eps_fg =
518 std::abs((a.fg(i, j) - b.fg(i, j)) / (a.fg(i, j) + b.fg(i, j)));
519 const auto eps_gf =
520 std::abs((a.gf(i, j) - b.gf(i, j)) / (a.gf(i, j) + b.gf(i, j)));
521 const auto eps_gg =
522 std::abs((a.gg(i, j) - b.gg(i, j)) / (a.gg(i, j) + b.gg(i, j)));
523 if (eps_ff > xfg)
524 xfg = eps_fg;
525 if (eps_gf > xgf)
526 xgf = eps_gf;
527 if (eps_gg > xgg)
528 xgg = eps_gg;
529 }
530 }
531 }
532 return std::max({xff, xfg, xgf, xgg});
533}
534
535//==============================================================================
536using GMatrix = RDMatrix<double>;
537using ComplexGMatrix = RDMatrix<std::complex<double>>;
538using ComplexDouble = std::complex<double>;
539
540} // namespace MBPT
Stores radial Dirac spinor: F_nk = (f, g)
Definition DiracSpinor.hpp:41
const std::vector< double > & f() const
Upper (large) radial component function, f.
Definition DiracSpinor.hpp:117
const Grid & grid() const
Resturns a const reference to the radial grid.
Definition DiracSpinor.hpp:107
const std::vector< double > & g() const
Lower (small) radial component function, g.
Definition DiracSpinor.hpp:124
const std::vector< double > & r() const
Grid points, r.
Definition Grid.hpp:75
Matrix class; row-major.
Definition Matrix.hpp:35
Matrix< T > & make_identity()
Constructs a diagonal unit matrix (identity), in place; only for square.
Definition Matrix.ipp:143
Matrix< T > & invert_in_place()
Inverts the matrix, in place. Uses GSL; via LU decomposition. Only works for double/complex<double>.
Definition Matrix.ipp:77
auto complex() const
Converts a real to complex matrix (changes type; returns a complex matrix)
Definition Matrix.ipp:205
auto imag() const
Returns imag part of complex matrix (changes type; returns a real matrix)
Definition Matrix.ipp:194
Matrix< T > & zero()
Sets all elements to zero, in place.
Definition Matrix.ipp:154
Matrix< T > & mult_elements_by(const Matrix< T > &a)
Muplitplies all the elements by those of matrix a, in place: M_ij *= a_ij.
Definition Matrix.ipp:270
auto real() const
Returns real part of complex matrix (changes type; returns a real matrix)
Definition Matrix.ipp:183
Definition RDMatrix.hpp:44
std::size_t index_to_fullgrid(std::size_t i) const
Converts an index on the sub-grid to the full grid.
Definition RDMatrix.hpp:368
RDMatrix< T > & operator+=(T aI)
Adition of identity: Matrix<T> += T : T assumed to be *Identity!
Definition RDMatrix.hpp:175
RDMatrix< std::complex< double > > complex() const
Converts a real to complex matrix (changes type; returns a complex matrix)
Definition RDMatrix.hpp:267
friend RDMatrix< T > operator+(RDMatrix< T > lhs, const RDMatrix< T > &rhs)
Matrix adition +,-.
Definition RDMatrix.hpp:160
RDMatrix< T > inverse() const
Returns inverse of matrix; original matrix unchanged.
Definition RDMatrix.hpp:297
double dr(std::size_t sub_index) const
returns dr at position along sub grid
Definition RDMatrix.hpp:361
friend RDMatrix< T > operator+(RDMatrix< T > M, T aI)
Adition of identity: Matrix<T> + T : T assumed to be *Identity!
Definition RDMatrix.hpp:188
friend RDMatrix< T > operator*(const RDMatrix< T > &a, const RDMatrix< T > &b)
Matrix multplication: C=A*B := Cij = \sum_k Aik*Bkj. Note: integration measure not included: call ....
Definition RDMatrix.hpp:200
RDMatrix< T > & invert_in_place()
Inversion (in place)
Definition RDMatrix.hpp:279
T & ff(std::size_t i, std::size_t j)
direct access to matrix elements
Definition RDMatrix.hpp:85
RDMatrix< T > & operator*=(const T x)
Scalar multiplication.
Definition RDMatrix.hpp:151
void zero()
Sets all matrix elements to zero.
Definition RDMatrix.hpp:112
friend RDMatrix< T > mult_elements(RDMatrix< T > lhs, const RDMatrix< T > &rhs)
Multiply elements (new matrix): Gij = Aij*Bij.
Definition RDMatrix.hpp:229
friend RDMatrix< T > operator*(const T x, RDMatrix< T > rhs)
Scalar multiplication.
Definition RDMatrix.hpp:170
RDMatrix< T > & operator-=(const RDMatrix< T > &rhs)
Matrix adition +,-.
Definition RDMatrix.hpp:143
friend RDMatrix< T > operator-(RDMatrix< T > lhs, const RDMatrix< T > &rhs)
Matrix adition +,-.
Definition RDMatrix.hpp:165
RDMatrix< T > & drj_in_place()
Multiplies by drj: Q_ij -> Q_ij*dr_j, in place.
Definition RDMatrix.hpp:304
RDMatrix< double > imag() const
Returns imag part of complex matrix (changes type; returns a real matrix)
Definition RDMatrix.hpp:257
friend RDMatrix< T > operator-(RDMatrix< T > M, T aI)
Adition of identity: Matrix<T> - T : T assumed to be *Identity!
Definition RDMatrix.hpp:192
void add(const DiracSpinor &ket, const DiracSpinor &bra, T k=T(1.0))
Adds k*|ket><bra| to matrix (used for building Green's functions)
Definition RDMatrix.hpp:374
RDMatrix< T > & dri_in_place()
Multiplies by dri: Q_ij -> Q_ij*dr_i, in place.
Definition RDMatrix.hpp:327
RDMatrix< double > real() const
Returns real part of complex matrix (changes type; returns a real matrix)
Definition RDMatrix.hpp:248
RDMatrix< T > dri() const
Multiplies by dri: Q_ij -> Q_ij*dr_i. Returns new matrix (orig unchanged)
Definition RDMatrix.hpp:355
RDMatrix< T > drj() const
Multiplies by drj: Q_ij -> Q_ij*dr_j. Returns new matrix (orig unchanged)
Definition RDMatrix.hpp:350
const LinAlg::Matrix< T > & ff() const
direct access to matrix's
Definition RDMatrix.hpp:95
RDMatrix< T > & mult_elements_by(const RDMatrix< T > &rhs)
Multiply elements (in place): Gij -> Gij*Bij.
Definition RDMatrix.hpp:221
DiracSpinor operator*(const DiracSpinor &Fn) const
Action of RDMatrix operator on DiracSpinor. Inludes Integration: G*F = Int[G(r,r')*F(r') dr'] = sum_j...
Definition RDMatrix.hpp:404
RDMatrix< T > & operator+=(const RDMatrix< T > &rhs)
Matrix adition +,-.
Definition RDMatrix.hpp:135
RDMatrix< T > conj() const
Returns conjugate of matrix.
Definition RDMatrix.hpp:238
RDMatrix< T > & operator-=(T aI)
Adition of identity: Matrix<T> -= T : T assumed to be *Identity!
Definition RDMatrix.hpp:181
std::vector< double > interpolate(const std::vector< double > &x_in, const std::vector< double > &y_in, const std::vector< double > &x_out, Method method=Method::cspline)
Performs interpolation using GSL (GNU Scientific Library)
Definition Interpolator.hpp:162
Many-body perturbation theory.
Definition CI_Integrals.hpp:9
bool equal(const RDMatrix< T > &lhs, const RDMatrix< T > &rhs)
Checks if two matrix's are equal (to within parts in 10^12)
Definition RDMatrix.hpp:458
double max_epsilon(const RDMatrix< T > &a, const RDMatrix< T > &b)
returns maximum relative diference [aij-bij/(aij+bij)] (abs) between two matrices
Definition RDMatrix.hpp:508
double max_delta(const RDMatrix< T > &a, const RDMatrix< T > &b)
returns maximum difference (abs) between two matrixs
Definition RDMatrix.hpp:486
double max_element(const RDMatrix< T > &a)
returns maximum element (by abs)
Definition RDMatrix.hpp:465