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
