lu.f90

使用WENO格式求解流体力学中的激波问题

!对矩阵A(N,N)进行LU分解
!无法分解者返回0,成功返回1
FUNCTION DECOMP(A,B,N)
IMPLICIT NONE
INTEGER R,I,K,N,DECOMP,IMAX,B(N),IP(N)
DOUBLE PRECISION FMAX,A(N,N),S,EPS,T(N)
PARAMETER(EPS=0.01)
DECOMP=1
R=1
DO WHILE(R<=N)
  FMAX=0
  IMAX=0
  DO I=R,N
    S=0
	DO K=1,R-1
	  S=S+A(I,K)*A(K,R)
	END DO
	A(I,R)=A(I,R)-S
	IF(ABS(A(I,R))>FMAX)THEN
	  FMAX=ABS(A(I,R))
	  IMAX=I
	END IF
  END DO
  IF(FMAX==0)THEN
    R=N+1
	DECOMP=0
  ELSE
    IP(R)=IMAX
    IF(IMAX/=R)THEN
      T(1:N)=A(R,1:N)
	  A(R,1:N)=A(IMAX,1:N)
	  A(IMAX,1:N)=T(1:N)
    END IF
    DO I=R+1,N
      A(I,R)=A(I,R)/A(R,R)
    END DO
    DO I=R+1,N
      S=0
	  DO K=1,R-1
	    S=S+A(R,K)*A(K,I)
	  END DO
	  A(R,I)=A(R,I)-S
    END DO
	R=R+1
  END IF
END DO
IF(DECOMP/=0)THEN
  DO R=1,N-1
    IF(IP(R)/=R)THEN
      K=B(R)
	  B(R)=B(IP(R))
	  B(IP(R))=K
    END IF
  END DO
END IF
END FUNCTION DECOMP