program test_matrix_multiplication
    
    ! Matlab code to test such a situation in Matlab
    ! a=rand(3600);
    ! tic
    ! b=a\speye(size(a,1));
    ! toc
    ! define parameters
    integer, parameter:: rkind=kind(1.D0)! precision of real 
    integer, parameter:: ikind=4 !precision of integers
    ! arrays that filled with random numbers for inverting matrices
    real(kind=rkind), allocatable, dimension(:,:):: a
    real(kind=rkind), allocatable, dimension(:,:):: b
    !dec$ attributes align: 32 :: a , b
    integer(kind=ikind),allocatable,dimension(:)::ipiv
    real(kind=rkind),allocatable,dimension(:)::work
    integer(kind=ikind)::a_rows ! first dimension of matrix
    integer(kind=ikind)::a_cols ! second dimension of matrix
    integer(kind=ikind)::b_rows ! first dimension of matrix
    integer(kind=ikind)::b_cols ! second dimension of matrix
    integer(kind=ikind)::nmax ! workspace query
    integer(kind=ikind)::seed
    ! for timing of the routines
    integer(kind=ikind):: clock_start 
    integer(kind=ikind):: clock_end
    integer(kind=ikind):: clock_rate
    character(len=3)::print_condition='yes' ! print only for small matrices
    
    !*******************fixed random seed*********************
    seed=4
    a_rows=3 ! change these numbers for larger matrices
    a_cols=3
    !***************Allocate Space on Heap*********************
    allocate(a(a_rows,a_cols))
    allocate(b(a_rows,a_cols))
    !**************Fill matrices with random numbers***********
    call random_number(a)
    !**************Create Positive Definite Matrix***********
    b=mat_mult(transpose(a),a) ! create a positive definite matrix to be inverted

    !**************Invert the matrix using dgetrf and dgetri****
    allocate(ipiv(a_rows),work(1))
    
    !**************start timing Using dgetrf/dgetri*************
    
    call system_clock(count_rate=clock_rate) 
    call system_clock(count=clock_start)
    
    ! do factorization
 
    call dgetrf(a_rows,a_rows,b,a_rows,ipiv,info)
   
    ! workspace query
    call dgetri(a_rows,b,a_rows,ipiv,work,-1,info)
    nmax=work(1)
   
    deallocate(work)
    allocate(work(nmax))
    ! invert matrix a
    
    call dgetri(a_rows,b,a_rows,ipiv,work,nmax,info)
    
    call system_clock(count=clock_end)
    write(*,"(a, f16.8)")'Total Time Taken For Inverting from dgetri Matrix = ',&
    real(clock_end-clock_start)/real(clock_rate), ' seconds'
    
    deallocate(ipiv,work)
    
    if (print_condition =='yes') then
        
        !****************print original matrix************
        call print_mat(a,a_rows,a_cols,'original_matrix')
        
        !****************print difference************
        call print_mat(b,a_rows,a_cols,'inverted_matrix')
    
    end if
    
    deallocate(a,b)
    
    contains

    function mat_mult(a,b) result (c)
    ! a wrapper function to multiply two matrices using mkl
    implicit none
    integer(kind=ikind):: a_rows ! no of rows in matrix a
    integer(kind=ikind):: a_cols ! no of cols in matrix a
    integer(kind=ikind):: b_cols ! no of cols in matrix b
    real(kind=rkind),intent(in),dimension(:,:) :: a, b
    real(kind=rkind),allocatable,dimension(:,:):: c
    real(kind=rkind):: alpha, beta ! two parameters required by dgemm subroutine 
    ! in atlas or mkl for multiplying matrices
    a_rows=size(a,1)
    a_cols=size(a,2)
    b_cols=size(b,2)
    alpha=1.0
    beta=0.0
    allocate(c(a_rows,b_cols))
    call dgemm('n', 'n', a_rows, b_cols, a_cols, alpha, a, a_rows, b, a_cols, beta, c, a_rows)
    end function mat_mult

    subroutine print_mat(matrix,rows,columns,name_matrix)
    ! function to print matrices
    ! variable assignment
    character(len=10):: col
    character(len=20):: formt
    character(len=*), intent(in):: name_matrix
    integer(kind=ikind):: n
    integer(kind=ikind), intent(in):: rows
    integer(kind=ikind), intent(in):: columns
    real(kind=rkind), intent(in), dimension(:,:):: matrix
    ! execution part
    write(col,"(i0)") columns
    col = trim(col)
    formt = '('//col//'f15.6))'
    write(*,"(a)"),' '
    write(*,"(a,a)"),'now printing: ', name_matrix
    write(*,"(a)"),' '
    do n = 1, rows
        write(*,formt) matrix(n,:)
    end do
    write(*,*),' '
    end subroutine print_mat

    end program test_matrix_multiplication