fortran中阶乘算法

因为，integer类型默认为四字节整数，最大握尘春值为2^31-1=2147483647。而13的阶乘为6227020800，已经超出兄拆integer上限。

可以改用real*8双精度类型，这样段耐n=33就不会溢出了。双精度类型的上限是10的308次方，最大可以计算到n=170，即170的阶乘。

求矩阵行列式是一个复杂的过程。有很虚晌多很多算法来做，但是各差弯锋有适用性。有的不适合病态矩阵等等。

以下是一个简单的全选主元高斯消去法。

摘自徐世良的闹铅《Fortran常用算法集》

Program Main

Implicit None

Real(8) :: rm(3,3) = reshape( (/1,2,4,5,7,3,13,5,7/) , (/3,3/) )

Real(8) :: rDet

call BSDet( rm , 3 , rDet )

write(*,*) rDet

End Program Main

SUBROUTINE BSDET(A,N,DET)

DIMENSION A(N,N)

DOUBLE PRECISION A,DET,F,D,Q

F=1.0

DET=1.0

DO 100 K=1,N-1

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))

IS=I

JS=J

END IF

10 CONTINUE

IF (Q+1.0.EQ.1.0) THEN

DET=0.0

RETURN

END IF

IF (IS.NE.K) THEN

F=-F

DO 20 J=K,N

D=A(K,J)

A(K,J)=A(IS,J)

A(IS,J)=D

20 CONTINUE

END IF

IF (JS.NE.K) THEN

F=-F

DO 30 I=K,N

D=A(I,JS)

A(I,JS)=A(I,K)

A(I,K)=D

30 CONTINUE

END IF

DET=DET*A(K,K)

DO 50 I=K+1,N

D=A(I,K)/A(K,K)

DO 40 J=K+1,N

40 A(I,J)=A(I,J)-D*A(K,J)

50 CONTINUE

100 CONTINUE

DET=F*DET*A(N,N)

RETURN

END

---