Line data    Source code

       1             : module m_unify_zmat
       2             :    use m_types
       3             : contains
       4             : 
       5           0 :    subroutine unify_zmat(eigval, z)
       6             :       implicit none
       7             :       real, intent(in) :: eigval(:)
       8             :       type(t_mat)      :: z 
       9             : 
      10             :       integer :: beg_group, end_group, n_g
      11           0 :       integer, allocatable :: groups(:)
      12             : 
      13           0 :       groups = make_groups(eigval)
      14             :    
      15           0 :       do n_g = 1, size(groups)
      16           0 :          beg_group = sum(groups(1:n_g-1)) + 1
      17           0 :          end_group = sum(groups(1:n_g))
      18           0 :          if(end_group <= z%matsize2) call unify_group(beg_group, end_group, z)
      19             :       enddo
      20           0 :    end subroutine unify_zmat
      21             : 
      22           0 :    subroutine unify_group(beg_group, end_group, z)
      23             :       use m_legendre_poly
      24             :       implicit none
      25             :       integer, intent(in)        :: beg_group, end_group
      26             :       type(t_mat), intent(inout) :: z 
      27           0 :       type(t_mat)   :: targ, lhs, new_basis
      28             :       real, allocatable :: x(:)
      29             :       integer :: nrhs
      30             : 
      31           0 :       nrhs = (end_group - beg_group) + 1
      32             : 
      33             :       ! leg poly can be symmetric and we don't want that
      34           0 :       x = linspace(-1.1, 1.0, z%matsize1)
      35             : 
      36           0 :       call targ%init(z%l_real, z%matsize1, nrhs)
      37           0 :       call set_sin_targ(targ)
      38             : 
      39           0 :       call lhs%init(z%l_real, z%matsize1, nrhs )
      40           0 :       if(lhs%l_real) then 
      41           0 :          lhs%data_r = z%data_r(:,beg_group:end_group)
      42             :       else 
      43           0 :          lhs%data_c = z%data_c(:,beg_group:end_group)
      44             :       endif
      45             : 
      46           0 :       call lhs%leastsq(targ)
      47             : 
      48           0 :       if(lhs%l_real) then 
      49           0 :          lhs%data_r = z%data_r(:,beg_group:end_group)
      50             :       else 
      51           0 :          lhs%data_c = z%data_c(:,beg_group:end_group)
      52             :       endif
      53             : 
      54           0 :       call new_basis%init(lhs) 
      55           0 :       call lhs%multiply(targ, res=new_basis)
      56           0 :       call mod_gram_schmidt(new_basis)
      57             : 
      58           0 :       if(lhs%l_real) then 
      59           0 :          z%data_r(:,beg_group:end_group) = new_basis%data_r
      60             :       else 
      61           0 :          z%data_c(:,beg_group:end_group) = new_basis%data_c
      62             :       endif
      63           0 :    end subroutine
      64             : 
      65           0 :    subroutine set_sin_targ(targ)
      66             :       use m_constants
      67             :       implicit none 
      68             :       type(t_mat), intent(inout) :: targ 
      69             :       integer :: i 
      70             :       real, allocatable :: x(:) 
      71             : 
      72           0 :       x = linspace(0.0,4.0, targ%matsize1)
      73           0 :       do i = 1, targ%matsize2 
      74           0 :          if(targ%l_real) then 
      75           0 :             targ%data_r(:,i) = sin(i*x)
      76             :          else 
      77           0 :             targ%data_c(:,i) = sin(i*x) + imagunit * cos(i*x)
      78             :          endif 
      79             :       enddo
      80             : 
      81           0 :    end subroutine set_sin_targ
      82             : 
      83           0 :    function proj_r(u,v) result(p)
      84             :       implicit none 
      85             :       real, intent(in) :: u(:), v(:)
      86             :       real :: p(size(u)) 
      87             : 
      88           0 :       p = dot_product(u,v) / dot_product(u,u) * u 
      89           0 :    end function
      90             : 
      91           0 :    function proj_c(u,v) result(p)
      92             :       implicit none 
      93             :       complex, intent(in) :: u(:), v(:)
      94             :       complex :: p(size(u)) 
      95             : 
      96           0 :       p = dot_product(u,v) / dot_product(u,u) * u 
      97           0 :    end function
      98             : 
      99           0 :    subroutine mod_gram_schmidt(M)
     100             :       implicit none 
     101             :       type(t_mat), intent(inout) :: M 
     102             :       
     103             :       integer :: k, i 
     104             : 
     105           0 :       do k = 1,M%matsize2 
     106           0 :          do i = 1,k-1 
     107           0 :             if(M%l_real) then 
     108           0 :                M%data_r(:,k) = M%data_r(:,k) - proj_r(M%data_r(:,i), M%data_r(:,k)) 
     109             :             else
     110           0 :                M%data_c(:,k) = M%data_c(:,k) - proj_c(M%data_c(:,i), M%data_c(:,k)) 
     111             :             endif 
     112             :          enddo 
     113             : 
     114             : 
     115           0 :          if(M%l_real) then
     116           0 :             M%data_r(:,k) = M%data_r(:,k) / norm2(M%data_r(:,k))
     117             :          else 
     118           0 :             M%data_c(:,k) = M%data_c(:,k) / norm2(abs(M%data_c(:,k)))
     119             :          endif
     120             :       enddo 
     121           0 :    end subroutine mod_gram_schmidt
     122             : 
     123           0 :    function linspace(beg, fin, n_steps) result(res)
     124             :       implicit none 
     125             :       real, intent(in)    :: beg, fin 
     126             :       integer, intent(in) :: n_steps
     127             :       real                :: res(n_steps), step
     128             : 
     129             :       integer :: i
     130             : 
     131           0 :       step = (fin - beg)  / (n_steps - 1.0)
     132             : 
     133           0 :       do i = 0,n_steps-1 
     134           0 :          res(i+1) = beg + step * i 
     135             :       enddo 
     136           0 :    end function linspace
     137             : 
     138           0 :    function make_groups(eigval) result(groups)
     139             :       implicit none
     140             :       real, intent(in)     :: eigval(:)
     141             :       integer, allocatable :: groups(:)
     142           0 :       integer              :: color(size(eigval))
     143             : 
     144             :       integer :: i, beg, g_cnt, n_groups
     145             : 
     146           0 :       g_cnt = 1
     147           0 :       beg  = 1 
     148             : 
     149           0 :       do while (beg <= size(eigval)) 
     150           0 :          i = beg
     151           0 :          do while (abs(eigval(i) - eigval(beg)) < 1e-8  ) 
     152           0 :             color(i) = g_cnt
     153           0 :             i = i + 1 
     154           0 :             if(i > size(eigval)) exit
     155             :          enddo 
     156           0 :          g_cnt = g_cnt + 1 
     157           0 :          beg = i 
     158             :       enddo
     159             : 
     160             : 
     161           0 :       n_groups = color(size(eigval))
     162           0 :       allocate(groups(n_groups))
     163             : 
     164           0 :       do i = 1,n_groups 
     165           0 :          groups(i) = count(i == color)
     166             :       enddo 
     167           0 :    end function make_groups
     168             : end module m_unify_zmat