1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132

use crate::number::{c64, Number};
use rayon::prelude::*;
use std::error::Error;

pub mod bd;
pub mod ci;
pub mod di;
pub mod ge;
pub mod gt;
pub mod kr;
pub mod operations;
pub mod operators;
pub mod po;
pub mod pt;
pub mod sp;
pub mod st;
pub mod sy;
pub mod to;
pub mod tr;

/// # Matrix
#[derive(Clone, Debug, Default, Hash, PartialEq)]
pub struct Matrix<T = f64>
where
    T: Number,
{
    rows: usize,
    cols: usize,
    elems: Vec<T>,
}

#[derive(thiserror::Error, Debug)]
pub enum MatrixError {
    #[error("Dimension mismatch.")]
    DimensionMismatch,
    #[error("BLAS routine error. routine: {routine}, info: {info}")]
    BlasRoutineError { routine: String, info: i32 },
    #[error("LAPACK routine error. routine: {routine}, info: {info}")]
    LapackRoutineError { routine: String, info: i32 },
    #[error("Others")]
    Others(Box<dyn Error + Send + Sync>),
}

impl From<Box<dyn Error + Send + Sync>> for MatrixError {
    fn from(e: Box<dyn Error + Send + Sync>) -> Self {
        MatrixError::Others(e)
    }
}

impl<T> Matrix<T>
where
    T: Number,
{
    pub fn new(rows: usize, cols: usize) -> Self {
        Self {
            rows,
            cols,
            elems: vec![T::default(); rows * cols],
        }
    }

    pub fn from(rows: usize, elems: Vec<T>) -> Self {
        Self {
            rows,
            cols: elems.len() / rows,
            elems,
        }
    }

    pub fn same_size(&self, rhs: &Matrix<T>) -> bool {
        self.rows == rhs.rows && self.cols == rhs.cols
    }

    pub fn rows(&self) -> usize {
        self.rows
    }

    pub fn cols(&self) -> usize {
        self.cols
    }

    pub fn vec(self) -> Vec<T> {
        self.elems
    }

    pub fn slice(&self) -> &[T] {
        &self.elems
    }

    pub fn reshape(mut self, rows: usize) -> Self {
        self.rows = rows;
        self.cols = self.elems.len() / rows;

        self
    }
}

impl Matrix<f64> {
    pub fn to_complex(&self) -> Matrix<c64> {
        Matrix::<c64>::from(
            self.rows,
            self.elems.par_iter().map(|&e| c64::new(e, 0.0)).collect(),
        )
    }
}

impl Matrix<c64> {
    pub fn to_real(&self) -> Matrix<f64> {
        Matrix::from(self.rows, self.elems.par_iter().map(|e| e.re).collect())
    }
}

pub trait Vector<T>
where
    T: Number,
{
    fn row_mat(self) -> Matrix<T>;
    fn col_mat(self) -> Matrix<T>;
}

impl<T> Vector<T> for Vec<T>
where
    T: Number,
{
    fn row_mat(self) -> Matrix<T> {
        Matrix::<T>::from(1, self)
    }

    fn col_mat(self) -> Matrix<T> {
        Matrix::<T>::from(self.len(), self)
    }
}