高斯消元法Fortran90程序

(三) 数组赋值 a(1:9,1:9) a(2:4,5:8)=a(3:5,1:4) 有效 a(1:4,1:3)=a(1:3,1:4) 无效 因形状不同 a(k+1:n,k)=0 给多个元素赋标量，有效 (四) where语句和where结构 where语句可用于为数组赋值，如a和b是两个形状相同的数组，将b中与 a的负元素对应位置上的元素置0，可写成 where(a<0) b=0 where结构 where(c/=0) a=b/c elsewhere a=0 c=1 end where
笫十五章 FORTRAN 90
FORTRAN 90 与FORTRAN 77 不同之处
一、说明语句
（一）固有数据类型 整型、实型、复型，其后都可以接一个下划线_和另一个 整型常量或有名整型常量用来指明种别参数。如 42_2 42_long 0.3_double (3.14_double,_7) ascii_’john Q.public’ 种别参数(kind parameter)：当它是一个整数1,2,3…时， 可能表示单精度双精度和四精度，在不同的系统中可能是 4,8,16。 （二）变量 integer (kind=2):: x,y character (len=20,kind=kanji) ::name kinji为一整型常量 real ::a=1.2,b=2.3, 初始化变量的值在执行过程中可变
六、数组
(六) 列表检索 1. 对无序列表的顺序检索 subroutine search_1(lost_card,card_number,found) integer,dimension(:),intent(in) :: lost_card integer,intent(in) ::card_number logcal,intent(out) ::found integer ::I found=.false. do I=1,size(lost_card) if(card_number==lost_card(I)) then found=.trun. exit end if end do end

程序设计举例
算法

算法是计算机软件中的一个基本概念，它是对解决实际问题的方法和步
骤的描述。设计算法是程序设计的核心，也是编写程序的基础。

正确的算法有以下几个显著特点：
(1)有穷性(flniteness)。 (2)确定性(certainty)。 (3)有效性(effcctivencss)。 (4)有零个或多个输入(input)。 (5)有一个或多个输出(output)。
简单的FORTRAN90程序分析

（6）空格并不是随处都可以使用的，一个像关键字、变量和常量名以 及操作符等的字符，其内部是不能使用空格的，空格会使字符失去其原 有的含义。 （7）FORTRAN90的注释语句是以感叹号为标志的，一行中感叹号后的 所有字符都被编译器作为注释语句而忽略。注释语句可以单独占一行， 也可在程序的其他语句后面出现。在FORTRAN90中，空行被作为注释 语句。 （8）如果遇到一条语句的长度超过了FORTRAN90所允许的行最大长度， 需要写成几行，FORTRAN90提供了一个续行符 ( & ) ，通过在语句末尾 添加续行符，编译器就会把下一行作为续行来处理。如果是把一个如关 键字这样的字符分成两行，那么需要在下行语句的开头再加一个续行符。
CHARACTER ( LEN = 10 ) :: STRING 子字符串是字符串中的若干个连续字符的集合。子字符
CHARACTER (LEN = 5 ) :: SUBSTRING 串的表示方法如下： CHARACTER (LEN = 1 ) :: CHAR STRING=’JANEJORDE’ 字符串变量名（[FIRST-POSITION] : [END-POSITION]） ! SUBSTRING的值为“JANE” SUBSTRING=STRING ( :5 ) 缺省的开始位置是１，缺省的结束位置是字符串的长度。 SUBSTRING=STRING ( 5:6 ) ! SUBSTRING的值为“J” 如果开始位置大于结束位置则子字符串为空，即它的长 SUBSTRING=STRING ( 3:7 ) ! SUBSTRING的值为“NEJO” N=7 度为０。 CHAR=’ABCDEFGHIJK’ ( N:N ) ! CHAR的值为“G”

字。以关键词end type结束。
派生类型的声明和使用
n 声明： n 赋值： n 引用结构成员：
例一：派生数据类型
模块(Module)
n 是一种不能直接执行的程序单元，但可以包含
数据声明和过程
n
使不同的过程实现过程和数据的共享 n 可以实现对自定义数据类型的封装
n 相对于program,function,subroutine,是一种特
例三：数据类型的范围
定义符号常量
n 符号常量在程序中不可以改变其值 n 定义符号常量可以增强程序的可读性和可维护
性 n 定义方法：加上parameter属性
real (8) ,parameter :: pi=3.1415926535_8
例四：定义符号常量
三、控制结构
n if语句:if then elseif else endif
n
在Compaq Visual Fortran中调用 CXML
n 第一种方法：在工程中加入cxml.lib
n 第二种方法：使用include
PROGRAM example INCLUDE 'CXML_INCLUDE.F90’
在Compaq Visual Fortran中调用 Intel MKL
n 添加库函数路径 n 在工程中加入mkl_s.lib
总结：Fortran90的新特性
n n n n n n n n
高级的表达能力：数组运算 (例如, X(1:N)=R(1:N)*COS(A(1:N))) 动态存储分配 (ALLOCATE, DEALLOCATE, ...) 导出类型(derived types)与操作符重载(operator overloading) 更好的变量声明方式(parameterized types),提高移植性 模块 更多的现代控制结构 (SELECT CASE，END SELECT,DO, END DO,EXIT...) 更多的内部函数 (date, precision, arrays, ...) 新的源程序书写格式－自由格式（free format） 指针

用类型说明语句定义变量的类型 ,可以改变隐含说 明语句和“I--N”规则的约定,有六种类型语句: INTEGER(整型说明语) REAL (实型说明语句) DOUBLE PRECISION(双精度型) COMPLEX (复型说明语句) LOGICAL(逻辑型说明语句) CHARACTER(字符型说明语句) 例:INTEGER A,K3,MAXW,CL REAL AREA,KEY,MAXL 类型说明语句优先级高于隐含说明语句 ,又高于隐 含规则 说明语句是非执行语句,应出现在所有执行语句前, 隐含说明语句还必须出现在类型说明语句前。
每条READ语句都从新的输入行开始读数。
输入数据类型要正确。 数据的之间用空格或逗号或回车符来分隔。 输入数据只有在按下回车键之后才有效。在上 机操作时,一定要记得在输入数据后按回车。 当输入数据中有 /符号时 ,/后面的输入数据没有 作用;如果/前面输入的数据不够用,其余要输入 的数据已经有值的保持不变,否则,数值(算术)型 与字符型数据取零值 ; 逻辑型数据取真值 ; 输入 多余的数据是没有作用的。 输入的数据,算术型数据必须是与对应的变量有 相同类型的常量 ; 字符型数据 , 当长度不够时以 空格填充 , 当长度超过时截掉多余部分 ; 逻辑型 可用以T或F或.T或.F打头的任意字符串输入。
例:求三个数(5,10,23)的平均值。 PROGRAM TEST ! 计算三个数的平均值 A=5;B=10;C=23 ave=(A+B+C)/3 write(*,*) 'AVE=',ave END PROGRA&
&M TEST !注意续行方法
§1.4 FORTRAN 90字符集 FORTRAN90字符集包括: 1) 大写与小写英文字母

fortran90整理1语句编译1.Build—Compile：编译；Build—Build：连接；Build—Exetuce：运⾏；或单击⼯具栏相应按钮。

注意：a、保存⽂件时将⾃动创建同名的project⽂件，形成*.dsp⽂件；b、同时还将⾃动创建同名的workspace，形成*.opt和*.dsw⽂件；c、编译连接后⾃动形成Debug⽬录，该⽬录中存放编译连接后⽂件。

如：*.obj,*.lnk,*.exe等2.Free Format（⾃由格式）1.！感叹号后⾯的⽂本都是注释。

2.每⾏可以编写132个字符。

3.⾏号放在每⾏程序的最前⾯。

4.⼀⾏程序代码的最后如果是符号&，代表下⼀⾏程序会和这⼀⾏连接。

如果⼀⾏程序代码的开头是符号&，代表它会和上⼀⾏程序连接。

3.书写格式⾏的书写(⾏的长度、分⾏、续⾏)1⼀⾏可以是0~132个字符，空格有意义，2语句最长不超过2640个字符3⼀⾏可以有多个语句，⽤“;”分隔4⼀个语句可分⾏写，读⾏标记为&（放在尾部），但如为关键字，⾸尾均加&。

5最多可有511个续⾏。

4.语句的分类注释语句：！后的所有字符都被编译器忽略1.可独占⼀⾏，可在其它语句之后，a)空⾏为注释⾏（固定格式⽤C和*）2.说明语句：⽤于说明变量的类型、属性等3.可执⾏语句：输⼊、赋值、输出……5.语句有位置规定：说明语句必须出现在可执⾏语句之前，格式说明语句（FORMAT语句）除外。

6.标志符⼩结注释标志符：1⾃由格式：！固定格式：C *2语句分隔符：分号；（仅⾃由格式可以使⽤）3续⾏符：⾃由格式：&4申明标号：1到5位⽆符号整数5空格：关键字、变量、常量内部不能⽤空格，但相邻两者之间须⽤空格6.FORTRA90源程序基本结构1、FORTRAN90程序是⼀种分块结构，由若⼲个程序单元块组成：主程序、外部⼦程序、模块、块数据单元⽆论是主程序单元，还是⼦程序单元，都是独⽴的程序单位，应该独⽴编写，它们的形式相似。

使 用 变 元 关 键 字 可 提 高 程 序 的 可 读 性 , 如 内 部 函 数 CMPLX ( X [ , Y] [ , KIND] ) 的X、Y和KIND为变元关键字,可使用 下面语句调用CMPLX函数：
CMPLX (10, 20, 4)
! 按序指定参数。
CMPLX (y=20, kind=4, x=10) ! 按变元关键字指定参数。
_wrong,U.S.A.,
下面是几个使用合法名称的语句例子。
INTEGER total
!total命名了一个整型变量
SUBROUTINE example !example命名了一个过程
PROGRAM area
!area命名了一个程序单元
Lable：DO I=1,N !Lable命名了一个DO循环
第三章 FORTRAN 90基础知识
3.3 关键字/示例
END FUNCTION mul(x,y)
mul=x*y END FUNCTION
第三章 FORTRAN 90基础知识
09/10秋学期
概述 分类 描述 示例
3.4 程序单元
一个F90程序可由多个程序单元组成(至少一个主程序单元)。 程序单元由数据声明和相关操作(语句)构成,必须用END语句结束。
3.2 名称
概述 语法描述 作用域 示例
<名称>→<英gt;∣“_”∣“$”}
说明： ① 名称只能由英文字母、数字、下划线符“_”和美元
符号“$”组成。 ② 名称第一个字符必须是英文字母。
3.2 名称/语法描述
③ 名称不能超过31个字符。
第三章 FORTRAN 90基础知识
模块单元：能被其它程序单元访问的一组定义(数据类型定义、 过程定义)、过程接口定义)所构成的程序单元，其中模块子程序 允许本模块或其它程序单元调用执行。由MODULE语句开始。

1.3 控制结构..................................................................................................................................... 1
1.4 数字处理........................................................................................................................................ 1 1.5 数组处理........................................................................................................................................ 1 1.6 指针................................................................................................................................................ 1 1.7 数据结构........................................................................................................................................ 2 1.8 用户定义的类型和操作符 ............................................................................................................ 2 1.9 过程................................................................................................................................................ 3 1.10 模块(modules).............................................................................................................................. 3 1.11 输入和输出.................................................................................................................................. 3 1.12 荒废的特性.................................................................................................................................. 4 2. Fortran 77, C, C++ 和 Fortran 90 的比较............................................................................................. 5 2.1 数值健壮性................................................................................................................................... 5 2.2 数据并行化部分........................................................................................................................... 7 2.3 数据抽象..................................................................................................................................... 12 2.4 面向对象编程............................................................................................................................. 12 2.5 函数型程序设计......................................................................................................................... 12 3. 数值健壮性........................................................................................................................................... 13 3.1 数值种类(Kind)参数化 .............................................................................................................. 13 3.2 精度选择...................................................................................................................................... 14 3.3 数值多态性................................................................................................................................. 14 3.4 数值近似模型............................................................................................................................. 15 3.5 环境查询..................................................................................................................................... 15 4. 数据并行性........................................................................................................................................... 17 4.1 数组操作...................................................................................................................................... 17 4.2 数组片段..................................................................................................................................... 21 4.3 动态数组...................................................................................................................................... 22 4.4 返回赋值数组的函数(array-valued functions) .......................................................................... 24 4.5 例子：高斯消元法 ..................................................................................................................... 26 参考文献.................................................................................................................................................... 28

例1.1 求
其中
2
x y u x y
2
a b x 2 2 a b
ab ab
a b a b y 4 a b
ab ab
这一题的算法并不难，可写成： (1)从键盘输入a、b的值。 ab x a b ,y (2)如果a<b，则 ab ， 4 x a b , y 否则 。 ab (3)计算u的值。 (4)输出u的值。
3. N-S图
由于传统流程图的缺点，1973年美国学者 I.Nassi和B.Shneiderman提出了一种新的流程图工 具─N-S图。N-S图以三种基本结构作为构成算法 的基本元素，每一种基本结构用一个矩形框来表 示，而且取消了流程线，各基本结构之间保持顺 序执行关系。N-S图可以保证程序具有良好的结 构，所以N-S图又叫做结构化流程图。
例1.3的算法：
(1)输入m和n的值。 (2)求m除以n的余数r。 (3) 若r=0 ，则转至第 (6) 步，否则执行第 (4) 步。 (4)n→m，r→n。 (5)转第(2)步。 (6)输出n。
算法的五个特征：
(1) 有穷性。算法中执行的步骤总是有限次数的， 不能无止境地执行下去。 (2) 确定性。算法中的每一步操作必须具有确切 的含义，不能有二义性。 (3) 有效性。算法中的每一步操作必须是可执行 的。 (4) 要有数据输入。算法中操作的对象是数据， 因此应提供有关数据。 (5) 要有结果输出。算法的目的是用来解决一个 给定的问题，因此应提供输出结果，否则算法 就2没有实际意义。
1.2.2 算法的描述
算法的描述有许多方法，常用的有：自 然语言、一般流程图、N－S图等。前面例 1.1~例1.3的算法是用自然语言──汉语描述 的，其优点是通俗易懂，但它不太直观， 描述不够简洁，且容易产生二义性。在实 际应用中常用流程图表示算法。

本文末给出Gauss-Jordan消去法的Fortran90源程序。

!/************************************************************* !程序：Gauss_Jordan消去法!过程：Gauss_Jordan(aa,b,n,sgn)!作用：aa为方阵，b为aa的逆，n为aa的阶! sgn为标识符，1表示求逆成功，0表示求逆失败!调用格式为：call Gauss_Jordan(aa,b,n,sgn)!*************************************************************/ subroutine Gauss_Jordan(aa,b,n,sgn)implicit noneinteger(4):: n,sgnreal(8):: aa(n,n),b(n,n)integer(4):: i,j,kreal(8),allocatable:: a(:,:)real(8):: tallocate(a(n,n))a=aa ! a代替aa进行运算sgn=1! 初始化b为单位阵do i=1,ndo j=1,nif(i==j) thenb(i,j)=1elseb(i,j)=0end ifend doend do! Gauss_Jordan消去法过程do k=1,nif(a(k,k)==0) thensgn=0;EXITend if! 化第k行使得a(k,k)为1t=1.0d0/a(k,k)do j=k,na(k,j)=a(k,j)*tend dodo j=1,nb(k,j)=b(k,j)*tend do! 完成第k列的计算do i=1,nif(i/=k)thent=a(i,k)do j=k,na(i,j)=a(i,j)-a(k,j)*tend dodo j=1,nb(i,j)=b(i,j)-b(k,j)*tend doend ifend doend doend subroutine Gauss_Jordanprogram equation_solveimplicit noneinteger n,sgnparameter(n=2)real tinteger i,j,kreal aa(n,n),b(n,n),c(n,1),x(n,1)real,allocatable:: a(:,:)read*,aa,callocate(a(n,n))a=aa ! a代替aa进行运算sgn=1! 初始化b为单位阵do i=1,ndo j=1,nif(i==j) thenb(i,j)=1elseb(i,j)=0end ifend doend do! Gauss_Jordan消去法过程do k=1,nif(a(k,k)==0) thensgn=0EXITend if! 化第k行使得a(k,k)为1t=1.0/a(k,k)do j=k,na(k,j)=a(k,j)*tend dodo j=1,nb(k,j)=b(k,j)*tend do ! 完成第k列的计算do i=1,nif(i/=k)thent=a(i,k)do j=k,na(i,j)=a(i,j)-a(k,j)*tend dodo j=1,nb(i,j)=b(i,j)-b(k,j)*tend doend ifend doend dox=matmul(b,c)write(*,20) ((b(i,j),j=1,n),i=1,n)20 format(5X,2F9.3)write(*,10) (x(i,1),i=1,n)10 format(5X,F5.3)end!正确的编程，但是需要能求得出逆矩阵的矩阵方程! write(*,20) ((b(i,j),j=1,n),i=1,n)! 20 format(5X,2F9.3) !print*," aa的逆矩阵 ",b !x=matmul(b,c)!write(*,10) (x(i,1),i=1,n)!10 format(5X,F5.3) !print*," aa的逆矩阵 ",b !end。

FORTRAN语言课程设计摘要：科技的日新月异使得计算机领域不断取得新的研究成果。

计算机在代替和延伸脑力劳动方面发挥越来越重要的作用，不仅在工业方面而且在日常生活和科研中也越来越离不开计算机。

特别是在天体运动方面需要运用到计算机处理大量的数据。

这次我选的实践课题是用Jacobi迭代和Gauss-Seidel迭代法求解线性方程组AX=B，这其中涉及的就是天体运动的轨迹问题，我利用从FORTRAN 90中学到的迭代、循环、子程序等知识设计程序，通过Fortran PowerStation 4.0进行运行、调试，不得不提的是QuickWin,它在绘制行星的运动轨迹上发挥出了相当大的贡献。

通过这次的实践我从中充分体会到了Fortran语言接近数学公式的自然描述，在计算机里具有很高的执行效率的最大特性。

同时我也看到了Fortran语言是一种极具发展潜力的语言，在数值计算中，Fortran语言仍然不可替代。

Fortran90标准引入了数组计算等非常利于矩阵运算的功能。

在数组运算时，Fortran能够自动进行并行运算，这是很多编程语言不具备的。

运用Fortran 语言，你能够运用很多现成的函数软件包，所以非常便利。

关键词：Fortran ；Jacobi迭代和Gauss-Seidel迭代；天体运动1设计思想这次的课程设计我选的是第三个课题，关于求解天体的运行轨道，原题如下：●用Jacobi迭代和Gauss-Seidel迭代法求解线性方程组AX=b。

一天文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立以太阳为原点的直角坐标系，在五个不同的点对小行星作了五次观察，测得轨道上五个点的坐标数据（单位：万公里）如下表所示：由开普勒第一定律知，小行星轨道为一椭圆，椭圆的一般方程可表示为：a1x2+2a2xy+a3y2+2a4x+2a5y+1=0分别将五个点的数据代入椭圆一般方程中，得到线性方程组，求出待定系数a1,a2,a3,a4,a5。

Fortran90均匀和正态（⾼斯）分布随机数⽣成程序module ran_mod! module contains four functions! ran1 and ran2 return a uniform random number between 0-1! normal1 and normal 2 return a normal distributioncontainsfunction ran1() !returns random number between 0 - 1implicit noneinteger :: flagdouble precision :: ran1save flagdata flag /0/if(flag==0) thencall random_seed(); flag = 1endifcall random_number(ran1) ! built in fortran 90 random number function returnend function ran1function ran2(idum)implicit noneinteger, parameter :: IM1=2147483563, IM2=2147483399integer, parameter :: IMM1=IM1-1, IA1=40014, IA2=40692integer, parameter :: IQ1=53668, IQ2=52774, IR1=12211, IR2=3791integer, parameter :: NTAB=32, NDIV=1+IMM1/NTABinteger, dimension (NTAB) :: ivinteger :: idum, idum2, iy, i, j, kdouble precision, parameter :: EPS=1.2D-7, RNMX=1.d0-EPS, AM=1.d0/IM1double precision :: ran2, RNMKsave iv, iy, idum2data idum2/123456789/, iv/NTAB*0/, iy/0/if(idum.le.0)thenidum = max(-idum,1)idum2 = idumdo j = NTAB+8,1,-1k = idum/IQ1idum = IA1*(idum-k*IQ1)-k*IR1if(idum.lt.0) idum = idum + IM1if(j.le.NTAB) iv(j) = idumend doiy = iv(1)end ifk = idum/IQ1idum = IA1*(idum-k*IQ1)-k*IR1if(idum.lt.0) idum = idum + IM1k = idum2/IQ2idum2 = IA2*(idum2-k*IQ2)-k*IR2if(idum2.lt.0) idum2 = idum2 + IM2j = 1+iy/NDIViy = iv(j)-idum2iv(j) = idumif(iy.lt.1) iy = iy + IMM1ran2 = min(AM*iy,RNMX)returnend functionfunction normal1(mean,sigma) !returns a normal distribution implicit noneinteger :: flagdouble precision :: normal1, tmp, mean, sigmadouble precision :: fac, gsave, rsq, r1, r2save flag, gsavedata flag /0/if (flag.eq.0) thenrsq=2.0d0do while(rsq.ge.1.0d0.or.rsq.eq.0.0d0) ! new from for dor1 = 2.0d0 * ran1() - 1.0d0r2 = 2.0d0 * ran1() - 1.0d0rsq = r1 * r1 + r2 * r2enddofac = sqrt(-2.0d0*log(rsq)/rsq)gsave = r1 * factmp = r2 * facflag = 1elsetmp = gsaveflag = 0endifnormal1 = tmp * sigma + meanreturnend function normal1function normal2(mean,sigma)implicit noneinteger :: flagdouble precision, parameter :: pi = 3.141592653589793239double precision :: u1, u2, y1, y2, normal2, mean, sigmasave flagdata flag /0/u1 = ran1(); u2 = ran1()if (flag.eq.0) theny1 = sqrt(-2.0d0*log(u1))*cos(2.0d0*pi*u2)normal2 = mean + sigma*y1flag = 1elsey2 = sqrt(-2.0d0*log(u1))*sin(2.0d0*pi*u2)normal2 = mean + sigma*y2flag = 0endifreturnend function normal2end module ran_modThe above codes are made in Fortran 90 language, if you have any question, you may write to sealin2008@/doc/d6ddff43336c1eb91a375daa.html。

!JACOBI迭代求解方程组program jacobiimplicit noneinteger :: i,jinteger :: k !迭代次数旗标save kreal,parameter :: e=0.001integer,parameter :: n=4real :: x(n),y(n),b(n)data b /1.38,-0.34,0.67,1.52/real :: Dreal :: a(n,n)!!!千万千万要注意Fortran的数据存放是按照列来存放的，如果变成还按照习惯的数组存放方式输入数据，会导致错误!!!data a /2.52,0.39,0.55,0.23,0.95,1.69,-1.25,-1.15,1.25,-0.45,1.96,-0.45,-0.85,0.49,-0.98,2.31/ write(*,*)"*********矩阵A的形式为*********"write(*,"(1x,4f6.2,/)")aforall(i=1:n)x(i)=0end forallk=1100 D=0do i=1,ny(i)=b(i)do j=1,nif(i/=j) y(i)=y(i)-a(i,j)*x(j)end doy(i)=y(i)/a(i,i)end dodo j=1,nD=abs(x(j)-y(j))end doforall(i=1:n)x(i)=y(i)end forallif(D>=e) thenk=k+1write(*,*)"迭代次数为：",kgoto 100elsegoto 200end if200 write(*,*)"*************************************"write(*,*)"用jacobi方法解得的结果X[t]为:"write(*,"(1x,4f6.2,/)")x(:)stopend program*******************************************************************************!GAUSS-SEIDEL迭代求解方程组program jacobiimplicit noneinteger :: i,jinteger :: k !迭代次数旗标save kreal,parameter :: e=0.001integer,parameter :: n=4real :: x(n),y(n),b(n)data b /1.38,-0.34,0.67,1.52/real :: Dreal :: a(n,n)!!!千万千万要注意Fortran的数据存放是按照列来存放的，如果变成还按照习惯的数组存放方式输入数据，会导致错误!!!data a /2.52,0.39,0.55,0.23,0.95,1.69,-1.25,-1.15,1.25,-0.45,1.96,-0.45,-0.85,0.49,-0.98,2.31/ write(*,*)"*********矩阵A的形式为*********"write(*,"(1x,4f6.2,/)")aforall(i=1:n)x(i)=0end forallk=1100 D=0do i=1,ny(i)=b(i)do j=1,nif(i<j) y(i)=y(i)-a(i,j)*x(j)if(i>j) y(i)=y(i)-a(i,j)*y(j)end doy(i)=y(i)/a(i,i)end dodo j=1,nD=abs(x(j)-y(j))end doforall(i=1:n)x(i)=y(i)end forallif(D>=e) thenk=k+1write(*,*)"迭代次数为：",kgoto 100elsegoto 200end if200 write(*,*)"*************************************"write(*,*)"用Gauss-seidel方法解得的结果X[t]为:"write(*,"(1x,4f6.2,/)")x(:)write(*,*)"由运行结果可以看出，运算量G-S比JACOBI少至少一倍"stopend programimplicit nonesubroutine gauss(a,b,n,x,l,js)dimension a(n,n),x(n),b(n),js(n) double precision a,b,x,t l=1do k=1,n-1 d=0.0do i=k,ndo j=k,nif(abs(a(i,j)).gt.d)then d=abs(a(i,j)) js(k)=j is=iEndifend doif(d+1.0.eq.1.0)then l=0else if(js(k).ne.k)thendo i=k,n t=a(i,k) a(i,k)=a(i,js(k)) a(i,js(k))=tend doendifif(is.ne.k)thendo j=k,n t=a(k,j) a(k,j)=a(is,j) a(is,j)=tend do t=b(k) b(k)=b(is) b(is)=tendifendifif(l.eq.0)then write(*,100)endifdo j=k+1,n a(k,j)=a(k,j)/a(k,k)end do b(k)=b(k)/a(k,k)do i=k+1,ndo j=k+1,n a(i,j)=a(i,j)-a(i,k)*a(j,k)end do b(i)=b(i)-a(i,k)*b(k)end doend do if(abs(a(n,n))+1.0.eq.1.0)then l=0 write(*,100)end if x(n)=b(n)/a(n,n)do i=n-1,1,-1 t=0.0do j=i+1,n t=t+a(i,j)*x(j)end do x(i)=b(i)-tend do 100 format(1x,'fall') js(n)=ndo k=n,1,-1if(js(k).ne.k)then t=x(k) x(k)=x(js(k)) x(js(k))=t end if end doend do。

以上程序如图 1 所示。
水研 2010 级
第3页共6页
张丽萍
扬州大学
解线性方程组
计算流体动力学
图 1 Fortran 解线性方程组程序
水研 2010 级
第4页共6页
张丽萍
扬州大学
解线性方程组
算例 1：
方程组如下：
⎧ x1 + x2 + x3 + 2x4 = 10
⎪⎪ ⎨ ⎪
3x1 3x1 +
+ 2x2 4x2 +
张丽萍
扬大学
解线性方程组
源程序如下：
program linear_equation implicit none integer::i,j,k,imax,t real::max,n ! 用矩阵(实型动态数组)将线性方程组表示出来 real,dimension(:,:),allocatable::a,m real,dimension(:),allocatable::x print*,'请输入线性方程组的元数 n：' read*,t !给动态数组分配内存 allocate(a(t,t+1),m(t,t+1),x(t)) print*,'请依次按行输入线性方程组的增广矩阵 a：' read*,((a(i,j),j=1,t+1),i=1,t) do k =1,t-1 !下面选取每列的最大列主元素 max=abs(a(k,k)) imax=k do i=k+1,t+1 if (abs(a(i,k))>max) then max=abs(a(i,k)) imax=i end if end do !将最大列元素所在的行与第 K 行进行交换 do j=k,t+1 m(k,j)=a(k,j) a(k,j)=a(imax,j)

高斯消元法程序一、引言高斯消元法是一种求解线性方程组的常用方法，通过将线性方程组转化为行阶梯形矩阵，进而求得方程组的解。

本文将介绍高斯消元法的原理和实现过程，并给出一个简单的高斯消元法的程序示例。

二、高斯消元法原理高斯消元法的核心思想是通过矩阵的初等行变换，将线性方程组转化为行阶梯形矩阵，进而求解方程组的解。

具体的步骤如下：1. 构造增广矩阵：将线性方程组的系数矩阵和常数向量合并成一个增广矩阵。

2. 第一步消元：通过初等行变换，使得增广矩阵的第一列除第一个元素外的其他元素都变为0。

3. 逐列消元：对于增广矩阵的每一列（除第一列外），重复以下步骤：a. 找到当前列中第一个非零元素所在的行，并将该行交换到当前列的第一行。

b. 通过初等行变换，使得当前列的第一个元素下方的所有元素都变为0。

4. 回代求解：从最后一行开始，逐步回代求解出方程组的解。

三、高斯消元法程序示例下面给出一个简单的高斯消元法的程序示例，该程序使用C语言编写：```c#include <stdio.h>#define N 3 // 方程组的未知数个数void gauss_elimination(double A[N][N+1], double x[N]) { int i, j, k;double factor, tmp;// 第一步消元for (i = 0; i < N; i++) {factor = A[i][i];for (j = i; j < N + 1; j++) {A[i][j] /= factor;}for (k = i + 1; k < N; k++) {factor = A[k][i];for (j = i; j < N + 1; j++) {A[k][j] -= factor * A[i][j];}}}// 回代求解for (i = N - 1; i >= 0; i--) { tmp = 0;for (j = i + 1; j < N; j++) { tmp += A[i][j] * x[j]; }x[i] = A[i][N] - tmp;}}int main() {double A[N][N+1] = {{2, 1, -1, 8},{-3, -1, 2, -11},{-2, 1, 2, -3}};double x[N];int i;gauss_elimination(A, x);printf("方程组的解为：\n"); for (i = 0; i < N; i++) {printf("x[%d] = %f\n", i, x[i]);}return 0;}```程序中，`gauss_elimination`函数实现了高斯消元法的过程，`main`函数中给出了一个包含3个未知数的线性方程组的示例，程序输出了方程组的解。