4
$\begingroup$

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.

Improve this question
$\endgroup$
9
  • $\begingroup$ 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. $\endgroup$
    – Tyler Olsen
    Sep 24, 2015 at 11:59
  • 2
    $\begingroup$ 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. $\endgroup$
    – hardmath ♦
    Sep 24, 2015 at 13:57
  • 1
    $\begingroup$ 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. $\endgroup$
    – Edmondo Giovannozzi
    Sep 24, 2015 at 16:45
  • 2
    $\begingroup$ 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? $\endgroup$
    – Brian Borchers
    Sep 25, 2015 at 2:37
  • 1
    $\begingroup$ 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. $\endgroup$
    – Edmondo Giovannozzi
    Sep 25, 2015 at 9:15

2 Answers 2

Reset to default
5
$\begingroup$

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
Improve this answer
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – Lysistrata
    Apr 1, 2018 at 8:06
  • 5
    $\begingroup$ @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. $\endgroup$
    – Anthony Scemama
    Apr 8, 2018 at 8:49
3
$\begingroup$

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.

Improve this answer
$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – Edmondo Giovannozzi
    Sep 25, 2015 at 9:17
  • $\begingroup$ For matrices up to 14x14, is it still worth to code it like that? $\endgroup$
    – Hydro Guy
    Sep 28, 2015 at 1:52
  • $\begingroup$ 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). $\endgroup$
    – Anthony Scemama
    Apr 8, 2018 at 8:54
  • $\begingroup$ @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. $\endgroup$
    – Anthony Scemama
    Apr 8, 2018 at 8:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.