Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
6534 Discussions

## Reduced Row Echelon Form of a matrix (rref) Beginner
259 Views

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.

2 Replies Employee
259 Views

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 Black Belt
259 Views

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
``` 