Bugs in LAPACK95 GESDD (was: who maintains lapack95?)
in [Fortran]
|
Prev: 100% without investment online part time jobs..(adsense,datawork,neobux..more jobs)
Next: C++ / Fortran interoperability
From: mecej4 on 5 Aug 2010 09:30 kamaraju wrote: >> What seems to be implied above is that the bug is in the Fortran 95 >> interface to dgesdd rather than in the Lapack 3.x routine DGESDD. Is that >> the case? <--CUT--> > ! However, The problem is with the la_dgesdd.f90 routine assigning > wrong value > ! to JOBZ variable around line 202. It assigns JOBZ='A' where as in > this case, > ! it should have assigned JOBZ='S'. > ! > program dummy > > use f95_lapack, only : LA_GETRF, LA_GESDD, LA_GESVD > implicit none > > ! This works with both gfortran-4.3, gfortran-4.4 > ! real (kind=8) :: lhs(20,10), U(20,10), S(10), VT(10,10) > ! > ! This gives a segmentation fault with la_gesdd but not with > la_gesvd. The > ! segmentation fault occurs only with gfortran-4.4. But the bug is > in > ! la_dgesdd.f90 > real (kind=8) :: lhs(21,10), U(21,10), S(10), VT(10,10) > > call random_number(lhs) > call random_number(U) > call random_number(S) > call random_number(VT) > write(*,*) 'calling la_gesvd' > call LA_GESVD( lhs, S, U, VT) > write(*,*) 'calling la_gesdd' > call LA_GESDD( lhs, S, U, VT) > end program dummy > > The problem is like this. When la_dgesdd.f90 function is called with a > rectangular matrix A(M, N) where M>N, U(M,N) is supplied as the > argument, the idea is to compute only the min(m,n) columns of U. > However the code does not check for this case (the changes should be > made around line 202). At least for my sample code, the JOBZ variable > should have been set to 'S' . But la_dgesdd.f90 sets it to 'A'. > > If you need any other information, please let me know. > I have spent some time looking at the Netlib sources for DGESDD (F77) and LA_DGESDD (F90), and I now think that the problem is more involved than a simple error in the JOBZ value as noted above by Raju. The problem is really in the incomplete mapping of optional arguments in the F90 wrapper to the F77 base routine. When using DGESDD, JOBZ has possible values of 'A', 'S','O' and 'N'. For 'N', the arguments U and VT are not to be computed and the corresponding arguments are not accessed. But, this routine has a F77 interface and all arguments must be present. The F9X wrapper LA_DGESDD makes many arguments optional, creates work arrays as needed and then calls DGESDD. LA_DGESDD has an argument JOB, with default value 'N', and other accepted values of and 'U', 'V'. Your example calls LA_GESDD with JOB not present. Therefore, JOB takes on the default value, and the arguments U and VT _should_not_ be present. The problem is three-fold. 1. The mapping between JOB (3 values) and JOBZ (4 values) is incomplete, and if I am not mistaken, this is a design error in LA_DGESDD. In particular, as noted below, no value of JOB maps to JOBZ='A'. 2. The wrapper does not properly check a wrong call such as yours. 3. The description of the arguments at http://www.netlib.org/lapack95/lug95/node325.html is incomplete. It says how the case JOB='U' and U not present is covered, and similarly for 'V'. It is not clear how the JOBZ value of 'O'. The case when the user wants the full SVD, i.e, both U and VT are present -- the case covered in the F77 routine by JOBZ='A'-- appears to have no counterpart at all. The solution that you gave (set JOBZ='S' in line 202 of LA_DGESDD) makes your incorrect call work correctly, but will surely break on many correct calls. I urge Lapack users in C.L.F. to look at this matter and give us the benefit of their advice. -- mecej4
From: juha ruokolainen on 6 Aug 2010 03:44 On Aug 5, 4:30 pm, mecej4 <mecej4.nyets...(a)opFeramail.com> wrote: > kamaraju wrote: > >> What seems to be implied above is that the bug is in the Fortran 95 > >> interface to dgesdd rather than in the Lapack 3.x routine DGESDD. Is that > >> the case? > <--CUT--> > > ! However, The problem is with the la_dgesdd.f90 routine assigning > > wrong value > > ! to JOBZ variable around line 202. It assigns JOBZ='A' where as in > > this case, > > ! it should have assigned JOBZ='S'. > > ! > > program dummy > > > use f95_lapack, only : LA_GETRF, LA_GESDD, LA_GESVD > > implicit none > > > ! This works with both gfortran-4.3, gfortran-4.4 > > ! real (kind=8) :: lhs(20,10), U(20,10), S(10), VT(10,10) > > ! > > ! This gives a segmentation fault with la_gesdd but not with > > la_gesvd. The > > ! segmentation fault occurs only with gfortran-4.4. But the bug is > > in > > ! la_dgesdd.f90 > > real (kind=8) :: lhs(21,10), U(21,10), S(10), VT(10,10) > > > call random_number(lhs) > > call random_number(U) > > call random_number(S) > > call random_number(VT) > > write(*,*) 'calling la_gesvd' > > call LA_GESVD( lhs, S, U, VT) > > write(*,*) 'calling la_gesdd' > > call LA_GESDD( lhs, S, U, VT) > > end program dummy > > > The problem is like this. When la_dgesdd.f90 function is called with a > > rectangular matrix A(M, N) where M>N, U(M,N) is supplied as the > > argument, the idea is to compute only the min(m,n) columns of U. > > However the code does not check for this case (the changes should be > > made around line 202). At least for my sample code, the JOBZ variable > > should have been set to 'S' . But la_dgesdd.f90 sets it to 'A'. > > > If you need any other information, please let me know. > > I have spent some time looking at the Netlib sources for DGESDD (F77) and > LA_DGESDD (F90), and I now think that the problem is more involved than a > simple error in the JOBZ value as noted above by Raju. The problem is > really in the incomplete mapping of optional arguments in the F90 wrapper > to the F77 base routine. > > When using DGESDD, JOBZ has possible values of 'A', 'S','O' and 'N'. > For 'N', the arguments U and VT are not to be computed and the > corresponding arguments are not accessed. But, this routine has a F77 > interface and all arguments must be present. > > The F9X wrapper LA_DGESDD makes many arguments optional, creates work arrays > as needed and then calls DGESDD. LA_DGESDD has an argument JOB, with > default value 'N', and other accepted values of and 'U', 'V'. Your example > calls LA_GESDD with JOB not present. Therefore, JOB takes on the default > value, and the arguments U and VT _should_not_ be present. > > The problem is three-fold. > > 1. The mapping between JOB (3 values) and JOBZ (4 values) is incomplete, > and if I am not mistaken, this is a design error in LA_DGESDD. In > particular, as noted below, no value of JOB maps to JOBZ='A'. > > 2. The wrapper does not properly check a wrong call such as yours. > > 3. The description of the arguments athttp://www.netlib.org/lapack95/lug95/node325.htmlis incomplete. It says > how the case JOB='U' and U not present is covered, and similarly for 'V'. > It is not clear how the JOBZ value of 'O'. The case when the user wants the > full SVD, i.e, both U and VT are present -- the case covered in the F77 > routine by JOBZ='A'-- appears to have no counterpart at all. > > The solution that you gave (set JOBZ='S' in line 202 of LA_DGESDD) makes > your incorrect call work correctly, but will surely break on many correct > calls. > > I urge Lapack users in C.L.F. to look at this matter and give us the benefit > of their advice. > > -- mecej4 I think the 'JOB' argument to LA_DGESDD is only intended to tell with what values to replace the input matrix 'A' with, if 'U' and/or 'VT' were not present. Thus the original test case is OK, given the argument documentation. Juha
From: juha ruokolainen on 6 Aug 2010 03:50 On Aug 5, 4:30 pm, mecej4 <mecej4.nyets...(a)opFeramail.com> wrote: > kamaraju wrote: > >> What seems to be implied above is that the bug is in the Fortran 95 > >> interface to dgesdd rather than in the Lapack 3.x routine DGESDD. Is that > >> the case? > <--CUT--> > > ! However, The problem is with the la_dgesdd.f90 routine assigning > > wrong value > > ! to JOBZ variable around line 202. It assigns JOBZ='A' where as in > > this case, > > ! it should have assigned JOBZ='S'. > > ! > > program dummy > > > use f95_lapack, only : LA_GETRF, LA_GESDD, LA_GESVD > > implicit none > > > ! This works with both gfortran-4.3, gfortran-4.4 > > ! real (kind=8) :: lhs(20,10), U(20,10), S(10), VT(10,10) > > ! > > ! This gives a segmentation fault with la_gesdd but not with > > la_gesvd. The > > ! segmentation fault occurs only with gfortran-4.4. But the bug is > > in > > ! la_dgesdd.f90 > > real (kind=8) :: lhs(21,10), U(21,10), S(10), VT(10,10) > > > call random_number(lhs) > > call random_number(U) > > call random_number(S) > > call random_number(VT) > > write(*,*) 'calling la_gesvd' > > call LA_GESVD( lhs, S, U, VT) > > write(*,*) 'calling la_gesdd' > > call LA_GESDD( lhs, S, U, VT) > > end program dummy > > > The problem is like this. When la_dgesdd.f90 function is called with a > > rectangular matrix A(M, N) where M>N, U(M,N) is supplied as the > > argument, the idea is to compute only the min(m,n) columns of U. > > However the code does not check for this case (the changes should be > > made around line 202). At least for my sample code, the JOBZ variable > > should have been set to 'S' . But la_dgesdd.f90 sets it to 'A'. > > > If you need any other information, please let me know. > > I have spent some time looking at the Netlib sources for DGESDD (F77) and > LA_DGESDD (F90), and I now think that the problem is more involved than a > simple error in the JOBZ value as noted above by Raju. The problem is > really in the incomplete mapping of optional arguments in the F90 wrapper > to the F77 base routine. > > When using DGESDD, JOBZ has possible values of 'A', 'S','O' and 'N'. > For 'N', the arguments U and VT are not to be computed and the > corresponding arguments are not accessed. But, this routine has a F77 > interface and all arguments must be present. > > The F9X wrapper LA_DGESDD makes many arguments optional, creates work arrays > as needed and then calls DGESDD. LA_DGESDD has an argument JOB, with > default value 'N', and other accepted values of and 'U', 'V'. Your example > calls LA_GESDD with JOB not present. Therefore, JOB takes on the default > value, and the arguments U and VT _should_not_ be present. > > The problem is three-fold. > > 1. The mapping between JOB (3 values) and JOBZ (4 values) is incomplete, > and if I am not mistaken, this is a design error in LA_DGESDD. In > particular, as noted below, no value of JOB maps to JOBZ='A'. > > 2. The wrapper does not properly check a wrong call such as yours. > > 3. The description of the arguments athttp://www.netlib.org/lapack95/lug95/node325.htmlis incomplete. It says > how the case JOB='U' and U not present is covered, and similarly for 'V'. > It is not clear how the JOBZ value of 'O'. The case when the user wants the > full SVD, i.e, both U and VT are present -- the case covered in the F77 > routine by JOBZ='A'-- appears to have no counterpart at all. > > The solution that you gave (set JOBZ='S' in line 202 of LA_DGESDD) makes > your incorrect call work correctly, but will surely break on many correct > calls. > > I urge Lapack users in C.L.F. to look at this matter and give us the benefit > of their advice. > > -- mecej4
From: juha ruokolainen on 6 Aug 2010 03:55
On Aug 5, 4:30 pm, mecej4 <mecej4.nyets...(a)opFeramail.com> wrote: > kamaraju wrote: > >> What seems to be implied above is that the bug is in the Fortran 95 > >> interface to dgesdd rather than in the Lapack 3.x routine DGESDD. Is that > >> the case? > <--CUT--> > > ! However, The problem is with the la_dgesdd.f90 routine assigning > > wrong value > > ! to JOBZ variable around line 202. It assigns JOBZ='A' where as in > > this case, > > ! it should have assigned JOBZ='S'. > > ! > > program dummy > > > use f95_lapack, only : LA_GETRF, LA_GESDD, LA_GESVD > > implicit none > > > ! This works with both gfortran-4.3, gfortran-4.4 > > ! real (kind=8) :: lhs(20,10), U(20,10), S(10), VT(10,10) > > ! > > ! This gives a segmentation fault with la_gesdd but not with > > la_gesvd. The > > ! segmentation fault occurs only with gfortran-4.4. But the bug is > > in > > ! la_dgesdd.f90 > > real (kind=8) :: lhs(21,10), U(21,10), S(10), VT(10,10) > > > call random_number(lhs) > > call random_number(U) > > call random_number(S) > > call random_number(VT) > > write(*,*) 'calling la_gesvd' > > call LA_GESVD( lhs, S, U, VT) > > write(*,*) 'calling la_gesdd' > > call LA_GESDD( lhs, S, U, VT) > > end program dummy > > > The problem is like this. When la_dgesdd.f90 function is called with a > > rectangular matrix A(M, N) where M>N, U(M,N) is supplied as the > > argument, the idea is to compute only the min(m,n) columns of U. > > However the code does not check for this case (the changes should be > > made around line 202). At least for my sample code, the JOBZ variable > > should have been set to 'S' . But la_dgesdd.f90 sets it to 'A'. > > > If you need any other information, please let me know. > > I have spent some time looking at the Netlib sources for DGESDD (F77) and > LA_DGESDD (F90), and I now think that the problem is more involved than a > simple error in the JOBZ value as noted above by Raju. The problem is > really in the incomplete mapping of optional arguments in the F90 wrapper > to the F77 base routine. > > When using DGESDD, JOBZ has possible values of 'A', 'S','O' and 'N'. > For 'N', the arguments U and VT are not to be computed and the > corresponding arguments are not accessed. But, this routine has a F77 > interface and all arguments must be present. > > The F9X wrapper LA_DGESDD makes many arguments optional, creates work arrays > as needed and then calls DGESDD. LA_DGESDD has an argument JOB, with > default value 'N', and other accepted values of and 'U', 'V'. Your example > calls LA_GESDD with JOB not present. Therefore, JOB takes on the default > value, and the arguments U and VT _should_not_ be present. > > The problem is three-fold. > > 1. The mapping between JOB (3 values) and JOBZ (4 values) is incomplete, > and if I am not mistaken, this is a design error in LA_DGESDD. In > particular, as noted below, no value of JOB maps to JOBZ='A'. > > 2. The wrapper does not properly check a wrong call such as yours. > > 3. The description of the arguments athttp://www.netlib.org/lapack95/lug95/node325.htmlis incomplete. It says > how the case JOB='U' and U not present is covered, and similarly for 'V'. > It is not clear how the JOBZ value of 'O'. The case when the user wants the > full SVD, i.e, both U and VT are present -- the case covered in the F77 > routine by JOBZ='A'-- appears to have no counterpart at all. > > The solution that you gave (set JOBZ='S' in line 202 of LA_DGESDD) makes > your incorrect call work correctly, but will surely break on many correct > calls. > > I urge Lapack users in C.L.F. to look at this matter and give us the benefit > of their advice. > > -- mecej4 I think the 'JOB' argument to LA_DGESDD is only intended to tell with what valuesto replace the input matrix 'A' with, if 'U' and/or 'VT' were not present. Thus the original test case is OK, given the argument documentation. Juha |