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 two very simple programs. They are 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 .

1) main program is in C. It allocates two arrays : a double array

rloc and a int array

iloc, then calls a fortran program :

debugdgesdd. This one calls a subroutine,

calcul, that loads the data and calls

dgesdd, using a sub-portion of the double array

rloc.

2) is identical to 1) except there is no C part : all is in Fortran.

In 1), depending on the fact that the sub-array begins 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 returned the singular value is -1 !!!

In 2), 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)C part

- Code: Select all
`#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);

}

Fortran part :

- Code: Select all
`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) - Code: Select all
`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