program main

  implicit none

  integer :: n,m,s,d,i,j,k
  real*8, allocatable :: a(:), q(:,:), l(:,:), p(:,:)

  !1)
  open(10,file='input')
  read(10,*) n,m
  close(10)

  if (n.eq.0.and.m.eq.0) then
     write(*,*) 'n and m are both zero, stop here'
     stop
  endif

  !2)
  if (n.gt.m) then
     s=n
  else
     s=m
  endif   
  write(*,*) 'Number of elements of a', s

  allocate(a(s))
  d=0
  do i=1,s
     a(i)=n+m+sqrt(dble(i))
     write(*,*) i, a(i)
     if (n+m.gt.sqrt(dble(i))) then
        d=d+1
     endif   
  enddo   
  write(*,*) 'Number of elements satisying (n+m)>sqrt(i)', d

  !3)
  if (d.ge.2) then
     allocate(q(d,d))
     allocate(l(d,d))
     allocate(p(d,d))
     j=0
     do i=1,s
        if (n+m.gt.sqrt(dble(i))) then
           j=j+1
           q(j,j)=a(i) 
           l(j,j)=a(i)
        endif  
     enddo   
     do j=1,d
        do k=j+1,d
           q(k,j)=k+j
           q(j,k)=q(k,j)
           l(k,j)=k-j
           l(j,k)=-l(k,j)
        enddo 
     enddo   
     call do_mat(q,l,d,p)
  else
     write(*,*) 'Matrix instructions are not done'
  endif

  !4)
  open(10,file='a.dat')
  do i=1,s
     write(10,100), i, a(i)
  enddo   
  close(10)

  if (d.ge.2) then
     open(11,file='q.dat')
     open(12,file='l.dat')
     open(13,file='p.dat')
     do i=1,d
        do j=1,d
           write(11,101) i,j,q(i,j)
           write(12,101) i,j,l(i,j)
           write(13,101) i,j,p(i,j)
        enddo
     enddo   
     close(11)
     close(12)
     close(13)
  endif 

100 format(I4,F10.4)
101 format(I4,I4,F10.4)

  stop

end program main  


subroutine do_mat(q,l,d,p)

 implicit none

 integer, intent(in) :: d
 real*8, intent(in)  :: q(d,d), l(d,d)
 real*8, intent(out) :: p(d,d)

 integer i,j,k

 p=0.d0
 do i=1,d
    do j=1,d
       do k=1,d
          p(i,j)=p(i,j)+q(i,k)*l(k,j)
       enddo   
    enddo
 enddo   

 return

end subroutine do_mat