!===============================================================================
! Program to demonstrate the use of MKL to invert a small matrix.
!
!  Link requirements:
!
!     Compiler: Qmkl 
!                 (in VS, choose Project > Properties > Fortran > Libraries
!                       > Use Intel MKL: Qmkl:Parallel) (parallel == generic)
!                 (in Command line: /Qmkl)
!
!     Linker:   mkl_Lapack95.lib
!                 (in VS, set as Linker > Input > Addl dependencies)
!                 (cmd line: mkl_lapack95.lib)
!     To run exe files outside the VS environment, must also include
!     libiomp5md.dll on the link line, possibly including its path. Or, locate
!     this file and put a copy in the Release and Debug folders.
!     c:\Program Files (x86)\Common Files\Intel\
!        Shared Libraries\redist\ia32\compiler
!
!  Specifying the path to these libraries appears to be unnecessary. Don't know
!  why.
!-------------------------------------------------------------------------------
   PROGRAM MKLinvert
      !include 'lapack.f90'
      !USE IFQWIN
      USE LAPACK95 !, only: getrf, getri, potrf, potri
      !use mkl95_lapack ! appears in Intel mkl samples; does same thing
      USE F95_PRECISION, ONLY: WP => SP ! or DP
      IMPLICIT NONE
      INTEGER ipiv(2) ! NOTE must have length >= A,B,C
      !REAL(WP)  :: A(2, 2), B(2, 2), C(2, 2)
      real(4) a(2,2), b(2,2), c(2,2) ! same result
     ! potrf,potri?
      character(1) uplo
      integer info
      
!---- Define an input matrix.      
      
      A(1,1) = 10.
      A(1,2) = 3.
      A(2,1) = A(1,2)
      A(2,2) = 5.
      print *, 'Input matrix:'
      print *, 'A0 =          ', A ! OK
      print *
      
!---- Call MKL functions to factor, then inverse.      

      ! Can invert in place with following calls,
      !CALL getrf (A, ipiv, info)
      !CALL getri (A, ipiv)
      ! But we will work on a copy to verify the result.
      ! Also demonstrate absence of optional parameters.
      !B = A_inv
      
      B = A
      print *, 'A1 =          ', A
      CALL getrf (B, ipiv)
      print *, 'A2 =          ', A
      CALL getri (B, ipiv)
      print *, 'A3 =          ', A 
      ! Note that ipiv is optional in getrf according to the documentation, but
      ! the result is wrong unless it is included. In getri it is (correctly)
      ! shown to be necessary. IPIV IS THE PIVOT INFORMATION, WHICH COMES FROM
      ! GETRF AND IS NEEDED BY GETRI. So it is optional when calling getrf only
      ! to factor the matrix, but is required when followed by getri to do the
      ! inversion.
      
!---- Show the result and check it.      

      PRINT *
      PRINT *, 'Matrix inversion using MKL functions getrf, getri'
      PRINT *, 'A =        ', A
      PRINT *, 'B = A_inv =', B
      PRINT *, 'A * B =    ', MATMUL (A, B)
      PRINT *, 'Identiy matrix result indicates inversion was correct.'
      
!---- Repeat the problem using routines potrf and potri. These are for a
!     positive-definite symmetric matrix. According to the documentation the no.
!     of floating-point operations for potri is about half of getri.

      C = A
      CALL potrf (C) ! Upper tri part of C is returned as Cholesky factor U
      CALL potri (C)
      !CALL potri (C, 'U', info) ! No help. info returns 0 (ok).
      
      PRINT *
      PRINT *, 'Result using potrf, potri'
      PRINT *, 'A =        ', A
      PRINT *, 'C = A_inv =', C
      PRINT *, 'A * C =    ', MATMUL (A, C)
      print *, 'Result using potrf, potri may be wrong. Don''t know why.'
      

   END PROGRAM MKLinvert