Issue with dgesdd in MKL

[bencteux]

Hello,
it happens that dgesdd returns negative singular values and "info=0".

To be more specific : I call dgesdd in a subroutine. Here is the piece of code.

[cpp]================================================         
         open(unit=12, file="c_3")
         l = 0
         do j = 1, 288
            do i = 1, 322
               l = l + 1
               read(12,'(e30.20)')rloc(l)
            enddo
         enddo
         close(12)

         ir1 = 1+322*288
         ir2 = ir1 + 322
         ir3 = ir2 + 322*322
         ir4 = ir3 + 288*288
         
         call dgesdd('A', 322, 288, rloc(1), 322,
     *        rloc(ir1),
     *        rloc(ir2), 322, rloc(ir3), 288,
     *        rloc(ir4), -1, iloc, info)

         mem = int(rloc(ir4))
         ir5 = ir4 + mem

         call dgesdd('A', 322, 288, rloc(1), 322,
     *        rloc(ir1),
     *        rloc(ir2), 322, rloc(ir3), 288,
     *        rloc(ir4), mem, iloc, info)
[/cpp]

===================================================

The 288 singular value which is returned is : -1, although info = 1, and no error is reported.

rloc and iloc are big size arrays, respectively : 7677071 and 6335, which is far more than needed. However, even when those arrays have exactly the needed size, it behaves exactly the same.

Of course, I suspect an error in my code and not in dgesdd : there must be something wrong with the rloc array. But I can not see how to track the bug.
Other precisions :
* the bug appears when I use MKL10.0, but not with a version of LAPACK 3.1.0 recompiled by myself with g77-3.3.5

Can someone give me any advice ? Thanks for help

[bencteux]

Hello,

I now have more details about this issue, but I am still puzzled with the behavior of dgesdd in MKL10.0. I would really appreciate some explanations.

Here are the new facts. I have made a very simple program : main program is in C. It allocates two arrays : a double array and a int array, then calls a fortran program. This one calls a subroutine that calls dgesdd, using a sub-portion of the double array. The complete code is reproduced at the end of the post. You can find the data file I use for "c_3" at : http://cermics.enpc.fr/~bencteux/c_3 .

Depending on the fact that the sub-array degins at a odd or an even integer (the variable "begin"), I get right or wrong singular values. By wrong singular values, I mean that one of the singular value is equal to -1 !!!

If I make all the program in fortran, it is just the other way around : I get wrong singular values when "begin" is even, and right ones when "begin" is odd !!

Does it mean dgesdd is not safe ? or MKL10.0 ? or both ?

Thanks for help,

Guy

================================

Code 1 : main program in C, calling fortran

#include
#include
#include
#include


void debugdgesdd_(double* rloc, int* nrloc, int* iloc, int* niloc);

int main(int argc, char * argv[]) {
double *rloc = 0;
int nrloc = 8432937, niloc = 2304, *iloc = 0;

rloc = (double*) calloc(nrloc, sizeof(double));
iloc = (int*) calloc(niloc, sizeof(int));
if (rloc == NULL || iloc == NULL) {
printf("Erreur a la creation du tableau rloc n");
exit (EXIT_FAILURE);
}
debugdgesdd_(rloc, &nrloc, iloc, &niloc);
}

Code 1 : the fortran part :

subroutine debugdgesdd(rloc, nrloc, iloc, niloc)
implicit none
integer nrloc, niloc, iloc(niloc)
double precision rloc(nrloc)
integer begin
c
begin = 2
call calcul(rloc(begin), 530501, iloc, niloc)
return
end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

subroutine calcul(rloc, nrloc, iloc, niloc)
implicit none
integer nrloc, niloc, iloc(niloc)
double precision rloc(nrloc)
c
integer l, i, j, ir1, ir2, ir3, ir4, ir5, mem, info
c
open(unit=12, file="c_3")
l = 0
do j = 1, 288
do i = 1, 322
l = l + 1
read(12,'(e30.20)')rloc(l)
enddo
enddo
close(12)
write(2,*) "Fichier c_3 relu. "
ir1 = 1+322*288
ir2 = ir1 + 288
ir3 = ir2 + 322*322
ir4 = ir3 + 288*288

call dgesdd('A', 322, 288, rloc(1), 322,
* rloc(ir1),
* rloc(ir2), 322, rloc(ir3), 288,
* rloc(ir4), -1, iloc, info)

mem = int(rloc(ir4))
ir5 = ir4 + mem

call dgesdd('A', 322, 288, rloc(1), 322,
* rloc(ir1),
* rloc(ir2), 322, rloc(ir3), 288,
* rloc(ir4), mem, iloc, info)

do i = 1, 288
write(2,'(i5,2x,e30.20)')i, rloc(ir1-1+i)
enddo
write(2,*) info, mem
write(2,*)"ir1, ir2, ir3, ir4, ir5, nrloc ",
* ir1, ir2, ir3, ir4, ir5, nrloc

return
end

=================================================

Code 2 : all in fortran

program debug_dgesdd
implicit none
integer nrloc, niloc
parameter(nrloc = 8432937, niloc = 2304)
integer iloc(niloc)
double precision rloc(nrloc)
c
call debugdgesdd(rloc, nrloc, iloc, niloc)
c
return
end

subroutine debugdgesdd(rloc, nrloc, iloc, niloc)
implicit none
integer nrloc, niloc, iloc(niloc)
double precision rloc(nrloc)
integer begin
c
begin = 2
call calcul(rloc(begin), 530501, iloc, niloc)
return
end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

subroutine calcul(rloc, nrloc, iloc, niloc)
implicit none
integer nrloc, niloc, iloc(niloc)
double precision rloc(nrloc)
c
integer l, i, j, ir1, ir2, ir3, ir4, ir5, mem, info
c
open(unit=12, file="c_3")
l = 0
do j = 1, 288
do i = 1, 322
l = l + 1
read(12,'(e30.20)')rloc(l)
enddo
enddo
close(12)
write(2,*) "Fichier c_3 relu. "
ir1 = 1+322*288
ir2 = ir1 + 288
ir3 = ir2 + 322*322
ir4 = ir3 + 288*288

call dgesdd('A', 322, 288, rloc(1), 322,
* rloc(ir1),
* rloc(ir2), 322, rloc(ir3), 288,
* rloc(ir4), -1, iloc, info)

mem = int(rloc(ir4))
ir5 = ir4 + mem

call dgesdd('A', 322, 288, rloc(1), 322,
* rloc(ir1),
* rloc(ir2), 322, rloc(ir3), 288,
* rloc(ir4), mem, iloc, info)

do i = 1, 288
write(2,'(i5,2x,e30.20)')i, rloc(ir1-1+i)
enddo
write(2,*) info, mem
write(2,*)"ir1, ir2, ir3, ir4, ir5, nrloc ",
* ir1, ir2, ir3, ir4, ir5, nrloc

return
end

[Michael_C_Intel4 · Accepted solution]

Hello,

the good idea is to submit an issue against dgesdd via Intel Premier Support.

Michael.