!----------------------------------------------------------------------
!     Copyright (c) 2002 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!$BI=Bj(B  CTPACK/ISPACK-F90 $B%F%9%H%W%m%0%i%`(B
!      1 $B<!856-3&CMLdBj(B
!
!$BMzNr(B  2002/01/24  $BC]9-??0l(B
!
program ctpack_matrix

  use ctpack
  use lumatrix
  implicit none
  integer, parameter :: im=8, km=5, nm=4
  real(8), parameter :: xmin=-1.0D0, xmax=1.0D0

  real(8), dimension(0:km,0:km) :: T

  real(8), dimension(0:km,0:im) :: G

  real(8), dimension(0:km,0:km) :: alu
  real(8), dimension(0:km) :: t_b, t_f, t_fsol
  integer, dimension(0:km) :: kp
  real(8), parameter                 :: pi=3.1415926535897932385D0

  real(8), dimension(0:km) :: factor

  integer :: k, i, m

  call ctinitial(im,km,xmin,xmax)

! ---- 2 $B3,HyJ,J}Dx<0(B d^2f/dx^2 = b(x), f=0 at x=xmin,xmax
!
!  $BJ}Dx<0$KBP1~$9$kO"N)J}Dx<0$r%A%'%S%7%'%U6u4V$G:n@.$9$k(B.
!      alu * t_f = t_b
!  alu $B$O(B d^2f/dx^2 $B$N%A%'%S%7%'%U6u4V$K$*$1$k1i;;;R(B

  t = 0.0D0
  do k=0,km
     t(k,k)=1.0                                 ! $B3F@.J,$,0l$D$E$DB8:_$9$k(B
  enddo

  alu = transpose( t_dx_t(t_dx_t(t)) )          ! $B78?t9TNs(B
                                                ! $B3F@.J,$KBP$9$kHyJ,$r7W;;(B
  t_b = t_g(cos(pi*g_x/2))                      ! $B1&JU(B

  g = g_t(t)
  alu(km-1,:) = g(:,0)                          ! $B6-3&>r7o(B(x=xmax)
  alu(km,:)   = g(:,im)                         ! $B6-3&>r7o(B(x=xmin)

  t_b(km-1:km)  = 0                             ! $B6-3&CM(B

  t_fsol = t_g(-cos(pi*g_x/2)/pi**2*4)          ! $B@52r(B

  call ludecomp(alu,kp)                         ! LU $BJ,2r(B
  t_f = lusolve(alu,kp,t_b)                     ! $B2r$r5a$a$k(B

  do k=0,km
     write(6,*) t_f(k), t_fsol(k)
  enddo

  call gropn(1)
  call grfrm
  call usgrph(im+1,real(g_x), real(g_t(t_f)))
  call uulin(im+1,real(g_x), real(g_t(t_fsol)))
  call grcls

! ---- 2 $B3,HyJ,J}Dx<0(B d^2f/dx^2 = b(x), df/dx=0 at x=xmin,xmax
!
!  $BN>C<(B Neumann $B$N>l9g$K$OGH?t(B 0 $B@.J,$,Dj$^$i$J$$(B. 
!  $BJLES>r7o$rM?$($k$Y$7(B. 
!  $B0J2<$ONN0hJ?6QCM$rM?$($k>l9g(B

  ! $BA4NN0h@QJ,MQ$N(B $B%A%'%S%7%'%U78?t$KBP$9$k(B factor
  ! \int_{-1}^{1}T_n(x)dx =  -2/(n**2-1), n=even
  factor=0
  do k=0,km,2
     factor(k)=-2.0D0/(k**2-1)
  enddo
  factor(0)=factor(0)/2                    ! 0 $B<!$NOB$K$O=E$_(B 1/2

  t = 0.0D0
  do k=0,km
     t(k,k)=1.0                                 ! $B3F@.J,$,0l$D$E$DB8:_$9$k(B
  enddo

  alu = transpose( t_dx_t(t_dx_t(t)) )          ! $B78?t9TNs(B
                                                ! $B3F@.J,$KBP$9$kHyJ,$r7W;;(B

  t_b = t_g(sin(pi*g_x/2))                      ! $B1&JU(B

  g = g_t(t_dx_t(t))
  alu(km-1,:) = g(:,0)                          ! $B6-3&>r7o(B(x=xmax)
  alu(km,:)   = g(:,im)                         ! $B6-3&>r7o(B(x=xmin)
  t_b(km-1:km)  = 0                             ! $B6-3&CM(B

  alu(0,:) = factor                             ! $BGH?t(B 0 $B@.J,$rDj$a$k<0(B
  t_b(0) = 0                                    ! $BJ?6QCM$r(B 0 $B$H$9$k(B

  t_fsol = t_g(-sin(pi*g_x/2)/pi**2*4)          ! $B@52r(B

  call ludecomp(alu,kp)                         ! LU $BJ,2r(B
  t_f = lusolve(alu,kp,t_b)                     ! $B2r$r5a$a$k(B

  do k=0,km
     write(6,*) t_f(k), t_fsol(k)
  enddo

  call gropn(1)
  call grfrm
  call usgrph(im+1,real(g_x), real(g_t(t_f)))
  call uulin(im+1,real(g_x), real(g_t(t_fsol)))
  call grcls

! ---- 1 $B3,HyJ,J}Dx<0(B x df/dx = b(x), f=0 at x=(xmin+xmax)/2
!
!  $B78?t$,(B x $B$K0MB8$9$k>l9g(B : $B<B6u4V$X$NJQ49$K$h$j9TNs$r:n@.(B
!  
  t = 0.0D0
  do k=0,km
     t(k,k)=1.0                                 ! $B3F@.J,$,0l$D$E$DB8:_$9$k(B
  enddo

  alu = transpose( t_g(spread(g_x,1,km+1)*g_t(t_dx_t(t))) )   ! $B78?t9TNs(B

  t_b = t_g(g_x**2)                             ! $B1&JU(B

  g = g_t(t)
  alu(km,:)   = g(:,im/2)                       ! $B6-3&>r7o(B(x=(xmin+xmax)/2)
  t_b(km)     = 0                               ! $B6-3&CM(B

  t_fsol = t_g(g_x**2/2)                        ! $B@52r(B

  call ludecomp(alu,kp)                         ! LU $BJ,2r(B
  t_f = lusolve(alu,kp,t_b)                     ! $B2r$r5a$a$k(B

  do k=0,km
     write(6,*) t_f(k), t_fsol(k)
  enddo

  call gropn(1)
  call grfrm
  call usgrph(im+1,real(g_x), real(g_t(t_f)))
  call uulin(im+1,real(g_x), real(g_t(t_fsol)))
  call grcls

end program ctpack_matrix
