精华内容
下载资源
问答
  • Python实现列主元高斯消去法与LU分解法 数值分析:Python实现列主元高斯消去法与LU分解法求解线性方程组 一、矩阵形式的线性代数方程组 二、高斯消去法 三、高斯列主元消去法 四、矩阵三角分解法(LU分解) 这里...

    Python实现列主元高斯消去法与LU分解法

    数值分析:Python实现列主元高斯消去法与LU分解法求解线性方程组

    一、矩阵形式的线性代数方程组

    在这里插入图片描述

    二、高斯消去法

    在这里插入图片描述
    在这里插入图片描述

    三、高斯列主元消去法

    在这里插入图片描述

    四、矩阵三角分解法(LU分解)

    这里只简单介绍Doolittle分解法。
    在这里插入图片描述
    在这里插入图片描述

    题目:编写列主元高斯消去法与LU分解法解线性方程组Ax=b。

    在这里插入图片描述
    列主元消去法代码实现:

    import math
    import numpy as np
    #目的:熟悉列主元消去法,以及三角分解法等直接求解线性方程组的算法
    #列主元消元法
    def CME(a,b,x):
    
        isdet0 = 0
        m, n = a.shape                  #矩阵a的行数和列数
        # j表示列
        for k in range(n - 1):          # k表示第一层循环,(0,n-1)行
                #在每次计算前,找到最大主元,进行换行
                ans = np.fabs(a[k][k])
                ik = k
                for i in range(k+1, n):
                    if ans < np.fabs(a[i][k]):  # fabs是绝对值,将a中绝
    对值最大的找出来
                        ik = i
                        ans = np.fabs(a[i][k])
    
                if np.fabs(ans) < 1e-10:
                    isdet0 = 1
                    break
                if ik != k :
                    for i in range(k,m):
                        temp = a[k][i]
                        a[k][i] = a[ik][i]
                        a[ik][i] = temp
                    temp = b[k]
                    b[k] = b[ik]
                    b[ik] = temp
    
                for i in range(k + 1, n):   # i表示第二层循环,(k+1,n)行,
    计算该行消元的系数
                    temp = a[i][k] / a[k][k]     #计算
    
                    for j in range(k,m):      # j表示列,对每一列进行运算
                        a[i][j] = a[i][j] - temp * a[k][j]
    
                    b[i] = b[i] - temp * b[k]
        # 回代求出方程解
        if np.fabs(a[n-1][n-1]) < 1e-10 :
            isdet0 = 1
        if isdet0 == 0:
               # x = np.zeros(n)
                x[n - 1] = b[n - 1] / a[n - 1][n - 1] #先算最后一位的x解
                for i in range(n - 2, -1, -1):    #依次回代倒着算每一个解
                    temp = 0
                    for j in range(n - 1, i,-1):
                        temp = temp + a[i][j]*x[j]
    
                    x[i] = (b[i]-temp) / a[i][i]
                for i in range(n):
                    print("x" + str(i + 1) + " = ", x[i])
                print("x" " = ", x)
    if __name__ == '__main__':      #当模块被直接运行时,以下代码块将被运行,当模块是被导入时,代码块不被运行。
        a = np.array([[3.01, 6.03, 1.99], [1.27, 4.16, -1.23], [0.987, -4.81, 9.34]])
        b = np.array([1.0, 1.0, 1.0])
        m,n = a.shape
        x = np.zeros(n)
        B = np.zeros((n, n))
        for i in range(n):
            for j in range(n):
                B[i][j] = a[i][j]
        CME(a,b,x)
        #验证
        for i in range(0, n):
            temp = 0
            for j in range(0, n):
                temp = temp + B[i][j] * x[j]
            print("%f ", temp)
    
    if __name__ == '__main__':
    main()
    
    

    LU分解法代码实现:

    import math
    import numpy as np
    #目的:熟悉列主元消去法,以及三角分解法等直接求解线性方程组的算法
    #列主元消元法
    def CME(a,b,x):
    
        isdet0 = 0
        m, n = a.shape                  #矩阵a的行数和列数
        # j表示列
        for k in range(n - 1):          # k表示第一层循环,(0,n-1)行
                #在每次计算前,找到最大主元,进行换行
                ans = np.fabs(a[k][k])
                ik = k
                for i in range(k+1, n):
                    if ans < np.fabs(a[i][k]):  # fabs是绝对值,将a中绝对值最大的找出来
                        ik = i
                        ans = np.fabs(a[i][k])
    
                if np.fabs(ans) < 1e-10:
                    isdet0 = 1
                    break
                if ik != k :
                    for i in range(k,m):
                        temp = a[k][i]
                        a[k][i] = a[ik][i]
                        a[ik][i] = temp
                    temp = b[k]
                    b[k] = b[ik]
                    b[ik] = temp
    
                for i in range(k + 1, n):   # i表示第二层循环,(k+1,n)行,计算该行消元的系数
                    temp = a[i][k] / a[k][k]     #计算
    
                    for j in range(k,m):      # j表示列,对每一列进行运算
                        a[i][j] = a[i][j] - temp * a[k][j]
    
                    b[i] = b[i] - temp * b[k]
        # 回代求出方程解
        if np.fabs(a[n-1][n-1]) < 1e-10 :
            isdet0 = 1
        if isdet0 == 0:
               # x = np.zeros(n)
                x[n - 1] = b[n - 1] / a[n - 1][n - 1] #先算最后一位的x解
                for i in range(n - 2, -1, -1):      #依次回代倒着算每一个解
                    temp = 0
                    for j in range(n - 1, i,-1):
                        temp = temp + a[i][j]*x[j]
    
                    x[i] = (b[i]-temp) / a[i][i]
                for i in range(n):
                    print("x" + str(i + 1) + " = ", x[i])
                print("x" " = ", x)
    #三角消元法
    def LU(a,b,x):
        m, n = a.shape  # 矩阵a的行数和列数
        y = np.array([0.0, 0.0, 0.0])
        for j in range(1,n):# L的第0列
            a[j][0] = a[j][0] / a[0][0]
        for i in range(1,n-1):# 求U的第i行 L的第i行
            for j in range(i,n):#求U的第i行的第j个元素
                sum = 0.0 #求和
                for s in range(0,i):
                    sum = sum +a[i][s] * a[s][j]
                a[i][j] = a[i][j] - sum
            #求L的第i列的第j个元素  在j行i列
            for j in range(i+1,n):
                sum = 0.0
                for s in range(0,i):
                    sum = sum + a[j][s] * a[s][i]
                a[j][i] = ( a[j][i] - sum ) / a[i][i]
        #求U[n-1][n-1]
        sum = 0.0 #求和
        for s in range(0,n-1):
            sum = sum + a[n-1][s] * a[s][n-1]
        a[n-1][n-1] = a[n-1][n-1] - sum
        y[0] = b[0]
        for i in range(1,n):
            sum = 0.0
            for j in range(0,i):
                sum = sum + a[i][j] * y[j]
            y[i] = b[i] - sum
        x[n-1] = y[n-1] / a[n-1][n-1]
        for i in range(n-2,-1,-1):#求x[i]
            sum = 0.0
            for j in range(n-1,i,-1):
                sum = sum + a[i][j] * x[j]
            x[i] = ( y[i] - sum ) / a[i][i]
        for i in range(n):
            print("x" + str(i + 1) + " = ", x[i])
        print("x" " = ", x)
    
    if __name__ == '__main__':      #当模块被直接运行时,以下代码块将被运行,当模块是被导入时,代码块不被运行。
        a = np.array([[3.01, 6.03, 1.99], [1.27, 4.16, -1.23], [0.987, -4.81, 9.34]])
        b = np.array([1.0, 1.0, 1.0])
        m,n = a.shape
        x = np.zeros(n)
        B = np.zeros((n, n))
        for i in range(n):
            for j in range(n):
                B[i][j] = a[i][j]
       # CME(a,b,x)
        LU(a,b,x)
        #验证
        for i in range(0, n):
            temp = 0
            for j in range(0, n):
                temp = temp + B[i][j] * x[j]
            print("%f ", temp)
    
    
    展开全文
  • 高斯消元法与LU分解

    千次阅读 2019-04-22 18:27:24
    1.高斯消去法 有一般方程:Ax=b; 图1 ​​ 正常思路消元到最后,可以求出,回带之前的方程,可以解出. 图2 但是一旦数据规模大了,那么便很麻烦,很难求解。 ...

    1.高斯消去法

            有一般方程:Ax=b;

          

    图1
    ​​

    正常思路消元到最后,可以求出 x_{4} ,回带之前的方程,可以解出 x_{3},x_{2},x_{1}.

                  

    图2

                                                                              

    但是一旦数据规模大了,那么便很麻烦,很难求解。

    Gauss消元法的本质是将 矩阵A分解为 L* U .L是一下三角矩阵,U是一上三角矩阵.

    A为未知数的初始L^{-1}*A*X=L^{-1}*b矩阵,消元的过程可以用矩阵L来替代,L_{n} *L_{n-1}*L_{n-1}*......*L_{1}*A=A_{n}  

    消到最后有 A_{n}*X=b_{n}

    最后的A_{n}就是图2.

    令  L=(L_{n} *L_{n-1}*L_{n-1}*......*L_{1})^{-1}.  则 L^{-1}*A*X=L^{-1}*b

    由于最后得到的矩阵 A_{n} 是上三角矩阵,把  A_{n} 记为  U , 可得 U*X=L^{-1}*b, 所以我们只需要求解L矩阵,U矩阵就可以了。

     

    2.LU分解

          将系数矩阵A转变成等价两个矩阵L和U的乘积 ,其中L和U分别是单位下三角矩阵和上三角矩阵。当A的所有顺序主子式都不为0时,矩阵A可以分解为A=LU(所有顺序主子式不为0,矩阵不一定不可以进行LU分解)。其中L是下三角矩阵,U是上三角矩阵。

       Doolittle直接分解法

    例题:

    求得的L矩阵,U矩阵

    LU分解的代码如下:

    % 输入矩阵x=[A,b]; Ax=b;
    A=[1 2 3;2 2 8;-3 -10 -2];b=[0;-4;-11];
    %x=[A,b];
    n=length(A);%A的行数
    %LU分解
    U=zeros(n);L=eye(n,n);
    % 共n次循环,每一次循环先求出U的这一行的全部数据,求出L的这一列的全部数据
    for i=1:n
        %求出U的这一行的数据
    	row=i;
        for col=row:n
    	    sum=0;
            %U第row行第col列
    		for k=1:row-1
    		    sum=sum+L(row,k)*U(k,col);
    		end	
    		U(row,col)=A(row,col)-sum;
        end
    	%求出L的这一列的全部数据
    	col=i;
        %第col列第row行
        for row=col+1:n
            sum=0;
            for k=1:col-1
                sum=sum+L(row,k)*U(k,col);
            end
            L(row,col)=(A(row,col)-sum)/U(col,col);
        end
    end
    L
    U

     

    3.高斯消元法的回带

               L矩阵,U矩阵已经求解完毕,

    令  UX=y

          Ly=b,求解X。

     

    •        解 Ly=b.   对 i=1来说,y_{1}=b_{1} . 当i>=2,    y_{i}=b_{i}-\prod ^{j=i-1} _{j=1}L_{i,j}*y_{i-1}
    •        解 Ux=y   U的求解不再赘述

     

    解Ly=b,Ux=y的代码如下: 

    for t=2:n   %解Lx=b
        b(t)=b(t)-L(t,1:t-1)*b(1:t-1);
    end
    b(n)=b(n)/U(n,n);   %解Ux=b
    for t=1:n-1;
        k=n-t;b(k)=(b(k)-U(k,k+1:n)*b(k+1:n))/U(k,k);
    end
    x=b;  %方程Ax=b的解即为x

     

    总的代码如下:

    clc;clear all;
    % 输入矩阵x=[A,b]; Ax=b;
    A=[1 2 3;2 2 8;-3 -10 -2];b=[0;-4;-11];
    %x=[A,b];
    n=length(A);%A的行数
    %LU分解
    U=zeros(n);L=eye(n,n);
    % 共n次循环,每一次循环先求出U的这一行的全部数据,求出L的这一列的全部数据
    for i=1:n
        %求出U的这一行的数据
    	row=i;
        for col=row:n
    	    sum=0;
            %U第row行第col列
    		for k=1:row-1
    		    sum=sum+L(row,k)*U(k,col);
    		end	
    		U(row,col)=A(row,col)-sum;
        end
    	%求出L的这一列的全部数据
    	col=i;
        %第col列第row行
        for row=col+1:n
            sum=0;
            for k=1:col-1
                sum=sum+L(row,k)*U(k,col);
            end
            L(row,col)=(A(row,col)-sum)/U(col,col);
        end
    end
    L
    U
    %LU分解完成
    %求U*x=(L^-1)*b
    y=zeros(n,1);
    y(1)=b(1);
    for i=2:n   %解Lx=b
        y(i)=b(i)-L(i,1:i-1)*y(1:i-1);
    end
    
    x=zeros(n,1);
    x(n)=y(n)/U(n,n);
    %b(n)=b(n)/U(n,n);   %解Ux=b
    for i=1:n-1;
        k=n-i;
    	x(k)=(y(k)-U(k,k+1:n)*x(k+1:n))/U(k,k);
    end
    x  %方程Ax=b的解即为x

     

    展开全文
  • 高斯消去法SSE并行化实验

    千次阅读 2018-04-20 21:09:57
    高斯消去法原理和伪代码: 高斯消去法LU分解),是线性代数中的一个算法,可用来求解线性方程组,并可以求出矩阵的秩,以及求出可逆方阵的逆矩阵。高斯消元法的原理是:若用初等行变换将增广矩阵化为 ,则AX = B...

    高斯消去法原理和伪代码:

     

    高斯消去法(LU分解),是线性代数中的一个算法,可用来求解线性方程组,并可以求出矩阵的秩,以及求出可逆方阵的逆矩阵。高斯消元法的原理是:若用初等行变换将增广矩阵化为 ,则AX = B与CX = D是同解方程组。所以我们可以用初等行变换把增广矩阵转换为行阶梯阵,然后回代求出方程的解。

     

    总结一套流程就是:

    原线性方程组——> 高斯消元法——> 下三角或上三角形式的线性方程组——>前向替换算法求解(对于上三角形式,采用后向替换算法)

     

    所以高斯消去法(LU分解)串行算法如下面伪代码所示:

    for k := 1 to n do
    
        for j := k to ndo
    
            A[k, j] := A[k, j]/A[k, k];
    
        for i := k + 1to n do
    
            for j := k + 1 to n do
    
                A[i, j] := A[i, j] - A[i, k] × A[k, j ];
    
            A[i, k] := 0;

    这其中,内嵌的第一个for循环的作用是把第k行的所有元素除以第一个非零元素,目的是第一个非零元为1

    而第二个内嵌的for循环(当然其中还内嵌了一个小的for循环)作用是从k+1行开始减去第k行乘以这一行行的第一个非零元,使得k+1行的第k列为0

     

     

    SSE/AVX介绍:

     

    Intel ICC和开源的GCC编译器支持的SSE/AVX指令的C接口声明在xmmintrin.hpmmintrin.h头文件中。其数据类型命名主要有__m128/__m256、__m128d/__m256i,默认为单精度(d表示双精度,i表示整型)。其函数的命名可大致分为3个使用“_”隔开的部分,3个部分的含义如下。

    第一个部分为_mm或_mm256。_mm表示其为SSE指令,操作的向量长度为64位或128位。_mm256表示AVX指令,操作的向量长度为256位。本节只介绍128位的SSE指令和256位的AVX指令。

    第二个部分为操作函数名称,如_add、_load、mul等,一些函数操作会增加修饰符,如loadu表示不对齐到向量长度的存储器访问。

    第三个部分为操作的对象名及数据类型,_ps表示操作向量中所有的单精度数据;_pd表示操作向量中所有的双精度数据;_pixx表示操作向量中所有的xx位的有符号整型数据,向量寄存器长度为64位;_epixx表示操作向量中所有的xx位的有符号整型数据,向量寄存器长度为128位;_epuxx表示操作向量中所有的xx位的无符号整型数据,向量寄存器长度为128位;_ss表示只操作向量中第一个单精度数据;si128表示操作向量寄存器中的第一个128位有符号整型。

    3个部分组合起来,就形成了一条向量函数,如_mm256_add_ps表示使用256位向量寄存器执行单精度浮点加法运算。由于使用指令级数据并行,因此其粒度非常小,需要使用细粒度的并行算法设计。SSE/AVX指令集对分支的处理能力非常差,而从向量中抽取某些元素数据的代价又非常大,因此不适合含有复杂逻辑的运算。

     

    现在对于接下来代码中要用到的几个SSE指令加以介绍:

     

    _mm_loadu_ps用于packed的加载(下面的都是用于packed的),不要求地址是16字节对齐,对应指令为movups。

     

    _mm_sub_ps(__m128_A,__m128_B);返回一个__m128的寄存器,仅将寄存器_A和寄存器_B最低对应位置的32bit单精度浮点数相乘,其余位置取寄存器_A中的数据,例如_A=(_A0,_A1,_A2,_A3),_B=(_B0,_B1,_B2,_B3),则返回寄存器为r=(_A0*_B0,_A1,_A2,_A3)

     

    _mm_storeu_ps(float *_V, __m128 _A);返回一个__m128的寄存器,Sets the low word to the single-precision,floating-pointValue of b,The upper 3 single-precision,floating-pointvalues are passed throughfrom a,r0=_B0,r1=_A1,r2=_A2,r3=_A3

     

    SSE算法设计与实现:

           通过分析高斯的程序可以发现,高斯消去法有两部分可以实现并行,分别是第一部分的除法和第二部分的减法。即:

    1.第一个内嵌的for循环里的A[k, j]:= A[k, j]/A[k, k]; 我们可以做除法并行

    2.第二个双层for循环里的A[i, j] := A[i, j] - A[i, k] × A[k, j ];我们可以做减法并行

     

     

    我们来看核心代码

    1.      首先没加上并行化的高斯消去法:

    float** normal_gaosi(float **matrix) //没加SSE串行的高斯消去法
    
    {
    
        for (int k = 0; k < N; k++)
    
        {
    
            float tmp =matrix[k][k];
    
            for (int j = k; j < N; j++)
    
            {
    
                matrix[k][j] = matrix[k][j] / tmp;
    
            }
    
            for (int i = k + 1; i < N; i++)
    
            {
    
                float tmp2 = matrix[i][k];
    
                for (int j = k + 1; j < N; j++)
    
                {
    
                    matrix[i][j] = matrix[i][j] - tmp2 * matrix[k][j];
    
                }
    
                matrix[i][k] = 0;
    
            }
    
        }
    
        return matrix;
    
    }

    2.      再来看加上并行化的高斯消去法:

    void SSE_gaosi(float **matrix) //加了SSE并行的高斯消去法
    
    {
    
        __m128 t1, t2, t3, t4;
    
        for (int k = 0; k < N; k++)
    
        {
    
            float tmp[4] = { matrix[k][k], matrix[k][k], matrix[k][k], matrix[k][k] };
    
            t1 = _mm_loadu_ps(tmp);
    
            for (int j = N - 4; j >=k; j -= 4) //从后向前每次取四个
    
            {
    
                t2 = _mm_loadu_ps(matrix[k] + j);
    
                t3 = _mm_div_ps(t2, t1);//除法
    
                _mm_storeu_ps(matrix[k] + j, t3);
    
            }
    
            if (k % 4 != (N % 4)) //处理不能被4整除的元素
    
            {
    
                for (int j = k; j % 4 != ( N% 4); j++)
    
                {
    
                    matrix[k][j] = matrix[k][j] / tmp[0];
    
                }
    
            }
    
            for (int j = (N % 4) - 1; j>= 0; j--)
    
            {   
    
                matrix[k][j] = matrix[k][j] / tmp[0];
    
            }
    
            for (int i = k + 1; i < N; i++)
    
            {
    
                float tmp[4] = { matrix[i][k], matrix[i][k], matrix[i][k], matrix[i][k] };
    
                t1 = _mm_loadu_ps(tmp);
    
                for (int j = N - 4; j >k;j -= 4)
    
                {
    
                    t2 = _mm_loadu_ps(matrix[i] + j);
    
                    t3 = _mm_loadu_ps(matrix[k] + j);
    
                    t4 = _mm_sub_ps(t2,_mm_mul_ps(t1, t3)); //减法
    
                    _mm_storeu_ps(matrix[i] + j, t4);
    
                }
    
                for (int j = k + 1; j % 4 !=(N % 4); j++)
    
                {
    
                    matrix[i][j] = matrix[i][j] - matrix[i][k] * matrix[k][j];
    
                }
    
                matrix[i][k] = 0;
    
            }
    
        }
    
    }

     

    实验结果分析:

    为了测试其性能,我们把矩阵的大小从8,64,512,1024,2048,4096的矩阵进行高斯消去,并将串行所花费的时间与并行所花费的时间进行对比。因为后边矩阵太大我们仅看时间对比即可

     

    1.       N=8时,由于数据量较小,所花时间差距并不大。


     

    2.       N=64,由于数据量较小,所花时间差距并不大。

     

    3.          N=512,这里开始我们看到时间上的变化了。之后随着数据量逐渐增加,并行的优势逐渐体现出来。

     

     

     

    4.      N=1024

     

     

    5.       N=2048

     

    6.      N=4096

     

     

     

     

    总的来说,优势并没有特别大,究其原因,我觉得是因为在做最后并行的步骤之前有很多固定的步骤是需要一定时间的,比如对齐,导致SSE并行的方法需要花时间代价在这上面,没有想象中得那么快。

     

     

    附上整个代码:

    #include<pmmintrin.h>
    
    #include<time.h>
    
    #include<xmmintrin.h>
    
    #include<iostream>
    
    #defineN 4096
    
     
    
    usingnamespace std;
    
     
    
    float** normal_gaosi(float **matrix) //没加SSE串行的高斯消去法
    
    {
    
        for (int k = 0; k < N; k++)
    
        {
    
            float tmp =matrix[k][k];
    
            for (int j = k; j < N; j++)
    
            {
    
                matrix[k][j] = matrix[k][j] / tmp;
    
            }
    
            for (int i = k + 1; i < N; i++)
    
            {
    
                float tmp2 = matrix[i][k];
    
                for (int j = k + 1; j < N; j++)
    
                {
    
                    matrix[i][j] = matrix[i][j] - tmp2 * matrix[k][j];
    
                }
    
                matrix[i][k] = 0;
    
            }
    
        }
    
        returnmatrix;
    
    }
    
     
    
    void SSE_gaosi(float **matrix) //加了SSE并行的高斯消去法
    
    {
    
        __m128 t1, t2, t3, t4;
    
        for (int k = 0; k < N; k++)
    
        {
    
            float tmp[4] = { matrix[k][k], matrix[k][k], matrix[k][k], matrix[k][k] };
    
            t1 = _mm_loadu_ps(tmp);
    
            for (int j = N - 4; j >=k; j -= 4) //从后向前每次取四个
    
            {
    
                t2 = _mm_loadu_ps(matrix[k] + j);
    
                t3 = _mm_div_ps(t2, t1);//除法
    
                _mm_storeu_ps(matrix[k] + j, t3);
    
            }
    
            if (k % 4 != (N % 4)) //处理不能被4整除的元素
    
            {
    
                for (int j = k; j % 4 != ( N% 4); j++)
    
                {
    
                    matrix[k][j] = matrix[k][j] / tmp[0];
    
                }
    
            }
    
            for (int j = (N % 4) - 1; j>= 0; j--)
    
            {   
    
                matrix[k][j] = matrix[k][j] / tmp[0];
    
            }
    
            for (int i = k + 1; i < N; i++)
    
            {
    
                float tmp[4] = { matrix[i][k], matrix[i][k], matrix[i][k], matrix[i][k] };
    
                t1 = _mm_loadu_ps(tmp);
    
                for (int j = N - 4; j >k;j -= 4)
    
                {
    
                    t2 = _mm_loadu_ps(matrix[i] + j);
    
                    t3 = _mm_loadu_ps(matrix[k] + j);
    
                    t4 = _mm_sub_ps(t2,_mm_mul_ps(t1, t3)); //减法
    
                    _mm_storeu_ps(matrix[i] + j, t4);
    
                }
    
                for (int j = k + 1; j % 4 !=(N % 4); j++)
    
                {
    
                    matrix[i][j] = matrix[i][j] - matrix[i][k] * matrix[k][j];
    
                }
    
                matrix[i][k] = 0;
    
            }
    
        }
    
    }
    
     
    
    void print(float **matrix) //输出
    
    {
    
        for (int i = 0; i < N; i++)
    
        {
    
            for (int j = 0; j < N; j++)
    
            {
    
                cout << matrix[i][j]<<" ";
    
            }
    
            cout << endl;
    
        }
    
    }
    
     
    
    int main()
    
    {
    
        srand((unsigned)time(NULL));
    
        float **matrix = newfloat*[N];
    
        float **matrix2 = newfloat*[N];
    
        for (int i = 0; i<N; i++)
    
        {
    
            matrix[i] = newfloat[N];
    
            matrix2[i] = matrix[i];
    
        }
    
     
    
        //cout << "我们生成了初始随机矩阵" << endl;
    
        for (int i = 0; i < N; i++)
    
        {
    
            for (int j = 0; j < N; j++)
    
            {
    
                matrix[i][j] = rand() % 100;
    
            }
    
        }
    
        //print(matrix);
    
     
    
        cout <<endl<<endl<<endl<<"不使用SSE串行的高斯消去法" << endl;
    
        clock_t  clockBegin,clockEnd;
    
        clockBegin = clock(); //开始计时
    
        float **B = normal_gaosi(matrix);
    
        clockEnd = clock();
    
        //print(matrix);
    
        cout << "总共耗时: " << clockEnd - clockBegin << "ms" << endl;
    
     
    
        cout <<endl<<endl<<endl<< "使用SSE并行的高斯消去法" << endl;  
    
        clockBegin = clock(); //开始计时
    
        SSE_gaosi(matrix2);
    
        clockEnd = clock();
    
        //print(matrix2);
    
        cout << "总共耗时: " << clockEnd - clockBegin << "ms" << endl;
    
        system("pause");
    
        return 0;
    
    }
    
     

     


    展开全文
  • 提出了基于righ-looking LU分解法的并行高斯消去算法,利用GPU(图形处理器)加速求解复系数回路阻抗方程组;采用GIS(地理信息系统)和虚拟现实技术,对潮流计算结果进行仿真可视化.仿真算例表明,该方法对节点编号无...
  • 一般方阵- 高斯消去法与矩阵的LU分解2.可逆矩阵- Doolittle/Crout分解3.分块方阵- 拟LU分解与拟LDU分解 文章目录《矩阵论》学习笔记(四)-1:4.1 矩阵的三角分解一、Gauss消去法的矩阵形式1.1. 引入1.2.矩阵论的...

    《矩阵论》学习笔记(四):4.1 矩阵的三角分解

    矩阵的三角分解
    1.一般方阵- 矩阵的LU/LDU分解
    2.可逆方阵- Doolittle/Crout/Gholesky分解
    3.分块方阵- 拟LU分解与拟LDU分解
    • 提出的目的:
      由于三角矩阵的计算,像行列式、逆矩阵、求解线性方程组等都是很方便的,所以考虑将普通矩阵分解成一些三角矩阵的乘积,从而简化运算。

    一、Gauss消去法的矩阵形式

    1.1. 引入

    对n元线性方程组Ax=wA\vec x=\vec w,求解通常采用高斯主元素消去法— 将增广矩阵进行行阶梯化简成三角矩阵。

    线性方程组Ax=wA\vec x=\vec w 增广矩阵 高斯主元素消去法
    {a11x1+a12x2+a13x3=w1a21x1+a22x2+a23x3=w2a31x1+a32x2+a33x3=w3\begin{cases}a_{11}x1+a_{12}x2+a_{13}x3=w1\\a_{21}x1+a_{22}x2+a_{23}x3=w2\\a_{31}x1+a_{32}x2+a_{33}x3=w3\end{cases} [a11a12a13w1a21a22a23w2a31a32a33w3]\left[ \begin{array}{cccc} a_{11} & a_{12} & a_{13}& w_1 \\ a_{21} & a_{22}& a_{23} & w_2 \\ a_{31}& a_{32}& a_{33} & w_3 \end{array} \right] [a11a12a13z10a22a23z200a33z3]\left[ \begin{array}{cccc} a_{11}'&a_{12}'&a_{13}' &z_1' \\0&a_{22}' &a_{23}'& z_2' \\ 0&0 &a_{33}' & z_3'\end{array} \right]

    1.2. 矩阵论中的三角分解理论

    为建立矩阵论的三角分解理论,使用矩阵理论描述高斯主元素消去法的消元过程。

    - 高斯消元过程:

    对矩阵A采用按自然顺序选主元素法进行消元。

    A(0)=AA^{(0)}=A,由于倍加初等变换不改变矩阵行列式的值,不断使用A的k阶顺序主子式k△k,构造Frobenius矩阵LkL_k,最终得到矩阵A(n1)A^{(n-1)}.
    A(n1)=[a11(0)a12(0)...a1n(0)a22(1)...a2n(1)...ann(n1)]A^{(n-1)}=\left[ \begin{array}{cccc} a_{11}^{(0)} & a_{12}^{(0)} &...& a_{1n}^{(0)} \\ & a_{22}^{(1)}& ... & a_{2n}^{(1)} \\ & &...\\ & & &a_{nn}^{(n-1)} \end{array} \right]

    -高斯消元过程 几个要点
    1- 过程中的特点 不使用行/列变换.
    2- 能进行到底的条件 前n-1个顺序主子式均不为0.
    3- 结果的特点 矩阵A高斯消元过程得到的LU分解是存在且唯一。
    其中,L是单位下三角阵,U是上三角阵。

    二、 矩阵的LU/LDU分解

    A=A(0)=L1A(1)=...=L1L2...Ln1A(n1)A=A^{(0)}=L_1A^{(1)}=...=L_1L_2...L_{n-1}A^{(n-1)}得到以下:

    2.1. LU/LDU分解公式

    – - 矩阵的LU分解公式:
    > A=LUA=LU.
    其中,L为下三角矩阵,U为下三角矩阵。L=L1L2...Ln1L=L_1L_2...L_{n-1}U=A(n1)=[a11(0)a12(0)...a1n(0)a22(1)...a2n(1)...ann(n1)]U=A^{(n-1)}=\left[\begin{array}{cccc} a_{11}^{(0)} & a_{12}^{(0)} &...& a_{1n}^{(0)} \\& a_{22}^{(1)}& ... & a_{2n}^{(1)} \\ & &...\\ & &&a_{nn}^{(n-1)} \end{array} \right](能够将一个矩阵表示成一个上三角矩阵×下三角矩阵的形式,简化计算)
    – - 矩阵的LDU分解公式:
    > A=LDUA=LDU.
    其中,L为单位下三角矩阵,D为对角矩阵,U为单位下三角矩阵。

    2.2. 三角分解的存在性与唯一性

    先要满足‘存在’,讨论 ‘是否是唯一’ 才有意义。

    - 存在性 唯一性
    奇异矩阵 1- 若满足前n-1个顺序主子式≠0 \Rightarrow 三角分解存在 [必要条件]
    2- 若三角分解存在,不一定满足前n-1个顺序主子式≠0。
    3-三角分解不存在
    - ( 奇异 det(A)=0\Rightarrow det(A)=0ann=0a_{nn}=0 )
    1- 若三角分解存在且满足前n-1个顺序主子式≠0
    \Rightarrow A=LU唯一;A=LDU唯一
    2- 若三角分解存在但不满足前n-1个顺序主子式≠0
    \Rightarrow A=LU/LDU不一定唯一
    3- 无
    非奇异矩阵 1- 若满足n-1个顺序主子式≠0 \Leftrightarrow 三角分解存在 [充要条件]。
    [实际上,非奇异矩阵的anna_{nn}≠0,n个主子式都≠0.]
    2- 若三角分解不存在,即不满足n个顺序主子式≠0
    \Leftrightarrow 存在P使PA的所有顺序主子式≠0。即A虽不存在,PA存在LU分解。 [条件更强(只满足前n-1个就行),但这是非奇异本身的性质]
    1- 三角分解存在(A的前(n-1)阶顺序主子式≠0)
    \Leftrightarrow A=LU唯一;A=LDU唯一
    2- 三角分解不存在
    \Leftrightarrow PA=LU唯一; PA=LDU唯一
    • 前(n-1)阶顺序主子式与唯一性
    前(n-1)阶顺序主子式 存在性与唯一性
    1- 均≠0 \Leftrightarrow A=LU/LDU存在且唯一,Doolittle/Crout/Gholesky分解唯一 [充要条件]
    2- 存在0 \Rightarrow A=LU/LUD不一定存在;
    若存在,也不一定唯一,可能有多种分解方式。在这里插入图片描述
    - 对以上的理解:
    这里所研究的方阵A的三角分解是从Gauss消元法得出的。根据该过程:
    1. 只要求前n-1个顺序主子式不为零,总可以得到LU分解(这里L是单位下三角阵,U是上三角阵),且这个分解唯一。
    而没要求最后第n阶顺序主子式为0(即det(A)可能为0,对应anna_{nn}=0),所以A可能是奇异的。若奇异,且LU存在,U的最右下角应为0。这个时候Gauss消元三角分解已经结束,即只进行到A(n-1)。
    2. 这里引出了三角分解的概念,同时引入一个问题— A=LU是否唯一。
    显然因为A=LDD1UA=LDD^{-1}U, LU不唯一。
    但自然顺序选主元下的Gauss消元法,执行过程要求前n-1个顺序主子式为零,即可导出唯一的LU分解(L是单位下三角阵,U是上三角阵),且这个分解可以写成LDU的形式(即把原U的对角元素取出作为对角阵D即可,相应的U变成单位上三角阵)。
    3. 从而引出定理-
    前n-1个顺序主子式不为0 当且仅当(“等价于”) A可唯一的分解为LDU】。
    但如果前n-1个顺序主子式中出现0,这样三角分解不存在或有多个。
    4. 推论:对于非奇异的情形,存在一般的LU分解即等价于所有顺序主子式为0。
    但要注意:要求顺序主子式为0的条件比较强。—> 对于非奇异的矩阵A,是存在P使PA的所有顺序主子式≠0。即A虽然没有,PA是有LU分解的。
    5. 对于奇异矩阵A也可以类似考虑置换矩阵来考察它是否满足前n-1个顺序主子式是否都为0。这其实对应着按列选主元(只有行置换),还是总体选主元(行列都可置换)。这样原矩阵A不能三角分解,但其置换后的矩阵可能可以进行LU分解。
    6. 一般的三角分解可能是不存在的,也可能是不唯一的。
    但按自然Gauss消元法,自然限定一边是单位三角矩阵,且因满足定理4.1,所以可以唯一确定。
    从这个意义上,上面探讨的三角分解条件可以推出一般的三角分解,反之则不然。
    但是结合初等变换与三角分解,仍可对一个秩为r的矩阵,进行类似计算。

    2.3. 三角分解的求解方法

    参见:[《线性代数及其应用》2.5]

    对可逆方阵A,若A可以化成阶梯型矩阵U,且简化过程仅使用行倍加变换:
    1. 存在单位下三角初等矩阵E1,...,EpE_1,...,E_p使得:Ep...E1A=UA=(Ep...E1)1UE_p...E_1*A=U \to A=(E_p...E_1)^{-1}U
    2. L=(Ep...E1)1L=(E_p...E_1)^{-1},满足(Ep...E1)L=I(E_p...E_1)*L=I
    - 即:AA左乘行变换U\to ULL左乘行变换I\to I
    - 三角分解的求解步骤 :
    1. AUAmnA \to U:A_{m*n} 经过初等行变换(无行交换)\toUmnU_{m*n},其中U是行阶梯矩阵.
    2. ULUmnU \to L:对U_{m*n} 确定主元列,每列元素除以列主元Lmm\to L_{m*m},其中L是单位下三角方阵.
    若主元列r<m,则用ImI_m补齐L的后(m-r)列.

    三、可逆矩阵的三角分解

    3.1. Doolittle/Crout分解的前提条件

    A是非奇异/可逆矩阵,且假定A的LDU分解存在。

    (这样A=LDU,D对角元素一定是非零的。如果按自然Gauss消元法,D最后一个元素可为0.)

    可逆矩阵A有特殊LU分解方法 分解公式
    1- Doolittle分解算法 A=L(DU)=LUA=L(DU)=LU'
    2- Crout分解算法 A=(LD)U=LUA=(LD)U=L'U
    3- Gholesky分解算法 A=GGTA=GG^T

    3.2. Doolittle分解

    A=L(DU)=LUA=L(DU)=LU',其中U=DUU'=DU
    L为对角元素为1的下三角矩阵,U’为上三角矩阵。
    L=[1l211...ln1ln2...1]U=[u11u12...l1nu22...u2n...unn]L'=\left[\begin{array}{cccc} 1& & & \\ l_{21}&1& & \\ &...\\ l_{n1} &l_{n2}&...&1\end{array} \right],U=\left[\begin{array}{cccc} u_{11} &u_{12} & ...&l_{1n} \\ & u_{22} &... &u_{2n} \\ &&...\\ & & &u_{nn} \end{array} \right],
    为了节省存储空间,将L和U’中元素写入矩阵A中相应位置上,得到如下矩阵:[l11u12...l1nl21l22...u2n...ln1ln2...lnn] \to \left[\begin{array}{cccc} l_{11}&u_{12} & ...&l_{1n} \\l_{21}& l_{22}&... &u_{2n} \\ &&...\\ l_{n1} &l_{n2}&...&l_{nn}\end{array} \right]

    3.3. Crout分解

    A=(LD)U=LUA=(LD)U=L'U,其中L=LDL'=LD
    L’为下三角矩阵,U为对角元素为1的上三角矩阵。
    L=[l11l21l22...ln1ln2...lnn]U=[1u12...l1n1...u2n...1]L'=\left[\begin{array}{cccc} l_{11}& & & \\ l_{21}& l_{22}& & \\ &...\\ l_{n1} &l_{n2}&...&l_{nn}\end{array} \right],U=\left[\begin{array}{cccc} 1&u_{12} & ...&l_{1n} \\ & 1&... &u_{2n} \\ &&...\\ & & &1\end{array} \right],
    为节省存储空间,将L’和U中元素写入矩阵A中相应位置上,得到如下矩阵:[l11u12...l1nl21l22...u2n...ln1ln2...lnn] \to \left[\begin{array}{cccc} l_{11}&u_{12} & ...&l_{1n} \\l_{21}& l_{22}&... &u_{2n} \\ &&...\\ l_{n1} &l_{n2}&...&l_{nn}\end{array} \right]

    • Crout分解算法步骤:
      L’的列 与 U的行 交替计算 (见课本p135)。

    3.4. Gholesky分解

    如果A不但是非奇异/可逆矩阵,A的LDU分解存在,而且是实对称正定矩阵,则有:

    A=ATLDU=UTDLTU=LTL=UTA=A^T、LDU=U^TDL^T \to U=L^T、L=U^T,从而:
    A=LDU=L(D)2U=(LD)(LD)T=GGTA=LDU=L(D')^2U=(LD')(LD')^T=GG^T,其中,G=LDG=LD'是下三角矩阵,D=(D)D'=\sqrt{(D)}.

    四、分块矩阵的拟LU分解与拟LDU分解

    • 提出的目的:
      对于高阶方阵,如果能将其分解成拟三角矩阵和拟对角矩阵,将会极大的减少计算量。
      对矩阵ARnnA∈R^{n*n},将A裂分成:
      A=[A11A12A21A22] A= \left[\begin{array}{cccc}A_{11}&A_{12} \\A_{21}& A_{22}\end{array} \right]
    1. 如果A11A_{11}是可逆的,可构造下三角阵,使得A能拟LU分解与拟LDU分解。
      det(A)=det(A11)det(A22A21A111A12)0\to det(A)=det(A_{11})det(A_{22}-A_{21}A_{11}^{-1}A_{12}) ≠0
    2. 如果A22A_{22}是可逆的,可构造下三角阵,使得A能拟LU分解与拟LDU分解。
      det(A)=det(A22)det(A11A12A221A21)0\to det(A)=det(A_{22})det(A_{11}-A_{12}A_{22}^{-1}A_{21}) ≠0
    • 推出结论:

    如果矩阵 ARmnBRnmA∈R^{m*n},B∈R^{n*m},则有:
    det(Im+ABT)=det(In+BTA)det(I_m+AB^T)=det(I_n+B^TA)

    五、三角分解的应用

    展开全文
  • 第6章 变治法预排序:用于下列场景,检验数组中元素的唯一性,...LU分解高斯消去法得到的上三角阵为U,在消去过程中行的乘数构成了下三角阵L(假设不存在行交换的情况下)。比如,在消去第i列的时候,把第i行乘以...
  • 高斯消去法 rref():将一个线性方程式的增广矩阵转化为除了对角线元素为1外,其余元素均为0 >> A=[1 2 1; 2 6 1; 1 1 4]; >> b=[2;7;3]; >> R=rref([A b]) R = 1 0 0 -3 0 1 0 2 0 0 1 1 *...
  • MATLAB LU函数

    2013-11-20 20:18:00
     高质量学习资源免费获取,专注但不限于【Linux】【C/C++/Qt... 高斯消元求解线性方程,包括把增广矩阵转换为三角矩阵形式的过程,消去阶段工作 步骤是把矩阵A分解成为下三角L和上三角U的乘积。这种...
  • 该教程包括LU分解和Cholesky分解的定义,Cholesky分解的条件,使用Numpy特征值函数测试正定性,Cholesky算法的派生和Python中的编码。 高斯-乔丹方法教程-逐步理论编码 在本教程中,使用符号和数字示例逐步说明了...
  • GaussElimination.m实现了高斯消去法求解方程组 L_GaussElimination.m实现了带列主元的高斯消去法求解方程组 LU.m实现了LU分解 LUP.m实现了带列主元的LU分解 reverse_and_det.m实现了利用LU分解求矩阵的逆和行列式 ...
  • 高斯消去法改写为紧凑形式,可以直接从矩阵A的元素计算出L,U元素的递推公式,因此可以使A分解为L和U。利用LU分解法,求解Ax=b等价于求解Ly=b和Ux=y。 我定义了一个创建希尔伯特矩阵的函数,通过输入行数(列数)来...
  • 高斯消去法rref() R = rref(A) 使用 Gauss-Jordan 消元法和部分主元消元法返回A的简化行阶梯形。 对增广矩阵 [A b] 使用rref()则可以求解 Ax=b 对应的线性方程组 LU因子化 [L,U,P] = lu(A) 将满矩阵或稀疏矩阵 A ...
  • 通常的逆矩阵可以用高斯消去法计算。十分有效。还可以使用LU分解,QR分解等。 二扩域中的逆矩阵则不同。看似简单,其实有别:它的所有元素定义在GF(2^m)中。从理论来看,似乎也可以用高斯消去法,只是计算规则 ...
  • 4.2 矩阵的LU分解法 4.3 特殊线性方程组的三角分解法 4.4 向量范数及矩阵范数 4.5 方程组的性态误差分析 习题四 第5章 解线性方程组的迭代法 5.1 常用迭代法 5.2 迭代法的收敛性 习题五 第6章 非线性方程求根 6.1 ...
  • 6.3 LU分解法 6.4 两类特殊矩阵方程 本章总结 习题6 计算实习6 第7章 线性方程组的迭代解法 7.1 迭代法的原理 7.2 古典迭代法及其收敛性 7.3 共轭梯度法 本章总结 习题7 计算实习7 第8章 矩阵特征值问题...
  • 算法设计分析基础

    2018-02-12 19:26:04
    6.2高斯消去法 6.2.1LU分解 6.2.2计算矩阵的逆 6.2.3计算矩阵的行列式 习题6.2 6.3平衡查找树 6.3.1AVL树 6.3.22—3树 习题6.3 6.4堆和堆排序 6.4.1堆的概念 6.4.2堆排序 习题6.4 6.5霍纳法则和二进制幂 6.5.1霍纳法...
  • 线性代数 | 复习笔记

    2020-07-26 18:36:49
    Schwarz 不等式三角不等式2 矩阵线性方程组对矩阵向量乘积的理解对线性方程组的理解可逆矩阵线性方程组的行图和列图3 高斯消元矩阵的初等行变换增广矩阵消去矩阵置换阵4 矩阵的运算矩阵乘法的性质分块矩阵矩阵...
  • 该部分就是针对线性方程组求解而设计的,内容包括:线性方程组的直接解法:Gauss消去法、Gauss列主元消去法、Gauss全主元消去法、列主元消去法应用『列主元求逆矩阵、列主元求行列式、矩阵的三角分解』、LU分解法、...
  • 数学算法原书光盘

    2012-05-11 17:38:03
    2.LU分解法 3.追赶法 4.五对角线性方程组解法 5.线性方程组解的迭代改善 6.范德蒙方程组解法 7.托伯利兹方程组解法 8.奇异值分解 9.线性方程组的共轭梯度法 10.对称方程组的乔列斯基分解法 11.矩阵的QR分解 12.松弛...
  • 2.LU分解法 3.追赶法 4.五对角线性方程组解法 5.线性方程组解的迭代改善 6.范德蒙方程组解法 7.托伯利兹方程组解法 8.奇异值分解 9.线性方程组的共轭梯度法 10.对称方程组的乔列斯基分解法 11.矩阵的QR分解 12.松弛...
  • 算法经典VC++数值分析

    热门讨论 2010-10-21 21:03:49
    1.2LU分解法 1.3追赶法 1.4五对角线性方程组解法 1.5线性方程组解的迭代改善 1.6范德蒙(Vandermonde)方程组解法 1.7托伯利兹(Toeplitz)方程组解法 1.8奇异值分解 1.9线性方程组的共轭梯度法 1.1对称方程组...
  • delphi算法集源码

    2012-07-24 09:46:04
    2.LU分解法 3.追赶法 4.五对角线性方程组解法 5.线性方程组解的迭代改善 6.范德蒙方程组解法 7.托伯利兹方程组解法 8.奇异值分解 9.线性方程组的共轭梯度法 10.对称方程组的乔列斯基分解法 11.矩阵的QR分解...
  • Visual C++ 常用数值算法集

    热门讨论 2012-03-19 11:57:59
    1.2LU分解法 1.3追赶法 1.4五对角线性方程组解法 1.5线性方程组解的迭代改善 1.6范德蒙(Vandermonde)方程组解法 1.7托伯利兹(Toeplitz)方程组解法 1.8奇异值分解 1.9线性方程组的共轭梯度法 1.1对称方程组...
  • 1.全主元高斯约当消去法2.LU分解法3.追赶法4.五对角线性方程组解法5.线性方程组解的迭代改善6.范德蒙方程组解法7.托伯利兹方程组解法8.奇异值分解9.线性方程组的共轭梯度法10.对称方程组的乔列斯基分解法11.矩阵的QR...
  • 2.LU分解法 3.追赶法 4.五对角线性方程组解法 5.线性方程组解的迭代改善 6.范德蒙方程组解法 7.托伯利兹方程组解法 8.奇异值分解 9.线性方程组的共轭梯度法 10.对称方程组的乔列斯基分解法 11.矩阵的QR分解...
  • 该部分就是针对线性方程组求解而设计的,内容包括:线性方程组的直接解法:Gauss消去法、Gauss列主元消去法、Gauss全主元消去法、列主元消去法应用『列主元求逆矩阵、列主元求行列式、矩阵的三角分解』、LU分解法、...

空空如也

空空如也

1 2
收藏数 33
精华内容 13
关键字:

lu分解法与高斯消去法