program test_matrix_inversion
    
    ! 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
    real(kind=rkind), allocatable, dimension(:,:):: c
    real(kind=rkind), allocatable, dimension(:,:):: dif
    !dec$ attributes align: 32 :: a , b, c, dif
    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
    ! input directory and filename information
    character(len=400)::inp_direc='E:\test_data_fortran_inversion\'
    character(len=400)::file_name='original_matrix.bin'
    character(len=400)::file_name2='inverted_matrix.bin'
    character(len=400)::file_name3='output_matrix'
    ! define variables for dimensions of the matrices
    character(len=3)::print_condition='yes'
 
    !*******************fixed random seed*********************
    seed=4
    !***************Enter dimensions of the matrix***********
    a_rows=1000
    b_rows=1000
    a_cols=1000
    b_cols=1000
    !**************allocate space on heap**********************
    allocate(a(a_rows,a_cols))
    allocate(b(b_rows,b_cols))
    allocate(c(a_rows,b_cols))
    allocate(dif(a_rows,b_cols))
    !**************Fill matrices with random numbers***********
    
    ! input_matrix
    a = read_binary_matrix(inp_direc,file_name,1000,1000)
    ! inverted_matrix_from_matlab
    b = read_binary_matrix(inp_direc,file_name2,1000,1000)
    c=a ! keep a copy of original matrix
    d=0.0d0 ! space for difference between two inverted matrices
    
    !**************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)
 
    call dgetrf(a_rows,a_rows,a,a_rows,ipiv,info)
    ! workspace query
    call dgetri(a_rows,a,a_rows,ipiv,work,-1,info)
    nmax=work(1)
    deallocate(work)
    allocate(work(nmax))
    ! invert matrix a
    call dgetri(a_rows,a,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'
    
    !**************calculate difference*************
    dif=a-b
    deallocate(ipiv,work)
    
    if (print_condition =='yes') then
        !****************print difference************
        call print_mat(dif(1:3,1:3),3,3,'difference between two inverted matrices')
    end if

    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
    
    function read_binary_matrix(inp_direc,file_name,matrix_rows,matrix_cols) result(matrix) 
    !**************************************************************************
    !purpose : read array written in binary format written in dimension rkind
    !***********specification part*********************************************
    implicit none
    character (len=*) :: inp_direc ! input directory from which to read file
    character (len=*) :: file_name ! name of the file to be read
    character (len=400) :: fname
    integer(kind=ikind),intent(in):: matrix_rows !number of cols in the file for sparse matrix
    integer(kind=ikind),intent(in):: matrix_cols
    integer(kind=ikind),parameter :: unitcons=200
    real(kind=rkind),dimension(:,:)::matrix(matrix_rows,matrix_cols)
    !**********execution part**************************************************
    fname=trim(adjustl(inp_direc))//trim(adjustl(file_name))
    open (unitcons,file=fname,form='unformatted',access='stream')
    read(unitcons) matrix
    close(unitcons)
    end function read_binary_matrix
    
    subroutine write_binary_matrix(inp_direc,file_name,matrix) 
    !**************************************************************************
    !written by : vineet yadav
    !date : 8/29/2010
    !purpose : read array writtent in binary format written in dimension rkind
    !***********specification part*********************************************
    implicit none
    character (len=*) :: inp_direc ! input directory from which to read file
    character (len=*) :: file_name ! name of the file to be read
    character (len=400) :: name
    integer (kind=ikind),parameter :: unitcons=200
    real (kind=rkind),intent(in), dimension(:,:)::matrix
    !**********execution part**************************************************
    name=trim(adjustl(inp_direc))//trim(adjustl(file_name))!file name with directory
    open (unitcons,file=name,form='unformatted',access='stream')
    write(unitcons) matrix
    close(unitcons)
    end subroutine write_binary_matrix


    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//'f25.16))'
    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_inversion