problem with dgetrf

[kon-konstantinov]

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

[mecej4]

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.

[Gennady_F_Intel]

Kon,

please try to relink this example with LP64 not iLP64 libraries and check the problem again.

--Gennady

[kon-konstantinov]

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

[kon-konstantinov]

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.

[mecej4]

The solution was given above, in the first sentence of Response #1.

[Konstantin_A_Intel]

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

[kon-konstantinov]

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

[mecej4]

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.

[Konstantin_A_Intel]

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

[kon-konstantinov]

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

[Konstantin_A_Intel]

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.

[kon-konstantinov]

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

[mecej4]

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?

[kon-konstantinov]

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?

[Konstantin_A_Intel]

Please try this one:

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

[kon-konstantinov]

I've used:
call dgemm('N', 'N', 2, 2, 3, alfa, a, 2, b, 3, beta, c, 2)

but still getting zero matrix.

[mecej4]

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