FORTRANのLU分解法テストプログラム
2003/08/18 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
実密行列を係数とする連立一次方程式Ax=bをガウス消去法でLU分解し、
LUx=bの数値解を前進後退代入で求めるプログラムのテストプログラム。
4種類のLU分解プログラムを一回にテストする。
1000次元までの乱数データの行列で、結果の精度、計算速度及びMFLOPSを
各LU分解法ごとに出力する。計算する次元数(N)を画面より与える。
2. プログラム
C=================================================================C
C Test Program of Solver for Dense-Real Linear Equation Ax=b C
C by Gauss Elimination with LU Decomposition C
C-----------------------------------------------------------------C
C Written by Yasunori Ushiro, 2003/08/17 C
C ( Hitachi Ltd. and Waseda University ) C
C=================================================================C
PARAMETER(ND=1000, NO=4)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(ND,ND), B(ND), IP(ND)
DIMENSION WA(ND,ND), WB(ND)
C Open File
OPEN(unit=1,FILE='F-GLU.txt')
CALL XCLOCK(T1,3)
EPS = 1.0E-14
C Size Inpuut
WRITE(6,*) 'Type In Matrix Size(N<=1000)'
READ(5,*) N
if(N.gt.ND) N=ND
WRITE(1,200) N
C Set WA,WB
SIG = 0.15
IDA = 1
IDB = 0
CALL DRMAT(WA,WB,N,ND,IDA,IDB,SIG)
C Loop for All Type LU Decomposition
do k=1,NO
C Set A,B
do j=1,N
B(j) = WB(j)
do i=1,N
A(i,j) = WA(i,j)
end do
end do
C LU Decompostion and Solve Ax=b
CALL XCLOCK(T1,5)
if(k.le.2) then
C No Pivoting
if(k.le.1) then
IER = 0
CALL GLU1(A,B,N,ND)
else
CALL GLU2(A,B,N,ND,EPS,IER)
end if
else if(k.le.3) then
C Partial Pivoting
CALL GLU3(A,B,N,ND,EPS,IP,IER)
else
CALL GLU4(A,B,N,ND,EPS,IP,IER)
end if
CALL XCLOCK(T2,5)
C Output
CPU = T2 - T1
FLOP = 0.0
IF(CPU.gt.0.0) FLOP = 2.0*N/3.0*N*N*1.0D-6/CPU
PRC = 0.0
SUM = 0.0
do i=1,N
PRC = PRC + (B(i) - 1.0)**2
SUM = SUM + B(i)**2
end do
PRC = DSQRT(PRC/SUM)
WRITE(1,300) k,IER,PRC,CPU,FLOP
WRITE(6,300) k,IER,PRC,CPU,FLOP
end do
C
stop
200 FORMAT(' Solve Ax=b by Gauss Elimination N=',I5,/
1 ,'Type. IER Precision Time(s) MFLOPS')
300 FORMAT(1H ,2I4,1X,E10.2,F9.1,F9.1)
end
C=================================================================C
SUBROUTINE XCLOCK(CPU,ID)
C=================================================================C
C CPU time Subroutine C
C CPU R*8 Out, CPU Time C
C ID I*4 In, Dummy ( Same Hitachi FORTRAN XCLOCK ) C
C-----------------------------------------------------------------C
C Written by Y. Ushiro ( Waseda and Hitachi ) 2002.10.29 C
C=================================================================C
REAL*8 CPU
INTEGER*2 I1, I2, I3, I4
C
IF(ID.GE.1) THEN
CALL GETTIM(I1,I2,I3,I4)
CPU = ( I1*60.0 + I2 )*60.0 + I3 + I4*0.01
END IF
C
RETURN
END