FORTRANのブロックLU分解法テストプログラム
2003/09/08 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
実密行列を係数とする連立一次方程式Ax=bをブロックガウス消去法でLU分解し、
LUx=bの数値解を前進後退代入で求めるプログラムのテストプログラム。
4種類のLU分解プログラムを一回にテストする。
1000次元までの乱数データの行列で、結果の精度、計算速度及びMFLOPSを
各LU分解法ごとに出力する。計算する次元数(N)及びブロック長(MB)を画面より与える。
2. プログラム
C=================================================================C
C Test Program of Solver for Dense-Real Linear Equation Ax=b C
C by Blocked 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=2000, 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-BLU.txt')
CALL XCLOCK(T1,3)
EPS = 1.0E-14
C Size Inpuut
WRITE(6,*) 'Type In N(<=2000) and MB(Block Size)'
READ(5,*) N,MB
if(N.gt.ND) N = ND
if(MB.gt.N) MB = N
WRITE(1,200) N,MB
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 Block-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 Block-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 BLU1(A,B,N,ND,MB)
else
CALL BLU2(A,B,N,ND,MB,EPS,IP,IER)
end if
else if(k.le.3) then
C Partial Pivoting
CALL BLU3(A,B,N,ND,MB,EPS,IP,IER)
else
CALL BLU4(A,B,N,ND,MB,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 Block-Gauss Elimination N,MB=',2I5,/
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