subroutine calculate_MATMUL                 &
  (NTOTAL,LDZ,W_single,Z_single,d_out)

!------------------------------
!  Eigenvalue decomposition
!
!  Employs the ESSL programs DSPEV  (double precision), SGEMUL (single precision)
!  or
!  Employs the LAPACK programs DSPEV  (double precision), SGEMM (single precision)
!-----------------------------------------------------------------------------
      integer :: NTOTAL,LDZ
      real W_single(1:NTOTAL)                  !! INPUT: Eigenvalues
      real Z_single(1:LDZ,1:NTOTAL)            !! INPUT: Eigenvectors
      real d_out(1:LDZ,1:NTOTAL)               !! OUTPUT: (I+A)**(-1/2)
!-----------------------------------------------------------------------------
      real*4,dimension(:,:),allocatable :: b,v,d
                                                                                                     
      EXTERNAL SGEMUL

!==============start calculation===================

      allocate(v(1:LDZ,1:NTOTAL))
      allocate(d(1:LDZ,1:NTOTAL))
      allocate(b(1:LDZ,1:NTOTAL))

      v(:,:)=Z_single(:,:)
      d(:,:)=0.0
      b(:,:)=0.0

      do j=1,NTOTAL
       d(j,j)=W_single(j)
        end do

!!--- LAPACK --------------------------------
!!    call SGEMM ('N','T',NTOTAL,NTOTAL,NTOTAL,1.0,d,NTOTAL,v,NTOTAL,0.0,b,NTOTAL)
!!--- ESSL --------------------------------
      call SGEMUL (d,NTOTAL,'N',v,NTOTAL,'T',b,NTOTAL,NTOTAL,NTOTAL,NTOTAL)
!!-----------------------------------------

      d(:,:)=0.0 

!!--- LAPACK --------------------------------
!!    call SGEMM ('N','N',NTOTAL,NTOTAL,NTOTAL,1.0,v,NTOTAL,b,NTOTAL,0.0,d,NTOTAL)
!!--- ESSL --------------------------------
      call SGEMUL (v,NTOTAL,'N',b,NTOTAL,'N',d,NTOTAL,NTOTAL,NTOTAL,NTOTAL)
!!-----------------------------------------

      d_out(:,:)=d(:,:)

      deallocate(d)
      deallocate(v)
      deallocate(b)

           END subroutine calculate_MATMUL