subroutine matrixsquare_up(m,n,A,B)
      implicit none

! compute B=AA^T. B is symmetrical so only its upper triangle is computed.
      integer m,n
      double precision A(m,n),B(m,m)
      
      integer i,j,k
      do i=1,m
        do j=i,m
          B(i,j)=0.0d0
          do k=1,n
            B(i,j)=B(i,j)+A(i,k)*A(j,k)
          enddo
        enddo
      enddo
      return
      end