golang/c++ Matrix | 青训营笔记

近日，作为三年级fw，参加了计算机能力测试，结果令人羞愧。某非计算机人士特此连夜观摩了线性代数和矩阵编写。痛哭涕泠(T_T)，不以言表。

由$m \times n$个数 $a_{ij}(i=1,2,\cdots,m;j=1,2,\cdots,n)$排成的$m$行$n$列的矩形数表:

称为一个$m \times n$矩阵.这$m\times n$个数称为该矩阵的元素, 其中元素$a_{ij}$位于该矩阵的第$i$行、第$j$列，称为该矩阵的$(i,j)$元.

package matrix

import (
    "math"
)

type SQ struct {
    //矩阵结构
    N, M int   // m是列数，n是行数
    Data [][]float64
}

// 矩阵定义
func (this *SQ) Set(m int, n int, data []float64) {
    // m是列数，n是行数，data是矩阵数据(由左到右由上到下填充)
    this.M = m
    this.N = n
    if len(data) != this.M * this.N {
        // 矩阵定义失败
        return
    }
    else
    {
        k := 0
        if this.M * this.N == len(data) {
            for i := 0; i < this.N; ++i {
                var tmpArr []float64
                for j := 0; j < this.M; ++j {
                    tmpArr = append(tmpArr, data\[k])
                    ++k
                }
                this.Data = append(this.Data, tmpArr)
            }
        }
        else
        {
            // ("矩阵定义失败")
            return
        }
    }
}

// a的列数和b的行数相等
// 矩阵乘法
func Mul(a SQ, b SQ) [][]float64 {
    if a.M == b.N {
        res := [][]float64{}
        // 以A行乘以B列
        for i := 0; i < a.N; ++i {
            t := []float64{}
            // a.data[i] 是第 i+1 行
            for j := 0; j < b.M; ++j {
                r := []float64(0)
                for k := 0; k < a.M; ++k {
                    r += a.Data[i][k] * b.Data[k][j]
                }
                t = append(t, r)
            }
            res = append(res, t)
        }
        return res
    }
    else
    {
        // 两矩阵无法进行相乘运算
        return [][]float64{}
    }
}

// 矩阵转置
func Transpose(a SQ) SQ {
    var b SQ
    b.M = a.N
    b.N = a.M
    var res = [][]float64{}

    for i := 0; i < a.M; ++i {
    	var t = []float64{}
    	for j := 0; j < a.N; ++j {
    		t = append(t, a.Data[j][i])
    	}
    	res = append(res, t)
    }
    b.Data = res
    return b

}

// 计算n阶行列式 (N = n - 1)
func Det(Matrix [][]float64, N int) float64 {
    var T0, T1, T2, Cha int
    var Num float64
    var B [][]float64

    if N > 0 {
    	Cha = 0
    	for i := 0; i < N; ++i {
    		var tmpArr []float64
    		for j := 0; j < N; ++j {
    			tmpArr = append(tmpArr, 0)
    		}
    		B = append(B, tmpArr)
    	}
    	Num = 0
    	for T0 = 0; T0 <= N; ++T0 {
    		for T1 = 1; T1 <= N; ++T1 {
    			for T2 = 0; T2 <= N - 1; ++T2 {
    				if T2 == T0 {
    					Cha = 1
    				}
    				B[T1 - 1][T2] = Matrix[T1][T2 + Cha]
    			} // T2循环
    			Cha = 0
    		}  // T1 循环
    		Num = Num + Matrix[0][T0] * Det(B, N - 1) * math.Pow(-1, float64(T0))
    	}   // T0 循环
    	return Num
    }
    else if N == 0 {
    	return Matrix[0][0]
    }
    return 0

}

// 矩阵求逆 (N = n - 1)
func Inverse(S1 SQ) (MatrixC [][]float64) {
    N := S1.N -1
    Matrix := S1.Data
    var T0, T1, T2, T3 int
    var B [][]float64
    for i := 0; i < N; ++i {
        var tmpArr []float64
        for j := 0; j < N; ++j {
            tmpArr = append(tmpArr, 0)
        }
        B = append(B, tmpArr)
    }

    for i := 0; i < N + 1; ++i {
    	var tmpArr []float64
    	for j := 0; j < N + 1; ++j {
    		tmpArr = append(tmpArr, 0)
    	}
    	MatrixC = append(MatrixC, tmpArr)
    }

    Chay := 0

    Chax := 0
    var add float64
    add = 1 / Det(Matrix, N)
    for T0 = 0; T0 <= N; ++T0 {
        for T3 = 0; T3 <= N; ++T3 {
            for T1 = 0; T1 <= N - 1; ++T1 {
                if T1 < T0 {
                Chax = 0
            }
            else
            {
                Chax = 1
            }
            for T2 = 0; T2 <= N - 1; ++T2 {
                if T2 < T3 {
                    Chay = 0
                }
                else
                {
                    Chay = 1
                }
                B[T1][T2] = Matrix[T1 + Chax][T2 + Chay]
            } // T2 循环
        }  // T1 循环
        Det(B, N - 1)
        MatrixC[T3][T0] = Det(B, N - 1) * add * (math.Pow(-1, float64(T0 + T3)))
    }

    return MatrixC
}

//
// Created by Xmx on 2022/7/26.
//

#include <iostream>
#include <cmath>
#include <fstream>
#include <initializer_list>

using namespace std;

template <class T>
class Matrix
{
public:
// 无参构造函数
Matrix() {}

    // 有参构造函数
    Matrix(int rows, int cols, T init_value)
    {
        set_rows(rows);
        set_cols(cols);
        new_mat();
        for (int i = 0; i < get_rows(); ++i)
            for (int j = 0; j < get_cols(); ++j)
                mat_[i][j] = init_value;
    }

    // 析构函数
    ~Matrix()
    {
        cout << "~Matrix" << endl;
        delete_mat();
    }

    // done:
    // 有参构造函数, 使用初始化列表
    Matrix(initializer_list<initializer_list<T>> listList)
    {
        rows_ = (listList.begin())->size();
        cols_ = listList.size();
        new_mat(); // 分配内存
        auto pr = listList.begin();
        auto pc = pr->begin();
        for (int i = 0; i < get_rows(); i++, pr++)
        {
            pc = pr->begin();
            for (int j = 0; j < get_cols(); j++, pc++)
                mat_[i][j] = *pc;
        }
    }

    // 重载拷贝构造函数, 实现深拷贝
    Matrix<T>(const Matrix<T> &matrix)
    {
        this->rows_ = matrix.get_rows();
        this->cols_ = matrix.get_cols();
        new_mat(); // 为拷贝接受对象开辟内存

        //拷贝数值
        for (int i = 0; i < matrix.get_rows(); ++i)
            for (int j = 0; j < matrix.get_cols(); ++j)
                mat_[i][j] = matrix.mat_[i][j];
    }

    //! 重载= 必须通过成员函数实现
    Matrix<T> &operator=(const Matrix<T> &matrix)
    {
        // 若其在堆区, 先释放mat_指针,再深拷贝即可
        if (mat_ != NULL)
        {
            delete mat_;
            mat_ = NULL;
        }

        // 深拷贝
        this->rows_ = matrix.get_rows();
        this->cols_ = matrix.get_cols();
        new_mat();
        for (int i = 0; i < matrix.get_rows(); ++i)
            for (int j = 0; j < matrix.get_cols(); ++j)
                mat_[i][j] = matrix.mat_[i][j];
        return *this;
    }

    bool operator==(const Matrix<T> &other)
    {
        if (this->rows_ == other.rows_ && this->cols_ == other.cols_)
        {
            if (this->mat_ != NULL && other.mat_ != NULL)
            {
                for (int i = 0; i < this->rows_; i++)
                    for (int j = 0; j < this->cols_; j++)
                        if (this->mat_[i][j] != other.mat_[i][j])
                            return false;
                return true;
            }
            else
                return false;
        }
        else
            return false;
    }

    bool operator!=(const Matrix<T> &other)
    {
        if (this->rows_ == other.rows_ && this->cols_ == other.cols_)
        {
            if (this->mat_ != NULL && other.mat_ != NULL)
            {
                for (int i = 0; i < this->rows_; i++)
                    for (int j = 0; j < this->cols_; j++)
                        if (mat_[i][j] != other.mat_[i][j])
                            return true;
                return true;
            }
            else
                return true;
        }
        else
            return true;
    }

    // 申请矩阵空间
    void new_mat()
    {
        this->mat_ = new T *[get_rows()];
        for (int i = 0; i < get_rows(); i++)
        {
            this->mat_[i] = new T[get_cols()];
        }
    }

    // 释放矩阵空间
    void delete_mat()
    {
        if (mat_ != NULL)
        {
            for (int i = 0; i < get_rows(); i++)
                delete[] this->mat_[i];
            delete[] mat_;
            mat_ = NULL;
        }
    }

    // *  ------------------------- 访问器与修改器 ---------------------------- *
    void set_rows(int rows)
    {
        this->rows_ = rows;
    }

    void set_cols(int cols)
    {
        this->cols_ = cols;
    }

    int get_rows() const
    {
        return this->rows_;
    }

    int get_cols() const
    {
        return this->cols_;
    }

    //* --------------------------- 文件和流 -------------------------------- *
    void load_matrix(string filename)
    {
        ifstream input;
        // input>>noskipws; //! 不忽略空格, 为了方便写入矩阵,不跳过空格
        input.open(filename);
        string read_contents = "";
        string temp;

        //! 读写模板
        // while (input>>temp)
        // {
        //     /* code */
        //     read_contents += temp;
        // }
        // cout << read_contents << endl;

        if (input.is_open())
        {
            int rows, cols;
            input >> rows >> cols;
            set_rows(rows);
            set_cols(cols);
            new_mat();
            for (int i = 0; i < get_rows(); i++)
                for (int j = 0; j < get_cols(); j++)
                    input >> mat_[i][j];

            cout << "read the txt successfully" << endl;
            cout << "rows: " << rows << " cols: " << cols << endl;
            input.close();
            cout << "The Matrix has been read successfully: " << endl;
            cout << *this << endl;
        }
        else
            cout << "File Opening Failed" << endl;
    }

    void write_matrix(string filename)
    {
        ofstream output(filename);
        if (output.is_open())
        {
            output << *this << endl;
            output.close();
            cout << "Matrix has been written" << endl;
        }
        else
            cout << "File Opening Failed" << endl;
    }

    // void show_mat()
    // {
    //     cout << endl;
    //     for (int i = 0; i < get_rows(); i++)
    //     {
    //         for (int j = 0; j < get_cols(); j++)
    //         {
    //             cout << this->mat_[i][j] << " ";
    //         }
    //         cout << endl;
    //     }
    //     cout << endl;
    // }

    // * ------------------------ 矩阵的基本运算函数 -------------------------- *
    void swap_element(T &element1, T &element2)
    {
        T temp = element1;
        element1 = element2;
        element2 = temp;
    }

    // 交换两行
    void swap_rows(int row1, int row2)
    {
        for (int i = 0; i < get_cols(); i++)
            swap_element(mat_[row1][i], mat_[row2][i]);
    }

    // 交换两列
    void swap_cols(int col1, int col2)
    {
        for (int i = 0; i < get_rows(); i++)
            swap_element(mat_[i][col1], mat_[i][col2]);
    }

    // *  ------------------------- 运算符重载的友元声明 ------------------------ *
    template <class U>
    friend istream &operator>>(istream &in, Matrix<U> &A); //! >>
    template <class U>
    friend ostream &operator<<(ostream &out, Matrix<U> &A); //! <<
    template <class U>
    friend Matrix<U> operator+(const Matrix<U> &A, const U &k); //! A + k
    template <class U>
    friend Matrix<U> operator+(const U &k, const Matrix<U> &A); //! k + A
    template <class U>
    friend Matrix<U> operator-(const Matrix<U> &A, const U &k); //! A - k
    template <class U>
    friend Matrix<U> operator-(const U &k, const Matrix<U> &A); //! k - A
    template <class U>
    friend Matrix<U> operator*(const U &k, const Matrix<U> &A); //! k * A
    template <class U>
    friend Matrix<U> operator*(const Matrix<U> &A, const U &k); //! A * k
    template <class U>
    friend Matrix<U> operator/(const U &k, const Matrix<U> &A); //! k / A
    template <class U>
    friend Matrix<U> operator/(const Matrix<U> &A, const U &k); //! A / k
    template <class U>
    friend Matrix<U> operator+(const Matrix<U> &A, const Matrix<U> &B); //! A + B
    template <class U>
    friend Matrix<U> operator-(const Matrix<U> &A, const Matrix<U> &B); //! A - B
    template <class U>
    friend Matrix<U> operator*(const Matrix<U> &A, const Matrix<U> &B); //! A * B
    template <class U>
    friend Matrix<U> operator/(const Matrix<U> &A, const Matrix<U> &B); //! A / B
    template <class U>
    friend Matrix<U> operator^(const Matrix<U> &A, const U &B); //! A^k

    //* ------------------------- 静态成员函数生成基本矩阵 --------------------------- */
    // 零矩阵
    static Matrix<T> zeros(int rows, int cols)
    {
        return Matrix<T>(rows, cols, 0);
    }

    // 单位矩阵
    static Matrix<T> units(int size)
    {
        Matrix<T> result = Matrix<T>::zeros(size, size); // 先生成一个size大小的 "方零矩阵"
        for (int i = 0; i < result.get_cols(); i++)      // 修改对角线元素
            result.mat_[i][i] = 1;
        return result;
    }

    //! done: 求逆矩阵(高斯消元实现)
    template <class U>
    static Matrix<U> inverse(const Matrix<U> &A)
    {
        double eps = 1e-6;
        // 返回double类型,防止精度丢失
        Matrix<double> temp = Matrix<double>::augmented_E(A);
        for (int i = 0; i < temp.get_rows(); i++)
        {
            if (fabs(temp.mat_[i][i]) < eps)
            {
                int j;
                for (j = i + 1; j < temp.get_rows(); j++)
                    if (fabs(temp.mat_[i][j]) > eps)
                        break;
            }

            // 初等变换实现消元
            for (int i = 0; i < temp.get_rows(); i++)
            {
                double t = temp.mat_[i][i];
                // 将temp.mat_[i][i] = 1;
                for (int j = i; j < temp.get_cols(); j++)
                    temp.mat_[i][j] /= t;

                for (int j = i + 1; j < temp.get_rows(); j++)
                {
                    double tt = -1 * (temp.mat_[j][i] / temp.mat_[i][i]);
                    for (int k = i; k < temp.get_cols(); k++)
                        temp.mat_[j][k] += tt * temp.mat_[i][k];
                }
            }
        }

        for (int i = temp.get_rows() - 1; i >= 0; i--)
        {
            for (int j = i - 1; j >= 0; j--)
            {
                double tt = -1 * (temp.mat_[j][i] / temp.mat_[i][i]);
                for (int k = i; k < temp.get_cols(); k++)
                    temp.mat_[j][k] += tt * temp.mat_[i][k];
            }
        }

        // 将temp左半部分传给result
        Matrix<double> result(A.get_rows(), A.get_cols(), 0);
        for (int i = 0; i < result.get_rows(); i++)
            for (int j = result.get_rows(); j < 2 * result.get_rows(); j++)
                result.mat_[i][j - result.get_rows()] = temp.mat_[i][j];

        // cout<<temp;

        return result;
    }

    //! done: 返回A的增广矩阵A|E
    static Matrix<T> augmented_E(const Matrix<T> &A)
    {
        int rows = A.get_rows();
        int cols = A.get_rows() + A.get_cols();
        Matrix<T> result = Matrix<T>(rows, cols, 0);
        for (int i = 0; i < result.get_rows(); i++)
        {
            for (int j = 0; j < result.get_cols(); j++)
            {
                if (j < A.get_cols())
                {
                    result.mat_[i][j] = A.mat_[i][j];
                }
                else
                {
                    if (j == A.get_cols() + i)
                        result.mat_[i][j] = 1;
                    else
                        result.mat_[i][j] = 0;
                }
            }
        }
        return result;
    }

    //! done: 增广矩阵A' , 返回A的增广矩阵
    static Matrix<T> augmented(const Matrix<T> &A, const T *input, int len)
    {
        int rows = A.get_rows();
        int cols = A.get_rows() + 1;
        for (int i = 0; i < len; i++)
            cout << input[i] << " ";
        // input数组长度为A的行数满足条件,可以生成增广矩阵
        if (len == rows)
        {
            Matrix<T> result = Matrix<T>(rows, cols, 0);
            for (int i = 0; i < result.get_rows(); i++)
            {
                for (int j = 0; j < result.get_cols(); j++)
                {
                    if (j < A.get_cols())
                        result.mat_[i][j] = A.mat_[i][j];
                    else
                        result.mat_[i][j] = input[i];
                }
            }
            return result;
        }
        else
            return Matrix<T>::zeros(1, 1);
    }

    //! done: 高斯消元(初等变换法)求解线性方程
    static string Gauss_eliminate(const Matrix<T> &A)
    {
        Matrix<T> temp = A;           // 用于计算的矩阵
        double ratio = 0, eps = 1e-6; // 默认最小主元素
        int n = temp.get_rows();      // 方程个数
        double result[n];             // 存储结果的数组
        string result_string = "";
        // 消元
        for (int k = 0; k < (n); k++)
        {
            for (int i = (k + 1); i < n; i++)
            {
                if (abs(temp.mat_[k][k]) < eps)
                    return "No solution!";

                double ratio = temp.mat_[i][k] / temp.mat_[k][k];
                for (int j = (k + 1); j < (n + 1); j++)
                    temp.mat_[i][j] -= ratio * temp.mat_[k][j];

                temp.mat_[i][k] = 0;
            }
        }
        result[n - 1] = temp.mat_[n - 1][n] / temp.mat_[n - 1][n - 1]; // 回代
        for (int i = (n - 2); i >= 0; i--)
        {
            double sum = 0;
            for (int j = (i + 1); j < n; j++)
            {
                sum += temp.mat_[i][j] * result[j];
                result[i] = (temp.mat_[i][n] - sum) / temp.mat_[i][i];
            }
        }

        // 生成结果字符串
        for (int i = 0; i < n; i++)
            result_string += "\nresult[" + to_string(i) + "]=" + to_string(result[i]) + "\n";

        return result_string;
    }

    //! done: 高斯核矩阵
    static Matrix<double> Gauss_kernel(int size, double sigma)
    {
        Matrix<double> result = Matrix<double>::zeros(size, size);
        const double PI = 4.0 * atan(1.0); // 圆周率赋值
        int center = size / 2;
        double sum = 0;
        for (int i = 0; i < size; i++)
        {
            for (int j = 0; j < size; j++)
            {
                result.mat_[i][j] = (1 / (2 * PI * sigma * sigma)) *
                                    exp(-((i - center) * (i - center) + (j - center) * (j - center)) /
                                        (2 * sigma * sigma));
                sum += result.mat_[i][j];
            }
        }
        return result;
    }

    //! done: H 为高斯滤波核矩阵, A 需要高斯滤波的矩阵, 返回滤波后的矩阵
    static Matrix<double> Gauss_filter(const Matrix<T> &A, const Matrix<double> &H)
    {
        Matrix<double> result = Matrix::zeros(A.get_rows(), A.get_cols());
        for (int i = 0; i < result.get_rows(); i++)
            for (int j = 0; j < result.get_cols(); j++)
                result.mat_[i][j] = A.mat_[i][j] * H.mat_[i][j];
        return result;
    }

private:
int rows\_;
int cols\_;
T \*\*mat\_;
};

//\*  --------------------------- 运算符重载 ------------------------------------- \*/

//! done: 重载>>
template <class U>
istream \&operator>>(istream \&in, Matrix<U> \&A)
{
cout << "Input the rows and columns for the matrix" << endl;
in >> A.rows\_ >> A.cols\_;
A.new\_mat();
for (int i = 0; i < A.get\_rows(); i++)
for (int j = 0; j < A.get\_cols(); j++)
in >> A.mat\_\[i]\[j];
return in;
}

//! done: 重载<<
template <class U>
ostream \&operator<<(ostream \&out, Matrix<U> \&A)
{
for (int i = 0; i < A.get\_rows(); i++)
{
out << " \[";
for (int j = 0; j < A.get\_cols(); j++)
out << " " << A.mat\_\[i]\[j];
out << " ]" << endl;
}
out << endl;
return out;
}

//! done A + k
template <class U>
Matrix<U> operator+(const Matrix<U> \&A, const U \&k)
{
if (A.get\_cols() == A.get\_rows()) // A 为方针,可以运算
{
Matrix<U> result = Matrix<U>(A);
for (int i = 0; i < result.get\_cols(); i++) // A + kE
result.mat\_\[i]\[i] += k;
return result;
}
else // 无法运算,返回1×1单位矩阵
return Matrix<U>::zeros(1, 1);
}

//! done k + A
template <class U>
Matrix<U> operator+(const U \&k, const Matrix<U> \&A)
{
if (A.get\_cols() == A.get\_rows()) // A 为方针,可以运算
{
Matrix<U> result = Matrix<U>(A);
for (int i = 0; i < result.get\_cols(); i++) // kE + A
result.mat\_\[i]\[i] += k;
return result;
}
else // 无法运算,返回1×1单位矩阵
return Matrix<U>::zeros(1);
}

//! done: A - k
template <class U>
Matrix<U> operator-(const Matrix<U> \&A, const U \&k)
{
if (A.get\_cols() == A.get\_rows()) // A 为方针,可以运算
{
Matrix<U> result = Matrix<U>(A);
for (int i = 0; i < result.get\_cols(); i++) // A - kE
result.mat\_\[i]\[i] -= k;
return result;
}
else // 无法运算,返回1×1单位矩阵
return Matrix<U>::zeros(1);
}

//! done k - A
template <class U>
Matrix<U> operator-(const U \&k, const Matrix<U> \&A)
{
if (A.get\_cols() == A.get\_rows()) // A 为方针,可以运算
{
Matrix<U> result = Matrix<U>(A);
for (int i = 0; i < result.get\_cols(); i++) // kE - A
for (int j = 0; j < result.get\_cols(); j++)
{
if (i == j)
result.mat\_\[i]\[j] = k - result.mat\_\[i]\[j];
else
result.mat\_\[i]\[j] = -result.mat\_\[i]\[j];
}
return result;
}
else // 无法运算,返回1×1单位矩阵
return Matrix<U>::zeros(1);
}

//! done: k \* A
template <class U>
Matrix<U> operator\*(const U \&k, const Matrix<U> \&A)
{
Matrix<U> result = Matrix<U>(A);
for (int i = 0; i < result.get\_rows(); i++)
for (int j = 0; j < result.get\_cols(); j++)
result.mat\_\[i]\[j] = result.mat\_\[i]\[j] \* k;
return result;
}

//! done: A \* k
template <class U>
Matrix<U> operator\*(const Matrix<U> \&A, const U \&k)
{
Matrix<U> result = Matrix<U>(A);
for (int i = 0; i < result.get\_rows(); i++)
for (int j = 0; j < result.get\_cols(); j++)
result.mat\_\[i]\[j] = result.mat\_\[i]\[j] \* k;
return result;
}

//! done: k / A
template <class U>
Matrix<U> operator/(const U \&k, const Matrix<U> \&A)
{
// k / A = kE \* A^-1
// A 为方针,可以运算
if (A.get\_rows() == A.get\_cols())
{
Matrix<U> result = Matrix<U>(A);
Matrix<U> kE = Matrix<U>::units(A.get\_rows());
kE = kE \* k; // 生成kE
return kE \* Matrix<U>::inverse(A);
}
else // 无法运算,返回零矩阵
return Matrix<U>::zeros(1, 1);
}

//! done: A / k
template <class U>
Matrix<U> operator/(const Matrix<U> \&A, const U \&k)
{
Matrix<U> result = Matrix<U>(A);
for (int i = 0; i < result.get\_rows(); i++)
for (int j = 0; j < result.get\_cols(); j++)
result.mat\_\[i]\[j] = result.mat\_\[i]\[j] / k;
return result;
}

//! done: A + B
template <class U>
Matrix<U> operator+(const Matrix<U> \&A, const Matrix<U> \&B)
{
// A 与 B 维度相同
if ((A.get\_cols() == B.get\_cols()) && (A.get\_rows() == B.get\_rows()))
{
Matrix<U> result = Matrix<U>(A);
for (int i = 0; i < result.get\_rows(); i++)
for (int j = 0; j < B.get\_cols(); j++)
result.mat\_\[i]\[j] += B.mat\_\[i]\[j];
return result;
}
else // 无法运算,返回零矩阵
return Matrix<U>::zeros(1, 1);
}

//! done: A - B
template <class U>
Matrix<U> operator-(const Matrix<U> \&A, const Matrix<U> \&B)
{
// A 与 B 维度相同
if ((A.get\_cols() == B.get\_cols()) && (A.get\_rows() == B.get\_rows()))
{
Matrix<U> result = Matrix<U>(A);
for (int i = 0; i < result.get\_rows(); i++)
for (int j = 0; j < B.get\_cols(); j++)
result.mat\_\[i]\[j] -= B.mat\_\[i]\[j];
return result;
}
else // 无法运算,返回零矩阵
return Matrix<U>::zeros(1, 1);
}

//! done: A \* B
template <class U>
Matrix<U> operator\*(const Matrix<U> \&A, const Matrix<U> \&B)
{
// 先判断能否运算
if (A.get\_cols() == B.get\_rows())
{
Matrix<U> result = Matrix<U>(A.get\_rows(), B.get\_cols(), 0);
for (int i = 0; i < result.get\_rows(); i++)
for (int j = 0; j < result.get\_cols(); j++)
for (int k = 0; k < result.get\_cols(); k++) // 所在行所在列元素相乘
result.mat\_\[i]\[j] += A.mat\_\[i]\[k] \* B.mat\_\[k]\[j];
return result;
}
else // 无法运算,返回零矩阵
return Matrix<U>::zeros(1, 1);
}

//! done: A / B
template <class U>
Matrix<U> operator/(const Matrix<U> \&A, const Matrix<U> \&B)
{
// A / B == A \* B^-1
// 先判断能否运算
if (A.get\_cols() == B.get\_rows())
return A \* (Matrix<U>::inverse(B));
else
return Matrix<U>::zeros(1, 1);
}

//! done: A^k
template <class U>
Matrix<U> operator^(const Matrix<U> \&A, int k)
{
Matrix<U> result = Matrix<U>::units(A.get\_cols());
Matrix<U> temp = A;

    // 快速幂进行运算
    while (k)
    {
        if (k % 2)
            result = result * temp;
        k /= 2;
        temp = temp * temp;
    }
    return result;

}

int main()
{
// Matrix<int> matrix1(3, 3, 3);
// Matrix<double> matrix1\_4(3, 4, 4);
// Matrix<int> matrix2(4, 4, 4);
// // Matrix<double> matrix3 = { {5,5,-1,7,54},{4,-9,20,12,-6},{9,-18,-3,1,21},{ 61,-8,-10,3,13 },{ 29,-28,-1,4,14 } };
// // Matrix<int> matrix4 = Matrix<int>(matrix1);
// Matrix<double> matrix5 = matrix1\_4 \* (0.3);
// Matrix<int> matrix6 = Matrix<int>::units(3);
// matrix1.show\_mat();
// matrix2.show\_mat();
// // matrix3.show\_mat();
// // matrix4.show\_mat();
// matrix5.show\_mat();
// matrix6.show\_mat();
// cout<\<matrix6;
// Matrix<double> matrix = Matrix<double>::augmented\_E(matrix1\_4);
// matrix.swap\_rows(0,1);
// cout << matrix;
// Matrix<double> matrix10;
// cin >> matrix10;
// Matrix<double> temp = Matrix<double>::inverse(matrix10);
// Matrix<double> temp1 = Matrix<double>::inverse(matrix10);
// cout << temp;
// temp = temp ^ 2;
// cout<\<temp;
// temp1 = temp1 \* temp1;
// cout<\<temp1;
// cout <<( matrix != matrix);

    //! 测试=, ==, != 重载
    Matrix<int> m1(3, 3, 3);
    Matrix<int> m2(3, 3, 4);
    cout << m1 << m2 << endl;
    cout << (m1 == m2);

    //     double input[] = {6,18,7};
    //     cout<<"len = "<< sizeof(input)/sizeof(input[0])<<endl;
    //      Matrix<double> Equations;
    //      cin >> Equations;
    //      Equations = Matrix<double>::augmented(Equations,input,sizeof(input)/sizeof(input[0]));
    //      cout<< Equations;
    //      cout<< Matrix<double>::Gauss_eliminate(Equations);

    //! 测试高斯滤波
    Matrix<double> temp(5, 5, 100);
    cout << temp << endl;
    temp = Matrix<double>::Gauss_filter(temp, Matrix<double>::Gauss_kernel(5, 5));
    cout << temp;

    //! 测试文件和流
    // Matrix<double> matrix_IO;
    // matrix_IO.load_matrix("input.txt");
    // matrix_IO.write_matrix("output.txt");

    // 测试初始化列表
    initializer_list<initializer_list<int>> initializerL_mat = {{1, 2, 3}, {2, 2, 2}, {3, 3, 3}};
    Matrix<int> matrix_initializer(initializerL_mat);
    cout << matrix_initializer;

}