program main

  implicit none

  integer                 :: n,nc
  complex(8), allocatable :: a(:,:),atc(:,:),b(:,:)
  complex(8), allocatable :: x(:),y(:)
  logical                 :: self
  integer,    allocatable :: xx(:)

  integer                 :: i,j,k
  real*8                  :: re,im

  type new
      integer  :: inew
      real*8   :: rnew
  end type new        

  type(new), allocatable :: xnew(:)

  !1)
  do
    write(*,*) 'Insert positive integer'
    read(*,*) n
    if (n.gt.1) exit
  enddo

  write(*,*) 'Consider only n>1'

  allocate(a(n,n),atc(n,n))
  do j=1,n
    do i=1,n
       read(*,*) re, im
       a(i,j)=dcmplx(re,im)
    enddo
  enddo


  !2)
  call do_atc(a,n,atc,self)
  if (self) then
     write(*,*) 'Matrix is self-adjoint'
  else
     write(*,*) 'Matrix is not self-adjoint'
  endif   
  write(*,*) 'Matrix A and its transpose conjugated'
  do j=1,n
     do i=1,n
        write(*,*) i,j,a(i,j),atc(i,j)
     enddo
  enddo   

  !3)
  allocate(x(n))
  call do_x(atc,n,x)

  allocate(y(n))
  call do_ax(a,n,x,y)

  !4)
  allocate(b(n,n))
  call do_b(x,y,n,b)

  nc=0
  do i=1,n
     do j=1,n
        if (mod(i+j,2).eq.0) nc=nc+1
     enddo
  enddo   

  allocate(xx(nc))
  k=0
  do i=1,n
     do j=1,n
        if (mod(i+j,2).eq.0) then 
           k=k+1
           xx(k) = int(aimag(b(i,j))) + i + j
        endif   
     enddo
  enddo

  call sort_xx(xx,nc)

  !5)
  allocate(xnew(nc))

  do i=1,nc
     xnew(i)%inew=xx(i)
     if (xx(i).gt.0) then
        xnew(i)%rnew=dsqrt(dble(xx(i)))
     else
        xnew(i)%rnew=0.d0     
     endif   
  enddo

  !6)
  open(10,file='out1.dat')
  open(11,file='out2.dat')
  open(12,file='out3.dat')

  do i=1,n
     write(10,90) i, dble(y(i)),aimag(y(i))
  enddo

  do i=1,nc
     write(11,100) i,xx(i)
     write(12,110) i,xnew(i)%inew,xnew(i)%rnew
  enddo

  close(10)
  close(11)
  close(12)

  deallocate(a,atc,b)
  deallocate(x,y,xx)
  deallocate(xnew)

 90 format(I5,F10.3,F10.3) 
100 format(2(I5))  
110 format(I5,I5,F10.3)


  stop

end program main      



subroutine do_atc(a,n,atc,self)

  implicit none

  integer,   intent(in)  :: n
  complex(8),intent(in)  :: a(n,n)
  complex(8),intent(out) :: atc(n,n)
  logical,   intent(out) :: self

  integer                :: i,j

  atc=transpose(a)
  atc=conjg(atc)

  self=.false.
  do j=1,n
     do i=1,n
        if (a(j,i).eq.atc(j,i)) then
           self=.true.
        else
           self=.false.
           return 
        endif        
     enddo   
  enddo

  return

end subroutine do_atc        


subroutine do_x(atc,n,x)

  implicit none

  integer,   intent(in)  :: n
  complex(8),intent(in)  :: atc(n,n)
  complex(8),intent(out) :: x(n)

  integer                :: i

  x=dcmplx(0.d0,0.d0)
  do i=1,n,2
     x=x+atc(:,i)
  enddo

  return

end subroutine do_x        


subroutine do_ax(a,n,x,y)

  implicit none

  integer,   intent(in)  :: n
  complex(8),intent(in)  :: a(n,n),x(n)
  complex(8),intent(out) :: y(n)

  integer                :: i,j

  y=dcmplx(0.d0,0.d0)
  do i=1,n
     do j=1,n
        y(i) = y(i) + a(i,j)*x(j)
     enddo
  enddo   

  return

end subroutine do_ax

subroutine do_b(x,y,n,b)

  implicit none

  integer,   intent(in)  :: n
  complex(8),intent(in)  :: x(n),y(n)
  complex(8),intent(out) :: b(n,n)


  integer                :: i 

  do i=1,n,2
     b(:,i)=y
  enddo

  do i=2,n,2
     b(:,i)=conjg(x)
  enddo

  return

end subroutine do_b        


subroutine sort_xx(xx,nc)

  implicit none

  integer,   intent(in)    :: nc
  integer,   intent(inout) :: xx(nc)
   
  integer                  :: i,j
  integer                  :: tmp

  do i=1,nc-1
    do j=2,nc-(i-1)
       if (xx(j-1).lt.xx(j)) then
          tmp=xx(j-1)
          xx(j-1)=xx(j)
          xx(j)=tmp    
       endif
    enddo
 enddo

 return

end subroutine sort_xx