* PACKAGE UMTLU   
*
*"  [HIS] 90/08/31(numaguti)
*
**********************************************************************
      SUBROUTINE LUMAKE    
     M         ( ALU   ,
     O           KP    ,
     D           JDIM  , NDIM   )
*
*"    ＬＵ行列は 入力行列に上書きされる.
*
*   [PARAM] 
      INTEGER    JDIM
      INTEGER    NDIM
*
*   [MODIFY] 
      ReaL*8       ALU ( JDIM, NDIM, NDIM )     
*
*   [OUTPUT] 
      INTEGER    KP  ( JDIM, NDIM )           
*
*   [INTERNAL WORK] 
      INTEGER    J, K, M, N
      ReaL*8       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*8       XV  ( IDIM, JDIM, NDIM )     
*
*   [INPUT] 
      ReaL*8       ALU ( JDIM, NDIM, NDIM )     
      INTEGER    KP  ( JDIM, NDIM )           
*
*   [INTERNAL WORK] 
      INTEGER    I, J, K, N, NN
      ReaL*8       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*8       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*8       XV  ( IDIM, JDIM, NDIM )     
*
*   [INPUT] 
      ReaL*8       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
