矩阵的加减法和数乘为线性运算，均为逐个元素进行运算。

形式化地，即 $AopB=A_{i,j}op{B}_{i,j}$，$op$ 为 $+/-$。

矩阵乘法要求两个矩阵的行数和列数相等，设 $A$ 为 $P\times M$ 的矩阵，$B$ 为 $M\times Q$ 的矩阵，则 $C=A\times B,C_{i,j}={\sum }_{k=1}^M{A}_{i,k}\times B_{k,j}$。即矩阵 $A$ 第 $i$ 行 $M$ 个数与矩阵 $B$ 第 $j$ 列 $M$ 个数分别相乘再相加得到的，如下图。
为了方便进行运算，我们会把矩阵设为二维数组，放在结构体里面，并重载运算符。
这是我的矩阵加减乘法的模版。

struct matrix{
	int a[N][N];
	matrix(){memset(a,0,sizeof(a));}
	matrix operator +(const matrix& T) const{
		matrix res;
		for(int i=1;i<N;i++)
			for(int j=1;j<N;j++)
				res.a[i][j]=(a[i][j]+T.a[i][j])%mod;
		return res;
	}
	matrix operator -(const matrix& T) const{
		matrix res;
		for(int i=1;i<N;i++)
			for(int j=1;j<N;j++)
				res.a[i][j]=(a[i][j]-T.a[i][j])%mod;
		return res;
	}
	matrix operator *(const matrix& T) const{
		matrix res;//这个写法常数会小一些
		for(int i=1;i<N;i++){
			for(int k=1;k<N;k++){
				int r=a[i][k];
				for(int j=1;j<N;j++)
					res.a[i][j]=(res.a[i][j]+1ll*r*T.a[k][j])%mod;
			}
		}
		return res;
	}
};

由于矩阵也有乘法，那自然也有矩阵的快速幂，写法与普通快速幂同理。
附洛谷矩阵快速幂模版代码：

#include<bits/stdc++.h>
using namespace std;
#define int long long
#define rd read()
const int N=110,mod=1e9+7;
int read(){
	int x=0,f=1;char ch=getchar();
	while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
	while(isdigit(ch)) x=x*10+ch-'0',ch=getchar();
	return x*f;
}
int n,k;
struct matrix{
	int a[N][N];
	matrix(){memset(a,0,sizeof(a));}
	matrix operator +(const matrix& T) const{
		matrix res;
		for(int i=1;i<N;i++)
			for(int j=1;j<N;j++)
				res.a[i][j]=(a[i][j]+T.a[i][j])%mod;
		return res;
	}
	matrix operator -(const matrix& T) const{
		matrix res;
		for(int i=1;i<N;i++)
			for(int j=1;j<N;j++)
				res.a[i][j]=(a[i][j]-T.a[i][j])%mod;
		return res;
	}
	matrix operator *(const matrix& T) const{
		matrix res;
		for(int i=1;i<N;i++){
			for(int k=1;k<N;k++){
				int r=a[i][k];
				for(int j=1;j<N;j++)
					res.a[i][j]=(res.a[i][j]+1ll*r*T.a[k][j])%mod;
			}
		}
		return res;
	}
}M;
signed main(){
	n=rd,k=rd;
	for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) M.a[i][j]=rd;
	matrix res;
	for(int i=1;i<=n;i++) res.a[i][i]=1;//res要初始化为单位矩阵，即单位矩阵*M^K
	while(k){
		if(k&1) res=res*M;
		M=M*M,k>>=1;
	}
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++) cout<<res.a[i][j]<<" ";
		cout<<"\n";
	}
	return 0;
}

而矩阵乘法的运用很多，最常用的便是加速递推，可以用其优化 dp（例如动态 dp）等等。最经典的便是斐波那契数列的应用，首先把原线性递推式转化为矩阵，然后用矩阵快速幂即可解决。

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define rd read()
ll read(){
	ll x=0,f=1;char ch=getchar();
	while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
	while(isdigit(ch)) x=x*10+ch-'0',ch=getchar();
	return x*f;
}
const int mod=1e9+7;ll n;
struct matrix{
	ll a[10][10];
	matrix(){memset(a,0,sizeof(a));}
	matrix operator *(const matrix& T)const{
		matrix res;
		for(int i=1;i<=2;i++){
			for(int k=1;k<=2;k++){
				ll r=a[i][k];
				for(int j=1;j<=2;j++)
					res.a[i][j]=(res.a[i][j]+1ll*r*T.a[k][j]%mod)%mod;
			}
		}
		return res;
	}
}bas,ans;
int main(){
	n=rd;ll x=n-2;
	if(n<=2){puts("1");return 0;}
	ans.a[1][1]=ans.a[1][2]=1;
	bas.a[1][1]=bas.a[1][2]=bas.a[2][1]=1;
	while(x){
		if(x&1) ans=ans*bas;
		bas=bas*bas,x>>=1;
	}
	printf("%lld\n",ans.a[1][1]%mod);
	return 0;
}

再来一道（P1939 矩阵加速（数列））

考虑$\left[\begin{array}{ccc}{a}_i& a_{i-1}& a_{i-2}\end{array}\right]$ 如何变为 $\left[\begin{array}{ccc}{a}_{i+1}& a_i& a_{i-1}\end{array}\right]$，这道题非常显然。

[ a i + 1 a i a i − 1 ] × [ 1 1 0 0 0 1 1 0 0 ] \begin{bmatrix} a_{i+1}& a_{i}&a_{i-1} \end{bmatrix}\times \begin{bmatrix} 1 & 1 &0 \\ 0& 0& 1\\ 1& 0&0 \end{bmatrix} [ai+1​​ai​​ai−1​​]× ​101​100​010​ ​。然后套个矩阵加速幂就可以了。