Reduced Row Echelon Form of a matrix (rref)

safa_o_ · ‎02-24-2015

Hi,

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

Victor_K_Intel1 · ‎02-26-2015

Hi,

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
              a(i,j)=b(j,i)
              d(i,j)=0D0
          end do
      end do
      
      do i=1,m
          do j=1,m
              c(i,j)=a(i,j)
          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(info.ne.0) then
          print *,"dgetrf returned info=",info
          stop
      end if
      
      call dgetri( m, c, m, ipiv, work, m, info )
      if(info.ne.0) then
          print *,"dgetri returned info=",info
          stop
      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    
      
      stop
   1  format(1x, 9(1x, D11.4)) 
      end

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

mecej4 · ‎03-03-2015

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 A₁ and multiplying that inverse by the remaining submatrix A₂, it would be more efficient to regard A₂ 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