传送门，这个讲得太好了吧！！！
矩阵乘法：
满足结合律和分配律，但是不满足交换律（因为乘出来的矩阵为左行右列）

1.简单的A*B矩阵

c(i,j) = sum a[i][k]*b[k][j];

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 110;
int a[N][N],b[N][N],c[N][N];
void solve(){
    
	int n,m,p;
	cin>>n>>m;
	for(int i=1;i<=n;i++){
    
		for(int j=1;j<=m;j++){
    
			cin>>a[i][j];
		}
	}
	cin>>p;
	for(int i=1;i<=m;i++){
    
		for(int j=1;j<=p;j++){
    
			cin>>b[i][j];
		}
	}
	for(int i=1;i<=n;i++){
    
		for(int j=1;j<=p;j++){
    
			for(int k=1;k<=m;k++){
    //共同的是左矩阵的列等于右矩阵的行，放最后 
				c[i][j]+=a[i][k]*b[k][j];//累加，最后得到的新矩阵的行和列分别等于左矩阵的行和右矩阵的列 
			}
		}
	}
	for(int i=1;i<=n;i++){
    
		for(int j=1;j<=p;j++){
    
			cout<<c[i][j]<<" ";
		}
		cout<<endl;
	}
}
int main(){
    
	int t;
	t=1;
	while(t--){
    
		solve();
	}
	return 0;
} 

2.斐波那契

//矩阵乘法加速递推 斐波那契数列 
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 10000;
int n;
void mul(int f[2],int a[2][2]){
    //对f操作 
	int c[2];
	memset(c,0,sizeof(c));
	for(int j=0;j<2;j++){
    
		for(int k=0;k<2;k++){
    
			c[j]=(c[j]+(ll)f[k]*a[k][j]) % mod;
		}
	}
	memcpy(f,c,sizeof(c));//复制 
}
void mulself(int a[2][2]){
    //矩阵乘法 
	int c[2][2];
	memset(c,0,sizeof(c));
	for(int i=0;i<2;i++){
    
		for(int j=0;j<2;j++){
    
			for(int k=0;k<2;k++){
    
				c[i][j]=(c[i][j]+(ll)a[i][k]*a[k][j])%mod;
			}
		}
	}
	memcpy(a,c,sizeof(c));//复制 
}
int main(){
    
	while(cin>>n&&n!=-1){
    
		int f[2]={
    0,1};//每次只需要保留2个数据 
		int a[2][2]={
    {
    0,1},{
    1,1}};//A矩阵 
		for(;n;n>>=1){
    //快速幂 
			if(n&1) mul(f,a);
			mulself(a);
		}
		cout<<f[0]<<endl;
	}
}

3.矩阵快速幂模板

详细讲解

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 105;
const int mod = 1e9+7;
struct matrix{
    
	ll m[N][N];
}; 
ll n,k;
matrix mul(matrix a,matrix b){
    //矩阵乘法 
	matrix ans;
	for(int i=1;i<=n;i++){
    
		for(int j=1;j<=n;j++){
    
			for(int k=1;k<=n;k++){
    
				ans.m[i][j]=(ans.m[i][j]+a.m[i][k]*b.m[k][j]%mod)%mod;
			}
		}
	}
	return ans;
}
matrix qk(matrix a,ll b){
    //矩阵快速幂 
	matrix ans;//答案矩阵
	//先初始化 
	for(int i=1;i<=n;i++){
    
		for(int j=1;j<=n;j++){
    
			if(i==j) ans.m[i][j]=1;
			else ans.m[i][j]=0;
		}
	}
	//快速幂的矩阵写法 
	while(b){
    
		if(b&1) ans = mul(a,ans); 
		a=mul(a,a);
		b>>=1;
	}
	return ans;
}
int main(){
    
	cin>>n>>k;
	matrix a;
	for(int i=1;i<=n;i++){
    
		for(int j=1;j<=n;j++){
    
			cin>>a.m[i][j];
		}
	}
	matrix ans = qk(a,k);//矩阵快速幂 
	for(int i=1;i<=n;i++){
    
		for(int j=1;j<=n;j++){
    
			cout<<ans.m[i][j]<<" ";
		}
		cout<<endl;
	}
	return 0;
}