matrix.h

#pragma once
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>
 
#include "complex_number.h"
#include "vector.h"
 
template<typename T>
class Matrix
{
private:
  std::vector<std::vector<T>>
      data;  
  int rows;  
  int cols;  
 
public:
  Matrix()
      : rows(0)
      , cols(0)
      , data {} {};
 
  Matrix(int r, int c)
      : rows(r)
      , cols(c)
      , data(r, std::vector<T>(c))
  {
  }
 
  Matrix(const std::vector<std::vector<T>>& mat);
 
  T& operator()(int row, int column);
 
  const T& operator()(int row, int column) const;
 
  std::vector<T>& operator[](int i) { return data[i]; }
 
  const std::vector<T>& operator[](int i) const { return data[i]; }
 
  // ==================== Basic Operations ====================
 
  Matrix<T> operator+(const Matrix<T>& other) const;
 
  Matrix<T> operator-(const Matrix<T>& other) const;
 
  Matrix<T> operator*(const Matrix<T>& other) const;
 
  Matrix<T> operator*(const T& scalar) const;
 
  Matrix<T> operator/=(const Matrix<T>& other) const;
 
  Matrix<T> operator-=(const Matrix<T>& other) const;
 
  friend std::ostream& operator<<(std::ostream& os, const Matrix<T>& m)
  {
    for (int i = 0; i < m.rows; ++i) {
      for (int j = 0; j < m.cols; ++j) {
        os << std::setw(12) << std::fixed << std::setprecision(6)
           << m.data[i][j];
      }
      if (i < m.rows - 1) {
        os << "\n";
      }
    }
    return os;
  }
 
  // ==================== Advanced Operations ====================
  T determinant() const;
 
  Matrix<T> inverse() const;
 
  Matrix<T> transpose() const;
 
  std::vector<T> eigenvalues() const;
 
  Matrix<T> eigenvectors() const;
 
  Matrix<T> adjoint() const;
 
  Matrix<T> submatrix(int row_start,
                      int row_end,
                      int col_start,
                      int col_end) const;
 
  T rank() const;
 
  bool isSymmetric() const;
 
  // ==================== Decompositions ====================
 
  std::pair<Matrix<T>, Matrix<T>> LUDecomposition() const;
 
  std::pair<Matrix<T>, Matrix<T>> QRDecomposition() const;
 
  // ==================== Vector Operations ====================
 
  std::vector<T> matVec(const std::vector<T>& v) const;
 
  T dot(const std::vector<T>& a, const std::vector<T>& b) const;
 
  void normalize(std::vector<T>& v) const;
 
  void lanczos(int m,
               std::vector<T>& alpha,
               std::vector<T>& beta,
               Matrix<T>& Q) const;
 
  Matrix<T> buildTridiagonal(const std::vector<T>& alpha,
                             const std::vector<T>& beta,
                             int m) const;
 
  // ==================== Matrix Norms ====================
 
  double frobeniusNorm() const;
 
  double spectralNorm() const;
 
  double norm(const std::vector<T>& v) const;
 
  // ==================== Utility ====================
 
  Matrix<T> getMinor(int row, int col) const;
 
  std::string toString() const;
 
  int getRows() const { return rows; }
 
  int getCols() const { return cols; }
 
  // ==================== Static Factory Methods ====================
 
  static Matrix<T> identity(int n);
 
  static Matrix<T> zeros(int r, int c);
 
  static Matrix<T> ones(int r, int c);
};