qr decomposition fails for certain matrices

[chrysn]

(i hope this doesn't get filed twice; i reported it from within octave,
but it seems to have missed the list)

--------
Bug report for Octave 3.2.3 configured for x86_64-pc-linux-gnu

Description:
-----------

  * for the attached matrix "bad", the qr decomposition creates NaN and Inf

Repeat-By:
---------

  * load "bad"
  * do `[Q,R] = qr(bad)` or `[Q,R,P] = qr(bad)`; you get several +/-Inf and
    NaN.

(for comparison, matlab can decompose this; moreover, the qr decomposition
should work with any given input)



Configuration (please do not edit this section):
-----------------------------------------------

uname output:     Linux hephaistos 2.6.30-2-amd64 #1 SMP Fri Sep 25 22:16:56 UTC 2009 x86_64 GNU/Linux
configure opts:   '--build=x86_64-linux-gnu' '--prefix=/usr' '--datadir=/usr/share' '--libdir=/usr/lib' '--libexecdir=/usr/lib' '--infodir=/usr/share/info' '--mandir=/usr/share/man' '--with-blas=-lblas-3gf' '--with-lapack=-llapackgf-3' '--with-hdf5' '--with-fftw' '--enable-shared' '--enable-rpath' '--disable-static' 'build_alias=x86_64-linux-gnu' 'CC=gcc' 'CFLAGS=-g -O2' 'LDFLAGS=' 'CPPFLAGS=' 'CXX=g++' 'CXXFLAGS=-g -O2' 'F77=gfortran' 'FFLAGS=-g -O2'
Fortran compiler: gfortran
FFLAGS:           -O2 -g
FLIBS:            -lgfortranbegin -lgfortran
CPPFLAGS:        
INCFLAGS:         -I. -I. -I./liboctave -I./src -I./libcruft/misc
C compiler:       gcc, version 4.3.4 (Debian 4.3.4-4)
CFLAGS:           -O2 -g
CPICFLAG:         -fPIC
C++ compiler:     g++, version 4.3.4
CXXFLAGS:         -O2 -g
CXXPICFLAG:       -fPIC
LD_CXX:           g++
LDFLAGS:          
LIBFLAGS:         -L.
RLD_FLAG:         -Wl,-rpath -Wl,/usr/lib/octave-3.2.3
BLAS_LIBS:        -llapackgf-3 -lblas-3gf
FFTW_LIBS:        -lfftw3 -lfftw3f
LIBS:             -lreadline  -lncurses -ldl -lhdf5 -lz -lm
LEXLIB:          
LIBGLOB:          
SED:              /bin/sed
DEFS:

  -DPACKAGE_NAME="" -DPACKAGE_TARNAME="" -DPACKAGE_VERSION=""
  -DPACKAGE_STRING="" -DPACKAGE_BUGREPORT="" -DPACKAGE_URL=""
  -DOCTAVE_SOURCE=1 -DSTDC_HEADERS=1 -DHAVE_SYS_TYPES_H=1
  -DHAVE_SYS_STAT_H=1 -DHAVE_STDLIB_H=1 -DHAVE_STRING_H=1 -DHAVE_MEMORY_H=1
  -DHAVE_STRINGS_H=1 -DHAVE_INTTYPES_H=1 -DHAVE_STDINT_H=1
  -DHAVE_UNISTD_H=1 -D__EXTENSIONS__=1 -D_ALL_SOURCE=1 -D_GNU_SOURCE=1
  -D_POSIX_PTHREAD_SEMANTICS=1 -D_TANDEM_SOURCE=1 -DSEPCHAR=':'
  -DSEPCHAR_STR=":" -D__NO_MATH_INLINES=1 -DCXX_NEW_FRIEND_TEMPLATE_DECL=1
  -DCXX_ISO_COMPLIANT_LIBRARY=1 -DHAVE_X_WINDOWS=1 -DHAVE_LIBM=1
  -DHAVE_QHULL=1 -DHAVE_PCRE_COMPILE=1 -DHAVE_PCRE=1 -DHAVE_REGEXEC=1
  -DHAVE_REGEX=1 -DHAVE_ZLIB_H=1 -DHAVE_ZLIB=1 -DHAVE_HDF5_H=1 -DHAVE_HDF5=1
  -DHAVE_H5GGET_NUM_OBJS=1 -DHAVE_FFTW3=1 -DHAVE_GLPK_H=1 -DHAVE_GLPK=1
  -DHAVE_CURL_CURL_H=1 -DHAVE_CURL=1 -DHAVE_MAGICK=1 -DHAVE_GL_GL_H=1
  -DHAVE_GL_GLU_H=1 -DHAVE_OPENGL=1 -DHAVE_FTGL_FTGL_H=1 -DHAVE_FTGL=1
  -DHAVE_FLTK=1 -DHAVE_IEEE754_DATA_FORMAT=1 -DF77_FUNC(name,NAME)=name
  ## _ -DF77_FUNC_(name,NAME)=name ## _ -DHAVE_BLAS=1 -DHAVE_QRUPDATE=1
  -DHAVE_SUITESPARSE_AMD_H=1 -DHAVE_AMD=1 -DHAVE_SUITESPARSE_UMFPACK_H=1
  -DHAVE_UMFPACK=1 -DUMFPACK_SEPARATE_SPLIT=1 -DHAVE_SUITESPARSE_COLAMD_H=1
  -DHAVE_COLAMD=1 -DHAVE_SUITESPARSE_CCOLAMD_H=1 -DHAVE_CCOLAMD=1
  -DHAVE_SUITESPARSE_CHOLMOD_H=1 -DHAVE_CHOLMOD=1 -DHAVE_SUITESPARSE_CS_H=1
  -DHAVE_CXSPARSE=1 -DHAVE_ARPACK=1 -DHAVE_GETHOSTNAME=1 -DHAVE_GETPWNAM=1
  -DHAVE_DEV_T=1 -DHAVE_INO_T=1 -DHAVE_NLINK_T=1 -DHAVE_NLINK_T=1
  -DHAVE_LONG_LONG_INT=1 -DHAVE_UNSIGNED_LONG_LONG_INT=1
  -DHAVE_SIGSET_T=1 -DHAVE_SIG_ATOMIC_T=1 -DSIZEOF_SHORT=2
  -DSIZEOF_INT=4 -DSIZEOF_LONG=8 -DSIZEOF_LONG_LONG=8 -DHAVE_ALLOCA_H=1
  -DHAVE_ALLOCA=1 -DHAVE_PLACEMENT_DELETE=1 -DHAVE_DYNAMIC_AUTO_ARRAYS=1
  -DHAVE_FAST_INT_OPS=1 -DSIZEOF_LONG_DOUBLE=16 -DSTDC_HEADERS=1
  -DHAVE_DIRENT_H=1 -DTIME_WITH_SYS_TIME=1 -DHAVE_SYS_WAIT_H=1
  -DHAVE_ASSERT_H=1 -DHAVE_CURSES_H=1 -DHAVE_DLFCN_H=1 -DHAVE_FCNTL_H=1
  -DHAVE_FLOAT_H=1 -DHAVE_GRP_H=1 -DHAVE_INTTYPES_H=1 -DHAVE_LIMITS_H=1
  -DHAVE_LOCALE_H=1 -DHAVE_MEMORY_H=1 -DHAVE_NCURSES_H=1 -DHAVE_POLL_H=1
  -DHAVE_PTHREAD_H=1 -DHAVE_PWD_H=1 -DHAVE_STDINT_H=1 -DHAVE_STDLIB_H=1
  -DHAVE_STRING_H=1 -DHAVE_SYS_IOCTL_H=1 -DHAVE_SYS_PARAM_H=1
  -DHAVE_SYS_POLL_H=1 -DHAVE_SYS_RESOURCE_H=1 -DHAVE_SYS_SELECT_H=1
  -DHAVE_SYS_STAT_H=1 -DHAVE_SYS_TIME_H=1 -DHAVE_SYS_TIMES_H=1
  -DHAVE_SYS_TYPES_H=1 -DHAVE_SYS_UTSNAME_H=1 -DHAVE_TERMCAP_H=1
  -DHAVE_UNISTD_H=1 -DHAVE_UTIME_H=1 -DHAVE_SSTREAM=1 -DHAVE_TERMIOS_H=1
  -DHAVE_TERMIO_H=1 -DHAVE_SGTTY_H=1 -DHAVE_GLOB_H=1 -DHAVE_FNMATCH_H=1
  -DHAVE_FNMATCH=1 -DHAVE_GLOB=1 -DHAVE_ATEXIT=1 -DHAVE_BASENAME=1
  -DHAVE_BCOPY=1 -DHAVE_BZERO=1 -DHAVE_CANONICALIZE_FILE_NAME=1
  -DHAVE_CHMOD=1 -DHAVE_DUP2=1 -DHAVE_ENDGRENT=1 -DHAVE_ENDPWENT=1
  -DHAVE_EXECVP=1 -DHAVE_EXPM1=1 -DHAVE_EXPM1F=1 -DHAVE_FCNTL=1
  -DHAVE_FORK=1 -DHAVE_FSTAT=1 -DHAVE_GETCWD=1 -DHAVE_GETEGID=1
  -DHAVE_GETEUID=1 -DHAVE_GETGID=1 -DHAVE_GETGRENT=1 -DHAVE_GETGRGID=1
  -DHAVE_GETGRNAM=1 -DHAVE_GETPGRP=1 -DHAVE_GETPID=1 -DHAVE_GETPPID=1
  -DHAVE_GETPWENT=1 -DHAVE_GETPWUID=1 -DHAVE_GETTIMEOFDAY=1 -DHAVE_GETUID=1
  -DHAVE_GETWD=1 -DHAVE_KILL=1 -DHAVE_LGAMMA=1 -DHAVE_LGAMMAF=1
  -DHAVE_LGAMMA_R=1 -DHAVE_LGAMMAF_R=1 -DHAVE_LINK=1 -DHAVE_LOCALTIME_R=1
  -DHAVE_LOG1P=1 -DHAVE_LOG1PF=1 -DHAVE_LSTAT=1 -DHAVE_MEMMOVE=1
  -DHAVE_MKDIR=1 -DHAVE_MKFIFO=1 -DHAVE_MKSTEMP=1 -DHAVE_ON_EXIT=1
  -DHAVE_PIPE=1 -DHAVE_POLL=1 -DHAVE_PUTENV=1 -DHAVE_RAISE=1
  -DHAVE_READLINK=1 -DHAVE_REALPATH=1 -DHAVE_RENAME=1 -DHAVE_RINDEX=1
  -DHAVE_RMDIR=1 -DHAVE_ROUNDL=1 -DHAVE_SELECT=1 -DHAVE_SETGRENT=1
  -DHAVE_SETLOCALE=1 -DHAVE_SETPWENT=1 -DHAVE_SETVBUF=1 -DHAVE_SIGACTION=1
  -DHAVE_SIGLONGJMP=1 -DHAVE_SIGPENDING=1 -DHAVE_SIGPROCMASK=1
  -DHAVE_SIGSUSPEND=1 -DHAVE_SNPRINTF=1 -DHAVE_STAT=1 -DHAVE_STRCASECMP=1
  -DHAVE_STRDUP=1 -DHAVE_STRERROR=1 -DHAVE_STRNCASECMP=1 -DHAVE_STRPTIME=1
  -DHAVE_STRSIGNAL=1 -DHAVE_SYMLINK=1 -DHAVE_TEMPNAM=1 -DHAVE_TGAMMAF=1
  -DHAVE_TRUNC=1 -DHAVE_UMASK=1 -DHAVE_UNAME=1 -DHAVE_UNLINK=1
  -DHAVE_USLEEP=1 -DHAVE_UTIME=1 -DHAVE_VFPRINTF=1 -DHAVE_VSPRINTF=1
  -DHAVE_VSNPRINTF=1 -DHAVE_WAITPID=1 -DHAVE_DECL_EXP2=1 -DHAVE_DECL_ROUND=1
  -DHAVE_DECL_TGAMMA=1 -DHAVE_EXP2=1 -DHAVE_ROUND=1 -DHAVE_TGAMMA=1
  -DHAVE_STRFTIME=1 -DHAVE_C99_VSNPRINTF=1 -DHAVE_LIBDL=1 -DHAVE_DLOPEN=1
  -DHAVE_DLSYM=1 -DHAVE_DLERROR=1 -DHAVE_DLCLOSE=1 -DHAVE_DLOPEN_API=1
  -DENABLE_DYNAMIC_LINKING=1 -DHAVE_TIMEVAL=1 -DHAVE_CMATH_ISNAN=1
  -DHAVE_CMATH_ISNANF=1 -DHAVE_CMATH_ISINF=1 -DHAVE_CMATH_ISINFF=1
  -DHAVE_CMATH_ISFINITE=1 -DHAVE_CMATH_ISFINITEF=1 -DHAVE_FINITE=1
  -DHAVE_ISNAN=1 -DHAVE_ISINF=1 -DHAVE_COPYSIGN=1 -DHAVE_DECL_SIGNBIT=1
  -DHAVE_ACOSH=1 -DHAVE_ACOSHF=1 -DHAVE_ASINH=1 -DHAVE_ASINHF=1
  -DHAVE_ATANH=1 -DHAVE_ATANHF=1 -DHAVE_ERF=1 -DHAVE_ERFF=1 -DHAVE_ERFC=1
  -DHAVE_ERFCF=1 -DHAVE_EXP2F=1 -DHAVE_LOG2=1 -DHAVE_LOG2F=1 -DHAVE_HYPOTF=1
  -DHAVE_STRUCT_STAT_ST_BLKSIZE=1 -DHAVE_STRUCT_STAT_ST_BLOCKS=1
  -DHAVE_STRUCT_STAT_ST_RDEV=1 -DHAVE_STRUCT_TM_TM_ZONE=1
  -DHAVE_TM_ZONE=1 -DUSE_READLINE=1 -DEXCEPTION_IN_MATH=1 -DRETSIGTYPE=void
  -DHAVE_DECL_SYS_SIGLIST=1 -DHAVE_POSIX_SIGNALS=1 -DRETSIGTYPE_IS_VOID=1
  -DHAVE_GETRUSAGE=1 -DHAVE_TIMES=1 -DYYTEXT_POINTER=1

User-preferences (please do not edit this section):
--------------------------------------------------

  EDITOR = vim
  EXEC_PATH = /usr/lib/octave/3.2.3/site/exec/x86_64-pc-linux-gnu:/usr/lib/octave/api-v37/site/exec/x86_64-pc-linux-gnu:/usr/lib/octave/site/exec/x86_64-pc-linux-gnu:/usr/lib/octave/3.2.3/exec/x86_64-pc-linux-gnu:/usr/bin:/usr/local/bin:/usr/bin:/bin:/usr/games
  IMAGE_PATH = .:/usr/share/octave/3.2.3/imagelib
  PAGER = pager
  PS1 = \s:\#>
  PS2 = >
  PS4 = +
  beep_on_error = 0
  completion_append_char =  
  crash_dumps_octave_core = 1
  echo_executing_commands = 0
  fixed_point_format = 0
  gnuplot_binary = gnuplot
# gnuplot_command_end = <no value or error in displaying it>
# gnuplot_command_plot = <no value or error in displaying it>
# gnuplot_command_replot = <no value or error in displaying it>
# gnuplot_command_splot = <no value or error in displaying it>
# gnuplot_command_title = <no value or error in displaying it>
# gnuplot_command_using = <no value or error in displaying it>
# gnuplot_command_with = <no value or error in displaying it>
  history_file = /home/chrysn/.octave_hist
  history_size = 1024
  ignore_function_time_stamp = system
  info_file = /usr/share/info/octave3.2.info
  info_program = info
  makeinfo_program = makeinfo
  max_recursion_depth = 256
  output_max_field_width = 5
  output_precision = 5
  page_output_immediately = 0
  page_screen_output = 1
# print_answer_id_name = <no value or error in displaying it>
  print_empty_dimensions = 1
  save_precision = 16
  saving_history = 1
  sighup_dumps_octave_core = 1
  sigterm_dumps_octave_core = 1
  silent_functions = 0
  split_long_rows = 1
  string_fill_char =  
  struct_levels_to_print = 2
  suppress_verbose_help_message = 0

# Created by Octave 3.2.3, Wed Oct 28 11:36:58 2009 CET <chrysn@hephaistos>
# name: bad
# type: matrix
# rows: 10
# columns: 10
 10.08281071474692 -6.365001220213275e-16 -1.095862284704128e-15 -2.696039612100385e-16 2.028106670699041e-16 7.87609959122284e-16 -8.69535116950435e-16 -4.752957577688955e-16 -7.999882714801519e-16 -9.526555401654391e-16
 3.072352834040245e-154 2.409834985187796 -7.080741419290825e-16 4.264640455734195e-16 -8.528245961335121e-16 1.662091446003039e-15 3.020979986170891e-16 -8.612567565067727e-16 -2.320203839266881e-16 4.274488504043562e-16
 1.284045590470485e-180 -6.001692936490212e-28 1.884356937399034 0.007710938642175293 -4.104785221720314e-16 -3.534126315668841e-17 1.729655694077542e-16 -3.49090542331505e-16 2.860087248915117e-17 2.068312017813794e-16
 1.990805540954218e-183 5.494494711313992e-30 0.007710938642175732 -1.831725592155371 8.958674039903884e-16 4.912516314330228e-16 1.729818694751086e-16 -3.080955978428914e-16 1.549335989744458e-16 2.953198444530853e-16
 2.689116890988173e-201 1.558664161541373e-47 4.215846637899941e-21 -2.127226215234033e-18 1.554273726563368 -6.572292182197986e-10 1.720608741233552e-17 -3.922726689584524e-16 6.419004949219602e-17 2.442336539058333e-16
 -7.597471931089994e-211 1.549041522751045e-57 -2.795771717229085e-31 4.89237461437074e-28 -6.572277398665318e-10 -1.421477099991831 -1.150360111462768e-16 2.244076894005266e-16 1.495432618687405e-16 3.055374795750419e-16
 1.962665636287987e-269 -3.018824720713485e-116 4.919718380222504e-89 -2.868212539135241e-86 1.124954879426486e-68 3.46216949595653e-59 0.8243629511918087 5.891479416258602e-12 4.214620209217503e-17 6.119273197461773e-17
 3.813852223721426e-281 1.000636734528609e-127 1.403004472579291e-100 -8.602656262752663e-98 1.599214745002984e-80 1.017409509385266e-70 5.890935017335692e-12 -0.7376758950027112 1.415418436186652e-16 -9.149250694666573e-17
 0 5.649230126363334e-175 5.167601812435397e-149 -5.991195123604499e-146 -2.604808135773483e-128 1.057036480301348e-118 1.972820712705903e-59 -4.574201222769327e-48 0.4756823873383532 3.113929142453984e-17
 0 6.208957057064851e-307 3.933340142416565e-281 -6.138020110399418e-278 -2.714866147239894e-260 1.129805731278581e-250 1.998330648384745e-191 -4.53409391582144e-180 5.171745669229912e-133 -0.1383051152355606


_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave

[John W. Eaton-3]

On  3-Nov-2009, chrysn wrote:

| (i hope this doesn't get filed twice; i reported it from within octave,
| but it seems to have missed the list)
|
| --------
| Bug report for Octave 3.2.3 configured for x86_64-pc-linux-gnu
|
| Description:
| -----------
|
|   * for the attached matrix "bad", the qr decomposition creates NaN and Inf
|
| Repeat-By:
| ---------
|
|   * load "bad"
|   * do `[Q,R] = qr(bad)` or `[Q,R,P] = qr(bad)`; you get several +/-Inf and
|     NaN.
|
| (for comparison, matlab can decompose this; moreover, the qr decomposition
| should work with any given input)

Octave just uses DGEQRF from LAPACK for computing the QR decomposition
in this case.  So unless we are calling that function incorrectly
(possible, but it seems unlikely, since it works for many other
cases...), I think the problem is with LAPACK.  I'm attaching a simple
program that demonstrates the problem with your data.  It makes the
same call to DGEQRF that Octave does.

If I use the LAPACK 3.2.1 library that is installed on my Debian
system, I see Inf and NaN values in the result.

If I use the older reference LAPACK 3.1 sources that are distributed
with Octave, I see the more reasonable result:

  required workspace =   320 we are using 560
   info =            0
   a =
     -0.10E+02    0.64E-15    0.11E-14    0.27E-15   -0.20E-15   -0.79E-15    0.87E-15    0.48E-15    0.80E-15    0.95E-15
      0.15-154   -0.24E+01    0.71E-15   -0.43E-15    0.85E-15   -0.17E-14   -0.30E-15    0.86E-15    0.23E-15   -0.43E-15
      0.64-181   -0.12E-27   -0.19E+01   -0.22E-03    0.41E-15    0.33E-16   -0.17E-15    0.35E-15   -0.29E-16   -0.21E-15
      0.99-184    0.11E-29    0.20E-02    0.18E+01   -0.90E-15   -0.49E-15   -0.17E-15    0.31E-15   -0.15E-15   -0.29E-15
      0.13-201    0.32E-47    0.11E-20    0.58E-18   -0.16E+01    0.56E-10   -0.17E-16    0.39E-15   -0.64E-16   -0.24E-15
     -0.38-211    0.32E-57   -0.74E-31   -0.13E-27   -0.21E-09    0.14E+01    0.12E-15   -0.22E-15   -0.15E-15   -0.31E-15
      0.97-270   -0.63-116    0.13E-88    0.78E-86    0.36E-68   -0.12E-58   -0.82E+00   -0.62E-12   -0.42E-16   -0.61E-16
      0.19-281    0.21-127    0.37-100    0.23E-97    0.51E-80   -0.36E-70    0.36E-11    0.74E+00   -0.14E-15    0.91E-16
      0.00E+00    0.12-174    0.14-148    0.16-145   -0.84-128   -0.37-118    0.12E-58    0.31E-47   -0.48E+00   -0.31E-16
      0.00E+00    0.13-306    0.10-280    0.17-277   -0.87-260   -0.40-250    0.12-190    0.31-179    0.54-132   -0.14E+00
   tau =
      0.20E+01
      0.20E+01
      0.20E+01
      0.20E+01
      0.20E+01
      0.20E+01
      0.20E+01
      0.20E+01
      0.20E+01
      0.00E+00


jwe


      program foo

      double precision a(10,10)
      double precision work(560)
      double precision tau(10)
      integer info

      a(1,1) = 10.08281071474692
      a(1,2) = -6.365001220213275D-16
      a(1,3) = -1.095862284704128D-15
      a(1,4) = -2.696039612100385D-16
      a(1,5) = 2.028106670699041D-16
      a(1,6) = 7.87609959122284D-16
      a(1,7) = -8.69535116950435D-16
      a(1,8) = -4.752957577688955D-16
      a(1,9) = -7.999882714801519D-16
      a(1,10) = -9.526555401654391D-16

      a(2,1) = 3.072352834040245D-154
      a(2,2) = 2.409834985187796
      a(2,3) = -7.080741419290825D-16
      a(2,4) = 4.264640455734195D-16
      a(2,5) = -8.528245961335121D-16
      a(2,6) = 1.662091446003039D-15
      a(2,7) = 3.020979986170891D-16
      a(2,8) = -8.612567565067727D-16
      a(2,9) = -2.320203839266881D-16
      a(2,10) = 4.274488504043562D-16

      a(3,1) = 1.284045590470485D-180
      a(3,2) = -6.001692936490212D-28
      a(3,3) = 1.884356937399034
      a(3,4) = 0.007710938642175293
      a(3,5) = -4.104785221720314D-16
      a(3,6) = -3.534126315668841D-17
      a(3,7) = 1.729655694077542D-16
      a(3,8) = -3.49090542331505D-16
      a(3,9) = 2.860087248915117D-17
      a(3,10) = 2.068312017813794D-16

      a(4,1) = 1.990805540954218D-183
      a(4,2) = 5.494494711313992D-30
      a(4,3) = 0.007710938642175732
      a(4,4) = -1.831725592155371
      a(4,5) = 8.958674039903884D-16
      a(4,6) = 4.912516314330228D-16
      a(4,7) = 1.729818694751086D-16
      a(4,8) = -3.080955978428914D-16
      a(4,9) = 1.549335989744458D-16
      a(4,10) = 2.953198444530853D-16

      a(5,1) = 2.689116890988173D-201
      a(5,2) = 1.558664161541373D-47
      a(5,3) = 4.215846637899941D-21
      a(5,4) = -2.127226215234033D-18
      a(5,5) = 1.554273726563368
      a(5,6) = -6.572292182197986D-10
      a(5,7) = 1.720608741233552D-17
      a(5,8) = -3.922726689584524D-16
      a(5,9) = 6.419004949219602D-17
      a(5,10) = 2.442336539058333D-16

      a(6,1) = -7.597471931089994D-211
      a(6,2) = 1.549041522751045D-57
      a(6,3) = -2.795771717229085D-31
      a(6,4) = 4.89237461437074D-28
      a(6,5) = -6.572277398665318D-10
      a(6,6) = -1.421477099991831
      a(6,7) = -1.150360111462768D-16
      a(6,8) = 2.244076894005266D-16
      a(6,9) = 1.495432618687405D-16
      a(6,10) = 3.055374795750419D-16

      a(7,1) = 1.962665636287987D-269
      a(7,2) = -3.018824720713485D-116
      a(7,3) = 4.919718380222504D-89
      a(7,4) = -2.868212539135241D-86
      a(7,5) = 1.124954879426486D-68
      a(7,6) = 3.46216949595653D-59
      a(7,7) = 0.8243629511918087
      a(7,8) = 5.891479416258602D-12
      a(7,9) = 4.214620209217503D-17
      a(7,10) = 6.119273197461773D-17

      a(8,1) = 3.813852223721426D-281
      a(8,2) = 1.000636734528609D-127
      a(8,3) = 1.403004472579291D-100
      a(8,4) = -8.602656262752663D-98
      a(8,5) = 1.599214745002984D-80
      a(8,6) = 1.017409509385266D-70
      a(8,7) = 5.890935017335692D-12
      a(8,8) = -0.7376758950027112
      a(8,9) = 1.415418436186652D-16
      a(8,10) = -9.149250694666573D-17

      a(9,1) = 0
      a(9,2) = 5.649230126363334D-175
      a(9,3) = 5.167601812435397D-149
      a(9,4) = -5.991195123604499D-146
      a(9,5) = -2.604808135773483D-128
      a(9,6) = 1.057036480301348D-118
      a(9,7) = 1.972820712705903D-59
      a(9,8) = -4.574201222769327D-48
      a(9,9) = 0.4756823873383532
      a(9,10) = 3.113929142453984D-17

      a(10,1) = 0
      a(10,2) = 6.208957057064851D-307
      a(10,3) = 3.933340142416565D-281
      a(10,4) = -6.138020110399418D-278
      a(10,5) = -2.714866147239894D-260
      a(10,6) = 1.129805731278581D-250
      a(10,7) = 1.998330648384745D-191
      a(10,8) = -4.53409391582144D-180
      a(10,9) = 5.171745669229912D-133
      a(10,10) = -0.1383051152355606

* query workspace requiremnt.      
      call dgeqrf (10, 10, a, 10, tau, work, -1, info)
      write (*,'(a,i5,a)')
     $     'required workspace = ', int (work(1)), ' we are using 560'

* compute decomposition
      call dgeqrf (10, 10, a, 10, tau, work, 560, info)

      print *, 'info = ', info
      print *, 'a = '
      do i = 1, 10
        write(*,'(10e12.2)') (a(i,j), j = 1, 10)
      enddo
      print *, 'tau = '
      do i = 1, 10
        write(*,'(e12.2)') tau(i)
      enddo

      end

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave

[Jaroslav Hajek-2]

On Tue, Nov 3, 2009 at 11:08 PM, John W. Eaton <jwe@...> wrote:

> On  3-Nov-2009, chrysn wrote:
>
> | (i hope this doesn't get filed twice; i reported it from within octave,
> | but it seems to have missed the list)
> |
> | --------
> | Bug report for Octave 3.2.3 configured for x86_64-pc-linux-gnu
> |
> | Description:
> | -----------
> |
> |   * for the attached matrix "bad", the qr decomposition creates NaN and Inf
> |
> | Repeat-By:
> | ---------
> |
> |   * load "bad"
> |   * do `[Q,R] = qr(bad)` or `[Q,R,P] = qr(bad)`; you get several +/-Inf and
> |     NaN.
> |
> | (for comparison, matlab can decompose this; moreover, the qr decomposition
> | should work with any given input)
>
> Octave just uses DGEQRF from LAPACK for computing the QR decomposition
> in this case.  So unless we are calling that function incorrectly
> (possible, but it seems unlikely, since it works for many other
> cases...), I think the problem is with LAPACK.  I'm attaching a simple
> program that demonstrates the problem with your data.  It makes the
> same call to DGEQRF that Octave does.
>
> If I use the LAPACK 3.2.1 library that is installed on my Debian
> system, I see Inf and NaN values in the result.
>
> If I use the older reference LAPACK 3.1 sources that are distributed
> with Octave, I see the more reasonable result:
>
>  required workspace =   320 we are using 560
>   info =            0
>   a =
>     -0.10E+02    0.64E-15    0.11E-14    0.27E-15   -0.20E-15   -0.79E-15    0.87E-15    0.48E-15    0.80E-15    0.95E-15
>      0.15-154   -0.24E+01    0.71E-15   -0.43E-15    0.85E-15   -0.17E-14   -0.30E-15    0.86E-15    0.23E-15   -0.43E-15
>      0.64-181   -0.12E-27   -0.19E+01   -0.22E-03    0.41E-15    0.33E-16   -0.17E-15    0.35E-15   -0.29E-16   -0.21E-15
>      0.99-184    0.11E-29    0.20E-02    0.18E+01   -0.90E-15   -0.49E-15   -0.17E-15    0.31E-15   -0.15E-15   -0.29E-15
>      0.13-201    0.32E-47    0.11E-20    0.58E-18   -0.16E+01    0.56E-10   -0.17E-16    0.39E-15   -0.64E-16   -0.24E-15
>     -0.38-211    0.32E-57   -0.74E-31   -0.13E-27   -0.21E-09    0.14E+01    0.12E-15   -0.22E-15   -0.15E-15   -0.31E-15
>      0.97-270   -0.63-116    0.13E-88    0.78E-86    0.36E-68   -0.12E-58   -0.82E+00   -0.62E-12   -0.42E-16   -0.61E-16
>      0.19-281    0.21-127    0.37-100    0.23E-97    0.51E-80   -0.36E-70    0.36E-11    0.74E+00   -0.14E-15    0.91E-16
>      0.00E+00    0.12-174    0.14-148    0.16-145   -0.84-128   -0.37-118    0.12E-58    0.31E-47   -0.48E+00   -0.31E-16
>      0.00E+00    0.13-306    0.10-280    0.17-277   -0.87-260   -0.40-250    0.12-190    0.31-179    0.54-132   -0.14E+00
>   tau =
>      0.20E+01
>      0.20E+01
>      0.20E+01
>      0.20E+01
>      0.20E+01
>      0.20E+01
>      0.20E+01
>      0.20E+01
>      0.20E+01
>      0.00E+00
>
>
> jwe
>
>
>      program foo
>
>      double precision a(10,10)
>      double precision work(560)
>      double precision tau(10)
>      integer info
>
>      a(1,1) = 10.08281071474692
>      a(1,2) = -6.365001220213275D-16
>      a(1,3) = -1.095862284704128D-15
>      a(1,4) = -2.696039612100385D-16
>      a(1,5) = 2.028106670699041D-16
>      a(1,6) = 7.87609959122284D-16
>      a(1,7) = -8.69535116950435D-16
>      a(1,8) = -4.752957577688955D-16
>      a(1,9) = -7.999882714801519D-16
>      a(1,10) = -9.526555401654391D-16
>
>      a(2,1) = 3.072352834040245D-154
>      a(2,2) = 2.409834985187796
>      a(2,3) = -7.080741419290825D-16
>      a(2,4) = 4.264640455734195D-16
>      a(2,5) = -8.528245961335121D-16
>      a(2,6) = 1.662091446003039D-15
>      a(2,7) = 3.020979986170891D-16
>      a(2,8) = -8.612567565067727D-16
>      a(2,9) = -2.320203839266881D-16
>      a(2,10) = 4.274488504043562D-16
>
>      a(3,1) = 1.284045590470485D-180
>      a(3,2) = -6.001692936490212D-28
>      a(3,3) = 1.884356937399034
>      a(3,4) = 0.007710938642175293
>      a(3,5) = -4.104785221720314D-16
>      a(3,6) = -3.534126315668841D-17
>      a(3,7) = 1.729655694077542D-16
>      a(3,8) = -3.49090542331505D-16
>      a(3,9) = 2.860087248915117D-17
>      a(3,10) = 2.068312017813794D-16
>
>      a(4,1) = 1.990805540954218D-183
>      a(4,2) = 5.494494711313992D-30
>      a(4,3) = 0.007710938642175732
>      a(4,4) = -1.831725592155371
>      a(4,5) = 8.958674039903884D-16
>      a(4,6) = 4.912516314330228D-16
>      a(4,7) = 1.729818694751086D-16
>      a(4,8) = -3.080955978428914D-16
>      a(4,9) = 1.549335989744458D-16
>      a(4,10) = 2.953198444530853D-16
>
>      a(5,1) = 2.689116890988173D-201
>      a(5,2) = 1.558664161541373D-47
>      a(5,3) = 4.215846637899941D-21
>      a(5,4) = -2.127226215234033D-18
>      a(5,5) = 1.554273726563368
>      a(5,6) = -6.572292182197986D-10
>      a(5,7) = 1.720608741233552D-17
>      a(5,8) = -3.922726689584524D-16
>      a(5,9) = 6.419004949219602D-17
>      a(5,10) = 2.442336539058333D-16
>
>      a(6,1) = -7.597471931089994D-211
>      a(6,2) = 1.549041522751045D-57
>      a(6,3) = -2.795771717229085D-31
>      a(6,4) = 4.89237461437074D-28
>      a(6,5) = -6.572277398665318D-10
>      a(6,6) = -1.421477099991831
>      a(6,7) = -1.150360111462768D-16
>      a(6,8) = 2.244076894005266D-16
>      a(6,9) = 1.495432618687405D-16
>      a(6,10) = 3.055374795750419D-16
>
>      a(7,1) = 1.962665636287987D-269
>      a(7,2) = -3.018824720713485D-116
>      a(7,3) = 4.919718380222504D-89
>      a(7,4) = -2.868212539135241D-86
>      a(7,5) = 1.124954879426486D-68
>      a(7,6) = 3.46216949595653D-59
>      a(7,7) = 0.8243629511918087
>      a(7,8) = 5.891479416258602D-12
>      a(7,9) = 4.214620209217503D-17
>      a(7,10) = 6.119273197461773D-17
>
>      a(8,1) = 3.813852223721426D-281
>      a(8,2) = 1.000636734528609D-127
>      a(8,3) = 1.403004472579291D-100
>      a(8,4) = -8.602656262752663D-98
>      a(8,5) = 1.599214745002984D-80
>      a(8,6) = 1.017409509385266D-70
>      a(8,7) = 5.890935017335692D-12
>      a(8,8) = -0.7376758950027112
>      a(8,9) = 1.415418436186652D-16
>      a(8,10) = -9.149250694666573D-17
>
>      a(9,1) = 0
>      a(9,2) = 5.649230126363334D-175
>      a(9,3) = 5.167601812435397D-149
>      a(9,4) = -5.991195123604499D-146
>      a(9,5) = -2.604808135773483D-128
>      a(9,6) = 1.057036480301348D-118
>      a(9,7) = 1.972820712705903D-59
>      a(9,8) = -4.574201222769327D-48
>      a(9,9) = 0.4756823873383532
>      a(9,10) = 3.113929142453984D-17
>
>      a(10,1) = 0
>      a(10,2) = 6.208957057064851D-307
>      a(10,3) = 3.933340142416565D-281
>      a(10,4) = -6.138020110399418D-278
>      a(10,5) = -2.714866147239894D-260
>      a(10,6) = 1.129805731278581D-250
>      a(10,7) = 1.998330648384745D-191
>      a(10,8) = -4.53409391582144D-180
>      a(10,9) = 5.171745669229912D-133
>      a(10,10) = -0.1383051152355606
>
> * query workspace requiremnt.
>      call dgeqrf (10, 10, a, 10, tau, work, -1, info)
>      write (*,'(a,i5,a)')
>     $     'required workspace = ', int (work(1)), ' we are using 560'
>
> * compute decomposition
>      call dgeqrf (10, 10, a, 10, tau, work, 560, info)
>
>      print *, 'info = ', info
>      print *, 'a = '
>      do i = 1, 10
>        write(*,'(10e12.2)') (a(i,j), j = 1, 10)
>      enddo
>      print *, 'tau = '
>      do i = 1, 10
>        write(*,'(e12.2)') tau(i)
>      enddo
>
>      end
>
> _______________________________________________
> Bug-octave mailing list
> Bug-octave@...
> https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave
>
>

I spent some time investigating this. Apparently the new xLARFP
routines are responsible for these problems.
These are described in the LAPACK working note 203:
http://www.netlib.org/lapack/lawnspdf/lawn203.pdf

Specifically, on page 4, the text says:

"Both imag(alpha)/gamma and norm(x)/gamma are less than one, so
abs(delta) <= abs(beta) and these computations introduce
no new overfl
ow possibilities into the existing routines."

It seems to me that this is utterly wrong. Apparently delta = norm(x)
* (norm(x) / beta) can't overflow, but it can easily underflow if
norm(x) is a tiny number; hence subsequent division asks for problems.
A possible remedy is to split into two scalings, one by norm(x) and
one by norm(x)/beta - still, even norm(x)/beta alone can underflow.

It seems to me that the goal of achieving beta >= 0 together with
keeping the leading one in the reflection vector is too ambitious.

Something like the following patch fixes the worst of the problem, but
still doesn't solve everything. Just raise a(:,1) to second power in
your test matrix and you're back in XXX. Hmmm...

diff --git a/dlarfp.f b/dlarfp.f
--- a/dlarfp.f
+++ b/dlarfp.f
@@ -134,8 +134,13 @@
             BETA = -BETA
             TAU = -ALPHA / BETA
          ELSE
-            ALPHA = XNORM * (XNORM/ALPHA)
-            TAU = ALPHA / BETA
+            ALPHA = XNORM/ALPHA
+            TAU = XNORM * ALPHA
+            IF (TAU >= SAFMIN) THEN
+               ALPHA = TAU
+            ELSE
+               CALL DSCAL ( N-1, ONE / XNORM, X, INCX )
+            END IF
             ALPHA = -ALPHA
          END IF
          CALL DSCAL( N-1, ONE / ALPHA, X, INCX )
diff --git a/slarfp.f b/slarfp.f
--- a/slarfp.f
+++ b/slarfp.f
@@ -134,8 +134,13 @@
             BETA = -BETA
             TAU = -ALPHA / BETA
          ELSE
-            ALPHA = XNORM * (XNORM/ALPHA)
-            TAU = ALPHA / BETA
+            ALPHA = XNORM/ALPHA
+            TAU = XNORM * ALPHA
+            IF (TAU >= SAFMIN) THEN
+               ALPHA = TAU
+            ELSE
+               CALL SSCAL ( N-1, ONE / XNORM, X, INCX )
+            END IF
             ALPHA = -ALPHA
          END IF
          CALL SSCAL( N-1, ONE / ALPHA, X, INCX )

--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave

[John W. Eaton-3]

On  4-Nov-2009, Jaroslav Hajek wrote:

| I spent some time investigating this. Apparently the new xLARFP
| routines are responsible for these problems.
| These are described in the LAPACK working note 203:
| http://www.netlib.org/lapack/lawnspdf/lawn203.pdf
|
| Specifically, on page 4, the text says:
|
| "Both imag(alpha)/gamma and norm(x)/gamma are less than one, so
| abs(delta) <= abs(beta) and these computations introduce
| no new overfl
| ow possibilities into the existing routines."
|
| It seems to me that this is utterly wrong. Apparently delta = norm(x)
| * (norm(x) / beta) can't overflow, but it can easily underflow if
| norm(x) is a tiny number; hence subsequent division asks for problems.
| A possible remedy is to split into two scalings, one by norm(x) and
| one by norm(x)/beta - still, even norm(x)/beta alone can underflow.
|
| It seems to me that the goal of achieving beta >= 0 together with
| keeping the leading one in the reflection vector is too ambitious.
|
| Something like the following patch fixes the worst of the problem, but
| still doesn't solve everything. Just raise a(:,1) to second power in
| your test matrix and you're back in XXX. Hmmm...

Is this patch for the LAPACK 3.2.x sources?

jwe
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave

[Jaroslav Hajek-2]

On Wed, Nov 4, 2009 at 5:09 PM, John W. Eaton <jwe@...> wrote:

> On  4-Nov-2009, Jaroslav Hajek wrote:
>
> | I spent some time investigating this. Apparently the new xLARFP
> | routines are responsible for these problems.
> | These are described in the LAPACK working note 203:
> | http://www.netlib.org/lapack/lawnspdf/lawn203.pdf
> |
> | Specifically, on page 4, the text says:
> |
> | "Both imag(alpha)/gamma and norm(x)/gamma are less than one, so
> | abs(delta) <= abs(beta) and these computations introduce
> | no new overfl
> | ow possibilities into the existing routines."
> |
> | It seems to me that this is utterly wrong. Apparently delta = norm(x)
> | * (norm(x) / beta) can't overflow, but it can easily underflow if
> | norm(x) is a tiny number; hence subsequent division asks for problems.
> | A possible remedy is to split into two scalings, one by norm(x) and
> | one by norm(x)/beta - still, even norm(x)/beta alone can underflow.
> |
> | It seems to me that the goal of achieving beta >= 0 together with
> | keeping the leading one in the reflection vector is too ambitious.
> |
> | Something like the following patch fixes the worst of the problem, but
> | still doesn't solve everything. Just raise a(:,1) to second power in
> | your test matrix and you're back in XXX. Hmmm...
>
> Is this patch for the LAPACK 3.2.x sources?
>
>

Yes. But it's not correct, so don't use it. I'll post a correction later.


--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave

[Marco atzeri-3]

--- Mer 4/11/09, Jaroslav Hajek <highegg@...> ha scritto:

[cut]

> I spent some time investigating this. Apparently the new
> xLARFP
> routines are responsible for these problems.
> These are described in the LAPACK working note 203:
> http://www.netlib.org/lapack/lawnspdf/lawn203.pdf
>
> Specifically, on page 4, the text says:
>
> "Both imag(alpha)/gamma and norm(x)/gamma are less than
> one, so
> abs(delta) <= abs(beta) and these computations
> introduce
> no new overfl
> ow possibilities into the existing routines."
>
> It seems to me that this is utterly wrong. Apparently delta
> = norm(x)
> * (norm(x) / beta) can't overflow, but it can easily
> underflow if
> norm(x) is a tiny number; hence subsequent division asks
> for problems.
> A possible remedy is to split into two scalings, one by
> norm(x) and
> one by norm(x)/beta - still, even norm(x)/beta alone can
> underflow.
>
> It seems to me that the goal of achieving beta >= 0
> together with
> keeping the leading one in the reflection vector is too
> ambitious.

looks like (*) bug0037 :: xLARFP and scaling
http://www.netlib.org/lapack/Errata/errata_3.2.1_20091001.html
or not ?

[cut]
>
> --
> RNDr. Jaroslav Hajek

Regards
Marco



     

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave

[Jaroslav Hajek-2]

On Wed, Nov 4, 2009 at 9:24 PM, Marco Atzeri <marco_atzeri@...> wrote:

>
>
> --- Mer 4/11/09, Jaroslav Hajek <highegg@...> ha scritto:
>
> [cut]
>
>> I spent some time investigating this. Apparently the new
>> xLARFP
>> routines are responsible for these problems.
>> These are described in the LAPACK working note 203:
>> http://www.netlib.org/lapack/lawnspdf/lawn203.pdf
>>
>> Specifically, on page 4, the text says:
>>
>> "Both imag(alpha)/gamma and norm(x)/gamma are less than
>> one, so
>> abs(delta) <= abs(beta) and these computations
>> introduce
>> no new overfl
>> ow possibilities into the existing routines."
>>
>> It seems to me that this is utterly wrong. Apparently delta
>> = norm(x)
>> * (norm(x) / beta) can't overflow, but it can easily
>> underflow if
>> norm(x) is a tiny number; hence subsequent division asks
>> for problems.
>> A possible remedy is to split into two scalings, one by
>> norm(x) and
>> one by norm(x)/beta - still, even norm(x)/beta alone can
>> underflow.
>>
>> It seems to me that the goal of achieving beta >= 0
>> together with
>> keeping the leading one in the reflection vector is too
>> ambitious.
>
> looks like (*) bug0037 :: xLARFP and scaling
> http://www.netlib.org/lapack/Errata/errata_3.2.1_20091001.html
> or not ?
>

Probably yes. It's marked as fixed in 3.2.2, so no need for me to try
(at least until I see their solution).

For anyone interested, here's a patch that simply forwards xLARFP to
the old xLARFG that can be used to patch 3.2.1.
Of course this brings back the old behavior of qr not yielding
nonnegative diagonals.

diff --git a/clarfp.f b/clarfp.f
--- a/clarfp.f
+++ b/clarfp.f
@@ -80,6 +80,8 @@
 *     ..
 *     .. Executable Statements ..
 *
+      CALL CLARFG( N, ALPHA, X, INCX, TAU )
+
       IF( N.LE.0 ) THEN
          TAU = ZERO
          RETURN
diff --git a/dlarfp.f b/dlarfp.f
--- a/dlarfp.f
+++ b/dlarfp.f
@@ -79,6 +79,8 @@
 *     ..
 *     .. Executable Statements ..
 *
+      CALL DLARFG( N, ALPHA, X, INCX, TAU )
+
       IF( N.LE.0 ) THEN
          TAU = ZERO
          RETURN
diff --git a/slarfp.f b/slarfp.f
--- a/slarfp.f
+++ b/slarfp.f
@@ -79,6 +79,8 @@
 *     ..
 *     .. Executable Statements ..
 *
+      CALL SLARFG( N, ALPHA, X, INCX, TAU )
+
       IF( N.LE.0 ) THEN
          TAU = ZERO
          RETURN
diff --git a/zlarfp.f b/zlarfp.f
--- a/zlarfp.f
+++ b/zlarfp.f
@@ -80,6 +80,8 @@
 *     ..
 *     .. Executable Statements ..
 *
+      CALL ZLARFG( N, ALPHA, X, INCX, TAU )
+
       IF( N.LE.0 ) THEN
          TAU = ZERO
          RETURN




--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave

[Marco atzeri-3]

--- Ven 6/11/09, Jaroslav Hajek <highegg@...> ha scritto:

> >
> > looks like (*) bug0037 :: xLARFP and scaling
> > http://www.netlib.org/lapack/Errata/errata_3.2.1_20091001.html
> > or not ?
> >
>
> Probably yes. It's marked as fixed in 3.2.2, so no need for
> me to try
> (at least until I see their solution).

in reality it is reported as

"Bug reports -- bug needs to be confirmed!"

so no solution is yet available.
So you are free to offer yours

>
> --
> RNDr. Jaroslav Hajek

Marco



     

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave