第六章 矩阵

一，矩阵相乘 A*B

1,矩阵乘法一般算法

前提：A的列数等于B的行数

原理：

代码：

//一般算法
struct mat()
{
    int r,c;
    double date[MAXN][MAXN];
};

int mul(mat& d,const nat& a,const mat& b)
{
    int i,j,k;
    if(a.c!=b.r)return 0;
    d.r=a.r;
    d.c=b.c;
    for(i=0;i<d.r;i++)
        for(j=0;j<d.c;j++)
            for(d.date[i][j]=k=0;k<a.c;k++)
                d.date[i][j]+=a.date[i][k]*b.date[k][j];
    return 1;
}

2，传统分治算法

C=AB

将A B C分成大小相等的方块

（分治算法的思想是，将一个难以解决的大问题，分割成一些规模较小的相同问题，分而治之）

如果n=2,则2个2阶方阵的乘积可以直接用(2)-(5)式计算出来，共需8次乘法和四次加法。

当矩阵的阶数大于2时，可以继续将子矩阵分块，直到子矩阵的阶数降为2.

这样就产生了一个分治降阶的递归算法。

3，Strassen算法

相较于分治算法，它减少了乘法运算的次数

做了七次乘法，再做若干次的加法减法

数据的存储结构：二维数组，二重指针

#include<iostream>
#include<cmath>
#define N 4
using namespace std;

//二阶矩阵相乘
void Matrix_Multiply(int A[][N],int B[][N],int C[][N])
{
    for(int i=0;i<2;i++)
        for(int j=0;j<2;j++)
        {
            C[i][j]=0;
            for(int t=0;t<2;t++)
                C[i][j]=C[i][j]+A[i][t]*B[t][j];
        }
}

//矩阵加法
void Add(int n,int A[][N],int B[][N],int R[][N])
{
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            R[i][j]=A[i][j]+B[i][j];
}

//矩阵减法
void Sub(int n,int A[][N],int B[][N],int R[][N])
{
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            R[i][j]=A[i][j]-B[i][j];
}

//Strassen算法
void Strassen(int n,int A[][N],int B[][N],int C[][N])
{
    //中间辅助矩阵
    int A11[N][N], A12[N][N], A21[N][N], A22[N][N];
    int B11[N][N], B12[N][N], B21[N][N], B22[N][N];
    int C11[N][N], C12[N][N], C21[N][N], C22[N][N];
    int AA[N][N], BB[N][N];
    int M1[N][N], M2[N][N], M3[N][N], M4[N][N], M5[N][N], M6[N][N], M7[N][N];

    if(n==2)
        Matrix_Multiply(A,B,C);
    else
    {
        for(int i=0;i<n/2;i++)
            for(int j=0;j<n/2;j++)
            {
                A11[i][j] = A[i][j];
                A12[i][j] = A[i][j + n / 2];
                A21[i][j] = A[i + n / 2][j];
                A22[i][j] = A[i + n / 2][j + n / 2];

                B11[i][j] = B[i][j];
                B12[i][j] = B[i][j + n / 2];
                B21[i][j] = B[i + n /2][j];
                B22[i][j] = B[i + n /2][j + n / 2];
            }
        Sub(n / 2, B12, B22, BB);
        Strassen(n / 2, A11, BB, M1);

        Add(n / 2, A11, A12, AA);
        Strassen(n / 2, AA, B22, M2);

        Add(n / 2, A21, A22, AA);
        Strassen(n / 2, AA, B11, M3);

        Sub(n / 2, B21, B11, BB);
        Strassen(n / 2, A22, BB, M4);

        Add(n / 2, A11, A22, AA);
        Add(n / 2, B11, B22, BB);
        Strassen(n / 2, AA, BB, M5);

        Sub(n / 2, A12, A22, AA);
        Add(n / 2, B21, B22, BB);
        Strassen(n / 2, AA, BB, M6);

        Sub(n / 2, A11, A21, AA);
        Add(n / 2, B11, B12, BB);
        Strassen(n / 2, AA, BB, M7);

        //C11 = M5 + M4 - M2 + M6
        Add(n / 2, M5, M4, AA);
        Sub(n / 2, M6, M2, BB);
        Add(n / 2, AA, BB, C11);

        //C12 = M1 + M2
        Add(n / 2, M1, M2, C12);

        //C21 = M3 + M4
        Add(n / 2, M3, M4, C21);

        //C22 = M5 + M1 - M3 - M7
        Sub(n / 2, M5, M3, AA);
        Sub(n / 2, M1, M7, BB);
        Add(n / 2, AA, BB, C22);

        for(int i=0;i<n/2;i++)
            for(int j=0;j<n/2;j++)
            {
                C[i][j]=C11[i][j];
                C[i][j+n/2]=C12[i][j];
                C[i+n/2][j]=C12[i][j];
                C[i+n/2][j+n/2]=C22[i][j];
            }
    }
}

int main()
{
    int A[N][N],B[N][N],C[N][N];
    
    cout<<"input A"<<endl;
    for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
            cin>>A[i][j];
            
    cout<<"input B"<<endl;
    for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
            cin>>B[i][j];
            
     Strassen(N,A,B,C);
     
     cout<<"C:"<<endl;
     for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
        {
            cout<<C[i][j]<<' ';
            if(j>0&&j%3==0)
                cout<<endl;
        }
        
    return 0;
}

二，数字方阵

1，模拟法生成N阶奇数幻方

#include<iostream>
using namespace std;
int n;
int f[20][20];
void solve(int step,int i,int j)
{
    if(step==n*n+1)return ;
    f[i][j]=step;
    if(f[(i-1+n)%n][(j+1+n)%n])
        solve(step+1,(i+1)%n,j);
    else
        solve(step+1,(i-1+n)%n,(j+1)%n);
}
int main()
{
    cin>>n;
    solve(1,0,n/2);
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++)
            cout<<f[i][j]<<' ';
        cout<<endl;
    }
    return 0;
}

2，N阶斜线方阵

#include<iostream>
using namespace std;
int n,step;
int f[21][21];
int main()
{
    cin>>n;
    
    for(int i=1;i<=n;i++)//左上
        for(int j=1;j<=i;j++)
            f[j][i+1-j]=++step;
            
    for(int i=2;i<=n;i++)//右下
        for(int j=i;j<=n;j++)
            f[j][i+n-j]=++step;
            
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)
            cout<<f[i][j]<<' ';
        cout<<endl;
    }
    
    return 0;
}

3，N阶数螺旋方阵

#include <iostream>
#include <string>
#define N 100
using namespace std;

int main()
{
    int n;
    cin>>n;
    int res[N][N];
    int k=1;
    for (int i=0;i<n/2;++i)
    {
        //上边
        for (int j=i;j<n-1-i;++j)
        {
            res[i][j]=k;
            ++k;
        }
        //右边
        for (int j=i;j<n-1-i;++j)
        {
            res[j][n-i-1]=k;
            ++k;
        }
        //下边
        for (int j=n-i-1;j>i;--j)
        {
            res[n-i-1][j]=k;
            ++k;
        }
        //左边
        for (int j=n-i-1;j>i;--j)
        {
            res[j][i]=k;
            ++k;
        }
    }
    if (1 == n%2)
    {
        res[n/2][n/2] = k;
    }
    for (int i=0;i<n;++i)
    {
        for (int j=0;j<n;++j)
        {
            cout<<res[i][j]<<' ';
        }
        cout<<endl;
    }
    return 0;
}

三，线性方程组及其解法

1，高斯消元法

用处：求解线性方程组，求解矩阵的秩，求矩阵的逆

步骤：消元，回代

高斯消元求解先行方程组模版：

#include<stdio.h>
    #include<algorithm>
    #include<iostream>
    #include<string.h>
    #include<math.h>
    using namespace std;
    const int MAXN=50;
    int a[MAXN][MAXN];//增广矩阵
    int x[MAXN];//解集
    bool free_x[MAXN];//标记是否是不确定的变元
    int gcd(int a,int b){
        if(b == 0) return a; else return gcd(b,a%b);
    }
    inline int lcm(int a,int b){
        return a/gcd(a,b)*b;//先除后乘防溢出
    }
    // 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解，但无整数解，
    //-1表示无解，0表示唯一解，大于0表示无穷解，并返回自由变元的个数)
    //有equ个方程，var个变元。增广矩阵行数为equ,分别为0到equ-1,列数为var+1,分别为0到var.
    int Gauss(int equ,int var){
        int i,j,k;
        int max_r;// 当前这列绝对值最大的行.
        int col;//当前处理的列
        int ta,tb;
        int LCM;
        int temp;
        int free_x_num;
        int free_index;

        for(int i=0;i<=var;i++){
            x[i]=0;
            free_x[i]=true;
        }

        //转换为阶梯阵.
        col=0; // 当前处理的列
        for(k = 0;k < equ && col < var;k++,col++){// 枚举当前处理的行.
        // 找到该col列元素绝对值最大的那行与第k行交换.(为了在除法时减小误差)
            max_r=k;
            for(i=k+1;i<equ;i++){
                if(abs(a[i][col])>abs(a[max_r][col])) max_r=i;
            }
            if(max_r!=k){// 与第k行交换.
                for(j=k;j<var+1;j++) swap(a[k][j],a[max_r][j]);
            }
            if(a[k][col]==0){// 说明该col列第k行以下全是0了，则处理当前行的下一列.
                k--;
                continue;
            }
            for(i=k+1;i<equ;i++){// 枚举要删去的行.
                if(a[i][col]!=0){
                    LCM = lcm(abs(a[i][col]),abs(a[k][col]));
                    ta = LCM/abs(a[i][col]);
                    tb = LCM/abs(a[k][col]);
                    if(a[i][col]*a[k][col]<0)tb=-tb;//异号的情况是相加
                    for(j=col;j<var+1;j++){
                        a[i][j] = a[i][j]*ta-a[k][j]*tb;
                    }
                }
            }
        }
        // 1. 无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0).
        for (i = k; i < equ; i++){ // 对于无穷解来说，如果要判断哪些是自由变元，那么初等行变换中的交换就会影响，则要记录交换.
            if (a[i][col] != 0) return -1;
        }
        // 2. 无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行，即说明没有形成严格的上三角阵.
        // 且出现的行数即为自由变元的个数.
        if (k < var){
            return var - k; // 自由变元有var - k个.
        }
        // 3. 唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵.
        // 计算出Xn-1, Xn-2 ... X0.
        for (i = var - 1; i >= 0; i--){
            temp = a[i][var];
            for (j = i + 1; j < var; j++){
                if (a[i][j] != 0) temp -= a[i][j] * x[j];
            }
            if (temp % a[i][i] != 0) return -2; // 说明有浮点数解，但无整数解.
            x[i] = temp / a[i][i];
        }
        return 0;
    }
    int main(void){
    //    freopen("in.txt", "r", stdin);
    //    freopen("out.txt","w",stdout);
        int i, j;
        int equ,var;
        while (scanf("%d %d", &equ, &var) != EOF){
            memset(a, 0, sizeof(a));
            for (i = 0; i < equ; i++){
                for (j = 0; j < var + 1; j++){
                    scanf("%d", &a[i][j]);
                }
            }
            int free_num = Gauss(equ,var);
            if (free_num == -1) printf("无解!\n");
            else if (free_num == -2) printf("有浮点数解，无整数解!\n");
            else if (free_num > 0){
                printf("无穷多解! 自由变元个数为%d\n", free_num);
                for (i = 0; i < var; i++){
                    if (free_x[i]) printf("x%d 是不确定的\n", i + 1);
                    else printf("x%d: %d\n", i + 1, x[i]);
                }
            }else{
                for (i = 0; i < var; i++){
                    printf("x%d: %d\n", i + 1, x[i]);
                }
            }
            printf("\n");
        }
        return 0;
    }

转自kuangbin大佬

更多高斯消元模版

2，LU分解法

该方法思想就是将一个矩阵分解为一个单位下三角矩阵L和一个上三角矩阵U，属于矩阵的三角分解法，又称杜利特尔（Doolittle）分解。其实高斯消元法的进行的每一步消元操作，都是等于左乘一个初等置换阵。

首先进行A=LU分解，按照分解的方法来，U的第一行和L的第一列算是初始化，可以直接用循环简单计算赋值，接下来计算剩下的行和列。这里有个计算的次序问题，必须先算U的一行后，才能计算L的一列，不然是错误的计算方法，还有一个注意点就是循环变量的值域问题。LU分解完后，就可以进行直接的求解了，先按LY=B，求解Y，然后按UX=Y，求解X。即求AX=B;先将A=LU,然后LUX=B,先LY=B,求出Y，然后UX=Y，求出X。

#include<iostream>
using namespace std;

const int n = 3;
//矩阵的ALU分解
void ALU(double a[n][n], double b[n])
{
    double l[n][n] = { 0 };
    double u[n][n] = { 0 };
    int i, r, k;
    //进行U的第一行的赋值
    for (i = 0; i<n; i++)
    {
        u[0][i] = a[0][i];
    }

    //进行L的第一列的赋值
    for (i = 1; i<n; i++)
    {
        l[i][0] = a[i][0] / u[0][0];
    }

    //计算U的剩下的行数和L的剩下的列数
    for (r = 1; r<n; r++)
    {
        for (i = r; i <n; i++)
        {
            double sum1 = 0;
            for (k = 0; k < r; k++)
            {
                sum1 += l[r][k] * u[k][i];
                //cout << "" << r << "" << sum1 << endl;
            }
            u[r][i] = a[r][i] - sum1;
        }


        if(r!=n)
        for(i=r+1;i<n;i++)
        {
            double sum2 = 0;
              for (k = 0; k<r; k++)
            {
                  sum2 += l[i][k] * u[k][r];
            }
                l[i][r] = (a[i][r] - sum2) / u[r][r];
        }

    }

    double y[n] = { 0 };
    y[0] = b[0];
    for (i = 1; i<n; i++)
    {
        double sum3 = 0;
        for (k = 0; k<i; k++)
            sum3 += l[i][k] * y[k];
        y[i] = b[i] - sum3;
    }

    double x[n] = { 0 };
    x[n - 1] = y[n - 1] / u[n - 1][n - 1];
    for (i = n - 2; i >= 0; i--)
    {
        double sum4 = 0;
        for (k = i + 1; k<n; k++)
            sum4 += u[i][k] * x[k];
        x[i] = (y[i] - sum4) / u[i][i];
    }
    for (i = 0; i<n; i++)
        cout << "x[" << i + 1 << "]=" << x[i] << endl;
    return;
}


int main()
{
    double a[3][3] = { 1,2,3,2,5,2,3,1,5 };
    double b[3] = { 14,18,20 };
    ALU(a, b);
    return 0;
}

 四，Matrix-Tree定理

https://blog.youkuaiyun.com/qq_33929112/article/details/53812196