数值分析与优化(一)——高斯消元法

基本原理

矩阵$A=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right]$
$A$做矩阵的行初等变换，使其成为上三角阵，之后即可得到与$A$等价(秩)的上三角阵$B$。

高斯消元法分为归一化和消元两部分，归一化是在每一遍消元前事先把需要用于消元的元素变成1，之后就对对应列进行消元。高斯消去法的流程是
对于正整数k=1,2,…,n，执行
(1)如果$a_{kk}^{(k)}=0$则算法失效结束计算，否则按(2)计算
(2)对于i=k+1,k+2,…,n计算
$$\begin{gathered} m_{ik}=a_{ik}^{(k)}/a_{kk}^{(k)} \\ {a}_{ij}^{(k+1)}=a_{ij}^{(k)}-m_{ik}{a}_{kj}^{(k)} \end{gathered}$$
这里存在两个问题，一个是就是如果矩阵不是满秩的，那么计算过程中可能会出现$a_{kk}^{(k)}=0$这时计算可能无法进行下去，如果这是求解线性方则程组或者逆阵这时就可以停止计算，但是如果要求矩阵的秩则需要让计算继续下去，用于确定秩的值。另外由于每次都是以对角线上的元素做归一化，但是当对角线上的元素的绝对值很小时很容易放大计算中的舍入误差。

为了保证计算误差能更小，并且计算可以进行下去需要对计算过程进行改进。每一次归一化之前都先将$a_{kk}^{(k)}$同一列下面的元素$a_{ik}^{(k)}$($i\ge k$)中找到绝对值最大的元素，将其和第k行交换之后再进行归一化和消去，这样就得到了列主元的高斯消去法。而对于求秩这类问题，在列主元的高斯消去时，如果遇到某一列中的主元为0，那么直接将k+1，继续计算，并统计下跳过的次数，直到最后完成整个消去过程，即可得到一个上三角阵，以及线性相关的行的数量。在列主元高斯消去的基础上还可以得到全主元的高斯消去，即比较时将$a_{kk}^{(k)}$右下角的元素全部比较一次，找到其中最大的，然后利用初等行变换和列变换将最大元素换到$a_{kk}^{(k)}$位置上，然后消去。这三种方法中最原始的顺序高斯消去法完全不需要任何数据间的比较，因此计算量最小，而且对应大规模矩阵很容易获得高的Cache利用率，但是稳定较差，当遇到对角线上元素很小时计算不稳定，会产生较大的误差。而列主元虽然多了比较的过程，但由于总是能找到较大的数作为对角线元素用来归一化，因此不容易产生大的误差，相应的计算量也会增加，尤其是矩阵规模增大之后，计算量的差距将会较为明显。而全选主元的高斯消去法由于总是寻找最大的可用值进行归一化稳定性最好，但是计算量也最大，而且高斯消元法作为矩阵计算中最基本的算法，在矩阵求逆、求解线性方程组、矩阵求秩等问题中都有应用，对于求解线性方程组这种问题列变换会使得计算结果发生一些变化，因此在使用列变换时需要记录下所有的列变换，最后在复原，而求等价上三角阵时列变换也不改变矩阵的秩，因此不需要记录列变换，这就需要根据要解决的具体问题进行分析设计。

这里举一个列主元消去的简单例子
$A=\left[\begin{array}{ccc}2& 4& 6\\ 1& 2& 7\\ 3& 6& 12\end{array}\right]$
找到第一列绝对值最大的元素，即第三行第一列，然后将第三行与第一行互换，并将互换后新的$a_{11}$归一化得到
$\left[\begin{array}{ccc}1& 2& 4\\ 1& 2& 7\\ 2& 4& 6\end{array}\right]$
然后利用第一行消去后两行的第一列元素得到
$\left[\begin{array}{ccc}1& 2& 4\\ 0& 0& 3\\ 0& 0& -2\end{array}\right]$
这是第一列完成了消元，然后开始对第二列进行消元，找到从第二行开始到最后一行中第二列绝对值最大的元素，是0，这时直接跳过，进行第三列消元。由于矩阵只有三列因此这一步可以不用继续消元，但是仍然要确定一次主元并且归一化一次，于是得到主元是-2，归一化后的上三角阵为
$\left[\begin{array}{ccc}1& 2& 4\\ 0& 0& 3\\ 0& 0& 1\end{array}\right]$
由于消去过程中仅出现了一次主元是0，因此矩阵的秩是$3-1=2$，这就是高斯消元法的基本原理

基本实现

首先要实现矩阵，这里因为是做示例，为了好解释，用一个自定义的matrix类来实现，行优先存储。

class matrix
{
    double* data;
    int _row,_col;
    public:
    matrix():data(nullptr){};
    matrix(int m,int n)
    {
        _row=m;
        _col=n;
        data=(double*)malloc(sizeof(double)*_row*_col);
    }
    matrix(const matrix& A)
    {
        _row=A._row;
        _col=A._col;
        data=(double*)malloc(sizeof(double)*_row*_col);
    }
    matrix(matrix&& A)
    {
        _row=A._row;
        _col=A._col;
        data=A.data;
        A.data=nullptr;
    }
    ~matrix()
    {
        free(data);
    }
    matrix& operator=(const matrix& A)
    {
        if(data!=nullptr) free(data);
        _row=A._row;
        _col=A._col;
        data=(double*)malloc(sizeof(double)*_row*_col);
        return *this;
    }
    matrix& operator=(matrix&& A)
    {
        data=A.data;
        A.data=nullptr;
        _row=A._row;
        _col=A._col;
        return *this;
    }
    double& operator()(int m,int n)
    {
        return data[m*_col+n];
    }
    const double& operator()(int m,int n) const
    {
        return data[m*_col+n];
    }
    double* operator[](int m)
    {
        return data+m*_col;
    }
    const double* operator[](int m) const
    {
        return data+m*_col;
    }
    int row()
    {
        return _row;
    }
    int col()
    {
        return _col;
    }
    int row() const
    {
        return _row;
    }
    int col() const
    {
        return _col;
    }
};
ostream& operator<<(ostream& out,const matrix& A)
{
    for(int i=0;i<A.row();i++)
    {
        for(int j=0;j<A.col();j++)
        {
            out<<A[i][j]<<'\t';
        }
        out<<'\n';
    }
    return out;
}

列主元高斯消去

由于顺序高斯消去的稳定性较差而且实现上较为简单，因此这里就不做了。这里先实现两种选主元的高斯消去中相对简单的列主元高斯消去。要完成列主元高斯消去需要实现矩阵的初等行变换，矩阵总共有三个初等行变换：整行乘/除同一个不为零常数、交换两行、将一行乘/除上一个常数加/减到另外一行上，分别记为op1、op2、op3。这里顺便加上全选主元需要的列交换op4。

void op1(matrix& A,int i,double k)
{
    if(k==0) return;
    for(int j=0;j<A.col();j++)
    {
        A[i][j]*=k;
    }
}
void op2(matrix& A,int i,int j)
{
    if(i==j) return;
    for(int k=0;k<A.col();k++)
    {
        double tmp=A[i][k];
        A[i][k]=A[j][k];
        A[j][k]=tmp;
    }
}
void op3(matrix& A,int i,int j,double k)
{
    if(k==0) return;
    for(int l=0;l<A.col();l++)
    {
        A[j][l]+=A[i][l]*k;
    }
}
void op4(matrix& A,int i,int j)
{
    if(i==j) return;
    for(int k=0;k<A.row();k++)
    {
        double tmp=A[k][i];
        A[k][i]=A[k][j];
        A[k][j]=tmp;
    }
}

然后就可以实现列主元高斯消去

void Gauss(matrix& A)
{
    for(int j=0;j<A.col();j++)
    {
        int max_r=j;
        double max=abs(A[j][j]);
        for(int i=j+1;i<A.row();i++)
        {
            if(max<abs(A[i][j])) 
            {
                max=abs(A[i][j]);
                max_r=i;
            }
        }
        op2(A,max_r,j);
        if(abs(max)<1e-12) 
        {
            continue;
        }
        op1(A,j,1/A[j][j]);
        for(int i=j+1;i<A.row();i++) op3(A,j,i,-A[i][j]);
    }
}

先用一个三阶方阵简单测试一下

int main()
{
    matrix A(3,3);
    for(int i=0;i<3;i++)
    {
        for(int j=0;j<3;j++) A(i,j)=i*3+j;
    }
    A(1,0)=0;
    A(1,1)=2;
    A(1,2)=4;
    cout<<A<<'\n';
    Gauss(A);
    cout<<A<<'\n';
}

计算结果如下

0	1	2	
0	2	4	
6	7	8	

1	1.16667	1.33333	
0	1	2	
0	0	0	

可以看到计算正确，得到了对应的上三角阵。

全选主元高斯消去

全选主元高斯消去需要列变换，因此需要在前面初等行变换op1、op2、op3的基础上在加上列交换op4，具体列交换代码已经在上面展示了这里直接展示高斯消去的代码。

void Gauss(matrix& A)
{
    for(int j=0;j<A.col();j++)
    {
        int max_r=j,max_c=j;
        double max=abs(A[j][j]);
        for(int i=j;i<A.row();i++)
        {
            for(int k=j;k<A.col();k++)
            {
                if(max<abs(A[i][k])) 
                {
                    max=abs(A[i][k]);
                    max_r=i;
                    max_c=k;
                }
            }
        }
        op2(A,max_r,j);
        op4(A,max_c,j);
        if(abs(max)<1e-12) 
        {
            return;
        }
        op1(A,j,1/A[j][j]);
        for(int i=j+1;i<A.row();i++) op3(A,j,i,-A[i][j]);
    }
}

测试一下

int main()
{
    matrix A(3,3);
    for(int i=0;i<3;i++)
    {
        for(int j=0;j<3;j++) A(i,j)=i*3+j;
    }
    A(1,0)=0;
    A(1,1)=2;
    A(1,2)=4;
    cout<<A<<'\n';
    Gauss(A);
    cout<<A<<'\n';
}

结果为

0	1	2	
0	2	4	
6	7	8	

1	0.875	0.75	
-0	1	2	
0	0	0

同样得到了等价的对角线元素为1的上三角阵。由于这里使用的列交换，因此对于求秩、求行列式值这类问题可以得到正确结果，但是对于求解线性方程组和求逆矩阵时就存在问题，需要记录下列交换的具体情况，然后再想办法还原回正确结果。

性能优化

这里对全选主元的高斯消去法做一些分析和优化。

在原始版本中按照计算顺序，首先寻找主元，然后将主元交换到$a_{kk}^{(k)}$的位置，再进行归一化和消元，然后再循环。这个过程中我们在寻找主元时需要将所有未完成消元的元素都搜索一遍，这个过程需要读取一次全部剩余元素，然后交换一次，交换过程只需要交换对应行和列，因此只用读写两行两列，最后再归一化，而归一化的过程中则有需要对所有剩余元素做一次计算，这就又要读写一次全部剩余元素。那么读写的数据量就是两次剩余元素和两行两列的数据。

这里我们尝试将两次全部剩余元素的读取合并为一次，可以看到在消元时，我们会读写全部剩余元素，以消去当前列，而消元完成后我们紧接着就会扫描全部剩余元素以找到主元，于是我们直接将两步合并，在消元的同时每计算完一个元素就比较一次从中找到最大的元素，这样我们就可以省掉一次读取的时间。

进一步由于列交换时，列读取是非连续访存，效率比较低，考虑是否能够将列交换如果我们找到主元素之后不交换列，而是就地用主元素所在的列归一化和消元，由于消元时是一行一行进行的，在完成一行消元时这一行的数据刚使用过，在这时进行列交换就能一定程度上减少列读取开销，但是另一方面由于消去过程仅对归一化那一行的下面行进行消去，因此这个过程做列交换只能交换对应列的部分元素，要完成完整的列交换还需要将归一化那一行上面的行进行列交换，这个过程中的非连续访存是不可避免的，不过这种合并操作仍能节约一定的IO开销。

最后我们观察每完成一次归一化、消元的过程，被归一化的行就不会再参与计算，而归一化对角线元素以下的值都会因为消元而变成零。换句话说一旦完成一次归一化和、消元剩余需要计算的元素就是其右下角的子矩阵，相当于计算元素减少了一行和一列，因此在做行初等变换时可以考虑到这一因素，对于交换前后均为0的元素就不再交换了，而对于会被消去的元素由于不再具有实际意义，因此可以不进行任何操作，这样可以节约一定的计算量。

通过修改矩阵初等变换的函数，以及高斯消去法的函数，得到的优化后的全选主元的高斯消去法。

void op1(matrix& A,int i,double k,int j,int c1,int c2)
{
    if(k==0) return;
    for(;j<A.col();j++)
    {
        A[i][j]*=k;
    }
    double tmp=A[i][c1];
    A[i][c1]=A[i][c2];
    A[i][c2]=tmp;
}
void op1(matrix& A,int i,double k,int j)
{
    if(k==0) return;
    for(;j<A.col();j++)
    {
        A[i][j]*=k;
    }
}
void op2(matrix& A,int i,int j,int k)
{
    if(i==j) return;
    for(;k<A.col();k++)
    {
        double tmp=A[i][k];
        A[i][k]=A[j][k];
        A[j][k]=tmp;
    }
}
void op3(matrix& A,int i,int j,double k,int l,int c1,int c2,std::tuple<int,int,double> &max)
{
    if(k==0) 
    {
        l++;
        double tmp=A[j][c1];
        A[j][c1]=A[j][c2];
        A[j][c2]=tmp;
        for(;l<A.col();l++)
        {
            if(abs(A[j][l])>std::get<2>(max))
            {
                std::get<0>(max)=j;
                std::get<1>(max)=l;
                std::get<2>(max)=abs(A[j][l]);
            }
        }
        return;
    }
    
    double tmp=A[j][c1];
    A[j][c1]=A[j][c2];
    A[j][c2]=tmp;

    //A[j][l]=0;
    l++;
    for(;l<A.col();l++)
    {
        A[j][l]+=A[i][l]*k;
        if(abs(A[j][l])>std::get<2>(max))
        {
            std::get<0>(max)=j;
            std::get<1>(max)=l;
            std::get<2>(max)=abs(A[j][l]);
        }
    }
    
}
void op3(matrix& A,int i,int j,double k,int l,std::tuple<int,int,double> &max)
{
    if(k==0) 
    {
        l++;
        for(;l<A.col();l++)
        {
            if(abs(A[j][l])>std::get<2>(max))
            {
                std::get<0>(max)=j;
                std::get<1>(max)=l;
                std::get<2>(max)=abs(A[j][l]);
            }
        }
        return;
    }
    //A[j][l]=0;
    l++;
    for(;l<A.col();l++)
    {
        A[j][l]+=A[i][l]*k;
        if(abs(A[j][l])>std::get<2>(max))
        {
            std::get<0>(max)=j;
            std::get<1>(max)=l;
            std::get<2>(max)=abs(A[j][l]);
        }
    }
    
}
void op4(matrix& A,int i,int j,int end)
{
    if(i==j) return;
    for(int k=0;k<end+1;k++)
    {
        double tmp=A[k][i];
        A[k][i]=A[k][j];
        A[k][j]=tmp;
    }
}
void Gauss(matrix& A)
{
    tuple<int,int,double> max(0,0,abs(A[0][0]));
    for(int i=0;i<A.row();i++)
    {
        for(int j=0;j<A.col();j++)
        {
            if(abs(A[i][j])>get<2>(max))
            {
                get<0>(max)=i;
                get<1>(max)=j;
                get<2>(max)=abs(A[i][j]);
            }
        }
    }
    for(int j=0;j<A.col();j++)
    {
        int max_r=get<0>(max),max_c=get<1>(max);
        op2(A,max_r,j,j);
        if(get<2>(max)<1e-12) 
        {
            return;
        }
        get<0>(max)=j+1;
        get<1>(max)=j+1;
        get<2>(max)=0;
        if(j!=max_c) 
        {
            op1(A,j,1/A[j][max_c],j,max_c,j);
            for(int i=j+1;i<A.row();i++) op3(A,j,i,-A[i][max_c],j,max_c,j,max);
            op4(A,max_c,j,j);
        }
        else
        {
            op1(A,j,1/A[j][max_c],j);
            for(int i=j+1;i<A.row();i++) op3(A,j,i,-A[i][max_c],j,max);
        } 
        
    }
}

这里用1200阶方阵测试一下

#define N 1200
int main()
{
    matrix A(N,N);
    for(int i=0;i<N;i++)
    {
        for(int j=0;j<N;j++) 
        {
            A(i,j)=rand()%N;
        }
    }
    clock_t s=clock();
    Gauss(A);
    cout<<"用时："<<double(clock()-s)/CLOCKS_PER_SEC<<"s\n";
    cout<<A[1][N-1]<<'\n';
}

使用gcc -O2优化后
原始的高斯消去法用时为

1.75671s

而改进后的高斯消去法用时为

1.2576s

可以看到计算效率得到了一定的提升，加速比大约是1.39，并没有想象中的提升明显。我在网上复制了一份其他人的全选主元的高斯消去代码，那份代码中并没有对操作进行合并，仅仅是在消去时那些被消为零的元素不做任何计算，其性能几乎和我改进后的代码相同，仅仅比我改进的代码慢了1%左右，几乎可以认为这个差距是测试的误差。这可能是因为高斯消去法都是顺序访存而且Cache利用率有限，所以并不能有像矩阵相乘那样通过一些操作顺序等方面的改进获得数倍的性能提升。

SIMD加速

没办法对访存进行优化，那只能使用SIMD加速高斯消去法了。利用SIMD来完成各种行变换和搜索主元完成如下的op函数和一些辅助函数。而高斯消去的函数则和之前改进的全选主元高斯消去函数完全相同。

void vec_max(const __m256d  &data,__m256d &max,const __m256i &pos,__m256i &max_pos)
{
    __m256d r=_mm256_cmp_pd(data,max,14);
    max=_mm256_max_pd(data,max);
    max_pos=_mm256_castpd_si256(_mm256_blendv_pd(_mm256_castsi256_pd(max_pos),_mm256_castsi256_pd(pos),r));
}
__m256d vec_abs(const __m256d& t1)
{
    __m256i t=_mm256_castpd_si256(t1);
    t=_mm256_and_si256(t,_mm256_set1_epi64x(0x7fffffffffffffff));
    return _mm256_castsi256_pd(t);
}
void op1(matrix& A,int i,double k,int j,int c1,int c2)
{
    if(k==0) return;
    double *p=A[i]+j;
    __m256d vk=_mm256_set1_pd(k);
    while(j+4<A.col())
    {
        __m256d t=_mm256_loadu_pd(p);
        t=_mm256_mul_pd(t,vk);
        _mm256_storeu_pd(p,t);
        p+=4;
        j+=4;
    }
    for(;j<A.col();j++)
    {
        A[i][j]*=k;
    }
    double tmp=A[i][c1];
    A[i][c1]=A[i][c2];
    A[i][c2]=tmp;
}
void op2(matrix& A,int i,int j,int k)
{
    if(i==j) return;
    for(;k<A.col();k++)
    {
        double tmp=A[i][k];
        A[i][k]=A[j][k];
        A[j][k]=tmp;
    }
}
void op3(matrix& A,int i,int j,double k,int l,int c1,int c2,std::tuple<int,int,double> &max)
{
    if(k==0) 
    {
        l++;
        double tmp=A[j][c1];
        A[j][c1]=A[j][c2];
        A[j][c2]=tmp;
        for(;l<A.col();l++)
        {
            if(abs(A[j][l])>std::get<2>(max))
            {
                std::get<0>(max)=j;
                std::get<1>(max)=l;
                std::get<2>(max)=abs(A[j][l]);
            }
        }
        return;
    }
    
    double tmp=A[j][c1];
    A[j][c1]=A[j][c2];
    A[j][c2]=tmp;
    A[j][l]=0;
    l+=1;
    __m256d vk=_mm256_set1_pd(k),M=_mm256_set1_pd(0);
    __m256i pos=_mm256_set_epi64x(l+3,l+2,l+1,l),max_pos=_mm256_set1_epi64x(0);
    double *p=A[j]+l,*q=A[i]+l;
    while(l+4<A.col())
    {
        __m256d t=_mm256_loadu_pd(p),a=_mm256_loadu_pd(q),rc;
        t=_mm256_fmadd_pd(a,vk,t);
        _mm256_storeu_pd(p,t);
        vec_max(vec_abs(t),M,pos,max_pos);
        pos=_mm256_add_epi64(pos,_mm256_set1_epi64x(4));
        p+=4;
        q+=4;
        l+=4;
    }
    double Max[4];
    long long int Max_pos[4];
    _mm256_storeu_pd(Max,M);
    _mm256_maskstore_epi64(Max_pos,_mm256_set1_epi64x(-1),max_pos);
    for(int i=0;i<4;i++)
    {
        if(Max[i]>std::get<2>(max))
        {
            std::get<0>(max)=j;
            std::get<1>(max)=Max_pos[i];
            std::get<2>(max)=Max[i];
        }
    }
    for(;l<A.col();l++)
    {
        A[j][l]+=A[i][l]*k;
        if(abs(A[j][l])>std::get<2>(max))
        {
            std::get<0>(max)=j;
            std::get<1>(max)=l;
            std::get<2>(max)=abs(A[j][l]);
        }
    }
}
void op4(matrix& A,int i,int j,int end)
{
    if(i==j) return;
    for(int k=0;k<end;k++)
    {
        double tmp=A[k][i];
        A[k][i]=A[k][j];
        A[k][j]=tmp;
    }
}
void Gauss(matrix& A)
{
    tuple<int,int,double> max(0,0,abs(A[0][0]));
    for(int i=0;i<A.row();i++)
    {
        for(int j=0;j<A.col();j++)
        {
            if(abs(A[i][j])>get<2>(max))
            {
                get<0>(max)=i;
                get<1>(max)=j;
                get<2>(max)=abs(A[i][j]);
            }
        }
    }
    for(int j=0;j<A.col();j++)
    {
        int max_r=get<0>(max),max_c=get<1>(max);
        op2(A,max_r,j,j);
        if(get<2>(max)<1e-12) 
        {
            return;
        }
        get<0>(max)=j+1;
        get<1>(max)=j+1;
        get<2>(max)=0;
        op1(A,j,1/A[j][max_c],j,max_c,j);
        for(int i=j+1;i<A.row();i++) op3(A,j,i,-A[i][max_c],j,max_c,j,max);
        op4(A,max_c,j,j);
    }
}

同样使用O2优化，1200阶方阵用时为

0.842909s

可以看到相较于前面优化后的1.25s的用时，加速比约为1.45，而与原始版本相比则达到2.07的加速比，快了约一倍。理论上使用AVX一次可以处理四个浮点数，理论加速比可以达到4。但是这里性能提升大约50%不到，我在调试SIMD程序时发现，其实所有的行变换操作的计算SIMD几乎没有任何加速，这可能是因为行变换都是顺序操作，因此编译器早已优化了。SIMD真正提高性能的是主元搜索过程，如果不使用SIMD加速主元搜索，性能不会有任何提升，这样可以推测搜索主元的操作过程在1200阶方阵1.25s的计算耗时中，可能大约占用了0.5s。

另外，由于高斯消去法的计算特点，这里使用SIMD读写数据时都用的是loadu和storeu两个操作，以保证内存不对齐时能正确读写数据，但是这两个操作的速度是慢于对齐内存读写的load和store操作的。可以通过先算一部分内存未对齐的数据然后从对齐位置开始使用AVX进行计算，也许这样可以再快一点，不过这个稍微有点麻烦，这里就不再做了。

紧凑存储

整个高斯消元的计算过程，在计算快要结束时每一行实际操作的数据已经非常少了，但是由于有大量的零元素在占据着内存，因此实际上消去时有用的数据并不是很紧凑。于是将每一列消去时将0元素从矩阵中删去，最终得到紧凑存储的上三角阵，即起始n个元素是矩阵右上角的第一行元素，紧接着n-1个元素是第二行的元素，以此类推，最后一行元素只占用一个存储空间。这样整个矩阵的数据更加紧凑，在一定程度上可以加速矩阵数据的读取，从而得到性能的提升。

void norm(double* row,int remain_cols,int max_col)
{
    double tmp=row[max_col];
    row[max_col]=row[0];
    row[0]=tmp;
    if(row[0]==0) return;
    double k=1.0/row[0];
    for(int i=0;i<remain_cols;i++)
    {
        row[i]*=k;
    }
    
}
void swap_row(double* row_i,double* row_j,int remain_cols)
{
    if(row_i==row_j) return;
    for(int i=0;i<remain_cols;i++)
    {
        double tmp=*row_i;
        *row_i=*row_j;
        *row_j=tmp;
        row_i++;
        row_j++;
    }
}
void elim(double* row_i,double* row_j,int row_num,int remain_cols,int move_dist,int max_col,std::tuple<int,int,double> &max)
{
    double k=-row_j[max_col];
    row_j[max_col]=row_j[0];
    if(k==0) 
    {
        for(int i=1;i<remain_cols;i++)
        {
            row_j[i-move_dist]=row_j[i];
            if(abs(row_j[i-move_dist])>std::get<2>(max))
            {
                std::get<0>(max)=row_num;
                std::get<1>(max)=i-1;
                std::get<2>(max)=abs(row_j[i-move_dist]);
            }
        }
        return;
    }
    for(int i=1;i<remain_cols;i++)
    {
        row_j[i-move_dist]=row_j[i]+row_i[i]*k;
        if(abs(row_j[i-move_dist])>std::get<2>(max))
        {
            std::get<0>(max)=row_num;
            std::get<1>(max)=i-1;
            std::get<2>(max)=abs(row_j[i-move_dist]);
        }
    }
}
void swap_col(double *p,int i,int j,int n,int end)
{
    if(i==j) return;
    int d=j-i;
    p=p+i;
    for(int k=0;k<end;k++)
    {
        double tmp=*p;
        *p=*(p+d);
        *(p+d)=tmp;
        p+=n-k-1;
    }
}
void up_show(const matrix &A,int k)
{
    const double *_A=A[0];
    for(int i=0;i<A.row();i++)
    {
        for(int j=0;j<A.col();j++)
        {
            if(i>=k) 
            {
                if(j<k) cout<<0<<'\t';
                else 
                {
                    if(fabs(*_A)<(1e-12)) 
                    {
                        cout<<0<<'\t';
                    }
                    else
                    {
                        cout<<*_A<<'\t';
                    }
                    _A++;
                }
            }
            else
            {
                if(j<i) cout<<0<<'\t';
                else 
                {
                    if(fabs(*_A)<(1e-12)) 
                    {
                        cout<<0<<'\t';
                    }
                    else
                    {
                        cout<<*_A<<'\t';
                    }
                    _A++;
                }
            }
        }
        cout<<'\n';
    }
}
void Gauss(matrix& A)
{
    double *norm_row=A[0],*elim_row=A[0];
    tuple<int,int,double> max(0,0,abs(A[0][0]));
    for(int i=0;i<A.row();i++)
    {
        for(int j=0;j<A.col();j++)
        {
            if(abs(A[i][j])>get<2>(max))
            {
                get<0>(max)=i;
                get<1>(max)=j;
                get<2>(max)=abs(A[i][j]);
            }
        }
    }
    int remain_cols=A.col();
    for(int i=0;i<A.row();i++)
    {
        int max_c=get<1>(max),max_dr=get<0>(max)-i;
        if(get<2>(max)<1e-12) 
        {
            return;
        }
        get<0>(max)=0;
        get<1>(max)=0;
        get<2>(max)=0;
        swap_row(norm_row,norm_row+max_dr*remain_cols,remain_cols);
        norm(norm_row,remain_cols,max_c);
        for(int row_num=i+1;row_num<A.row();row_num++)
        {
            elim_row+=remain_cols;
            elim(norm_row,elim_row,row_num,remain_cols,row_num-i,max_c,max);
        }
        swap_col(A[0],max_c+i,i,N,i);
        norm_row+=remain_cols;
        elim_row=norm_row;
        remain_cols--;
    }
}

同样使用1200阶方程O2优化测试一下，用时为

0.914547s

这里暂时没有用SIMD加速，前面没有SIMD加速的程序用时约为1.25s，有一定的性能提升。

再利用SIMD加速计算和比较过程

void vec_max(const __m256d  &data,__m256d &max,const __m256i &pos,__m256i &max_pos)
{
    __m256d r=_mm256_cmp_pd(data,max,14);
    max=_mm256_max_pd(data,max);
    max_pos=_mm256_castpd_si256(_mm256_blendv_pd(_mm256_castsi256_pd(max_pos),_mm256_castsi256_pd(pos),r));
}
__m256d vec_abs(const __m256d& t1)
{
    __m256i t=_mm256_castpd_si256(t1);
    t=_mm256_and_si256(t,_mm256_set1_epi64x(0x7fffffffffffffff));
    return _mm256_castsi256_pd(t);
}
void get_max(matrix& A,std::tuple<int,int,double> &max)
{
    double* p=A[0];
    for(int i=0;i<A.row();i++)
    {
        __m256d M=_mm256_set1_pd(0);
        __m256i pos=_mm256_set_epi64x(3,2,1,0),max_pos=_mm256_set1_epi64x(0);
        int j=0;
        for(;j+4<A.col();j+=4)
        {
            __m256d t=_mm256_loadu_pd(p);
            vec_max(vec_abs(t),M,pos,max_pos);
            pos=_mm256_add_epi64(pos,_mm256_set1_epi64x(4));
            p+=4;
        }
        double Max[4];
        long long int Max_pos[4];
        _mm256_storeu_pd(Max,M);
        _mm256_maskstore_epi64(Max_pos,_mm256_set1_epi64x(-1),max_pos);
        for(int k=0;k<4;k++)
        {
            if(Max[k]>std::get<2>(max))
            {
                std::get<0>(max)=i;
                std::get<1>(max)=Max_pos[k];
                std::get<2>(max)=Max[k];
            }
        }
        for(;j<A.col();j++)
        {
            if(abs(*p)>std::get<2>(max))
            {
                std::get<0>(max)=i;
                std::get<1>(max)=j;
                std::get<2>(max)=abs(*p);
            }
            p++;
        }
    }
}
void norm(double* row,int remain_cols,int max_col)
{
    double tmp=row[max_col];
    row[max_col]=row[0];
    row[0]=tmp;
    if(row[0]==0) return;
    double k=1.0/row[0];
    for(int i=0;i<remain_cols;i++)
    {
        row[i]*=k;
    }
    
}
void swap_row(double* row_i,double* row_j,int remain_cols)
{
    if(row_i==row_j) return;
    for(int i=0;i<remain_cols;i++)
    {
        double tmp=*row_i;
        *row_i=*row_j;
        *row_j=tmp;
        row_i++;
        row_j++;
    }
}
void elim(double* row_i,double* row_j,int row_num,int remain_cols,int move_dist,int max_col,std::tuple<int,int,double> &max)
{
    double k=-row_j[max_col];
    row_j[max_col]=row_j[0];
    if(k==0) 
    {
        for(int i=1;i<remain_cols;i++)
        {
            row_j[i-move_dist]=row_j[i];
            if(abs(row_j[i-move_dist])>std::get<2>(max))
            {
                std::get<0>(max)=row_num;
                std::get<1>(max)=i-1;
                std::get<2>(max)=abs(row_j[i-move_dist]);
            }
        }
        return;
    }
    double *ri=row_i+1,*rj=row_j+1;
    __m256d vk=_mm256_set1_pd(k),M=_mm256_set1_pd(0);
    __m256i pos=_mm256_set_epi64x(3,2,1,0),max_pos=_mm256_set1_epi64x(0);
    int i=1;
    for(;i+4<remain_cols;i+=4)
    {
        __m256d t=_mm256_loadu_pd(rj),a=_mm256_loadu_pd(ri);
        t=_mm256_fmadd_pd(a,vk,t);
        _mm256_storeu_pd(rj-move_dist,t);
        vec_max(vec_abs(t),M,pos,max_pos);
        pos=_mm256_add_epi64(pos,_mm256_set1_epi64x(4));
        ri+=4;
        rj+=4;
        /*
        row_j[i-move_dist]=row_j[i]+row_i[i]*k;
        if(abs(row_j[i-move_dist])>std::get<2>(max))
        {
            std::get<0>(max)=row_num;
            std::get<1>(max)=i-1;
            std::get<2>(max)=abs(row_j[i-move_dist]);
        }
        */
    }
    double Max[4];
    long long int Max_pos[4];
    _mm256_storeu_pd(Max,M);
    _mm256_maskstore_epi64(Max_pos,_mm256_set1_epi64x(-1),max_pos);
    for(int l=0;l<4;l++)
    {
        if(Max[l]>std::get<2>(max))
        {
            std::get<0>(max)=row_num;
            std::get<1>(max)=Max_pos[l];
            std::get<2>(max)=Max[l];
        }
    }
    for(;i<remain_cols;i++)
    {
        row_j[i-move_dist]=row_j[i]+row_i[i]*k;
        if(abs(row_j[i-move_dist])>std::get<2>(max))
        {
            std::get<0>(max)=row_num;
            std::get<1>(max)=i-1;
            std::get<2>(max)=abs(row_j[i-move_dist]);
        }
    }
}
void swap_col(double *p,int i,int j,int n,int end)
{
    if(i==j) return;
    int d=j-i;
    p=p+i;
    for(int k=0;k<end;k++)
    {
        double tmp=*p;
        *p=*(p+d);
        *(p+d)=tmp;
        p+=n-k-1;
    }
}
void Gauss(matrix& A)
{
    double *norm_row=A[0],*elim_row=A[0];
    tuple<int,int,double> max(0,0,0);
    get_max(A,max);
    int remain_cols=A.col();
    for(int i=0;i<A.row();i++)
    {
        int max_c=get<1>(max),max_dr=get<0>(max)-i;
        if(get<2>(max)<1e-12) 
        {
            return;
        }
        get<0>(max)=0;
        get<1>(max)=0;
        get<2>(max)=0;
        swap_row(norm_row,norm_row+max_dr*remain_cols,remain_cols);
        norm(norm_row,remain_cols,max_c);
        for(int row_num=i+1;row_num<A.row();row_num++)
        {
            elim_row+=remain_cols;
            elim(norm_row,elim_row,row_num,remain_cols,row_num-i,max_c,max);
        }
        swap_col(A[0],max_c+i,i,N,i);
        norm_row+=remain_cols;
        elim_row=norm_row;
        remain_cols--;
    }
}

同样O2优化，1200阶方阵，用时为

0.784726s

对比前面SIMD加速的计算代码，同样有一定的性能提升。

列优先存储

高斯消元法以行变换为主，列变换仅在交换主元位置时才需要用到，因此这个算法对于矩阵本身是行优先存储还是列优先存储性能上必然存在差异。这里测试一些列优先存储时高斯消元法的性能，首先是不考虑紧凑存储，不使用SIMD情况下列优先存储的全选主元高斯消去法的测试。

void copy_scale(double* col,double* dst,int n)
{
    for(int i=0;i<n;i++) *(dst++)=*(col++);
}
void swap_scale(double *scale,int r1,int r2)
{
    double tmp=scale[r1];
    scale[r1]=scale[r2];
    scale[r2]=tmp;
}
void elim(double* src,double *dst,double* scale,int n,int rows,int col_num,int r1,int r2,std::tuple<int,int,double> &max)
{
    int i=n;
    double tmp=src[r1];
    src[r1]=src[r2];
    src[r2]=tmp;
    tmp=src[i]=src[i]/scale[i];
    i+=1;
    for(;i<rows;i++)
    {
        dst[i]=src[i]+scale[i]*tmp;
        if(abs(dst[i])>std::get<2>(max))
        {
            std::get<0>(max)=i;
            std::get<1>(max)=col_num;
            std::get<2>(max)=abs(dst[i]);
        }
    }
}
void swap_col(double *col_i,double *col_j,int rows)
{
    if(col_i==col_j) return;
    for(int i=0;i<rows;i++)
    {
        double tmp=*col_i;
        *col_i=*col_j;
        *col_j=tmp;
        col_i++;
        col_j++;
    }
}
void Gauss(matrix& A)
{
    tuple<int,int,double> max(0,0,abs(A[0][0]));
    for(int i=0;i<A.col();i++)
    {
        for(int j=0;j<A.row();j++)
        {
            if(abs(A[i][j])>get<2>(max))
            {
                get<0>(max)=j;
                get<1>(max)=i;
                get<2>(max)=abs(A[i][j]);
            }
        }
    }
    for(int j=0;j<A.col();j++)
    {
        int max_r=get<0>(max),max_c=get<1>(max);
        if(get<2>(max)<1e-12) 
        {
            return;
        }
        get<0>(max)=j+1;
        get<1>(max)=j+1;
        get<2>(max)=0;
        swap_col(A[j],A[max_c],A.row());
        swap_scale(A[j],max_r,j);
        for(int i=j+1;i<A.col();i++) elim(A[i],A[i],A[j],j,A.row(),i,max_r,j,max);
        A[j][j]=1;
        for(int i=j+1;i<A.col();i++) A[j][i]=0;
    }
}

同样使用1200阶方阵，O2优化测试用时为

0.930971s

可以看到在没有紧凑存储和SIMD加速时列优先存储的高斯消元法就已经具有较高的性能了。考虑消去过程顺序操作为主，而交换过程中列交换和行交换必然一个是顺序操作，一个是非连续内存操作。随着高斯消去的进行后续需要参与消去的行中，零元素逐渐变多，需要进行行交换的数据越来越少，行交换的操作次数会逐渐减少，而列交换总是整列进行的，零元素的增加并不会减少列交换的操作次数。对于行优先存储而言，列交换是非连续内存操作，这样计算后期消去过程的执行速度越来越快，而列交换的执行速度没有发生改变，大部分时间都消耗在了列交换过程上。列优先存储消去过程是顺序操作，而列交换变成了顺序操作，行交换变成了非连续内存操作，这样非连续内存操作的次数随着计算的进行会越来越少，而列交换虽然操作量不变但是从原来的行优先的非连续操作变成现在的顺序操作，性能得到提升。