Another issue with dgesdd

Open discussion regarding features, bugs, issues, vendors, etc.

Another issue with dgesdd

Postby bencteux » Sat Dec 06, 2008 5:59 am

Hello,
it also arrives 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.

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


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
Guy
bencteux
 
Posts: 13
Joined: Wed Apr 16, 2008 11:52 am

Re: Another issue with dgesdd

Postby bencteux » Sun Dec 07, 2008 5:54 am

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
bencteux
 
Posts: 13
Joined: Wed Apr 16, 2008 11:52 am

Re: Another issue with dgesdd

Postby Julien Langou » Fri Dec 12, 2008 3:56 pm

Bonjour Guy,

1- if INFO=1 in output of DGESDD, I am not surprised that you get garbage in your singular value array (e.g. negative singular values). INFO=0 is the flag to tell you everything went fine.

2- It seems that MKL has a problem while LAPACK-3.1 is working fine. Is this the case? In that case, maybe Intel fellows would be interested to here from you and they probably be more helpful than us.

Best wishes,
Julien.
Julien Langou
 
Posts: 727
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Another issue with dgesdd

Postby bencteux » Fri Dec 12, 2008 4:11 pm

Bonjour Julien,
in those tests, dgesdd always retruns "info = 0".

In fact, further tests have shown me that it is not only the problem of MKL.
I have some matrices for which MKL fails and not the debian package of lapack/blas.
But for some others, it is just the opposite : the debian package fails and not the MKL.

Are there bug reports on dgesdd ?

Thanks,
Guy
bencteux
 
Posts: 13
Joined: Wed Apr 16, 2008 11:52 am

Re: Another issue with dgesdd

Postby Julien Langou » Fri Dec 12, 2008 4:16 pm

Well then yes, if DGESDD (reference LAPACK) returns with INFO = 0 and negative singular values that's pretty bad. So now the first thing to do is to look with a more careful eyes to your codes. I do not have time now. I hope someone else can have a look. Julien.
Julien Langou
 
Posts: 727
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Another issue with dgesdd

Postby bencteux » Sat Dec 13, 2008 3:22 am

I also do not have much time before Christmas, but I will come back on the forum with this issue and with simple test cases in January.
Thanks
Guy
bencteux
 
Posts: 13
Joined: Wed Apr 16, 2008 11:52 am

Re: Another issue with dgesdd

Postby joshtanga » Thu Jan 29, 2009 2:31 am

I am noticing a related problem. I am using the ACML BLAS and calling dgesdd through MATLAB (thanks for the previous help). However, successive runs of my code on the SAME matrix produces singular vectors whose values differ by a factor of -1. This occurs even though info is returned as 0. An portion of output from two successive runs is below.

Code: Select all
vv =
   -0.5428   -0.8297   -0.1300
   -0.8380    0.5248    0.1494
   -0.0558    0.1901   -0.9802


Code: Select all
vv =
    0.5428    0.8297   -0.1300
   -0.8380    0.5248   -0.1494
    0.0558   -0.1901   -0.9802
joshtanga
 
Posts: 6
Joined: Tue Oct 21, 2008 10:24 pm

Re: Another issue with dgesdd

Postby bencteux » Thu Jan 29, 2009 9:04 am

Singular vectors are defined up to a "-1" factor (or, if you are working with complex matrix, up to a multiplicative constant of module 1),
since they have to verify : A u_i = sigma_i v_i.
Here sigma_i is the singular value, u_i and v_i are singular vectors.
Thus, there is no problem with your results.

In this post, I was telling about a singular value that was equal to -1. This revealed an issue in the computation since the singular values must be non negative.

Guy
bencteux
 
Posts: 13
Joined: Wed Apr 16, 2008 11:52 am

Re: Another issue with dgesdd

Postby Julie » Fri Sep 25, 2009 4:12 pm

Guy,
I was trying to reproduce your problem with our latest version of LAPACK, but I cannot reproduce your bug.
The C and Fortran version with begin sets to 1 or 2 give all the same results.
Julie
Julie
 
Posts: 299
Joined: Wed Feb 23, 2005 12:32 am
Location: ICL, Denver. Colorado

Re: Another issue with dgesdd

Postby bencteux » Mon Sep 28, 2009 2:58 am

Thank you Julie.It probably means that the bug has been fixed.
I can go back to dgesdd now !
Guy
bencteux
 
Posts: 13
Joined: Wed Apr 16, 2008 11:52 am


Return to User Discussion

Who is online

Users browsing this forum: No registered users and 2 guests