!c Description:
!c   係数行列の作成
!c
!c Current Code Owner: 
!c   sugiyama@gfd-dennou.org
!c
!c Histry: 
!c   Version    Date          Comment
!c   -------    ----------    --------
!c   1.0        2003-02-25    作成
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2003, All rights reserved

module interface_matrix_coeff
  interface matrix_coeff
     
     subroutine matrix_coeff_gen(mol, phase, struct, A)
       real(8), intent(in)     :: mol(:)      !化学種のモル数
       real(8), intent(in)     :: struct(:,:) !化学種に含まれる元素数
       real(8), intent(out)    :: A(:,:)      !係数行列
       integer, intent(in)     :: phase(:)    !化学種の相番号
     end subroutine matrix_coeff_gen
     
     subroutine matrix_coeff_sym(mol, phase, struct, A)
       real(8), intent(in)     :: mol(:)      !化学種のモル数
       real(8), intent(in)     :: struct(:,:) !化学種に含まれる元素数
       real(8), intent(out)    :: A(:)        !係数行列
       integer, intent(in)     :: phase(:)    !化学種の相番号
     end subroutine matrix_coeff_sym
     
  end interface
end module interface_matrix_coeff




subroutine matrix_coeff_gen(mol, phase, struct, A)
  use interface_phase_count
  use interface_phase_pack
  
  implicit none
  real(8), intent(in)  :: mol(:)      !化学種のモル数
  real(8), intent(in)  :: struct(:,:) !化学種に含まれる元素数
  real(8), intent(out) :: A(:,:)      !係数行列
  integer, intent(in)  :: phase(:)    !化学種の相番号

  integer, allocatable :: phase_tmp(:)  
  integer              :: smax        !存在する化学種の数
  integer              :: pmax        !存在する相の数
  integer              :: emax        !存在する元素の数
  integer              :: e
  integer              :: f
  integer              :: p
  integer              :: max
  integer              :: min
  
  
!!!
!!!変数の初期化
!!!
  smax = size(struct, 1)
  emax = size(struct, 2)
  call phase_count(phase, pmax)
  A = 0.0D0
  
  allocate(phase_tmp(smax))
  phase_tmp = phase
  call phase_pack(phase_tmp)

!!!
!!!行列要素の計算
!!!
  do f = 1, emax
     do e = 1, emax
        A(f, e) = sum(struct(1:smax, e) * struct(1:smax, f) * mol(1:smax))
     end do
  end do
  
  do p = 1, pmax  
     min = minloc(phase_tmp, 1, phase_tmp == p)
     max = min + count(phase_tmp == p) - 1
     
     do e = 1, emax        
        A(emax + p, e) = sum(struct(min:max, e) * mol(min:max))     
     end do
     A(:, emax + p) = A(emax + p, :)
  end do

!!!
!!!後片付け
!!!
  deallocate(phase_tmp)

end subroutine matrix_coeff_gen



subroutine matrix_coeff_sym(mol, phase, struct, A)
  use interface_phase_count
  use interface_phase_pack
  
  implicit none
  real(8), intent(in)  :: mol(:)      !化学種のモル数
  real(8), intent(in)  :: struct(:,:) !化学種に含まれる元素数
  real(8), intent(out) :: A(:)        !係数行列
  integer, intent(in)  :: phase(:)    !化学種の相番号

  real(8), allocatable :: FA(:,:)      !係数行列
  integer, allocatable :: phase_tmp(:)  
  integer              :: smax        !存在する化学種の数
  integer              :: pmax        !存在する相の数
  integer              :: emax        !存在する元素の数
  integer              :: e
  integer              :: f
  integer              :: p
  integer              :: max
  integer              :: min
  
  
!!!
!!!変数の初期化
!!!
  smax = size(struct, 1)
  emax = size(struct, 2)
  call phase_count(phase, pmax)
  allocate(FA(emax+pmax, emax+pmax))
  A = 0.0D0; FA = 0.0d0
  
  allocate(phase_tmp(smax))
  phase_tmp = phase
  call phase_pack(phase_tmp)

!!!
!!!行列要素の計算
!!!
  do f = 1, emax
     do e = 1, emax
        FA(f, e) = sum(struct(1:smax, e) * struct(1:smax, f) * mol(1:smax))
     end do
  end do
  
  do p = 1, pmax  
     min = minloc(phase_tmp, 1, phase_tmp == p)
     max = min + count(phase_tmp == p) - 1
     
     do e = 1, emax        
        FA(emax + p, e) = sum(struct(min:max, e) * mol(min:max))     
     end do
     FA(:, emax + p) = FA(emax + p, :)
  end do


!!!
!!!対象行列圧縮モードに変換 (SSL2)
!!!
  e = 1; f = 1
  do p = 1, emax+pmax
     A(e:f) = FA(p, 1:p)
     e = f + 1
     f = f + p + 1     
  end do


!!!
!!!後片付け
!!!
  deallocate(phase_tmp)

end subroutine matrix_coeff_sym

