Program Listing for File matrix.hpp

↰ Return to documentation for file (src_util/matrix.hpp)

#ifndef matrix_h
#define matrix_h

#include <iostream>

#include <complex>
#include <cstring>
#include <ctype.h>

#include "parser.hpp"
#include "vector_num.hpp"

using namespace std;

#define matrix_SMALL_VALUE 1e-10
//#define BOUND_CHECK

template <typename T> struct matrix {
  size_t r;
  size_t c;
  vector<T> v;

  auto data() ->decltype(v.data()) { return v.data(); }
  auto size() ->decltype(v.size()) { return v.size(); }

  matrix() : r(0), c(0) {}

  matrix(size_t _r) : r(_r), c(_r) { v.resize(r * c); }

  matrix(size_t _r, size_t _c) : r(_r), c(_c) { v.resize(r * c); }

  matrix(const size_t _r, const size_t _c, const vector<T> &_v) : r(_r), c(_c), v(_v) {}

  matrix(const size_t _r, const vector<T> &_v) : r(_r), c(_r), v(_v) {}

  matrix(const matrix<T> &A) : r(A.r), c(A.c), v(A.v) {}

  matrix(const matrix<T> &A, const matrix<T> &B) : r(A.r), c(A.c + B.c) {
    v.resize(A.v.size() + B.v.size());
    copy(A.v.data(), A.v.data() + A.v.size(), v.data());
    copy(B.v.data(), B.v.data() + B.v.size(), &v[A.v.size()]);
  }


  template <typename S> matrix<T> &assign(const matrix<S> &A) {
    if ((void *)this == (void *)&A)
      return *this;
    else if (v.size() == 0) {
      set_size(A.r, A.c);
      for (size_t i = 0; i < A.v.size(); ++i) v[i] = A.v[i];
    } else if (r == A.r and c == A.c) {
      for (size_t i = 0; i < A.v.size(); ++i) v[i] = A.v[i];
    } else if (v.size() > A.v.size()) {
      for (size_t i = 0; i < A.r; ++i) {
        for (size_t j = 0; j < A.c; ++j) { v[i + r * j] = A.v[i + A.r * j]; }
      }
    } else {
      for (size_t i = 0; i < r; ++i) {
        for (size_t j = 0; j < c; ++j) { v[i + r * j] = A.v[i + A.r * j]; }
      }
    }
    return *this;
  }



  inline const T &operator()(const size_t &i, const size_t &j) const { return (v[i + r * j]); }

  T &operator()(const size_t &i, const size_t &j) { return (v[i + r * j]); }

  void concatenate(const matrix<T> &A, const matrix<T> &B) {
    set_size(A.r, A.c + B.c);
    copy(&A.v[0], &A.v[0] + A.v.size(), &v[0]);
    copy(&B.v[0], &B.v[0] + B.v.size(), &v[A.v.size()]);
  }


  void sub_matrix(size_t i, matrix<T> &A) {
    QCM_ASSERT(r == c and A.r == A.c and i + A.r <= r);
    for (size_t k = 0; k < A.r; ++k) {
      for (size_t l = 0; l < A.r; ++l) { A(k, l) = v[i + k + r * (i + l)]; }
    }
  }



  template <typename S>
  void move_sub_matrix(size_t R, size_t C, size_t i, size_t j, size_t I, size_t J, matrix<S> &A, double z = 1.0) const {
    QCM_ASSERT(j + C <= c and i + R <= r and I + R <= A.r and J + C <= A.c);
    for (size_t k = 0; k < R; ++k) {
      for (size_t l = 0; l < C; ++l) { A(I + k, J + l) += z * v[i + k + r * (j + l)]; }
    }
  }

  template <typename S>
  void move_sub_matrix_conjugate(size_t R, size_t C, size_t i, size_t j, size_t I, size_t J, matrix<S> &A, double z = 1.0) const {
    QCM_ASSERT(j + C <= c and i + R <= r and I + R <= A.r and J + C <= A.c);
    for (size_t k = 0; k < R; ++k) {
      for (size_t l = 0; l < C; ++l) { A(I + k, J + l) += z * conjugate(v[i + k + r * (j + l)]); }
    }
  }


  template <typename S>
  void move_sub_matrix_HC(size_t R, size_t C, size_t i, size_t j, size_t I, size_t J, matrix<S> &A, double z = 1.0) const {
    QCM_ASSERT(j + C <= c and i + R <= r and I + C <= A.r and J + R <= A.c);
    for (size_t k = 0; k < R; ++k) {
      for (size_t l = 0; l < C; ++l) { A(I + l, J + k) += z * conjugate(v[i + k + r * (j + l)]); }
    }
  }


  template <typename S>
  void move_sub_matrix_transpose(size_t R, size_t C, size_t i, size_t j, size_t I, size_t J, matrix<S> &A, double z = 1.0) const {
    QCM_ASSERT(j + C <= c and i + R <= r and I + C <= A.r and J + R <= A.c);
    for (size_t k = 0; k < R; ++k) {
      for (size_t l = 0; l < C; ++l) { A(I + l, J + k) += z * v[i + k + r * (j + l)]; }
    }
  }



  void insert_row(size_t i, const vector<T> &A) {
    for (size_t j = 0; j < c; j++) v[i + r * j] = A[j];
  }



  void insert_column(size_t i, const vector<T> &A) {
    for (size_t j = 0; j < r; j++) v[j + r * i] = A[j];
  }



  void extract_row(size_t i, vector<T> &A) {
    for (size_t j = 0; j < c; j++) A[j] = v[i + r * j];
  }

  vector<T> extract_row(size_t i) {
    vector<T> A(c);
    for (size_t j = 0; j < c; j++) A[j] = v[i + r * j];
    return A;
  }


  void extract_column(size_t i, vector<T> &A) { copy(&v[i * r], &v[i * r] + r, &A[0]); }

  vector<T> extract_column(size_t i) {
    vector<T> A(r);
    copy(&v[i * r], &v[i * r] + r, &A[0]);
    return A;
  }


  inline void zero() { to_zero(v); }


  void apply_mask(const vector<size_t> &mask) {
    for (auto &i : mask) v[i] = 0.0;
  }


  inline void set_size(size_t _r, size_t _c) {
    r = _r;
    c = _c;
    v.clear();
    v.resize((size_t)(r * c));
  }
  inline void set_size(size_t _r) { set_size(_r, _r); }
  template <typename S> inline void set_size(matrix<S> mat) { set_size(mat.r, mat.c); }


  inline void import(T *x) { copy(x, v.size(), &v[0]); }



  void real_part(const matrix<Complex> &A) {
    for (size_t i = 0; i < v.size(); i++) v[i] = real(A.v[i]);
  }


  void identity() {
    zero();
    for (size_t i = 0; i < r; i++) v[i + i * r] = T(1.0);
  }


  template <typename S> inline void diagonal(const vector<S> &x) {
    for (size_t i = 0; i < r; ++i) v[i + i * r] = x[i];
  }


  template <typename S> inline void add_to_diagonal(const S x) {
    for (size_t i = 0; i < r; ++i) v[i + i * r] += x;
  }

  template <typename S> inline void add_to_diagonal(const vector<S> &x) {
    for (size_t i = 0; i < r; ++i) v[i + i * r] += x[i];
  }


  void transpose() {
    vector<T> tmp(v);
    swap(r, c);
    for (size_t i = 0; i < r; ++i)
      for (int j = 0; j < c; ++j) v[i + r * j] = tmp[j + c * i];
  }


  void hermitian_conjugate() {
    vector<T> tmp(v);
    swap(r, c);
    for (size_t i = 0; i < r; ++i)
      for (size_t j = 0; j < c; ++j) v[i + r * j] = conjugate(tmp[j + c * i]);
  }


  T trace() {
    T z(0.0);
    for (size_t i = 0; i < r; i++) z += v[i + r * i];
    return z;
  }


  T trace(size_t start_row, size_t end_row) {
    T z(0.0);
    for (size_t i = start_row; i < end_row; i++) z += v[i + r * i];
    return z;
  }


  template <typename S> T trace_product(matrix<S> A, bool transpose = false) {
    T z(0.0);
    if (transpose) {
      QCM_ASSERT(A.r == r and A.c == c);
      for (size_t i = 0; i < r; i++) {
        for (size_t j = 0; j < c; j++) { z += v[i + r * j] * A.v[i + r * j]; }
      }
    } else {
      QCM_ASSERT(A.r == c and A.c == r);
      for (size_t i = 0; i < r; i++) {
        for (size_t j = 0; j < c; j++) { z += v[i + r * j] * A.v[j + c * i]; }
      }
    }
    return z;
  }



  template <typename S> matrix<T> &operator-=(const S &a) {
    for (size_t i = 0; i < r; i++) v[i + r * i] -= a;
    return *this;
  }



  template <typename S> matrix<T> &operator-=(const matrix<S> &A) {
    for (size_t i = 0; i < v.size(); i++) v[i] = v[i] - A.v[i];
    return *this;
  }

  template <typename S> matrix<T> &operator+=(const S &a) {
    for (size_t i = 0; i < r; i++) v[i + r * i] += a;
    return *this;
  }



  template <typename S> matrix<T> &operator+=(const matrix<S> &A) {
    for (size_t i = 0; i < v.size(); i++) v[i] = v[i] + A.v[i];
    return *this;
  }



  double diff_sq(const matrix<T> &A, double a = 1.0) {
#ifdef BOUND_CHECK
    QCM_ASSERT(r == A.r and c == A.c);
#endif
    double z = 0.0;
    for (size_t i = 0; i < v.size(); ++i) {
      double zz = abs(v[i] - A.v[i] * a);
      z += zz * zz;
    }
    return (z);
  }



  template <typename S, typename U> void product(const matrix<S> &A, const matrix<U> &B) {
#ifdef BOUND_CHECK
    QCM_ASSERT(A.r == r and B.c == c and A.c == B.r);
#endif
    QCM_ASSERT(&A != this and &B != this);
    for (size_t j = 0; j < c; j++) {
      size_t jj = j * B.r;
      for (size_t i = 0; i < r; ++i) {
        T z = 0.0;
        for (size_t k = 0; k < A.c; ++k) z += A.v[i + k * r] * B.v[k + jj];
        v[i + j * r] = z;
      }
    }
  }


  template <typename S, typename U> void product(const matrix<S> &A, const matrix<U> &B, size_t M) {
#ifdef BOUND_CHECK
    QCM_ASSERT(A.r == r and B.c == c and A.c == B.r);
#endif
    if(M>A.c) M = A.c;
    QCM_ASSERT(&A != this and &B != this);
    for (size_t j = 0; j < c; j++) {
      size_t jj = j * B.r;
      for (size_t i = 0; i < r; ++i) {
        T z = 0.0;
        for (size_t k = 0; k < M; ++k) z += A.v[i + k * r] * B.v[k + jj];
        v[i + j * r] = z;
      }
    }
  }





  template <typename S, typename U> void product_transpose(const matrix<S> &A, const matrix<U> &B) {
#ifdef BOUND_CHECK
    QCM_ASSERT(A.r == r and B.r == c and A.c == B.c);
#endif
    QCM_ASSERT(&A != this and &B != this);
    for (size_t j = 0; j < c; j++) {
      for (size_t i = 0; i < r; ++i) {
        T z = 0.0;
        for (size_t k = 0; k < A.c; ++k) z += A.v[i + k * r] * B.v[j + k * r];
        v[i + j * r] = z;
      }
    }
  }



  template <typename S, typename U> void product_hermitian_conjugate(const matrix<S> &A, const matrix<U> &B) {
#ifdef BOUND_CHECK
    QCM_ASSERT(A.r == r and B.r == c and A.c == B.c);
#endif
    QCM_ASSERT(&A != this and &B != this);
    for (size_t j = 0; j < c; j++) {
      for (size_t i = 0; i < r; ++i) {
        T z = 0.0;
        for (size_t k = 0; k < A.c; ++k) z += A.v[i + k * r] * conjugate(B.v[j + k * r]);
        v[i + j * r] = z;
      }
    }
  }



  void add(T d) {
    for (size_t i = 0; i < r; i++) v[i + i * r] += d;
  }

  template <typename S, typename B> void add(const matrix<S> &A, B d) {
    for (size_t i = 0; i < v.size(); i++) v[i] += d * A.v[i];
  }

  template <typename S> void simil(const matrix<S> &A) {
#ifdef BOUND_CHECK
    QCM_ASSERT(A.r == c and r == c);
#endif
    T *w;
    w = new T[r * c];
    memcpy(w, v.data(), r * c * sizeof(T));

    for (size_t j = 0; j < c; j++) {
      for (size_t i = 0; i < r; ++i) {
        T z(0);
        for (size_t k = 0; k < c; ++k) {
          T zz(0);
          for (size_t l = 0; l < c; ++l) zz += A(l, j) * w[k + r * l];
          z += zz * conjugate(A(k, i));
        }
        v[i + j * r] = z;
      }
    }
    delete[] w;
  }

  template <typename S> void simil_inv(const matrix<S> &A) {
#ifdef BOUND_CHECK
    QCM_ASSERT(A.c == r and r == c);
#endif
    T *w;
    w = new T[r * c];
    memcpy(w, v.data(), r * c * sizeof(T));

    for (size_t j = 0; j < c; j++) {
      for (size_t i = 0; i < r; ++i) {
        T z(0);
        for (size_t k = 0; k < c; ++k) {
          T zz(0);
          for (size_t l = 0; l < c; ++l) zz += conjugate(A(j, l)) * w[k + r * l];
          z += zz * A(i, k);
        }
        v[i + j * r] = z;
      }
    }
    delete[] w;
  }

  template <typename S> void simil_transpose(const matrix<S> &A) {
#ifdef BOUND_CHECK
    QCM_ASSERT(A.r == r and A.c == c);
#endif
    T *w;
    w = new T[r * c];
    memcpy(w, v.v, r * c * sizeof(T));

    for (size_t j = 0; j < c; j++) {
      for (size_t i = 0; i < r; ++i) {
        T z(0);
        for (size_t k = 0; k < c; ++k) {
          T zz(0);
          for (size_t l = 0; l < c; ++l) zz += A(l, j) * w[k + r * l];
          z += zz * A(k, i);
        }
        v[i + j * r] = z;
      }
    }
    delete[] w;
  }

  template <typename S> void simil_conjugate(const matrix<S> &A) {
#ifdef BOUND_CHECK
    QCM_ASSERT(A.r == r and A.c == c);
#endif
    T *w;
    w = new T[r * c];
    memcpy(w, v.v, r * c * sizeof(T));

    for (size_t j = 0; j < c; j++) {
      for (size_t i = 0; i < r; ++i) {
        T z(0);
        for (size_t k = 0; k < c; ++k) {
          T zz(0);
          for (size_t l = 0; l < c; ++l) zz += conjugate(A(l, j)) * w[k + r * l];
          z += zz * conjugate(A(k, i));
        }
        v[i + j * r] = z;
      }
    }
    delete[] w;
  }

  template <typename S1, typename S2> void simil(const matrix<S1> &A, const matrix<S2> &B) {
#ifdef BOUND_CHECK
    QCM_ASSERT(A.r == B.r and B.c == B.r and A.c == r);
#endif
    for (size_t j = 0; j < c; j++) {
      for (size_t i = 0; i < r; ++i) {
        T z = 0.0;
        for (size_t k = 0; k < B.c; ++k) {
          T zz = 0.0;
          for (size_t l = 0; l < B.c; ++l) zz += A(l, j) * B(k, l);
          z += zz * conjugate(A(k, i));
        }
        v[i + j * r] += z;
      }
    }
  }


  void apply_add(T *x, T *y) {
    for (size_t j = 0; j < c; ++j) {
      for (size_t i = 0; i < r; ++i) y[i] += v[i + j * r] * x[j];
    }
  }

  void apply_add(const vector<T> &x, vector<T> &y) const {
    for (size_t j = 0; j < c; ++j) {
      for (size_t i = 0; i < r; ++i) y[i] += v[i + j * r] * x[j];
    }
  }


  void apply(T *x, T *y) {
    memset(y, 0, r * sizeof(*y));
    apply_add(x, y);
  }

  void apply(const vector<T> &x, vector<T> &y) {
    to_zero(y);
    apply_add(x, y);
  }


  void left_apply_add(const vector<T> &x, vector<T> &y) const {
    for (size_t j = 0; j < c; ++j) {
      for (size_t i = 0; i < r; ++i) y[j] += v[i + j * r] * conjugate(x[i]);
    }
  }


  bool is_diagonal(double accuracy = 1.0e-8) {
    QCM_ASSERT(r == c);
    for (size_t i = 0; i < r; ++i) {
      for (size_t j = 0; j < c; ++j) {
        if (abs(v[i + r * j]) > accuracy and i != j) return false;
      }
    }
    return true;
  }

  inline void cconjugate() { conjugate(v); }


  T qform(vector<T> &x, vector<T> &y) {
    vector<T> w(r);
    apply(y, w);
    return x * w;
  }

  bool printable() {
    size_t max_dim_print = global_int("max_dim_print");
    if (r > max_dim_print or c > max_dim_print) return false;
    return true;
  }

  bool is_unitary(double threshold) {
    if (r != c) return false;
    for (size_t i = 0; i < r; ++i) {
      for (size_t j = 0; j <= i; ++j) {
        T z(0);
        for (size_t k = 0; k < r; ++k) { z += conjugate(v[k + i * r]) * v[k + j * r]; }
        if (i == j and abs(z - 1.0) > threshold) return false;
        if (i != j and abs(z) > threshold) return false;
      }
    }
    return true;
  }



  bool is_hermitian(double threshold=1e-10) {
    if (r != c) return false;
    for (size_t i = 0; i < r; ++i) {
      if (imag(v[i + i * r]) > threshold) return false;
      for (size_t j = 0; j < i; ++j) {
        if (abs(v[i + j * r] - conjugate(v[j + i * r])) > threshold) return false;
      }
    }
    return true;
  }



  bool Gram_Schmidt(const vector<T> &x, size_t col) {
    QCM_ASSERT(r == c and x.size() == r);

    vector<vector<T>> A(r);
    A[0] = x;
    if (!normalize(A[0])) return false;
    insert_column(col, A[0]);

    size_t k = 1;
    for (size_t i = 0; i < r; i++) {
      if (i == col) continue;
      A[k].resize(r);
      A[k][i] = 1.0;
      for (size_t j = 0; j < k; j++) {
        T z = A[j][i]; // A[j]*A[k]
        mult_add(-z, A[j], A[k]);
      }
      if (!normalize(A[k])) return false;
      insert_column(i, A[k]);
      k++;
    }
    if (!is_unitary(1.0e-6)) return false;
    return true;
  }



  void Cayley(const matrix<T> &A) {
    zero();
    QCM_ASSERT(r == c and A.r == A.c and r == A.r);
    matrix<T> X(A);
    X.add(-1.0);
    X.v *= -1.0;
    X.inverse();
    matrix<T> Y(A);
    Y.add(1.0);
    product(Y, X);
  }



  void Cayley_inverse(const matrix<T> &M) {
    zero();
    QCM_ASSERT(r == c and M.r == M.c and r == M.r);
    matrix<T> X(M);
    X.add(1.0);
    X.inverse();
    matrix<T> Y(M);
    Y.add(-1.0);
    product(X, Y);
  }



  void antisymmetrize() {
    for (size_t i = 0; i < r; i++) {
      for (size_t j = 0; j < i; j++) {
        (*this)(i, j) = ((*this)(i, j) - conjugate((*this)(j, i))) * 0.5;
        (*this)(j, i) = -conjugate((*this)(i, j));
      }
      T z = (*this)(i, i);
      (*this)(i, i) = 0.5 * (z - conjugate(z));
    }
  }


  void symmetrize() {
    for (size_t i = 0; i < r; i++) {
      for (size_t j = 0; j < i; j++) {
        (*this)(i, j) = ((*this)(i, j) + conjugate((*this)(j, i))) * 0.5;
        (*this)(j, i) = conjugate((*this)(i, j));
      }
      T z = (*this)(i, i);
      (*this)(i, i) = 0.5 * (z + conjugate(z));
    }
  }



  friend std::istream &operator>>(std::istream &flux, matrix<T> &A) {
    for (size_t i = 0; i < A.r; i++) {
      for (size_t j = 0; j < A.c; j++) flux >> A(i, j);
    }
    return flux;
  }



  friend std::ostream &operator<<(std::ostream &flux, const matrix<T> &A) {
    size_t max_dim_print = global_int("max_dim_print");
    if (A.r > max_dim_print or A.c > max_dim_print) return flux;
    for (size_t i = 0; i < A.r; ++i) {
      for (size_t j = 0; j < A.c; ++j) flux << chop(A(i, j)) << '\t';
      flux << '\n';
    }
    flux << "\n";
    return flux;
  }


  // below are declarations of routines that are type specific
  double norm();
  double norm2();
  void inverse();
  T determinant();
  void eigensystem(vector<double> &d, matrix<T> &U) const; // finds the eigenvalues and eigenvectors of real symmetric matrix
  void eigenvalues(vector<double> &d);                     // finds the eigenvalues of hermitian matrix
  bool is_real(double accuracy = 1e-6);
  bool is_orthogonal(double accuracy = 1e-6);
  bool cholesky(matrix<double> &A);
  bool triangular_inverse();


  //------------------------------------------------------------------------------
};

matrix<Complex> to_complex_matrix(const matrix<double> &x);
inline matrix<complex<double>> to_complex_matrix(const matrix<complex<double>> &v)
{
  return v;
}

inline matrix<double> to_real_matrix(const matrix<double> &v)
{
  return v;
}
inline matrix<double> to_real_matrix(const matrix<complex<double>> &v)
{
  matrix<double> x(v.r, v.c);
  for(size_t k = 0; k < v.v.size(); ++k) x.v[k] = v.v[k].real();
  return x;
}

matrix<Complex> hermitian_matrix_from_real_vector(size_t d, const vector<double> &x);
void hermitian_matrix_to_real_vector(const matrix<Complex> &M, double *x);

#endif