用fortran语言进行高斯消去法的openmp的并行编写,但是结果却不正确,求解

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了用fortran语言进行高斯消去法的openmp的并行编写,但是结果却不正确,求解相关的知识,希望对你有一定的参考价值。

本人正在自学openmp的并行编程,现在编写高斯消去法求解系数的增广矩阵的求解程序,并行的程序和串行程序求解的结果并不一样,求指点问题出在哪里,谢谢!!!
下面程序中A矩阵式m行n列的矩阵:
SUBROUTINE GSS(A,MM,NN,M,N)
USE OMP_LIB

IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(MM,NN)

INTEGER I,J,K,NTH,L,NT,T,ID
REAL*8 BUFFER

NT=0
T=0
K=0

PRINT *,'Matrix A is already'

C 开始并行计算
!$OMP PARALLEL SHARED(A,NTH) PRIVATE(ID,I,J,L,BUFFER)
ID=OMP_GET_THREAD_NUM()
!$OMP MASTER
NTH=OMP_GET_NUM_THREADS()
PRINT *,'Compute the coefficient matrix with',NTH,'threads'
!$OMP END MASTER
PRINT *,'Thread',ID,'is already'
!$OMP END PARALLEL

c 开始进行计算,首先消去下三角

do k=1,m-1

!$OMP PARALLEL SHARED(A,NTH,K,M,N) PRIVATE(I,J,L,BUFFER)
!$OMP SINGLE
L=K
DO I=K+1,M
IF(ABS(A(I,K)).GT. ABS(BUFFER))THEN
BUFFER=A(I,K)
L=I
END IF
END DO
IF(L .NE. K)THEN
DO J=K,N
BUFFER=A(K,J)
A(K,J)=A(L,J)
A(L,J)=BUFFER
END DO
END IF
IF(ABS(A(K,K)) .LT. 1.0E-20 )THEN
STOP
PRINT *,'除数为0'
END IF
!$OMP END SINGLE
!$OMP DO SCHEDULE(DYNAMIC,1)
do i=k+1,m
buffer=A(i,k)/A(k,k)
do j=k,n
A(i,j)=A(i,j)-buffer*A(k,j)
end do
end do
!$OMP END DO
!$OMP END PARALLEL
end do

print *,'**************消去下三角后********************'

c 接下来消去上三角

do k=m,2,-1
!$OMP PARALLEL SHARED(A,NTH,K,M,N) PRIVATE(I,J,L,BUFFER)
!$OMP DO SCHEDULE(DYNAMIC,1)
do i=k-1,1,-1
buffer=A(i,k)/A(k,k)
do j=i+1,n
A(i,j)=A(i,j)-buffer*A(k,j)
end do
end do
!$OMP END DO
!$OMP END PARALLEL
end do

print *,'**************消去上三角后********************'

c 系数矩阵运算完成,求解矩阵x
!$OMP PARALLEL SHARED(A,NTH,K) PRIVATE(I,J,L,BUFFER)
!$OMP DO SCHEDULE(DYNAMIC,1)
do i=1,m
a(i,n)=a(i,n)/a(i,i)
end do
!$OMP END DO
!$OMP END PARALLEL
c 运算结束,关闭并行区域,处理结果
c 初步检查运算结果
print *,'**************************************************'
print *,'Done'

end

参考技术A 并行程序段中,要求前后左右计算没有关联。
比如你要对屏幕上的任何一个点与字符A做异或运算,这些点的运算之间是没有关联的,这时你可以安排进行并行计算。追问

是没有关联的 我只是在矩阵消去某一行下面的行的时候进行了并行

高斯消去法与矩阵三角分解法(LU分解)


引言

  1. 研究数值解法的必要性

  2. 线性代数方程组的常用解法
  • 直接法
    只包含有限次四则运算。在计算过程中不发生舍入误差的假定下,计算结果就是原方程组的精确解。(适用于小规模的n阶稠密线性方程组)
  • 迭代法
    从解的某个近似值出发,通过构造一个无穷序列去逼近精确解的方法(一般有限步内得不到精确解)。实现这一极限过程每一步的结果是把前一步所得的结果施行相同的演算步骤得到的。(适用于大规模的n阶稀疏线性方程组)
    Remark:由于运算过程中舍入误差的存在,实际上
    直接方法得到的解也是方程组的近似解。

一、预备知识

1.1 矩阵基本运算和和行列式

1.2 矩阵的特征值和谱半径


设A是n阶方阵,如果数λ和n维非零列向量x使关系式Ax=λx成立,那么这样的数λ称为矩阵A特征值,非零向量x称为A的对应于特征值λ的特征向量。

  • 特征值性质

    t r ( A ) tr(A) tr(A)为A 的

  • 注意

  • 谱半径

    习题

1.3 特殊矩阵

  • 对称正定矩阵

二、Gauss消去法

  • 方程组的矩阵表示

    高斯消元的基本思想:将矩阵A的下三角部分消为零,即化为上三角形

2.1 不消主元的Gauss消去法

首先将方程组Ax=b化为上三角方程组,此过程称为消去过程,再求解上三角方程组,此过程称为回代过程.

例.


总结:Gauss消去法解Ax=b分两步:消元回代.

  • 不消主元的Gauss消去法
    将方程组Ax=b的系数矩阵与右端项合并为:

    第一步:

    第二步:

    第k步:

    第n-1步:


    回代:
  • 定理

2.2 高斯主元素消去法

在计算过程中的舍入误差增大能得到控制,该方法就是稳定的。
小主元是不稳定的根源,这就需要采用“选主元素”技术,即选取绝对值最大的元素作为主元

1. Gauss列主元消去法

a k k ( k ) a_kk^(k) akk(k) 称为第k步的主元

  • Gauss列主元消去法:
    第一步:在第一列中选取最大值的元素,并将该行 a i , 1 a_i,1 ai,1与第一行 a 1 , 1 a_1,1 a1,1进行交换。(若该列中第一行 a 1 , 1 a_1,1 a1,1最大,无需交换)


    第k步:设第 步将增广矩阵(A,b)已化为


    结论:当 k k k 1 1 1变化到 n − 1 n-1 n1时, ( A , b ) (A,b) (A,b)成为上三角形式,再用回代公式求 x x x.
    例. 用Gauss列主元解方程组

2. Gauss全主元消去法

概念不好理解,直接上例题+解析。
例.

通常用列主元消去法。


三、矩阵的三角分解法

3.1 用直接三角分解法求解方程组

  1. 高斯消元与LU分解的等价性





  • L与U ⭐⭐
  1. 矩阵三角分解的基本定理

3.2 Doolittle分解法(杜利特尔)



题目

  • 例题

  1. 谱半径,就是特征值绝对值(复数取模)中的最大值,先求特征值。再取模,分别得到√5,√5,因此谱半径是√5。

  • 例题


以上是关于用fortran语言进行高斯消去法的openmp的并行编写,但是结果却不正确,求解的主要内容,如果未能解决你的问题,请参考以下文章

高斯消去法与矩阵三角分解法(LU分解)

《数值分析》-- 高斯消去法与矩阵三角分解法(LU分解)

13.高斯消去法——三角矩阵

OpenMP 的全局变量

intel fortran如何实现单机多核并行运算

如何判断矩阵是不是能够进行LU分解