|
原帖由 家有小丫头 于 2006-7-4 21:56 发表
丫头碰到麻烦老师了,出个鬼题目.请大家帮帮忙
"求线性方程组的通解"要用什么函数做.
请这里的高手门能帮帮忙给个答案,很急很急.
丫头在这里先谢过,谢谢了Sample Text
这是用fortran编的子程序,稍改一下,到c就行了
应该可以看懂吧~~顺便说一下,这是用的高斯-约当消去方法!
说明:A双精度实型二维数组,n*n,输入参数,存放方程组的系数矩阵,返回时被破坏
B双精度实型二维数组,n*M,输入、输出参数,调用时存放M组常数向量,返回M组解向量
N整形输入参数,阶数
M整形输入参数,右端常数向量的组数
L整形输出参数,若L=0,说明系数矩阵奇异,求解失败,若L!=0,正常返回求解正确
SUBROUTINE AGJDN(A,B,N,M,L,JS)
DIMENSION A(N,N),B(N,M),JS(N)
DOUBLE PRECISION A,B,D
L=1
DO 100 K=1,N
Q=0.0
DO 10 I=K,N
DO 10 J=K,N
IF (ABS(A(I,J)).GT.Q) THEN
Q=ABS(A(I,J))
JS(K)=J
IS=I
END IF
10 CONTINUE
IF (Q+1.0.EQ.1.0) THEN
WRITE(*,20)
L=0
RETURN
END IF
20 FORMAT(1X,' FAIL ')
DO 30 J=K,N
D=A(K,J)
A(K,J)=A(IS,J)
A(IS,J)=D
30 CONTINUE
DO 40 J=1,M
D=B(K,J)
B(K,J)=B(IS,J)
B(IS,J)=D
40 CONTINUE
DO 50 I=1,N
D=A(I,K)
A(I,K)=A(I,JS(K))
A(I,JS(K))=D
50 CONTINUE
DO 60 J=K+1,N
60 A(K,J)=A(K,J)/A(K,K)
DO 70 J=1,M
70 B(K,J)=B(K,J)/A(K,K)
DO 90 I=1,N
IF (I.NE.K) THEN
DO 80 J=K+1,N
80 A(I,J)=A(I,J)-A(I,K)*A(K,J)
DO 85 J=1,M
85 B(I,J)=B(I,J)-A(I,K)*B(K,J)
END IF
90 CONTINUE
100 CONTINUE
DO 110 K=N,1,-1
DO 110 J=1,M
D=B(K,J)
B(K,J)=B(JS(K),J)
B(JS(K),J)=D
110 CONTINUE
RETURN
END
[ 本帖最后由 flybaly 于 2006-7-6 13:27 编辑 ] |
|