* PACKAGE UMTLU !" ＬＵ分解による連立１次方程式の解 * *" [HIS] 90/08/31(numaguti) * ********************************************************************** SUBROUTINE LUMAKE !" 行列のＬＵ分解（部分ピボット選択） M ( ALU , O KP , D JDIM , NDIM ) * *" ＬＵ行列は 入力行列に上書きされる. * * [PARAM] INTEGER JDIM INTEGER NDIM * * [MODIFY] REAL ALU ( JDIM, NDIM, NDIM ) !" 入力／ＬＵ行列 * * [OUTPUT] INTEGER KP ( JDIM, NDIM ) !" ピボット * * [INTERNAL WORK] INTEGER J, K, M, N REAL PIVOT, TEMP * * DO 3000 K = 1, NDIM-1 DO 3000 J = 1, JDIM * *" < 1. 部分ピボット選択 > * * PIVOT = ALU ( J,K,K ) KP ( J,K ) = K DO 1100 M = K+1, NDIM IF ( ABS( ALU ( J,M,K ) ) .GT. ABS( PIVOT ) ) THEN PIVOT = ALU ( J,M,K ) KP ( J,K ) = M ENDIF 1100 CONTINUE * IF ( KP ( J,K ) .NE. K ) THEN DO 1200 N = 1, NDIM TEMP = ALU ( J,K,N ) ALU ( J,K,N ) = ALU ( J,KP(J,K),N ) ALU ( J,KP(J,K),N ) = TEMP 1200 CONTINUE ENDIF * *" < 2. ＬＵ分解 > * DO 2100 N = K+1, NDIM *" U[kj], j>k+1 ALU ( J,K,N ) = ALU( J,K,N ) / PIVOT * DO 2100 M = K+1, NDIM *" L[i,k+1], i>=k+1 ALU ( J,M,N ) = ALU( J,M,N ) & - ALU( J,M,K ) * ALU( J,K,N ) * 2100 CONTINUE 2200 CONTINUE * 3000 CONTINUE * RETURN END *********************************************************************** SUBROUTINE LUSOLV !" ＬＵ分解による解の計算 M ( XV , I ALU , KP , D IDIM , JDIM , NDIM ) * *" 解は右辺の入力ベクトルに上書きされる * * [PARAM] INTEGER IDIM INTEGER JDIM INTEGER NDIM * * [MODIFY] REAL XV ( IDIM, JDIM, NDIM ) !" 右辺ベクトル／解 * * [INPUT] REAL ALU ( JDIM, NDIM, NDIM ) !" ＬＵ行列 INTEGER KP ( JDIM, NDIM ) !" ピボット * * [INTERNAL WORK] INTEGER I, J, K, N, NN REAL TEMP * * *" < 1. ピボット選択による並び換え > * DO 1100 K = 1, NDIM-1 DO 1100 I = 1, IDIM * DO 1110 J = 1, JDIM IF ( KP ( J,K ) .NE. K ) THEN TEMP = XV ( I,J,K ) XV ( I,J,K ) = XV ( I,J,KP(J,K) ) XV ( I,J,KP(J,K) ) = TEMP ENDIF 1110 CONTINUE * 1100 CONTINUE * *" < 2. 前進代入 > * DO 2100 N = 1, NDIM *" Y[i] DO 2110 I = 1, IDIM DO 2110 J = 1, JDIM XV ( I,J,N ) = XV ( I,J,N ) / ALU ( J,N,N ) 2110 CONTINUE * DO 2130 NN = N+1, NDIM DO 2120 I = 1, IDIM DO 2120 J = 1, JDIM XV ( I,J,NN ) = XV ( I,J,NN ) & - XV ( I,J,N ) * ALU ( J,NN,N ) 2120 CONTINUE 2130 CONTINUE * 2100 CONTINUE * *" < 3. 後退代入 > * DO 3100 K = NDIM-1, 1, -1 DO 3100 N = K+1, NDIM *" X[k] DO 3110 I = 1, IDIM DO 3110 J = 1, JDIM XV ( I,J,K ) = XV ( I,J,K ) & - XV ( I,J,N ) * ALU ( J,K,N ) 3110 CONTINUE * 3100 CONTINUE * RETURN END ********************************************************************** SUBROUTINE LUMAK3 !" 行列のＬＵ分解 [ ３重対角行列 ] M ( ALU , D JDIM , NDIM ) * *" ＬＵ行列は 入力行列に上書きされる * * [PARAM] INTEGER JDIM INTEGER NDIM * * [MODIFY] REAL ALU ( JDIM, NDIM, -1:1 ) !" 入力／ＬＵ行列 * * [INTERNAL WORK] INTEGER J, K * * DO 1000 J = 1, JDIM ALU ( J,1,1 ) = ALU ( J,1,1 ) / ALU ( J,1,0 ) 1000 CONTINUE * DO 1100 K = 2, NDIM-1 DO 1100 J = 1, JDIM ALU ( J,K,0 ) = ALU ( J,K,0 ) & - ALU ( J,K,-1 ) * ALU ( J,K-1,1 ) ALU ( J,K,1 ) = ALU ( J,K,1 ) / ALU ( J,K,0 ) 1100 CONTINUE * DO 1200 J = 1, JDIM ALU ( J,NDIM,0 ) = ALU ( J,NDIM,0 ) & - ALU ( J,NDIM,-1 ) * ALU ( J,NDIM-1,1 ) 1200 CONTINUE * RETURN END *********************************************************************** SUBROUTINE LUSOL3 !" ＬＵ分解による解の計算 [ ３重対角行列 ] M ( XV , I ALU , D IDIM , JDIM , NDIM ) * *" 解は右辺の入力ベクトルに上書きされる * * [PARAM] INTEGER IDIM INTEGER JDIM INTEGER NDIM * * [MODIFY] REAL XV ( IDIM, JDIM, NDIM ) !" 右辺ベクトル／解 * * [INPUT] REAL ALU ( JDIM, NDIM, -1:1 ) !" ＬＵ行列 * * [INTERNAL WORK] INTEGER I, J, K, L * *" < 1. 前進代入 > * DO 1000 I = 1, IDIM DO 1000 J = 1, JDIM XV ( I,J,1 ) = XV ( I,J,1 ) / ALU ( J,1,0 ) 1000 CONTINUE * DO 1100 L = 2, NDIM *" Y[l] DO 1110 I = 1, IDIM DO 1110 J = 1, JDIM XV ( I,J,L ) = ( XV ( I,J,L ) & - XV ( I,J,L-1 ) * ALU ( J,L,-1 ) ) & / ALU ( J,L,0 ) 1110 CONTINUE 1100 CONTINUE * *" < 2. 後退代入 > * DO 2000 K = NDIM-1, 1, -1 *" X[k] DO 2010 I = 1, IDIM DO 2010 J = 1, JDIM XV ( I,J,K ) = XV ( I,J,K ) & - XV ( I,J,K+1 ) * ALU ( J,K,1 ) 2010 CONTINUE 2000 CONTINUE * RETURN END