! 拡散方程式の陰解法(backword 法)
! \DP{u}{t} = \kappa \DP[2]{u}{x}
! \kappa=0.01
!
      program back

      implicit none

      real(8), parameter :: pi = 4.0d0*atan(1.0d0)
      real(8), parameter :: kappa = 0.01d0
      real(8), parameter :: xmin = 0.0d0
      real(8), parameter :: xmax = 1.0d0
      real(8), parameter :: tmin = 0.0d0
      real(8), parameter :: tmax = 10.0d0

      real(8), allocatable :: U_a(:) !U^{n+1}
      real(8), allocatable :: U_b(:) !U^{n}
      real(8), allocatable :: X_cal(:) !xの値
      real(8), allocatable :: L_mat(:,:)
      real(8), allocatable :: U_mat(:,:)
      real(8), allocatable :: A_mat(:,:)

      real, allocatable :: U_sing(:)
      real, allocatable :: X_sing(:)
      real, allocatable :: genmitu_sing(:)

      real(8) :: dx
      real(8) :: dt
      real(8) :: f_i
      real(8) :: LU

      integer :: x_step
      integer :: t_step
      integer :: x_lupe
      integer :: t_lupe
      integer :: i_lupe
      integer :: k_lupe
      integer :: j_lupe
      integer :: IWS
      integer :: scheme

      character(10) chara

      write(*,*) "Enter in time steps. & 
        & if you want to enter more 10000000000,&
        &  enter in 'a'. and rekey. default is 800."! 500.
                                             
      read(*,'(a)') chara
      !時間分割数の仮入力
                                    
      call empty(chara, 200, t_step)
      !仮文字列が空白の場合デフォルト値200を入れる
                                                 
      write(*,*) "You chose ", t_step
 
      write(*,*) "Enter in X steps. & 
        & if you want to enter more 10000000000,&
        &  enter in 'a'. and rekey. default is 10."!20 "
                       
      read(*,'(a)') chara
      !空間分割数の仮文字列入力
                              
      call empty(chara, 50, x_step)
!      call empty(chara, 3, x_step)
      !仮文字列が空白の場合デフォルト値10を入れる
                                           
      write(*,*) "You chose", x_step
      open(10,file='test.dat', status='replace')

      allocate(U_a(0:x_step))
      allocate(U_b(0:x_step))
      allocate(X_cal(0:x_step))
      allocate(U_sing(0:x_step))
      allocate(X_sing(0:x_step))
      allocate(A_mat(1:x_step-1,1:x_step-1))
      allocate(L_mat(1:x_step-1,1:x_step-1))
      allocate(U_mat(1:x_step-1,1:x_step-1))
      allocate(genmitu_sing(0:x_step))


      dx = (abs(xmax) - abs(xmin))/dble(x_step)
      dt = (abs(tmax) - abs(tmin))/dble(t_step)
      !------------------------------------------

      call swpset('LDUMP',.true.)
      call swpset('LWAIT',.false.)
      call swpset('lwait1',.false.)

      call sgpwsn
      read(*,'(a)') chara
      call empty(chara,1,IWS)
      call gropn(IWS)

!-----初期条件----------------------------------
      do x_lupe = 0, x_step
      
        X_cal(x_lupe) = xmin + dx*dble(x_lupe)

        U_b(x_lupe) = sin(pi*X_cal(x_lupe))

      end do

      X_sing=real(X_cal)
      U_sing=real(U_b)

      call grfrm
      call grswnd(0.0,1.0,0.0,1.0)

      call uspfit
      call grstrf
      call usdaxs
      call uulin(x_step+1, X_sing(0:x_step), U_sing(0:x_step))

      A_mat(:,:)=0.0d0
      U_mat(:,:)=0.0d0
      L_mat(:,:)=0.0d0
      LU = 0.0d0
      U_b(0) = U_b(x_step) 
      U_a = U_b

      do x_lupe = 1, x_step-1

         A_mat(x_lupe,x_lupe) = 1.0d0 + 2.0d0*(kappa*dt/(dx**2))
         A_mat(x_lupe,x_lupe+1) = - (kappa*dt/(dx**2))
         A_mat(x_lupe+1,x_lupe) = - (kappa*dt/(dx**2))

      end do

!      A_mat(x_step-1,x_step-1) = 1.0d0 + 2.0d0*(kappa*dt/(dx**2))
!        A_mat(1,1) = 1.0d0
!        A_mat(1,2) = 2.0d0
!        A_mat(1,3) = 3.0d0
!        A_mat(1,4) = 4.0d0
!        A_mat(2,1) = 2.0d0
!        A_mat(2,2) = 3.0d0
!        A_mat(2,3) = 4.0d0
!        A_mat(2,4) = 1.0d0
!        A_mat(3,1) = 3.0d0
!        A_mat(3,2) = 4.0d0
!        A_mat(3,3) = 1.0d0
!        A_mat(3,4) = 2.0d0
!        A_mat(4,1) = 4.0d0
!        A_mat(4,2) = 1.0d0
!        A_mat(4,3) = 2.0d0
!        A_mat(4,4) = 3.0d0
!        A_mat(1,1) = 1.0d0
!        A_mat(1,2) = 2.0d0
!        A_mat(1,3) = 3.0d0
!        A_mat(2,1) = 2.0d0
!        A_mat(2,2) = 1.0d0
!        A_mat(2,3) = 2.0d0
!        A_mat(3,1) = 3.0d0
!        A_mat(3,2) = 2.0d0
!        A_mat(3,3) = 1.0d0

!-----LU分解--------------------------------
      ! 1行目の計算
!      do i_lupe = 1,x_step
!        
!        L_mat(i_lupe,1) = A_mat(i_lupe, 1)
!
!        U_mat(1,i_lupe) = A_mat(1,i_lupe)/L_mat(1,1)
!      end do
!
!      !2行目以降
!      do k_lupe = 2, x_step - 1
!
!        L_mat(k_lupe, 2) = &
!            & A_mat(k_lupe,2) - L_mat(k_lupe,1) * U_mat(1,2) 
!
!        do i_lupe = 3, k_lupe
!
!          do j_lupe = 2, i_lupe-1
!
!            LU = LU + L_mat(k_lupe,j_lupe)*U_mat(j_lupe,i_lupe)
!
!          end do
!
!          L_mat(k_lupe,i_lupe) = & 
!            & A_mat(k_lupe,i_lupe) - L_mat(k_lupe,1) * U_mat(1,i_lupe) &
!            & - LU
!
!          LU = 0.0d0
!
!        end do
!
!
!        do i_lupe = k_lupe + 1, x_step
!
!          do j_lupe = 2, k_lupe -1
!
!            LU = LU + L_mat(k_lupe,j_lupe)*U_mat(j_lupe,i_lupe)
!          end do
!
!          U_mat(k_lupe,i_lupe) = &
!            & (A_mat(k_lupe,i_lupe) - L_mat(k_lupe,1) * U_mat(1,i_lupe)&
!            & - LU)/L_mat(k_lupe,k_lupe)
!
!        end do
!
!      end do
!
!      LU = 0.0d0
!      !x_step行目の計算
!      do j_lupe = 2, x_step-1
!        
!        LU = LU + L_mat(x_step,j_lupe)*U_mat(j_lupe,x_step)
!
!      end do        
!        L_mat(x_step,x_step) = &
!          & A_mat(x_step,x_step) - L_mat(x_step,1)*U_mat(1,x_step) &
!          & - LU
!      
!      do x_lupe = 1, x_step
!        do t_lupe = 1, x_step
!        write(10,*) (L_mat(x_lupe,t_lupe)),U_mat(x_lupe,t_lupe)
!        end do
!        end do 
      call LUdecom(x_step-1,A_mat,L_mat,U_mat)

!      do x_lupe = 1, x_step
!        do t_lupe = 1, x_step
!        write(10,*) L_mat(x_lupe,t_lupe),U_mat(x_lupe,t_lupe)
!        end do
!        end do 


!      write(10,*) "(1,1)", L_mat(1,1)*U_mat(1,1), & 
!        & "(1,2)", L_mat(1,1)*U_mat(1,2), &
!        & "(1,3)", L_mat(1,1)*U_mat(1,3)
!
!      write(10,*) "(2,1)", L_mat(2,1)*U_mat(1,1), & 
!        & "(2,2)", L_mat(2,1)*U_mat(1,2) + L_mat(2,2), &
!        & "(2,3)", L_mat(2,1)*U_mat(1,3) + L_mat(2,2)*U_mat(2,3)
!
!      write(10,*) "(3,1)", L_mat(3,1)*U_mat(1,1), & 
!        & "(3,2)", L_mat(3,1)*U_mat(1,2) + L_mat(3,2), &
!        & "(3,3)", L_mat(3,1)*U_mat(1,3) + L_mat(3,2)*U_mat(2,3) +& 
!        & L_mat(3,3)
!      write(10,*) "(1,1)", L_mat(1,1)*U_mat(1,1), & 
!        & "(1,2)", L_mat(1,1)*U_mat(1,2), &
!        & "(1,3)", L_mat(1,1)*U_mat(1,3), &
!        & "(1,4)", L_mat(1,1)*U_mat(1,4)
!      write(10,*) "(2,1)", L_mat(2,1)*U_mat(1,1), & 
!        & "(2,2)", L_mat(2,1)*U_mat(1,2) + L_mat(2,2), &
!        & "(2,3)", L_mat(2,1)*U_mat(1,3) + L_mat(2,2)*U_mat(2,3), &
!        & "(2,4)", L_mat(2,1)*U_mat(1,4) + L_mat(2,2)*U_mat(2,4)
!      write(10,*) "(3,1)", L_mat(3,1)*U_mat(1,1), & 
!        & "(3,2)", L_mat(3,1)*U_mat(1,2) + L_mat(3,2), &
!        & "(3,3)", L_mat(3,1)*U_mat(1,3) + L_mat(3,2)*U_mat(2,3) +& 
!        & L_mat(3,3), &
!        & "(3,4)", L_mat(3,1)*U_mat(1,4) + L_mat(3,2)*U_mat(2,4) +&
!        & L_mat(3,3)*U_mat(3,4)
!      write(10,*) "(4,1)", L_mat(4,1)*U_mat(1,1), & 
!        & "(4,2)", L_mat(4,1)*U_mat(1,2) + L_mat(4,2), &
!        & "(4,3)", L_mat(4,1)*U_mat(1,3) + L_mat(4,2)*U_mat(2,3) +& 
!        & L_mat(4,3), &
!        & "(4,4)", L_mat(4,1)*U_mat(1,4) + L_mat(4,2)*U_mat(2,4) +&
!        & L_mat(4,3)*U_mat(3,4) + L_mat(4,4)
!      do x_lupe = 1,3
!
!       write(10,*) "(",x_lupe,"1)", L_mat(x_lupe,1), &
!         & "(",x_lupe,"2)", L_mat(x_lupe,2), &
!         & "(",x_lupe,"3)", L_mat(x_lupe,3)
!
!      end do
!      write(10,*) ""
!      do x_lupe = 1,3
!
!       write(10,*) "(",x_lupe,"1)", U_mat(x_lupe,1), &
!         & "(",x_lupe,"2)", U_mat(x_lupe,2), &
!         & "(",x_lupe,"3)", U_mat(x_lupe,3)
!
!      end do
      do t_lupe = 1, t_step
      !xyの計算

        U_a(1) = U_b(1)/L_mat(1,1)
        do x_lupe = 2, x_step-1

          LU = 0.0d0
    
          do j_lupe = 1, x_lupe - 1

            LU = LU + L_mat(x_lupe,j_lupe) * U_a(j_lupe)
        
          end do

          U_a(x_lupe) = (U_b(x_lupe) - LU)/L_mat(x_lupe,x_lupe)

        end do

!        U_b(1:x_step-1) = U_a(x_step-1:1:-1)
!        do x_lupe = 2, x_step-1
!
!          LU = 0.0d0
!      
!          do j_lupe = 1, x_lupe-1 
!
!            LU = U_mat(x_lupe,j_lupe) * U_b(j_lupe)
!
!          end do
!
!          U_b(x_lupe) = U_a(x_lupe) - LU
!
!        end do

        do x_lupe = x_step-1, 1, -1

          LU = 0.0d0
      
          do j_lupe = x_step-1, x_lupe+1, -1 

            LU = U_mat(x_lupe,j_lupe) * U_b(j_lupe)

          end do

          U_b(x_lupe) = U_a(x_lupe) - LU

        end do
        U_b(0) = 0.0d0
        U_b(x_step) = 0.0d0
!        U_a(0) = U_b(0)
!        U_a(x_step) = U_b(x_step)
!        U_a(1:x_step-1) = U_b(x_step-1:1:-1)
        
!        U_sing = real(U_a)
        U_sing = real(U_b)
!        U_b = U_a
        do x_lupe = 0, x_step
          genmitu_sing(x_lupe) = &
            & sin(pi*X_sing(x_lupe))*exp(-(pi**2)*kappa*dt*t_lupe)
        end do

        call grfrm
        call grswnd(0.0,1.0,0.0,1.0)
        call usspnt(x_step+1,X_sing(0:x_step),U_sing(0:x_step))
        call usspnt(x_step+1,X_sing(0:x_step),genmitu_sing(0:x_step))


        call uspfit
        call grstrf
        call usdaxs
        call uuslnt(1)
        call uulin(x_step+1, X_sing(0:x_step), U_sing(0:x_step))
        call uuslnt(2)
        call uulin(x_step+1, X_sing(0:x_step), genmitu_sing(0:x_step))

      end do

      call grcls

      close(10)
      deallocate(U_a,U_b,X_sing,X_cal,U_sing,L_mat,A_mat,U_mat)
      stop

      end program back

!-----LU分解-------------------------------------
      subroutine LUdecom(x_step, A_mat,L_mat,U_mat)

      implicit none

      integer k_lupe
      integer i_lupe
      integer j_lupe
      integer, intent(in) :: x_step

      real(8), intent(in) :: A_mat(1:x_step,1:x_step)
      real(8), intent(out) :: L_mat(1:x_step,1:x_step)
      real(8), intent(out) :: U_mat(1:x_step,1:x_step)
      real(8) :: LU

      LU = 0.0d0

!      do k_lupe = 1, x_step
!
!        U_mat(k_lupe,k_lupe) = 1.0d0
!
!      end do

      do i_lupe = 1,x_step
        
        L_mat(i_lupe,1) = A_mat(i_lupe, 1)

        U_mat(1,i_lupe) = A_mat(1,i_lupe)/L_mat(1,1)
      end do

      !2行目以降
      do k_lupe = 2, x_step

        L_mat(k_lupe, 2) = &
          & A_mat(k_lupe,2) - L_mat(k_lupe,1) * U_mat(1,2) 

!       U_mat(k_lupe, 3) = &
!          & (A_mat(k_lupe, 3) - L_mat(k_lupe,1)*U_mat(1,3))/ l_mat(
        do i_lupe = 3, k_lupe

          do j_lupe = 2, i_lupe-1

            LU = LU + L_mat(k_lupe,j_lupe)*U_mat(j_lupe,i_lupe)

          end do

          L_mat(k_lupe,i_lupe) = & 
            & A_mat(k_lupe,i_lupe) - L_mat(k_lupe,1) * U_mat(1,i_lupe) &
            & - LU

          LU = 0.0d0

        end do

        do i_lupe = k_lupe + 1, x_step

          do j_lupe = 2, k_lupe -1

            LU = LU + L_mat(k_lupe,j_lupe)*U_mat(j_lupe,i_lupe)

          end do

          U_mat(k_lupe,i_lupe) = &
            & (A_mat(k_lupe,i_lupe) - L_mat(k_lupe,1) * U_mat(1,i_lupe)&
            & - LU)/L_mat(k_lupe,k_lupe)

          LU = 0.0d0

        end do

      end do

!      LU = 0.0d0
      !x_step行目の計算
!      do j_lupe = 2, x_step -1
!        
!        LU = LU + L_mat(x_step,j_lupe)*U_mat(j_lupe,x_step)
!
!      end do
!
!      L_mat(x_step,x_step) = &
!        & A_mat(x_step,x_step) - L_mat(x_step,1)*U_mat(1,x_step) &
!        & - LU
!
      return
      end subroutine LUdecom
!-----入力判断サブルーチン-------------------------
!入力診断のサブルーチン.
!キーボード入力が空白のときはデフォルト値を入れる
!
      subroutine empty(chara,default_step,step)

      !------入力診断の宣言文----------------------
      implicit none

      character(10) chara

      integer default_step
      integer step

      !------入力診断------------------------------
      if(len_trim(chara) == 0)then

        step = default_step
        !空白の場合デフォルト値を代入
      else

        if(chara == "a")then

          read(*,*) step
          !aの場合さらにもう一度代入
        else

          read(chara,*) step
          !仮文字列から実際の変数への代入

        end if     
             
      end if
      return
      end subroutine empty
