Efficently invert tiny matrix in Fortran

I have a piece of code in Fortran90 in which I have to solve both a non-linear (with the Newton-Raphson method, for which I have to invert the Jacobian matrix) and a linear system of equations. When I say very small I mean n unknowns for both operations, with $n\leq4$. Unfortunately, $n$ is not known a priori. What do you think is the fastest option? I thought of writing explicit formulas for cases with $n=1,2$ and using other methods for $n=3,4$ (e.g. some functions of the Intel MKL libraries), for the sake of performance. Is this sensible or should I write explicit formulas for the inverse matrix also for $n=3,4$?

The code is going to be called very many times in a Finite Element Method analysis, for this reason I was looking for the fastest solution.

• The explicit inverses for matrices up to 3x3 are very easy to write, so you should at least write those. As for 4x4, I don't use them much, so I just drop to an LU decomp to invert. It's still fast enough for my needs, but you may find that you require more speed, in which case writing the explicit inverse may be worth the (1-time) effort. Sep 24 '15 at 11:59
• One should also question the need to invert a matrix, rather than to solve linear systems. Unless the same system has to be solved a large number of times, solving a system by elimination/factorization is more efficient than computing a matrix inverse and using it to multiply the right hand side. Sep 24 '15 at 13:57
• I made a test, for what it is worth, between a hand calculated inverse of a 4x4 matrix (I used Maxima to symbolically calculate it, of course) and a standard lapack DGESV. The solution of the system with the hand calculated inverse was from five to ten times faster then the Lapack one. Of course I didn't check if the matrix was invertible, etc. To take with extreme care. Sep 24 '15 at 16:45
• Have you convinced yourself that your code is currently spending most of its time solving these systems of equations when you simply call the appropriate LAPACK routine? If not, then there's no point in optimizing this because you won't get a significant speedup. What portion of your run time is currently taken up by this? Sep 25 '15 at 2:37
• Few words about that simple test. I should say that in a real program, I would use lapack (I will be worried by numerical stability, etc.). I made tests in single and double precision, with ifort using standard -O2 optimization as well with -fast optimization (which includes ipa). The difference between single and double precision gives me an hint on the effect of cache miss. The problem is that I'm not controlling the Lapack, while I could compile that hand made routine together with the main program and the compiler could do all the optimizations it liked. But I agree with Brian. Sep 25 '15 at 9:15

Here are the routines we use in QMC=Chem for small matrices. QMC=Chem is under GPL license, so if your code is also under GPL, you can use them.

EDIT: To get the inverse, you need to multiply the result by 1.d0/det_l.

subroutine invert2(a,LDA,na,det_l)
implicit none
double precision :: a (LDA,na)
integer          :: LDA
integer          :: na
double precision :: det_l
double precision :: b(2,2)

double precision :: f
b(1,1) = a(1,1)
b(2,1) = a(2,1)
b(1,2) = a(1,2)
b(2,2) = a(2,2)
det_l = a(1,1)*a(2,2) - a(1,2)*a(2,1)
a(1,1) = b(2,2)
a(2,1) = -b(2,1)
a(1,2) = -b(1,2)
a(2,2) = b(1,1)
end

subroutine invert3(a,LDA,na,det_l)
implicit none
double precision, intent(inout) :: a (LDA,na)
integer, intent(in)             :: LDA
integer, intent(in)             :: na
double precision, intent(inout) :: det_l
double precision :: b(4,3)
integer :: i
double precision :: f
det_l = a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2)) &
-a(1,2)*(a(2,1)*a(3,3)-a(2,3)*a(3,1)) &
+a(1,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1))
do i=1,4
b(i,1) = a(i,1)
b(i,2) = a(i,2)
b(i,3) = a(i,3)
enddo
a(1,1) =  b(2,2)*b(3,3) - b(2,3)*b(3,2)
a(2,1) =  b(2,3)*b(3,1) - b(2,1)*b(3,3)
a(3,1) =  b(2,1)*b(3,2) - b(2,2)*b(3,1)

a(1,2) =  b(1,3)*b(3,2) - b(1,2)*b(3,3)
a(2,2) =  b(1,1)*b(3,3) - b(1,3)*b(3,1)
a(3,2) =  b(1,2)*b(3,1) - b(1,1)*b(3,2)

a(1,3) =  b(1,2)*b(2,3) - b(1,3)*b(2,2)
a(2,3) =  b(1,3)*b(2,1) - b(1,1)*b(2,3)
a(3,3) =  b(1,1)*b(2,2) - b(1,2)*b(2,1)

end

subroutine invert4(a,LDA,na,det_l)
implicit none
double precision, intent(inout) :: a (LDA,na)
integer, intent(in)             :: LDA
integer, intent(in)             :: na
double precision, intent(inout) :: det_l
double precision :: b(4,4)
integer :: i,j
double precision :: f
det_l =  a(1,1)*(a(2,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))  &
-a(2,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))  &
+a(2,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))) &
-a(1,2)*(a(2,1)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))  &
-a(2,3)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))  &
+a(2,4)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))) &
+a(1,3)*(a(2,1)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))  &
-a(2,2)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))  &
+a(2,4)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))) &
-a(1,4)*(a(2,1)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))  &
-a(2,2)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))  &
+a(2,3)*(a(3,1)*a(4,2)-a(3,2)*a(4,1)))
do i=1,4
b(1,i) = a(1,i)
b(2,i) = a(2,i)
b(3,i) = a(3,i)
b(4,i) = a(4,i)
enddo

a(1,1) =  b(2,2)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(2,3)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))+b(2,4)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))
a(2,1) = -b(2,1)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))+b(2,3)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))-b(2,4)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))
a(3,1) =  b(2,1)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))-b(2,2)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(2,4)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))
a(4,1) = -b(2,1)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))+b(2,2)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))-b(2,3)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))

a(1,2) = -b(1,2)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))+b(1,3)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))-b(1,4)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))
a(2,2) =  b(1,1)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(1,3)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(1,4)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))
a(3,2) = -b(1,1)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))+b(1,2)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))-b(1,4)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))
a(4,2) =  b(1,1)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))-b(1,2)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))+b(1,3)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))

a(1,3) =  b(1,2)*(b(2,3)*b(4,4)-b(2,4)*b(4,3))-b(1,3)*(b(2,2)*b(4,4)-b(2,4)*b(4,2))+b(1,4)*(b(2,2)*b(4,3)-b(2,3)*b(4,2))
a(2,3) = -b(1,1)*(b(2,3)*b(4,4)-b(2,4)*b(4,3))+b(1,3)*(b(2,1)*b(4,4)-b(2,4)*b(4,1))-b(1,4)*(b(2,1)*b(4,3)-b(2,3)*b(4,1))
a(3,3) =  b(1,1)*(b(2,2)*b(4,4)-b(2,4)*b(4,2))-b(1,2)*(b(2,1)*b(4,4)-b(2,4)*b(4,1))+b(1,4)*(b(2,1)*b(4,2)-b(2,2)*b(4,1))
a(4,3) = -b(1,1)*(b(2,2)*b(4,3)-b(2,3)*b(4,2))+b(1,2)*(b(2,1)*b(4,3)-b(2,3)*b(4,1))-b(1,3)*(b(2,1)*b(4,2)-b(2,2)*b(4,1))

a(1,4) = -b(1,2)*(b(2,3)*b(3,4)-b(2,4)*b(3,3))+b(1,3)*(b(2,2)*b(3,4)-b(2,4)*b(3,2))-b(1,4)*(b(2,2)*b(3,3)-b(2,3)*b(3,2))
a(2,4) =  b(1,1)*(b(2,3)*b(3,4)-b(2,4)*b(3,3))-b(1,3)*(b(2,1)*b(3,4)-b(2,4)*b(3,1))+b(1,4)*(b(2,1)*b(3,3)-b(2,3)*b(3,1))
a(3,4) = -b(1,1)*(b(2,2)*b(3,4)-b(2,4)*b(3,2))+b(1,2)*(b(2,1)*b(3,4)-b(2,4)*b(3,1))-b(1,4)*(b(2,1)*b(3,2)-b(2,2)*b(3,1))
a(4,4) =  b(1,1)*(b(2,2)*b(3,3)-b(2,3)*b(3,2))-b(1,2)*(b(2,1)*b(3,3)-b(2,3)*b(3,1))+b(1,3)*(b(2,1)*b(3,2)-b(2,2)*b(3,1))

end

• Downvoted by me. The program does not work and it should have been tested before posting here. The line "do i=1,4" in subroutine Invert3 is clearly wrong. Apr 1 '18 at 8:06
• @Lysistrata : No, the do i=1,4 is not wrong. It is to enforce automatic vectorization with AVX. The dimension of b is (4,3). This code is correct, but I just forgot to mention you need to divide by det_l to get the inverse, my bad. Apr 8 '18 at 8:49

For such small matrices I simply call LAPACK ...

! Computes inverse of a matrix
subroutine MatInv(A,Ainv)
implicit none
real(8) :: A(:,:),Ainv(:,:)
real(8),allocatable :: work(:)
integer :: n,info
integer,allocatable :: ipiv(:)
Ainv=A
n=size(Ainv,1)
allocate(ipiv(n),work(n))
call dgetrf(n,n,Ainv,n,ipiv,info)
call dgetri(n,Ainv,n,ipiv,work,n,info)
deallocate(ipiv,work)
end subroutine MatInv


Not worth it to roll your own stuff because of potential bugs and code bloat.

• For a small matrix I wouldn't use an allocatable variable. Declaring work, and ipiv as real(wp) :: work(size(A,1)) and integer :: ipiv(size(A,1)) will be better. Sep 25 '15 at 9:17
• For matrices up to 14x14, is it still worth to code it like that? Sep 28 '15 at 1:52
• I would not call Lapack for small matrices because the overhead of the function call can be of the order of the computation time. And for matrices less by 5x5, the O(N!) algorithm is faster than the O(N^3). Apr 8 '18 at 8:54
• @HydroGuy : It may be preferable to declare ipiv(64) and work(64) hard-coded with a static max size for small problems (<64x64), and use the allocate for larger problems. See my edit. Apr 8 '18 at 8:58