ブロックLU分解FORTRANプログラムNO.4
2003/09/08 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
実密行列を係数とする連立一次方程式Ax=bをブロックガウス消去法でLU分解し、
LUx=bの数値解を前進後退代入で求める。
部分軸交換付で、軸交換は枢軸の右側だけを交換する。
ブロックサイズ(MB×MB)の単位に主力部の計算をするキャッシュ対策が組み込まれている
2. プログラム
C=================================================================C
SUBROUTINE BLU4(A,B,N,ND,MB,EPS,IP,IER)
C=================================================================C
C Real-Dense LU Decomposition by Block Gauss Elimination C
C and Solve Ax=b by Substitution C
C A ---> LU Decomposition with Partial Pivoting C
C Type-2 Partial Pivoting (Changing partial rows) C
C C=C=A*B is used by Blocked Mult C
C-----------------------------------------------------------------C
C A(ND,N) R*8, I/O, Input and Output Matrix C
C B(N) R*8, I/O, Right vector and Solution vector C
C N I*4, In, Matrix Order of A C
C ND I*4, In, Array Size of A ( ND >= N ) C
C MB I*4, In. Block Size C
C EPS R*8, In, Value for Singularity check C
C IP(N) I*4, Out, Pivot Number C
C IER I*4, Out, 0 : Normal Execution C
C 1 : Singular Stop C
C 2 ; Parameter Error C
C-----------------------------------------------------------------C
C Written by Yasunori Ushiro , 2003/09/04 C
C ( Hitachi Ltd. and Waseda University ) C
C Ver.2 Tuning by Blocked Multiplication C
C=================================================================C
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(ND,N), B(N), IP(N)
C----- Block Gauss Elimination Step ------------
C Check Parameter
IER = 0
if(ND.lt.N .or. MB.lt.1) then
IER = 2
go to 100
end if
do kk=1,N,MB
KK1 = kk + MB
NLG = N - KK1 + 1
MBK = MIN(kk+MB-1,N)
C Block Process-1
do k=kk,MBK
C Search Maximum Value in k's column
KPIV = k
PIV = abs(A(k,k))
do i=k+1,N
if(abs(A(i,k)).gt.PIV) then
KPIV = i
PIV = abs(A(i,k))
end if
end do
C Check Singularity
if(PIV.lt.EPS) then
IER = 1
go to 100
end if
IP(k) = KPIV
C Change A(k,*) <--> A(KPIV,*)
if(KPIV.ne.k) then
do j=kk,N
AKJ = A(k,j)
A(k,j) = A(KPIV,j)
A(KPIV,j) = AKJ
end do
end if
C Pivot Value
DPIV = 1.0/A(k,k)
A(k,k) = DPIV
C A(*,k)=A(*,k)/A(k,k)
do i=k+1,N
A(i,k) = A(i,k)*DPIV
end do
C Lower-Block Elimination
do j=k+1,MBK
AKJ = A(k,j)
do i=k+1,N
A(i,j) = A(i,j) - AKJ*A(i,k)
end do
end do
end do
C Block Process-2 ( Upper-Block Elimination )
do k=MBK+1,N
do j=kk,MBK-1
AJK = A(j,k)
do i=j+1,MBK
A(i,k) = A(i,k) - AJK*A(i,j)
end do
end do
end do
C Block Process-3 ( Main Elimination by C=C-A*B) )
if(NLG.ge.1) then
CALL BMULT(A(KK1,kk),A(kk,KK1),A(KK1,KK1),NLG,MB,ND)
end if
end do
C---- Solve LUx=b by Substitution -------------
C Interchange Entries of B
do jj=1,N-1,MB
ME = MIN(jj+MB-1,N-1)
do j=jj,ME
k = IP(j)
BW = B(k)
B(k) = B(j)
B(j) = BW
end do
C Forward Substitution
do j=jj,ME
BJ = B(j)
do i=j+1,N
B(i) = B(i) - BJ*A(i,j)
end do
end do
end do
C Backword Substitution
do j=N,1,-1
BJ = B(j)*A(j,j)
B(j) = BJ
do i=1,j-1
B(i) = B(i) - BJ*A(i,j)
end do
end do
C
100 continue
C
RETURN
END
C=================================================================C
SUBROUTINE BMULT(A,B,C,N,MB,ND)
C=================================================================C
C C=C-A*B by Block Multiplication with MB Block Size C
C-----------------------------------------------------------------C
C A(ND,MB) R*8, In, 1st Input Matrix, Size=N*MB C
C B(ND,N) R*8, In, 2nd Input Matrix, Size=MB*N C
C C(ND,N) R*8, I/O, Input/Output Matrix, Size=N*N C
C N I*4, In, Column Size of Matrix A C
C MB I*4, In, Block Size and Row Size of Matrix A C
C ND I*4, In, Array Size of A,B,C ( ND >= N ) C
C-----------------------------------------------------------------C
C Written by Yasunori Ushiro and Shouko Sakuma, 2003/09/04 C
C ( Hitachi Ltd. and Waseda University ) C
C=================================================================C
IMPLICIT REAL*8(A-H,O-z)
DIMENSION A(ND,MB), B(ND,N), C(ND,N)
C Outer Loop
do ii=1,N,MB
MC = MIN(MB,N-ii+1)
C Main Loop
do jj=1,N,MB
MR = MIN(MB,N-jj+1)
C C=C-A*B
do j=1,MR
jp = j+jj-1
do k=1,MB
BKJ = B(k,jp)
do i=1,MC
C(i+ii-1,jp) = C(i+ii-1,jp) - BKJ*A(i+ii-1,k)
end do
end do
end do
end do
end do
C
RETURN
END