Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
The Intel sign-in experience has changed to support enhanced security controls. If you sign in, click here for more information.
6743 Discussions

Reduced Row Echelon Form of a matrix (rref)



I'm newbie to the Intel  MKL library and I'm trying to convert a code from Matlab to C using the C interface of Intel MKL routines.

I haven't found a function to achieve the RREF of a matrix in Intel MKL. In fact, the RREF in Matlab performs gaussian elimination with partial pivoting and I want to apply it to a 5x9 matrix.

Here is an example of RREF applied to a 5x9 matrix A

    0.8147    0.0975    0.1576    0.1419    0.6557    0.7577    0.7060    0.8235    0.4387
    0.9058    0.2785    0.9706    0.4218    0.0357    0.7431    0.0318    0.6948    0.3816
    0.1270    0.5469    0.9572    0.9157    0.8491    0.3922    0.2769    0.3171    0.7655
    0.9134    0.9575    0.4854    0.7922    0.9340    0.6555    0.0462    0.9502    0.7952
    0.6324    0.9649    0.8003    0.9595    0.6787    0.1712    0.0971    0.0344    0.1869

The RREF of A is composed of an 5x5 identity matrix and a 5x4 matrix

    1.0000         0         0         0         0   -0.8409    1.9764   -2.4048   -2.7594
         0    1.0000         0         0         0    3.7698   -5.9264    8.8249    8.0095
         0         0    1.0000         0         0    3.8571   -4.0608    7.6411    7.0677
         0         0         0    1.0000         0   -8.0047    9.2175  -17.0036  -15.1016
         0         0         0         0    1.0000    2.4444   -1.5157    4.7735    4.4750

I tried to apply the LU factorisation with dgetrf and also the LAPACKE_dgels function but I dindn't get the wanted result.

Thank you for your help


0 Kudos
2 Replies


Assume the matrix under consideration is represented as a compound matrix A=(A1|A2). What you want can be described by the following formula A1^{-1}*(A1|A2). See the code in Fortran below. Please don't be surprised with a small difference appered due to low accuracy of data you provided (only 4 digits).

      program test
      implicit none
      integer m,n,i,j
      parameter (m=5,n=9)
      double precision a(m,n),b(n,m),c(m,m), work(m), d(m,n),alpha,beta
      integer info, ipiv(m)
      data alpha/1D0/, beta/0D0/
      data b/ 0.8147D0, 0.0975D0, 0.1576D0, 0.1419D0, 0.6557D0, 0.7577D0, 0.7060D0, 0.8235D0, 0.4387D0,
     &        0.9058D0, 0.2785D0, 0.9706D0, 0.4218D0, 0.0357D0, 0.7431D0, 0.0318D0, 0.6948D0, 0.3816D0,
     &        0.1270D0, 0.5469D0, 0.9572D0, 0.9157D0, 0.8491D0, 0.3922D0, 0.2769D0, 0.3171D0, 0.7655D0,
     &        0.9134D0, 0.9575D0, 0.4854D0, 0.7922D0, 0.9340D0, 0.6555D0, 0.0462D0, 0.9502D0, 0.7952D0,
     &        0.6324D0, 0.9649D0, 0.8003D0, 0.9595D0, 0.6787D0, 0.1712D0, 0.0971D0, 0.0344D0, 0.1869D0/

      do i=1,m
          do j=1,n
          end do
      end do
      do i=1,m
          do j=1,m
          end do
      end do

      print *,"initial matrix (A1|A2)"
      do i=1,m
          print 1,(a(i,j), j=1,n)
      end do    
      call dgetrf(m,m,c,m,ipiv,info)
      if( then
          print *,"dgetrf returned info=",info
      end if
      call dgetri( m, c, m, ipiv, work, m, info )
      if( then
          print *,"dgetri returned info=",info
      end if
      call dgemm('N', 'N', m, n, m, alpha, c, m, a, m, beta, d, m)
      print *,"result A^{-1}*(A1|A2)"
      do i=1,m
          print 1,(d(i,j), j=1,n)
      end do    
   1  format(1x, 9(1x, D11.4)) 


The output of the code is as follows:

  initial matrix (A1|A2)
   0.8147D+00  0.9750D-01  0.1576D+00  0.1419D+00  0.6557D+00  0.7577D+00  0.7060D+00  0.8235D+00  0.4387D+00
   0.9058D+00  0.2785D+00  0.9706D+00  0.4218D+00  0.3570D-01  0.7431D+00  0.3180D-01  0.6948D+00  0.3816D+00
   0.1270D+00  0.5469D+00  0.9572D+00  0.9157D+00  0.8491D+00  0.3922D+00  0.2769D+00  0.3171D+00  0.7655D+00
   0.9134D+00  0.9575D+00  0.4854D+00  0.7922D+00  0.9340D+00  0.6555D+00  0.4620D-01  0.9502D+00  0.7952D+00
   0.6324D+00  0.9649D+00  0.8003D+00  0.9595D+00  0.6787D+00  0.1712D+00  0.9710D-01  0.3440D-01  0.1869D+00
 result A^{-1}*(A1|A2)
   0.1000D+01  0.0000D+00  0.0000D+00  0.0000D+00  0.4441D-15 -0.8384D+00  0.1973D+01 -0.2399D+01 -0.2755D+01
   0.0000D+00  0.1000D+01 -0.1776D-14  0.0000D+00  0.0000D+00  0.3764D+01 -0.5919D+01  0.8812D+01  0.7998D+01
   0.8882D-15  0.1776D-14  0.1000D+01  0.0000D+00  0.1776D-14  0.3852D+01 -0.4055D+01  0.7631D+01  0.7059D+01
   0.3553D-14  0.0000D+00  0.0000D+00  0.1000D+01  0.0000D+00 -0.7995D+01  0.9205D+01 -0.1698D+02 -0.1508D+02
  -0.8882D-15  0.0000D+00 -0.4441D-15 -0.8882D-15  0.1000D+01  0.2442D+01 -0.1512D+01  0.4768D+01  0.4470D+01

Black Belt

For small example problems efficiency may not be considered important, but note that if you are going to perform the same calculation with larger problems or do the calculation a large number of times, instead of forming the inverse of the leading m X m submatrix A1 and multiplying that inverse by the remaining submatrix A2, it would be more efficient to regard A2 as consisting of multiple right-hand sides, and use a single call to DGETRS instead of calling DGETRI and DGEMM as was done in #2. The corresponding lines of code:

      call dgetrs('N',m,n-m,c,m,ipiv,a(1,m+1),m,info)
      print *,"result A^{-1}*A2)"
      do i=1,m
          print 1,(a(i,j), j=m+1,n)
      end do