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 <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>


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