program strassen
  ! Algoritmo ricorsivo di Strassen per il calcolo del prodotto
  ! tra matrici con un costo O(n^{log7/log2})
  implicit none
  real(8),dimension(:,:),allocatable :: a,b,c,mul
  integer :: n,k,null,fmm

  write(*,'(a)',advance='no') "Inserisci la potenza di 2 (max 8) "
  read(*,*) k
  n=2**k

  allocate(a(n,n),b(n,n),c(n,n),mul(n,n))
 
  call random_number(a)
  call random_number(b)

  write(*,'(a)') "Calcolo la moltiplicazione con Strassen"
  null=fmm(a,b,c,n)
  write(*,'(a)') "Calcolo la moltiplicazione con matmul"
  mul=matmul(a,b)

  write(*,*) "La differenza assoluta e` ",maxval(abs(c-mul))
end program strassen

recursive function fmm (a,b,c,n) result(null)
  implicit none
  
  integer :: n,null
  real(8),dimension(n,n) :: a,b,c
  real(8),dimension(n/2,n/2) :: a11,a12,a21,a22,b11,b12,b21,b22
  real(8),dimension(n/2,n/2) :: p1,p2,p3,p4,p5,p6,p7,c11,c12,c21,c22

  if (n==1) then
     c=a*b;null=1
  else
     a11=a(1:n/2,1:n/2);a12=a(1:n/2,n/2+1:n)
     a21=a(n/2+1:n,1:n);a22=a(n/2+1:n,n/2+1:n)
     b11=b(1:n/2,1:n/2);b12=b(1:n/2,n/2+1:n)
     b21=b(n/2+1:n,1:n);b22=b(n/2+1:n,n/2+1:n)
     
     null=fmm(a11+a22,b11+b22,p1,n/2)
     null=fmm(a21+a22,b11,p2,n/2)
     null=fmm(a11,b12-b22,p3,n/2)
     null=fmm(a22,b21-b11,p4,n/2)
     null=fmm(a11+a12,b22,p5,n/2)
     null=fmm(a21-a11,b11+b12,p6,n/2)
     null=fmm(a12-a22,b21+b22,p7,n/2)
     
     c11=p1+p4-p5+p7
     c21=p2+p4
     c12=p3+p5
     c22=p1+p3-p2+p6
     
     c(1:n/2,1:n/2)=c11;c(1:n/2,n/2+1:n)=c12
     c(n/2+1:n,1:n)=c21;c(n/2+1:n,n/2+1:n)=c22

  end if
end function fmm
