Complex left-divide and SVD give ireproducible results

Dirk Laurie dpl at sun.ac.za
Mon Jul 28 08:14:25 CDT 2008


--------
Bug report for Octave 2.1.73 configured for i486-pc-linux-gnu    

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

The following operations give irreproducible results:
    z=A\b  when A is a nonsquare complex matrix
    [U,S,V]=svd(A) when A is a complex matrix

By 'irreproducible' is meant that the results are not identical
to the last bit when the operation is run more than once on
identical input, although the results in both cases are numerically
correct in the sense that b-A*z and A-U*S*V respectively are at
roundoff error level.  The two bugs are reported as one because 
it is thought probable that they share a common cause.

A possible cause might be that an uninitialized value is used in 
a comparison of which the result would normally not critically 
affect the output, e.g. whether columns in a full-rank orthogonal 
factorization are to be interchanged. 


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

The left-divide operation is so capricious that it suffices to 
generate a random test case.  The following code exposes the bug:
 
A=rand(200,40)+1i*rand(200,40); b=rand(200,1); x0=A\b;
for iter=1:10, x=A\b; err(iter)=norm(x0-x,Inf); end
err

The otput should be all zeros.  If I save the above code in an M-file,
I get nonzero results most of the time, whether running the code in 
an interactive shell or from the command line.

The SVD operation can also be exposed by a random test case, 
albeit one in which the matrix is scaled down, as follows:

A=eps*(rand(40,8)+1i*rand(40,8)); [U0,S0,V0]=svd(A);
for iter=1:10,
    [U,S,V]=svd(A);
    err(iter)=max([norm(U-U0,Inf),norm(S-S0,Inf),norm(V-V0,Inf)]);
    end
err

The same comments as before apply to the output.

The same errors are also present, in fact manifesting more often, on
all more recent Octave releases that I have access to, including
Octave 3.0.0 as supplied with Ubuntu Hardy.


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

uname output:     Linux gawie 2.6.22-14-generic #1 SMP Tue Feb 12 07:42:25 UTC 2008 i686 GNU/Linux
configure opts:   '--prefix=/usr' '--datadir=/usr/share' '--libdir=/usr/lib' '--libexecdir=/usr/lib' '--infodir=/usr/share/info' '--mandir=/usr/share/man' '--with-blas=-lblas-3' '--with-lapack=-llapack-3' '--with-hdf5' '--with-fftw' '--with-f77=/usr/bin/gfortran' '--enable-shared' '--enable-rpath' '--disable-static' '--build' 'i486-linux-gnu' 'build_alias=i486-linux-gnu' 'CC=/usr/bin/gcc' 'CXX=/usr/bin/g++' 'F77=/usr/bin/gfortran'
Fortran compiler: /usr/bin/gfortran
FFLAGS:           -O2
F2C:              
F2CFLAGS:         
FLIBS:            -L/usr/lib/gcc/i486-linux-gnu/4.1.2 -L/usr/lib/gcc/i486-linux-gnu/4.1.2/../../../../lib -L/lib/../lib -L/usr/lib/../lib -lhdf5 -lz -lgfortranbegin -lgfortran -lm
CPPFLAGS:         
INCFLAGS:         -I. -I. -I./liboctave -I./src -I./libcruft/misc 
C compiler:       /usr/bin/gcc, version 4.1.2 20061103 (prerelease) (Ubuntu 4.1.1-18ubuntu2)
CFLAGS:           -O2
CPICFLAG:         -fPIC
C++ compiler:     /usr/bin/g++, version 4.1.2
CXXFLAGS:         -O2
CXXPICFLAG:       -fPIC
LD_CXX:           /usr/bin/g++
LDFLAGS:          -s
LIBFLAGS:         -L.
RLD_FLAG:         -Wl,-rpath -Wl,/usr/lib/octave-2.1.73
BLAS_LIBS:        -llapack-3 -lblas-3
FFTW_LIBS:        -lfftw3
LIBS:             -lreadline  -lncurses -ldl -lhdf5 -lz -lm 
LEXLIB:           
LIBPLPLOT:        
LIBDLFCN:         
LIBGLOB:          %OCTAVE_CONF_LIBGLOB%
SED:              /bin/sed
DEFS:

  -DPACKAGE_NAME="" -DPACKAGE_TARNAME="" -DPACKAGE_VERSION=""
  -DPACKAGE_STRING="" -DPACKAGE_BUGREPORT="" -DOCTAVE_SOURCE=1
  -D_GNU_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 -DSEPCHAR=':' -DSEPCHAR_STR=":" -D__NO_MATH_INLINES=1
  -DCXX_NEW_FRIEND_TEMPLATE_DECL=1 -DCXX_ISO_COMPLIANT_LIBRARY=1
  -DCXX_ABI=gnu_v3 -DHAVE_LIBM=1 -DHAVE_HDF5_H=1 -DHAVE_HDF5=1
  -DHAVE_H5GGET_NUM_OBJS=1 -DHAVE_FFTW3=1 -DHAVE_IEEE754_DATA_FORMAT=1
  -DF77_FUNC(name,NAME)=name ## _ -DF77_FUNC_(name,NAME)=name ## _
  -DHAVE_BLAS=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=4 -DSIZEOF_LONG_LONG=8
  -DHAVE_ALLOCA_H=1 -DHAVE_ALLOCA=1 -DNPOS=std::string::npos
  -DHAVE_PLACEMENT_DELETE=1 -DHAVE_DYNAMIC_AUTO_ARRAYS=1 -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_LIMITS_H=1 -DHAVE_MEMORY_H=1
  -DHAVE_NCURSES_H=1 -DHAVE_POLL_H=1 -DHAVE_PWD_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_SSTREAM=1 -DHAVE_TERMIO_H=1
  -DHAVE_SGTTY_H=1 -DHAVE_ATEXIT=1 -DHAVE_BASENAME=1 -DHAVE_BCOPY=1
  -DHAVE_BZERO=1 -DHAVE_CANONICALIZE_FILE_NAME=1 -DHAVE_DUP2=1
  -DHAVE_ENDGRENT=1 -DHAVE_ENDPWENT=1 -DHAVE_EXECVP=1 -DHAVE_FCNTL=1
  -DHAVE_FORK=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_LINK=1 -DHAVE_LOCALTIME_R=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_RENAME=1 -DHAVE_RINDEX=1
  -DHAVE_RMDIR=1 -DHAVE_ROUND=1 -DHAVE_SELECT=1 -DHAVE_SETGRENT=1
  -DHAVE_SETPWENT=1 -DHAVE_SETVBUF=1 -DHAVE_SIGACTION=1 -DHAVE_SIGLONGJMP=1
  -DHAVE_SIGPENDING=1 -DHAVE_SIGPROCMASK=1 -DHAVE_SIGSUSPEND=1 -DHAVE_STAT=1
  -DHAVE_STRCASECMP=1 -DHAVE_STRDUP=1 -DHAVE_STRERROR=1 -DHAVE_STRFTIME=1
  -DHAVE_STRNCASECMP=1 -DHAVE_STRPTIME=1 -DHAVE_SYMLINK=1 -DHAVE_TEMPNAM=1
  -DHAVE_UMASK=1 -DHAVE_UNLINK=1 -DHAVE_USLEEP=1 -DHAVE_VFPRINTF=1
  -DHAVE_VSPRINTF=1 -DHAVE_VSNPRINTF=1 -DHAVE_WAITPID=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_FINITE=1 -DHAVE_ISNAN=1 -DHAVE_ISINF=1 -DHAVE_COPYSIGN=1
  -DHAVE_DECL_SIGNBIT=1 -DHAVE_ACOSH=1 -DHAVE_ASINH=1 -DHAVE_ATANH=1
  -DHAVE_ERF=1 -DHAVE_ERFC=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 -DSYS_SIGLIST_DECLARED=1
  -DHAVE_SYS_SIGLIST=1 -DHAVE_POSIX_SIGNALS=1 -DHAVE_GETRUSAGE=1
  -DHAVE_TIMES=1 -DGNUPLOT_BINARY="gnuplot" -DGNUPLOT_HAS_FRAMES=

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

  DEFAULT_EXEC_PATH = "/usr/lib/octave/2.1.73/site/exec/i486-pc-linux-gnu:/usr/lib/octave/site/exec/i486-pc-linux-gnu:/usr/lib/octave/2.1.73/exec/i486-pc-linux-gnu:/usr/bin"
  DEFAULT_LOADPATH = ".:/usr/lib/octave/2.1.73/site/oct/i486-pc-linux-gnu//:/usr/lib/octave/site/oct/api-v13/i486-pc-linux-gnu//:/usr/lib/octave/site/oct/i486-pc-linux-gnu//:/usr/share/octave/2.1.73/site/m//:/usr/share/octave/site/api-v13/m//:/usr/share/octave/site/m//:/usr/lib/octave/2.1.73/oct/i486-pc-linux-gnu//:/usr/share/octave/2.1.73/m//"
  EDITOR = "emacs"
  EXEC_PATH = ":/home/dirk/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/bin/X11:/usr/games"
  IMAGEPATH = ".:/usr/share/octave/2.1.73/imagelib//"
  INFO_FILE = "/usr/share/info/octave2.1.info"
  INFO_PROGRAM = "info"
  LOADPATH = "/home/dirk/octave/lib//::/usr/local/share/octave/site-m//:"
  PAGER = "pager"
  PS1 = "> "
  PS2 = 
  PS4 = "+ "
  automatic_replot = 1
  beep_on_error = 1
  completion_append_char = " "
  crash_dumps_octave_core = 1
  default_save_format = "ascii"
  echo_executing_commands = 0
  fixed_point_format = 0
  gnuplot_binary = "gnuplot"
  gnuplot_command_end = "\n"
  gnuplot_command_plot = "pl"
  gnuplot_command_replot = "rep"
  gnuplot_command_splot = "sp"
  gnuplot_command_title = "t"
  gnuplot_command_using = "u"
  gnuplot_command_with = "w"
  gnuplot_has_frames = 1
  history_file = "/home/dirk/.octave_hist"
  history_size = 1024
  ignore_function_time_stamp = "system"
  max_recursion_depth = 256
  output_max_field_width = 20
  output_precision = 15
  page_output_immediately = 0
  page_screen_output = 1
  print_answer_id_name = 1
  print_empty_dimensions = 1
  print_rhs_assign_val = 0
# return_last_computed_value = <no value or error in displaying it>
  save_precision = 15
  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 = 1
  warn_assign_as_truth_value = 1
  warn_divide_by_zero = 1
  warn_empty_list_elements = 0
  warn_fortran_indexing = 0
  warn_function_name_clash = 0
  warn_future_time_stamp = 1
  warn_imag_to_real = 0
  warn_missing_semicolon = 0
  warn_neg_dim_as_zero = 0
  warn_num_to_str = 1
  warn_resize_on_range_error = 0
  warn_separator_insert = 0
  warn_single_quote_string = 0
  warn_str_to_num = 0
  warn_undefined_return_values = 1
  warn_variable_switch_label = 0


More information about the Bug-octave mailing list