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