! ---------------------------------------------------
! Example program
! Matrix manipulation: intrinsic routines
!
! Notice how the Fortran internal order of elements
! is the transpose of their usual mathematical order! 
! ---------------------------------------------------
program test
  implicit none
  !
  ! interface needed for the assumed-shape 
  ! matrix output routine 
  interface
     subroutine writematrix(A)
       real(kind(1.d0)),intent(in):: A(:,:)
     end subroutine writematrix
  end interface

  integer, parameter:: dbl=kind(1.d0),n=3,m=5
  integer:: i,j,k
  real(dbl):: A(n,n),B(n,n),C(n,m),v(n*n)

  k = 0
  do i = 1, n
     do j = 1, m
        C(i,j) = dble(i+j)+1.d-15 ! notice the subtraction, see (*)
        if(j>n) cycle
        A(i,j) = 10*i+j
        B(i,j) = dble(i/j)
        k = k + 1
        v(k) = dble(k)
     end do
  end do
  print*,'A'
  write(*,99) A
  print*,'B'
  write(*,99) B
  print*,'transpose(A)'
  write(*,99) transpose(A)
  !
  ! matrix product
  ! --------------
  print*,'matmul(A,B)'
  write(*,99) matmul(A,B)
  print*,'transpose(matmul(A,B))'
  write(*,99) transpose(matmul(A,B))
  print*,'manually computed product'
  print*,'-------------------------'
  write(*,'(2a3,a10)') 'i','j','AB(i,j)'
  !
  do i = 1, n
     do j = 1, n
        !write(*,'(2i3,f10.2)') i,j,sum(A(i,:)*B(:,j))
        ! same thing using intrinsic dot_product
        write(*,'(2i3,f10.2)') i,j,dot_product(A(i,:),B(:,j))
     end do
  end do

  print*,'C'
  call writematrix(C)


  print*
  print*
  print*,'------------------------------------------------------------'
  ! (*) WARNING: comparison of real numbers 
  print*,'comparison of real numbers may cause unecpected results'
  print*,'C(2,1)=',C(2,1)
  print*,'output of C>3.d0 logical test:'
  write(*,'(3(5l3,/))') transpose(C>3.d0)
  print*,'because there was a tiny mistake of the size 1.d-15'
  print*,'------------------------------------------------------------'
  print*
  print*
  
  print*,'------------------------------------------------------------'
  print*,'Some vector manipulations:'
  print*,'v'
  write(*,'(10f7.2)') v
  print*,'maxloc(v)'
  write(*,'(10i2)') maxloc(v)
  print*,"why setting i = maxloc(v) won't work :"
  print*,'shape(i) = ',shape(i)
  print*,'shape(maxloc(v)) = ',shape(maxloc(v))
  print*,'circular shift left: cshift(v,2)'
  write(*,'(10f7.2)') cshift(v,2) 
  print*,'shift right with zero boundary, eoshift(v,-1)'
  write(*,'(10f7.2)') eoshift(v,-1)
  print*,'shift right two steps with boundary 5.d0, eoshift(v,-2,5.d0)'
  write(*,'(10f7.2)') eoshift(v,-2,5.d0)
  A= reshape(v,(/n,n/))
  print*,'reshape(v,(/3,3/))'
  write(*,99) A
  print*,'------------------------------------------------------------'



98 format(5(3f7.2,/))
99 format(3(3f7.2,/))
end program test

! ------------------------
! Matrix output
! of sssumed-shape array A
! ------------------------
subroutine writematrix(A)
  implicit none
  integer::i,j,n,m
  real(kind(1.d0)),intent(in):: A(:,:)
  n = ubound(A,1)
  m = ubound(A,2)
  print*,'writematrix: output in mathematical order'
  do i = 1, n
     write(*,'(10f7.2)') (A(i,j) ,j=1,m)
  end do
end subroutine writematrix