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

problem with dgetrf

kon-konstantinov
Beginner
‎10-04-2011 10:41 PM
1,399 Views
Hi,
I have problems calling dgetrf. A small test program is compiled as follows:
/opt/intel/bin/ifort test_intel.f90 -L/opt/intel/composerxe-2011.2.137/mkl/lib/intel64 -I/opt/intel/composerxe-2011.2.137/mkl/include/intel64/ilp64 -lmkl_blas95_ilp64 -lmkl_lapack95_ilp64 -Wl,--start-group /opt/intel/composerxe-2011.2.137/mkl/lib/intel64/libmkl_intel_ilp64.a /opt/intel/composerxe-2011.2.137/mkl/lib/intel64/libmkl_intel_thread.a /opt/intel/composerxe-2011.2.137/mkl/lib/intel64/libmkl_core.a -Wl,--end-group -openmp -lpthread -o test_intel
When I use:

call dgetrf( A, IPIV, INFO )

I get the following:
The matrix has dimensions: n = m = lda = 4
Original matrix A
5.00 7.00 6.00 5.00
7.00 10.00 8.00 7.00
6.00 8.00 10.00 9.00
5.00 7.00 9.00 10.00
forrtl: severe (174): SIGSEGV, segmentation fault occurred
Image PC Routine Line Source
test_intel 0000000000408253 Unknown Unknown Unknown
test_intel 000000000040817C Unknown Unknown Unknown
test_intel 0000000000407DCD Unknown Unknown Unknown
test_intel 000000000040799C Unknown Unknown Unknown
libc.so.6 000000374FA1D994 Unknown Unknown Unknown
test_intel 00000000004078A9 Unknown Unknown Unknown

If I use fortran 77 call:
call dgetrf( m, n, a, lda, ipiv, info )

I get :
MKL ERROR: Parameter 4 was incorrect on entry to DGETRF
parameter INFO is : -4

I would greatly appreciate if somebody can help me to solve my problem. Thanks in advance.

Kon
0 Kudos

Link Copied

17 Replies
mecej4
Black Belt
‎10-04-2011 11:39 PM
1,399 Views
If you want to use the F95 interface with optional arguments, use the generic name GETRF, and USE mkl95_lapack to pick up the interface.

[fortran]program getrfEx
use mkl95_lapack,only : getrf
implicit none
integer, parameter :: N = 4
integer :: i,j
real(8), dimension(N,N) :: A = & 
  (/ 5,7,6,5,7,10,8,7,6,8,10,9,5,7,9,10/)
integer, dimension(N) :: piv
integer :: info

call getrf( a, piv, info )
write(*,10)((a(i,j),j=1,N),i=1,N)
stop

10 format(1x,4F12.6)  
end program getrfEx
[/fortran]
If you call the F77 routine, you must make sure that your variables are declared with the correct types and contain proper values. Since you showed no declarations or initializations, I cannot say what was wrong.
0 Kudos
Copy link
Gennady_F_Intel
Moderator
‎10-05-2011 06:34 AM
1,399 Views
Kon,
please try to relink this example with LP64 not iLP64 libraries and check the problem again.
--Gennady
0 Kudos
Copy link
kon-konstantinov
Beginner
‎10-05-2011 02:47 PM
1,399 Views

Hi,

Here is my code compiled with:

/opt/intel/bin/ifort test_intel.f90 -L/opt/intel/composerxe-2011.2.137/mkl/lib/intel64 -I/opt/intel/composerxe-2011.2.137/mkl/include/intel64/ilp64 -lmkl_blas95_ilp64 -lmkl_lapack95_ilp64 -Wl,--start-group /opt/intel/composerxe-2011.2.137/mkl/lib/intel64/libmkl_intel_ilp64.a /opt/intel/composerxe-2011.2.137/mkl/lib/intel64/libmkl_intel_thread.a /opt/intel/composerxe-2011.2.137/mkl/lib/intel64/libmkl_core.a -Wl,--end-group -openmp -lpthread -o test_intel


program test_lapack
USE mkl95_LAPACK, ONLY: GETRF
!implicit none
integer :: i, j, k, l, m, n, LDA, INFO, LWORK
real(8), dimension(:,:), allocatable :: A
real(8), dimension(:), allocatable :: WORK
integer, dimension(:), allocatable :: IPIV
n = 4; m = 4
LDA = n
allocate(A(LDA,n), IPIV(n), WORK(n) )
! define lower half of matrix
A(1,1)= 5;
A(2,1)=7; A(2,2)=10
A(3,1)=6; A(3,2)=8; A(3,3)= 10
A(4,1)=5; A(4,2)=7; A(4,3)=9; A(4,4)= 10
! define upper half by symmetry
do i = 1, n
do j = i+1, n
A(i,j)=A(j,i)
end do
end do
print*, " Original matrix A "
do i = 1, n
write(*,fmt='(4f10.2)') A(i,:)
enddo
! call dsytrf( a, ipiv, info )
! call dsytri( a, ipiv, 'L', info )
! call dpotrf( A, 'L', INFO )
! call dpotrf( 'L', n, a, lda, info )
call dgetrf( A, IPIV, INFO )
! call dgetrf( m, n, a, lda, ipiv, info )
print*, " parameter INFO is : ", INFO
print*, " LU decomposition using LAPACK DGETRF routine "
!print*, " Cholesky decomposition using LAPACK DPOTRF routine "
do i = 1, n
write(*,fmt='(4f10.2)') A(i,:)
enddo
LWORK = n

! call dgetri( a, ipiv ,info )
!print*, " parameter INFO is : ", INFO
!print*, " Matrix inverse using DGETRI after LU decomposition using DGETRF routines "
!!print*, " Matrix inverse using DPOTRI after Cholesky decomposition using DPOTRF routines "
!do i = 1, n
!write(*,fmt='(4f10.2)') A(i,:)
!enddo

end program test_lapack

This is the result:
Original matrix A
5.00 7.00 6.00 5.00
7.00 10.00 8.00 7.00
6.00 8.00 10.00 9.00
5.00 7.00 9.00 10.00
forrtl: severe (174): SIGSEGV, segmentation fault occurred
Image PC Routine Line Source
test_intel 0000000000408253 Unknown Unknown Unknown
test_intel 000000000040817C Unknown Unknown Unknown
test_intel 0000000000407DCD Unknown Unknown Unknown
test_intel 000000000040799C Unknown Unknown Unknown
libc.so.6 000000368461D994 Unknown Unknown Unknown
test_intel 00000000004078A9 Unknown Unknown Unknown


0 Kudos
Copy link
kon-konstantinov
Beginner
‎10-05-2011 02:57 PM
1,399 Views
Hi Gennady,

Thanks for your reply. I tryed to compile with LP64 instead of iLP64 but I had the same problem. Here is my code:

program test_lapack
USE mkl95_LAPACK, ONLY: GETRF
!implicit none
integer :: i, j, k, l, m, n, LDA, INFO, LWORK
real(8), dimension(:,:), allocatable :: A
real(8), dimension(:), allocatable :: WORK
integer, dimension(:), allocatable :: IPIV
n = 4; m = 4
LDA = n
allocate(A(LDA,n), IPIV(n), WORK(n) )
! define lower half of matrix
A(1,1)= 5;
A(2,1)=7; A(2,2)=10
A(3,1)=6; A(3,2)=8; A(3,3)= 10
A(4,1)=5; A(4,2)=7; A(4,3)=9; A(4,4)= 10
! define upper half by symmetry
do i = 1, n
do j = i+1, n
A(i,j)=A(j,i)
end do
end do
print*, " Original matrix A "
do i = 1, n
write(*,fmt='(4f10.2)') A(i,:)
enddo
! call dsytrf( a, ipiv, info )
! call dsytri( a, ipiv, 'L', info )
! call dpotrf( A, 'L', INFO )
! call dpotrf( 'L', n, a, lda, info )
call dgetrf( A, IPIV, INFO )
! call dgetrf( m, n, a, lda, ipiv, info )
print*, " parameter INFO is : ", INFO
print*, " LU decomposition using LAPACK DGETRF routine "
!print*, " Cholesky decomposition using LAPACK DPOTRF routine "
do i = 1, n
write(*,fmt='(4f10.2)') A(i,:)
enddo
LWORK = n

! call dgetri( a, ipiv ,info )
!print*, " parameter INFO is : ", INFO
!print*, " Matrix inverse using DGETRI after LU decomposition using DGETRF routines "
!!print*, " Matrix inverse using DPOTRI after Cholesky decomposition using DPOTRF routines "
!do i = 1, n
!write(*,fmt='(4f10.2)') A(i,:)
!enddo

end program test_lapack

This is the result:

Original matrix A
5.00 7.00 6.00 5.00
7.00 10.00 8.00 7.00
6.00 8.00 10.00 9.00
5.00 7.00 9.00 10.00
forrtl: severe (174): SIGSEGV, segmentation fault occurred
Image PC Routine Line Source
test_intel 0000000000408163 Unknown Unknown Unknown
test_intel 0000000000407DCD Unknown Unknown Unknown
test_intel 000000000040799C Unknown Unknown Unknown
libc.so.6 000000368461D994 Unknown Unknown Unknown
test_intel 00000000004078A9 Unknown Unknown Unknown

I am stuck, I don't know what to do.
0 Kudos
Copy link
mecej4
Black Belt
‎10-05-2011 04:24 PM
1,399 Views
The solution was given above, in the first sentence of Response #1.
0 Kudos
Copy link
Konstantin_A_Intel
Employee
‎10-05-2011 08:19 PM
1,399 Views
Hi Kosi,

mecej4 is right, you used f95 interface incorrectly. Correct interfaces are as follows:

FORTRAN 77:

call sgetrf(m,n,a,lda,ipiv,info)

call dgetrf(m,n,a,lda,ipiv,info)

call cgetrf(m,n,a,lda,ipiv,info)

call zgetrf(m,n,a,lda,ipiv,info)

Fortran 95:

call getrf(a[,ipiv][,info])


But in your program you used F77 name (dgetrf) with the arguments of F95 interface. Please use either
call getrf( A, IPIV, INFO )
or
call dgetrf( m, n, a, lda, ipiv, info )

Also make sure you've changed ALL "ilp64" substrings in your compilation link with "lp64". In my program it looked like:

ifort test.f90 -L$MKLROOT/lib/intel64 -I$MKLROOT/include/intel64 -lmkl_blas95_lp64 -lmkl_lapack95_lp64 -Wl,--start-group $MKLROOT/lib/intel64/libmkl_intel_lp64.a $MKLROOT/lib/intel64/libmkl_intel_thread.a $MKLROOT/lib/intel64/libmkl_core.a -Wl,--end-group -openmp -lpthread -o test_intel


bash-4.1$ ./run.sh
Original matrix A
5.00 7.00 6.00 5.00
7.00 10.00 8.00 7.00
6.00 8.00 10.00 9.00
5.00 7.00 9.00 10.00
parameter INFO is : 0
LU decomposition using LAPACK DGETRF routine
7.00 10.00 8.00 7.00
0.86 -0.57 3.14 3.00
0.71 0.25 2.50 4.25
0.71 0.25 -0.20 0.10

0 Kudos
Copy link
kon-konstantinov
Beginner
‎10-06-2011 05:04 PM
1,399 Views
Hi Konstantin,

Thank you very much for your help. I followed your instructions and it worked, I got the same answers.
However, my final goal is to inver the matrix. For this purpose I used GETRI subroutine.
I've added to my program:

USE mkl95_LAPACK, ONLY: GETRF, GETRI
call getri( A, IPIV ,INFO )
The result is the same as from GETRF:
Original matrix A
5.00 7.00 6.00 5.00
7.00 10.00 8.00 7.00
6.00 8.00 10.00 9.00
5.00 7.00 9.00 10.00
parameter INFO is : 0
Cholesky decomposition using LAPACK DPOTRF routine
7.00 10.00 8.00 7.00
0.86 -0.57 3.14 3.00
0.71 0.25 2.50 4.25
0.71 0.25 -0.20 0.10
parameter INFO is : 0
Matrix inverse using GETRI after using GETRF routines
7.00 10.00 8.00 7.00
0.86 -0.57 3.14 3.00
0.71 0.25 2.50 4.25
0.71 0.25 -0.20 0.10
Again I am doing something wrong. Could you help me again please.
One more question. I need a routine that calculates a generalised inverse of non-positive definite matrix.
Do you know such routine from LAPACK. If yes can you tell me how to use it. Thanks very much again.
By the way my first name is also Konstantin.

Cheers,

Kon
0 Kudos
Copy link
mecej4
Black Belt
‎10-06-2011 08:23 PM
1,399 Views
You probably forgot to CALL GETRI. It is not possible to comment on what could be wrong if you do not show us your code.
0 Kudos
Copy link
Konstantin_A_Intel
Employee
‎10-06-2011 08:46 PM
1,399 Views
In order to get inverse of symmetrical indefinite matrix please use a pair of sytrf/sytri function:

call sytrf( a, uplo, ipiv,info )
call sytri( a, ipiv, uplo, info )

uplo should be either 'U' or 'L' depending on upper or lower part of your matrix is meaninful.

Please refer to online MKL documentation for more details of any functionality.

http://software.intel.com/sites/products/documentation/hpc/mkl/updates/10.3.5/mklman/index.htm

Regards,
Konstantin

0 Kudos
Copy link
kon-konstantinov
Beginner
‎10-06-2011 09:16 PM
1,399 Views
Hi,
I didn't forget to call GETRI. Here is my code compiled with:
/opt/intel/bin/ifort test_intel.f90 -L/opt/intel/composerxe-2011.2.137/mkl/lib/intel64 -I/opt/intel/composerxe-2011.2.137/mkl/include/intel64/lp64 -lmkl_blas95_ilp64 -lmkl_lapack95_ilp64 -Wl,--start-group /opt/intel/composerxe-2011.2.137/mkl/lib/intel64/libmkl_intel_lp64.a /opt/intel/composerxe-2011.2.137/mkl/lib/intel64/libmkl_intel_thread.a /opt/intel/composerxe-2011.2.137/mkl/lib/intel64/libmkl_core.a -Wl,--end-group -openmp -lpthread -o test_intel




program test_lapack
!Serial number : CZMG-JX2MSWXW
!use deeba
! USE mkl95_PRECISION, ONLY: WP => SP
USE mkl95_LAPACK, ONLY: GETRF, GETRI
! USE mkl95_LAPACK, ONLY: GETRI
!implicit none
integer :: i, j, k, l, m, n, LDA, INFO, LWORK

real(8), dimension(:,:), allocatable :: A
real(8), dimension(:), allocatable :: WORK
integer, dimension(:), allocatable :: IPIV
CHARACTER(LEN=250) :: fin, fout

!PRINT*, ' Enter input file name i : '
!READ(*,*) fin(i)
!PRINT*, ' Enter output file name: '
!READ(*,*) fout
n = 4; m = 4
LDA = n
allocate(A(LDA,n), IPIV(n), WORK(n) )
! define lower half of matrix
A(1,1)= 5;
A(2,1)=7; A(2,2)=10
A(3,1)=6; A(3,2)=8; A(3,3)= 10
A(4,1)=5; A(4,2)=7; A(4,3)=9; A(4,4)= 10


! define upper half by symmetry
do i = 1, n
do j = i+1, n
A(i,j)=A(j,i)
end do
end do
print*, " Original matrix A "
do i = 1, n
write(*,fmt='(4f10.2)') A(i,:)
enddo



! call dsytrf( a, ipiv, info )
! call dsytri( a, ipiv, 'L', info )
! call dpotrf( A, 'L', INFO )
! call dpotrf( 'L', n, a, lda, info )

call getrf( A, IPIV, INFO )
! call dgetrf( m, n, a, lda, ipiv, info )



print*, " parameter INFO is : ", INFO
!print*, " LU decomposition using LAPACK DGETRF routine "
print*, " Cholesky decomposition using LAPACK DPOTRF routine "
do i = 1, n
write(*,fmt='(4f10.2)') A(i,:)
enddo
LWORK = n

call getri( A, IPIV ,INFO )
print*, " parameter INFO is : ", INFO
print*, " Matrix inverse using GETRI after LU decomposition using GETRF routines "
!print*, " Matrix inverse using DPOTRI after Cholesky decomposition using DPOTRF routines "
do i = 1, n
write(*,fmt='(4f10.2)') A(i,:)
enddo




end program test_lapack

This is the result:

Original matrix A
5.00 7.00 6.00 5.00
7.00 10.00 8.00 7.00
6.00 8.00 10.00 9.00
5.00 7.00 9.00 10.00
parameter INFO is : 0
Cholesky decomposition using LAPACK DPOTRF routine
7.00 10.00 8.00 7.00
0.86 -0.57 3.14 3.00
0.71 0.25 2.50 4.25
0.71 0.25 -0.20 0.10
parameter INFO is : 0
Matrix inverse using DGETRI after using DGETRF routines
7.00 10.00 8.00 7.00
0.86 -0.57 3.14 3.00
0.71 0.25 2.50 4.25
0.71 0.25 -0.20 0.10


0 Kudos
Copy link
Konstantin_A_Intel
Employee
‎10-06-2011 09:38 PM
1,399 Views
Hi Konstantin,

It seems you're still compiling your code with ilp64 modules that leads to incorrect results. Please replace-lmkl_blas95_ilp64 -lmkl_lapack95_ilp64with-lmkl_blas95_lp64 -lmkl_lapack95_lp64. It should work.
0 Kudos
Copy link
kon-konstantinov
Beginner
‎10-13-2011 09:50 PM
1,399 Views
Hi Konstantin,

I've done that and it worked. However, I have another problem. I want to use dgemm from BLAS.
Here is my code:

program test_blas

USE mkl95_BLAS, ONLY: gemm
!implicit none
integer :: i, j, k, l, m, n, lda
real(8), DIMENSION(:,:), ALLOCATABLE :: A, B, C

real(kind=8)::alpha,beta


alpha=1.0;beta=1.0
l
allocate (A(2,3), B(3,2), C(2,2) )

A(1,1) = 1; A(1,2) = 2; A(1,3) = 3
A(2,1) = 4; A(2,2) = 5; A(2,3) = 6
!
B(1,1) = 2; B(1,2) = 3
B(2,1) = 4; B(2,2) = 5
B(3,1) = 6; B(3,2) = 7
Print*, " Original matriz A : "
do i = 1, 2
write(*,fmt='(13f10.3)') a(i,:)
enddo
Print*, " Original matriz B : "
do i = 1, 3
write(*,fmt='(13f10.3)') b(i,:)
enddo


call dgemm('N', 'N', 2, 2, 3, alfa, a, 2, b, 3, 1.0, c, 2)
print*, " Product matrix C : "

do i = 1, 2
write(*,fmt='(13f10.3)') c(i,:)
enddo

end program test_blas

The result is:
Product matrix C :
0.000 0.000
0.000 0.000

The program is compiled with:
/opt/intel/bin/ifort test_blas.f90 -L/opt/intel/composerxe-2011.2.137/mkl/lib/intel64 -I/opt/intel/composerxe-2011.2.137/mkl/include/intel64/lp64 -lmkl_blas95_lp64 -lmkl_lapack95_lp64 -Wl,--start-group /opt/intel/composerxe-2011.2.137/mkl/lib/intel64/libmkl_intel_lp64.a /opt/intel/composerxe-2011.2.137/mkl/lib/intel64/libmkl_intel_thread.a /opt/intel/composerxe-2011.2.137/mkl/lib/intel64/libmkl_core.a -Wl,--end-group -openmp -lpthread -fast -o test_blas

What is wrong again? Thanks.

Cheers,

Kon
0 Kudos
Copy link
mecej4
Black Belt
‎10-13-2011 09:57 PM
1,399 Views
Explain why you wrote

call dgemm('N', 'N', 2, 2, 2, alfa, a, 2, b, 2, beta, c, 2)

How is the subroutine to know that matrix a is 2 X 3 and matrix b is 3 X 2, when there is no '3' or variable with value 3 anywhere in the argument list?
0 Kudos
Copy link
kon-konstantinov
Beginner
‎10-13-2011 10:28 PM
1,399 Views
OK that was a mistake.
The new call is:
call dgemm('N', 'N', 2, 3, 3, alfa, a, 2, b, 3, beta, c, 2)
and the result is the same. What is wrong?
0 Kudos
Copy link
Konstantin_A_Intel
Employee
‎10-13-2011 11:26 PM
1,399 Views
Please try this one:

call dgemm('N', 'N', 2, 2, 3, alfa, a, 2, b, 3, beta, c, 2)
0 Kudos
Copy link
kon-konstantinov
Beginner
‎10-13-2011 11:40 PM
1,399 Views
I've used:
call dgemm('N', 'N', 2, 2, 3, alfa, a, 2, b, 3, beta, c, 2)

but still getting zero matrix.
0 Kudos
Copy link
mecej4
Black Belt
‎10-13-2011 11:57 PM
1,399 Views
What is the value of alfa? Why don't you try using IMPLICIT NONE? Or include interfaces for MKL routines:

!DEC$NOFREEFORM
include 'mkl_blas.fi'
!DEC$FREEFORM


0 Kudos
Copy link
Reply

Community support is provided during standard business hours (Monday to Friday 7AM - 5PM PST). Other contact methods are available here.

Intel does not verify all solutions, including but not limited to any file transfers that may appear in this community. Accordingly, Intel disclaims all express and implied warranties, including without limitation, the implied warranties of merchantability, fitness for a particular purpose, and non-infringement, as well as any warranty arising from course of performance, course of dealing, or usage in trade.

For more complete information about compiler optimizations, see our Optimization Notice.