#! /bin/sh # # To apply this patch, cd to the top level Octave source directory and # run this file through /bin/sh. It will first remove any files that # have been deleted from the source distribution since the last # release and then update the sources with patch(1). # # Diffs for updating *.ps, *.dvi, and *.info* files are not included # because they can be recreated from the Texinfo files using TeX and # makeinfo. # # Diffs for updating parse.cc and y.tab.h are not included because # they can be recreated from the file parse.y using bison. # # Diffs for updating lex.cc are not included because it can be # recreated from lex.l using flex. # # Diffs for updating the configure script are not included because # it can be recreated from configure.in using autoconf. # # John W. Eaton # jwe@bevo.che.wisc.edu # University of Wisconsin-Madison # Department of Chemical Engineering if test -f src/octave.cc ; then true else echo '***********************************************************' 1>&2 echo 'You must run this script in the top-level octave directory!' 1>&2 echo '***********************************************************' 1>&2 exit 1 fi ### ### Special commands should go here. ### rm -f make/ChangeLog make/README make/common.make make/config.make rm -f make/dist.make make/library.make make/makevars.make make/misc.make rm -f make/paths.make make/programs.make make/rdepend.make make/texi.make rm -f make/tkpathsea.make rmdir make echo 'patching existing files' patch -p1 << \PATCH_EOF diff -cNr octave-2.0.7/ChangeLog octave-2.0.8/ChangeLog *** octave-2.0.7/ChangeLog Wed Jun 4 12:44:21 1997 --- octave-2.0.8/ChangeLog Mon Jun 23 14:27:12 1997 *************** *** 1,3 **** --- 1,41 ---- + Mon Jun 23 14:26:56 1997 John W. Eaton + + * Version 2.0.8 released. + + Mon Jun 23 09:16:56 1997 John W. Eaton + + * configure.in (BOUNDS_CHECKING): Fix comment, allow bounds + checking to be enabled. + + Fri Jun 20 14:26:17 1997 John W. Eaton + + * configure.in: Define SH_LD, SH_LDFLAGS, and RLD_FLAG for + sparc-sun-sunos4*. + + Fri Jun 20 14:26:17 1997 John W. Eaton + + * mkoctfile.in: Handle --strip. + + Thu Jun 12 12:29:09 1997 John W. Eaton + + * make: Delete directory. + * octMakefile.in: Don't distribute files from make directory. + + Wed Jun 11 16:28:36 1997 John W. Eaton + + * mkoctfile.in: Allow more options, support for C and Fortran + source, existing object files, etc. Based on patch from Guido + Dietz . + + Fri Jun 6 21:55:46 1997 John W. Eaton + + * mkoctfile.in: Use eval to invoke compilation and linking + commands. + + Thu Jun 5 11:44:04 1997 John W. Eaton + + * mkoctfile.in (SH_LDFLAGS): Fix previous change. + Wed Jun 4 12:43:30 1997 John W. Eaton * Version 2.0.7 released. diff -cNr octave-2.0.7/NEWS octave-2.0.8/NEWS *** octave-2.0.7/NEWS Wed Jun 4 12:48:35 1997 --- octave-2.0.8/NEWS Mon Jun 23 14:02:57 1997 *************** *** 1,3 **** --- 1,16 ---- + Summary of changes for version 2.0.8: + ------------------------------------ + + This is a bug-fixing release. There are only a few new user-visible + features. + + * If the argument to eig() is symmetric, Octave uses the specialized + Lapack subroutine for symmetric matrices for a significant + increase in performance. + + * It is now possible to use the mkoctfile script to create .oct + files from multiple source and object files. + Summary of changes for version 2.0.7: ------------------------------------ diff -cNr octave-2.0.7/configure.in octave-2.0.8/configure.in *** octave-2.0.7/configure.in Mon Jun 2 13:49:53 1997 --- octave-2.0.8/configure.in Mon Jun 23 09:17:34 1997 *************** *** 179,186 **** BOUNDS_CHECKING=false AC_ARG_ENABLE(bounds-check, ! [ --enable-bounds-check for internal array classes (default is yes)], ! [if test "$enableval" = no; then BOUNDS_CHECKING=false; fi], []) if $BOUNDS_CHECKING; then AC_DEFINE(BOUNDS_CHECKING, 1) fi --- 179,186 ---- BOUNDS_CHECKING=false AC_ARG_ENABLE(bounds-check, ! [ --enable-bounds-check for internal array classes (default is no)], ! [if test "$enableval" = yes; then BOUNDS_CHECKING=true; fi], []) if $BOUNDS_CHECKING; then AC_DEFINE(BOUNDS_CHECKING, 1) fi *************** *** 605,610 **** --- 605,613 ---- else FPICFLAG=-PIC fi + SH_LD=ld + SH_LDFLAGS="-assert nodefinitions" + RLD_FLAG='-L$(libdir)' ;; sparc-sun-solaris2*) if $f77_is_g77; then diff -cNr octave-2.0.7/doc/ChangeLog octave-2.0.8/doc/ChangeLog *** octave-2.0.7/doc/ChangeLog Wed Jun 4 12:44:11 1997 --- octave-2.0.8/doc/ChangeLog Mon Jun 23 14:27:07 1997 *************** *** 1,3 **** --- 1,7 ---- + Mon Jun 23 14:26:56 1997 John W. Eaton + + * Version 2.0.8 released. + Wed Jun 4 12:43:30 1997 John W. Eaton * Version 2.0.7 released. diff -cNr octave-2.0.7/kpathsea/ChangeLog octave-2.0.8/kpathsea/ChangeLog *** octave-2.0.7/kpathsea/ChangeLog Wed Jun 4 02:33:42 1997 --- octave-2.0.8/kpathsea/ChangeLog Thu Jun 12 12:28:15 1997 *************** *** 1,3 **** --- 1,12 ---- + Thu Jun 12 12:24:09 1997 John W. Eaton + + * Makefile.in: Replace all ac_include commands with the files + that they reference. + (Makefile): Don't depend on files in ../make. + (depend depend.make): Depend on Makefile, not ../make/rdepend.make. + (MakeTeXPK): Depend on Makefile, not ../make/paths.make. + (texmf.sed): Depend on Makefile, not ../make/makevars.make + Wed Jun 4 02:33:28 1997 John W. Eaton * common.ac (program_invocation_name): Fix check. diff -cNr octave-2.0.7/kpathsea/Makefile octave-2.0.8/kpathsea/Makefile *** octave-2.0.7/kpathsea/Makefile Wed Jun 4 12:46:10 1997 --- octave-2.0.8/kpathsea/Makefile Thu Jun 12 12:31:36 1997 *************** *** 97,102 **** --- 97,103 ---- default_texsizes = 300:600 # End of paths.make. + # makevars.make -- the directory names we pass. # It's important that none of these values contain [ @%], for the sake # of kpathsea/texmf.sed. *************** *** 170,180 **** --- 171,183 ---- ##CFLAGS = -g $(XCFLAGS) ##endif # End of common.make. + # library.make -- stuff only useful for libraries. AR = ar ARFLAGS = cq RANLIB = ranlib # End of library.make. + # programs.make -- used by Makefiles for executables only. # Linking. Don't include $(CFLAGS), since ld -g under Linux forces # static libraries, including libc.a and libX*.a *************** *** 188,193 **** --- 191,197 ---- CCLD = $(CC) link_command = $(CCLD) -o $@ $(LDFLAGS) # End of programs.make. + # texi.make -- making .dvi and .info from .texi. MAKEINFO = makeinfo *************** *** 269,275 **** kpsewhich: libkpathsea.a kpsewhich.o $(link_command) kpsewhich.o $(LOADLIBES) ! MakeTeXPK: MakeTeXPK.in $(top_srcdir)/../make/paths.make sed -e 's,${prefix}/share,$(datadir),g' \ -e 's,@psheaderdir@,$(psheaderdir),g' \ -e 's,@fontnamedir@,$(fontnamedir),g' \ --- 273,279 ---- kpsewhich: libkpathsea.a kpsewhich.o $(link_command) kpsewhich.o $(LOADLIBES) ! MakeTeXPK: MakeTeXPK.in Makefile sed -e 's,${prefix}/share,$(datadir),g' \ -e 's,@psheaderdir@,$(psheaderdir),g' \ -e 's,@fontnamedir@,$(fontnamedir),g' \ *************** *** 289,295 **** # The idea is to turn each var=value into s%@var@%value%g. Seems simpler # to put the substitutions in a file than to play shell games. ! texmf.sed: $(top_srcdir)/../make/makevars.make echo $(makevars) \ | tr ' ' '\012' \ | sed -e 's/^/s%@/' -e 's/=/@%/' -e 's/$$/%/' -e 's/$$/g/' \ --- 293,299 ---- # The idea is to turn each var=value into s%@var@%value%g. Seems simpler # to put the substitutions in a file than to play shell games. ! texmf.sed: Makefile echo $(makevars) \ | tr ' ' '\012' \ | sed -e 's/^/s%@/' -e 's/=/@%/' -e 's/$$/%/' -e 's/$$/g/' \ *************** *** 444,450 **** config.status: $(srcdir)/configure $(SHELL) $(srcdir)/configure --no-create --verbose ! Makefile: $(srcdir)/Makefile.in config.status $(top_srcdir)/../make/*.make $(SHELL) config.status # This rule isn't used for web2c or the top-level Makefile, but it --- 448,454 ---- config.status: $(srcdir)/configure $(SHELL) $(srcdir)/configure --no-create --verbose ! Makefile: $(srcdir)/Makefile.in config.status $(SHELL) config.status # This rule isn't used for web2c or the top-level Makefile, but it *************** *** 517,523 **** # "kpathsea/..." in the sources. But that means we have to remove the # directory prefixes and all the system include files. # And is generated, not part of the distribution. ! depend depend.make:: c-auto.h $(top_srcdir)/../make/rdepend.make $(CC) -M $(CPPFLAGS) *.c \ | sed -e 's,\.\./kpathsea/,$$(kpathsea_srcdir)/,g' \ -e 's,$$(kpathsea_srcdir)/paths.h,paths.h,g' \ --- 521,527 ---- # "kpathsea/..." in the sources. But that means we have to remove the # directory prefixes and all the system include files. # And is generated, not part of the distribution. ! depend depend.make:: c-auto.h Makefile $(CC) -M $(CPPFLAGS) *.c \ | sed -e 's,\.\./kpathsea/,$$(kpathsea_srcdir)/,g' \ -e 's,$$(kpathsea_srcdir)/paths.h,paths.h,g' \ *************** *** 527,532 **** --- 531,537 ---- | grep -v '^ *\\$$' \ >depend.make # End of rdepend.make. + absolute.o: absolute.c $(kpathsea_srcdir)/config.h ./c-auto.h $(kpathsea_srcdir)/c-std.h \ $(kpathsea_srcdir)/c-unistd.h $(kpathsea_srcdir)/systypes.h \ $(kpathsea_srcdir)/c-memstr.h \ diff -cNr octave-2.0.7/kpathsea/Makefile.in octave-2.0.8/kpathsea/Makefile.in *** octave-2.0.7/kpathsea/Makefile.in Tue May 20 14:33:34 1997 --- octave-2.0.8/kpathsea/Makefile.in Thu Jun 12 12:28:09 1997 *************** *** 1,16 **** # Makefile for kpathsea --kb@cs.umb.edu. version = 2.52 ! ac_include ../make/paths.make ! ac_include ../make/makevars.make # Add -DNO_DEBUG to disable debugging, for unnoticeably better performance. DEFS = $(XDEFS) ! ac_include ../make/common.make ! ac_include ../make/library.make ! ac_include ../make/programs.make ! ac_include ../make/texi.make # Install these header files (except c-auto.h). install_headers = *.h --- 1,208 ---- # Makefile for kpathsea --kb@cs.umb.edu. version = 2.52 ! # paths.make -- installation directories. ! # ! # The compile-time paths are defined in kpathsea/paths.h, which is built ! # from kpathsea/paths.h.in and these definitions. See kpathsea/INSTALL ! # for a description of how the various path-related files are used and ! # created. ! ! # Do not change prefix and exec_prefix in Makefile.in! ! # configure doesn't propagate the change to the other Makefiles. ! # Instead, give the -prefix/-exec-prefix options to configure. ! # (See kpathsea/INSTALL for more details.) This is arguably ! # a bug, but it's not likely to change soon. ! prefix = @prefix@ ! exec_prefix = @exec_prefix@ ! ! # Architecture-dependent executables. ! bindir = $(exec_prefix)/bin ! ! # Architecture-independent executables. ! scriptdir = $(bindir) ! ! # Architecture-dependent files, such as lib*.a files. ! libdir = $(exec_prefix)/lib ! ! # Architecture-independent files. ! datadir = $(prefix)/lib ! ! # Header files. ! includedir = $(prefix)/include ! ! # GNU .info* files. ! infodir = $(prefix)/info ! ! # Unix man pages. ! manext = 1 ! mandir = $(prefix)/man/man$(manext) ! ! # TeX & MF-specific directories. Not all of the following are relevant ! # for all programs, but it seems cleaner to collect everything in one place. ! ! # The default paths are now in kpathsea/paths.h.in. Passing all the ! # paths to sub-makes can make the arg list too long on system V. ! ! # The root of the tree. ! texmf = $(datadir)/texmf ! ! # TeX and MF source files. ! texinputdir = $(texmf)/tex ! mfinputdir = $(texmf)/mf ! ! # MakeTeXPK.site, texmf.cnf, etc. ! web2cdir = $(texmf)/web2c ! ! # The top-level font directory. ! fontdir = $(texmf)/fonts ! ! # Memory dumps (.fmt and .base). ! fmtdir = $(texmf)/ini ! basedir = $(fmtdir) ! ! # Pool files. ! texpooldir = $(texmf)/ini ! mfpooldir = $(texpooldir) ! ! # If install_fonts=true, the PostScript/LaserJet TFM and VF files for ! # the builtin fonts get installed in subdirectories of this directory, ! # named for the typeface families of these directories. If you don't ! # have the default directory setup, you will want to set ! # install_fonts=false. Ditto for install_macros. ! install_fonts = true ! install_macros = true ! ! # Where the .map files from fontname are installed. ! fontnamedir = $(texmf)/fontname ! ! # Where the dvips configuration files get installed, and where ! # psfonts.map is. ! dvipsdir = $(texmf)/dvips ! psheaderdir = $(dvipsdir) ! ! # MakeTeXPK will go here to create dc*. ! dcfontdir = $(fontdir)/public/dc ! ! # MakeTeXPK will go here if it exists to create nonstandard CM fonts, ! # e.g., cmr11. See ftp.cs.umb.edu:pub/tex/sauter.tar.gz. The Sauter ! # files must be in your regular MFINPUTS. ! sauterdir = $(fontdir)/public/sauter ! ! # If a font can't be found close enough to its stated size, we look for ! # each of these sizes in the order given. This colon-separated list is ! # overridden by the envvar TEXSIZES, and by a program-specific variable ! # (e.g., XDVISIZES), and perhaps by a config file (e.g., in dvips). ! default_texsizes = 300:600 ! ! # End of paths.make. ! ! # makevars.make -- the directory names we pass. ! # It's important that none of these values contain [ @%], for the sake ! # of kpathsea/texmf.sed. ! makevars = prefix=$(prefix) exec_prefix=$(exec_prefix) \ ! bindir=$(bindir) scriptdir=$(scriptdir) libdir=$(libdir) \ ! datadir=$(datadir) infodir=$(infodir) includedir=$(includedir) \ ! manext=$(manext) mandir=$(mandir) \ ! texmf=$(texmf) web2cdir=$(web2cdir) \ ! texinputdir=$(texinputdir) mfinputdir=$(mfinputdir) \ ! fontdir=$(fontdir) \ ! fmtdir=$(fmtdir) basedir=$(basedir) \ ! texpooldir=$(texpooldir) mfpooldir=$(mfpooldir) \ ! install_fonts=$(install_fonts) \ ! dvipsdir=$(dvipsdir) psheaderdir=$(psheaderdir) \ ! default_texsizes='$(default_texsizes)' ! # End of makevars.make. # Add -DNO_DEBUG to disable debugging, for unnoticeably better performance. DEFS = $(XDEFS) ! # common.make -- used by all Makefiles. ! SHELL = /bin/sh ! @SET_MAKE@ ! top_srcdir = @top_srcdir@ ! srcdir = @srcdir@ ! VPATH = @srcdir@ ! ! CPICFLAG = @CPICFLAG@ ! ! SHLEXT = @SHLEXT@ ! ! SH_LD = @SH_LD@ ! SH_LDFLAGS = @SH_LDFLAGS@ ! ! SHARED_LIBS = @SHARED_LIBS@ ! ! CC = @CC@ ! # CFLAGS is used for both compilation and linking. ! CFLAGS = @CFLAGS@ $(XCFLAGS) ! ! # Do not override CPPFLAGS; change XCPPFLAGS, CFLAGS, XCFLAGS, or DEFS instead. ! CPPFLAGS = $(XCPPFLAGS) -I. -I$(srcdir) \ ! -I$(kpathsea_parent) -I$(kpathsea_srcdir_parent) \ ! $(prog_cflags) @CPPFLAGS@ $(DEFS) ! ! %.o : %.c ! $(CC) $(CPPFLAGS) $(CFLAGS) -c $< ! ! pic/%.o : %.c ! $(CC) $(CPPFLAGS) $(CPICFLAG) $(CFLAGS) -c $< -o $@ ! ! # Installation. ! INSTALL = @INSTALL@ ! INSTALL_PROGRAM = @INSTALL_PROGRAM@ ! INSTALL_DATA = @INSTALL_DATA@ ! ! # This is used to recursively copy a fonts/ or tex/ directory to ! # $(fontdir) or $(texinputdir). ! # The first arg is `.', and the second is the target directory. ! CP_R = cp -r ! ! # This is so kpathsea will get remade automatically if you change ! # something in it and recompile from the package directory. ! kpathsea_parent = .. ! kpathsea_dir = $(kpathsea_parent)/kpathsea ! kpathsea_srcdir_parent = $(top_srcdir)/.. ! kpathsea_srcdir = $(kpathsea_srcdir_parent)/kpathsea ! kpathsea = $(kpathsea_dir)/libkpathsea.a ! ! ##ifeq ($(CC), gcc) ! ##XDEFS = -Wall -Wpointer-arith $(warn_more) ! ##CFLAGS = -g $(XCFLAGS) ! ##endif ! # End of common.make. ! ! # library.make -- stuff only useful for libraries. ! AR = ar ! ARFLAGS = cq ! RANLIB = @RANLIB@ ! # End of library.make. ! ! # programs.make -- used by Makefiles for executables only. ! # Linking. Don't include $(CFLAGS), since ld -g under Linux forces ! # static libraries, including libc.a and libX*.a ! LDFLAGS = @LDFLAGS@ $(XLDFLAGS) ! LIBS = @LIBS@ ! # proglib is for web2c; ! # XLOADLIBES is for the installer. ! LOADLIBES= $(proglib) $(kpathsea) $(LIBS) -lm $(XLOADLIBES) ! ! # Why separate CCLD from CC? No particular reason. ! CCLD = $(CC) ! link_command = $(CCLD) -o $@ $(LDFLAGS) ! # End of programs.make. ! ! # texi.make -- making .dvi and .info from .texi. ! ! MAKEINFO = makeinfo ! MAKEINFO_FLAGS = --paragraph-indent=2 -I$(srcdir) ! TEXI2DVI = texi2dvi ! ! .SUFFIXES: .info .dvi .texi ! .texi.info: ! -$(MAKEINFO) $(MAKEINFO_FLAGS) $< -o $@ ! .texi.dvi: ! -$(TEXI2DVI) $(TEXI2DVI_FLAGS) $< # Install these header files (except c-auto.h). install_headers = *.h *************** *** 81,87 **** kpsewhich: libkpathsea.a kpsewhich.o $(link_command) kpsewhich.o $(LOADLIBES) ! MakeTeXPK: MakeTeXPK.in $(top_srcdir)/../make/paths.make sed -e 's,@datadir@,$(datadir),g' \ -e 's,@psheaderdir@,$(psheaderdir),g' \ -e 's,@fontnamedir@,$(fontnamedir),g' \ --- 273,279 ---- kpsewhich: libkpathsea.a kpsewhich.o $(link_command) kpsewhich.o $(LOADLIBES) ! MakeTeXPK: MakeTeXPK.in Makefile sed -e 's,@datadir@,$(datadir),g' \ -e 's,@psheaderdir@,$(psheaderdir),g' \ -e 's,@fontnamedir@,$(fontnamedir),g' \ *************** *** 101,107 **** # The idea is to turn each var=value into s%@var@%value%g. Seems simpler # to put the substitutions in a file than to play shell games. ! texmf.sed: $(top_srcdir)/../make/makevars.make echo $(makevars) \ | tr ' ' '\012' \ | sed -e 's/^/s%@/' -e 's/=/@%/' -e 's/$$/%/' -e 's/$$/g/' \ --- 293,299 ---- # The idea is to turn each var=value into s%@var@%value%g. Seems simpler # to put the substitutions in a file than to play shell games. ! texmf.sed: Makefile echo $(makevars) \ | tr ' ' '\012' \ | sed -e 's/^/s%@/' -e 's/=/@%/' -e 's/$$/%/' -e 's/$$/g/' \ *************** *** 236,242 **** fi .PHONY: bin-dist ! ac_include ../make/config.make info: kpathsea.info dvi: kpathsea.dvi --- 428,472 ---- fi .PHONY: bin-dist ! # config.make -- autoconf rules to remake the Makefile, c-auto.h, etc. ! ! ##ifdef HOSTNAME ! ##ac_dir = $(gnu)/share/autoconf ! ##autoconf = $(ac_dir)/acspecific.m4 $(ac_dir)/acgeneral.m4 $(ac_dir)/acsite.m4 ! ##autoheader = $(ac_dir)/acconfig.h ! ## ! ### I define $(autoconf) to acgeneral.m4 and the other Autoconf files, so ! ### configure automatically gets remade in the sources with a new Autoconf ! ### release. But it would be bad for installers with Autoconf to remake ! ### configure (not to mention require Autoconf), so I take out the variable ! ### $(autoconf) definition before release. ! ##configure_in = $(srcdir)/configure.in $(kpathsea_srcdir)/common.ac ! ##$(srcdir)/configure: $(configure_in) $(autoconf) ! ## cd $(srcdir) && autoconf ! ##endif ! ! config.status: $(srcdir)/configure ! $(SHELL) $(srcdir)/configure --no-create --verbose ! ! Makefile: $(srcdir)/Makefile.in config.status ! $(SHELL) config.status ! ! # This rule isn't used for web2c or the top-level Makefile, but it ! # doesn't hurt. We don't depend on config.status because configure ! # always rewrites config.status, even when it doesn't change. Thus it ! # might be newer than c-auto.h when we don't need to remake the latter. ! c-auto.h: $(srcdir)/stamp-auto ! $(srcdir)/stamp-auto: $(srcdir)/c-auto.h.in ! $(SHELL) config.status ! touch $(srcdir)/stamp-auto ! ! ##ifdef HOSTNAME ! ### autoheader reads acconfig.h (and c-auto.h.top) automatically. ! ##$(srcdir)/c-auto.h.in: $(srcdir)/stamp-auto.in ! ##$(srcdir)/stamp-auto.in: $(configure_in) $(autoheader) $(srcdir)/acconfig.h ! ## cd $(srcdir) && autoheader ! ## touch $(srcdir)/stamp-auto.in ! ##endif info: kpathsea.info dvi: kpathsea.dvi *************** *** 245,251 **** | sed -n -e '/^Installation/,/wrong fonts/'p >$@ add-info-toc $@ ! ac_include ../make/misc.make mostlyclean:: rm -f kpsewhich --- 475,514 ---- | sed -n -e '/^Installation/,/wrong fonts/'p >$@ add-info-toc $@ ! # misc.make -- cleaning, etc. ! TAGS: *.c *.h ! if pwd | grep kpathsea >/dev/null; then \ ! etags *.c *.h; else etags -i $(kpathsea_dir)/TAGS *.c *.h; fi ! ! mostlyclean:: ! rm -f *.o $(program) $(programs) squeeze $(library).a pic/*.o ! if $(SHARED_LIBS); then rm -f *.$(SHLEXT); fi ! ! clean:: mostlyclean ! rm -f *.dvi *.lj ! ! distclean:: clean ! rm -f Makefile MakeTeXPK *.pool ! rm -f config.status config.log config.cache c-auto.h ! rm -f stamp-picdir ! -rmdir pic ! ! # Although we can remake configure and c-auto.h.in, we don't remove ! # them, since many people may lack Autoconf. Use configclean for that. ! maintainer-clean:: distclean ! rm -f TAGS *.info* configure stamp-auto ! ! extraclean:: ! rm -f *.aux *.bak *.bbl *.blg *.dvi *.log *.orig *.pl *.rej ! rm -f *.i *.s *.tfm *.vf *.vpl *\#* *gf *pk *~ ! rm -f CONTENTS.tex a.out core mfput.* texput.* ! ! configclean: ! rm -f configure c-auto.h.in c-auto.h ! ! # Prevent GNU make 3.[59,63) from overflowing arg limit on system V. ! .NOEXPORT: ! # End of misc.make. mostlyclean:: rm -f kpsewhich *************** *** 253,259 **** rm -f paths.h texmf.cnf texmf.sed so_locations rm -f libkpathsea.a libkpathsea.$(SHLEXT) ! ac_include ../make/rdepend.make absolute.o: absolute.c $(kpathsea_srcdir)/config.h ./c-auto.h $(kpathsea_srcdir)/c-std.h \ $(kpathsea_srcdir)/c-unistd.h $(kpathsea_srcdir)/systypes.h \ $(kpathsea_srcdir)/c-memstr.h \ --- 516,537 ---- rm -f paths.h texmf.cnf texmf.sed so_locations rm -f libkpathsea.a libkpathsea.$(SHLEXT) ! # rdepend.make -- rules for remaking the dependencies. ! # Have to use -M, not -MM, since we use instead of ! # "kpathsea/..." in the sources. But that means we have to remove the ! # directory prefixes and all the system include files. ! # And is generated, not part of the distribution. ! depend depend.make:: c-auto.h Makefile ! $(CC) -M $(CPPFLAGS) *.c \ ! | sed -e 's,\.\./kpathsea/,$$(kpathsea_srcdir)/,g' \ ! -e 's,$$(kpathsea_srcdir)/paths.h,paths.h,g' \ ! -e 's,/usr[^ ]* ,,g' \ ! -e 's,/usr[^ ]*$$,,g' \ ! -e 's,dvi2xx.o,dvilj.o dvilj2p.o dvilj4.o dvilj4l.o,' \ ! | grep -v '^ *\\$$' \ ! >depend.make ! # End of rdepend.make. ! absolute.o: absolute.c $(kpathsea_srcdir)/config.h ./c-auto.h $(kpathsea_srcdir)/c-std.h \ $(kpathsea_srcdir)/c-unistd.h $(kpathsea_srcdir)/systypes.h \ $(kpathsea_srcdir)/c-memstr.h \ diff -cNr octave-2.0.7/libcruft/ChangeLog octave-2.0.8/libcruft/ChangeLog *** octave-2.0.7/libcruft/ChangeLog Wed Jun 4 12:44:19 1997 --- octave-2.0.8/libcruft/ChangeLog Mon Jun 23 14:27:11 1997 *************** *** 1,3 **** --- 1,11 ---- + Mon Jun 23 14:26:56 1997 John W. Eaton + + * Version 2.0.8 released. + + Fri Jun 6 21:53:17 1997 John W. Eaton + + * slatec-fn/xdgamma.f: New file. + Wed Jun 4 12:43:30 1997 John W. Eaton * Version 2.0.7 released. diff -cNr octave-2.0.7/libcruft/blas/dsymv.f octave-2.0.8/libcruft/blas/dsymv.f *** octave-2.0.7/libcruft/blas/dsymv.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/blas/dsymv.f Mon Jun 23 14:10:50 1997 *************** *** 0 **** --- 1,265 ---- + * + ************************************************************************ + * + SUBROUTINE DSYMV ( UPLO, N, ALPHA, A, LDA, X, INCX, + $ BETA, Y, INCY ) + * .. Scalar Arguments .. + DOUBLE PRECISION ALPHA, BETA + INTEGER INCX, INCY, LDA, N + CHARACTER*1 UPLO + * .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) + * .. + * + * Purpose + * ======= + * + * DSYMV performs the matrix-vector operation + * + * y := alpha*A*x + beta*y, + * + * where alpha and beta are scalars, x and y are n element vectors and + * A is an n by n symmetric matrix. + * + * Parameters + * ========== + * + * UPLO - CHARACTER*1. + * On entry, UPLO specifies whether the upper or lower + * triangular part of the array A is to be referenced as + * follows: + * + * UPLO = 'U' or 'u' Only the upper triangular part of A + * is to be referenced. + * + * UPLO = 'L' or 'l' Only the lower triangular part of A + * is to be referenced. + * + * Unchanged on exit. + * + * N - INTEGER. + * On entry, N specifies the order of the matrix A. + * N must be at least zero. + * Unchanged on exit. + * + * ALPHA - DOUBLE PRECISION. + * On entry, ALPHA specifies the scalar alpha. + * Unchanged on exit. + * + * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). + * Before entry with UPLO = 'U' or 'u', the leading n by n + * upper triangular part of the array A must contain the upper + * triangular part of the symmetric matrix and the strictly + * lower triangular part of A is not referenced. + * Before entry with UPLO = 'L' or 'l', the leading n by n + * lower triangular part of the array A must contain the lower + * triangular part of the symmetric matrix and the strictly + * upper triangular part of A is not referenced. + * Unchanged on exit. + * + * LDA - INTEGER. + * On entry, LDA specifies the first dimension of A as declared + * in the calling (sub) program. LDA must be at least + * max( 1, n ). + * Unchanged on exit. + * + * X - DOUBLE PRECISION array of dimension at least + * ( 1 + ( n - 1 )*abs( INCX ) ). + * Before entry, the incremented array X must contain the n + * element vector x. + * Unchanged on exit. + * + * INCX - INTEGER. + * On entry, INCX specifies the increment for the elements of + * X. INCX must not be zero. + * Unchanged on exit. + * + * BETA - DOUBLE PRECISION. + * On entry, BETA specifies the scalar beta. When BETA is + * supplied as zero then Y need not be set on input. + * Unchanged on exit. + * + * Y - DOUBLE PRECISION array of dimension at least + * ( 1 + ( n - 1 )*abs( INCY ) ). + * Before entry, the incremented array Y must contain the n + * element vector y. On exit, Y is overwritten by the updated + * vector y. + * + * INCY - INTEGER. + * On entry, INCY specifies the increment for the elements of + * Y. INCY must not be zero. + * Unchanged on exit. + * + * + * Level 2 Blas routine. + * + * -- Written on 22-October-1986. + * Jack Dongarra, Argonne National Lab. + * Jeremy Du Croz, Nag Central Office. + * Sven Hammarling, Nag Central Office. + * Richard Hanson, Sandia National Labs. + * + * + * .. Parameters .. + DOUBLE PRECISION ONE , ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) + * .. Local Scalars .. + DOUBLE PRECISION TEMP1, TEMP2 + INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY + * .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME + * .. External Subroutines .. + EXTERNAL XERBLA + * .. Intrinsic Functions .. + INTRINSIC MAX + * .. + * .. Executable Statements .. + * + * Test the input parameters. + * + INFO = 0 + IF ( .NOT.LSAME( UPLO, 'U' ).AND. + $ .NOT.LSAME( UPLO, 'L' ) )THEN + INFO = 1 + ELSE IF( N.LT.0 )THEN + INFO = 2 + ELSE IF( LDA.LT.MAX( 1, N ) )THEN + INFO = 5 + ELSE IF( INCX.EQ.0 )THEN + INFO = 7 + ELSE IF( INCY.EQ.0 )THEN + INFO = 10 + END IF + IF( INFO.NE.0 )THEN + CALL XERBLA( 'DSYMV ', INFO ) + RETURN + END IF + * + * Quick return if possible. + * + IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) + $ RETURN + * + * Set up the start points in X and Y. + * + IF( INCX.GT.0 )THEN + KX = 1 + ELSE + KX = 1 - ( N - 1 )*INCX + END IF + IF( INCY.GT.0 )THEN + KY = 1 + ELSE + KY = 1 - ( N - 1 )*INCY + END IF + * + * Start the operations. In this version the elements of A are + * accessed sequentially with one pass through the triangular part + * of A. + * + * First form y := beta*y. + * + IF( BETA.NE.ONE )THEN + IF( INCY.EQ.1 )THEN + IF( BETA.EQ.ZERO )THEN + DO 10, I = 1, N + Y( I ) = ZERO + 10 CONTINUE + ELSE + DO 20, I = 1, N + Y( I ) = BETA*Y( I ) + 20 CONTINUE + END IF + ELSE + IY = KY + IF( BETA.EQ.ZERO )THEN + DO 30, I = 1, N + Y( IY ) = ZERO + IY = IY + INCY + 30 CONTINUE + ELSE + DO 40, I = 1, N + Y( IY ) = BETA*Y( IY ) + IY = IY + INCY + 40 CONTINUE + END IF + END IF + END IF + IF( ALPHA.EQ.ZERO ) + $ RETURN + IF( LSAME( UPLO, 'U' ) )THEN + * + * Form y when A is stored in upper triangle. + * + IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN + DO 60, J = 1, N + TEMP1 = ALPHA*X( J ) + TEMP2 = ZERO + DO 50, I = 1, J - 1 + Y( I ) = Y( I ) + TEMP1*A( I, J ) + TEMP2 = TEMP2 + A( I, J )*X( I ) + 50 CONTINUE + Y( J ) = Y( J ) + TEMP1*A( J, J ) + ALPHA*TEMP2 + 60 CONTINUE + ELSE + JX = KX + JY = KY + DO 80, J = 1, N + TEMP1 = ALPHA*X( JX ) + TEMP2 = ZERO + IX = KX + IY = KY + DO 70, I = 1, J - 1 + Y( IY ) = Y( IY ) + TEMP1*A( I, J ) + TEMP2 = TEMP2 + A( I, J )*X( IX ) + IX = IX + INCX + IY = IY + INCY + 70 CONTINUE + Y( JY ) = Y( JY ) + TEMP1*A( J, J ) + ALPHA*TEMP2 + JX = JX + INCX + JY = JY + INCY + 80 CONTINUE + END IF + ELSE + * + * Form y when A is stored in lower triangle. + * + IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN + DO 100, J = 1, N + TEMP1 = ALPHA*X( J ) + TEMP2 = ZERO + Y( J ) = Y( J ) + TEMP1*A( J, J ) + DO 90, I = J + 1, N + Y( I ) = Y( I ) + TEMP1*A( I, J ) + TEMP2 = TEMP2 + A( I, J )*X( I ) + 90 CONTINUE + Y( J ) = Y( J ) + ALPHA*TEMP2 + 100 CONTINUE + ELSE + JX = KX + JY = KY + DO 120, J = 1, N + TEMP1 = ALPHA*X( JX ) + TEMP2 = ZERO + Y( JY ) = Y( JY ) + TEMP1*A( J, J ) + IX = JX + IY = JY + DO 110, I = J + 1, N + IX = IX + INCX + IY = IY + INCY + Y( IY ) = Y( IY ) + TEMP1*A( I, J ) + TEMP2 = TEMP2 + A( I, J )*X( IX ) + 110 CONTINUE + Y( JY ) = Y( JY ) + ALPHA*TEMP2 + JX = JX + INCX + JY = JY + INCY + 120 CONTINUE + END IF + END IF + * + RETURN + * + * End of DSYMV . + * + END diff -cNr octave-2.0.7/libcruft/blas/dsyr2.f octave-2.0.8/libcruft/blas/dsyr2.f *** octave-2.0.7/libcruft/blas/dsyr2.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/blas/dsyr2.f Mon Jun 23 14:10:50 1997 *************** *** 0 **** --- 1,233 ---- + * + ************************************************************************ + * + SUBROUTINE DSYR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA ) + * .. Scalar Arguments .. + DOUBLE PRECISION ALPHA + INTEGER INCX, INCY, LDA, N + CHARACTER*1 UPLO + * .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) + * .. + * + * Purpose + * ======= + * + * DSYR2 performs the symmetric rank 2 operation + * + * A := alpha*x*y' + alpha*y*x' + A, + * + * where alpha is a scalar, x and y are n element vectors and A is an n + * by n symmetric matrix. + * + * Parameters + * ========== + * + * UPLO - CHARACTER*1. + * On entry, UPLO specifies whether the upper or lower + * triangular part of the array A is to be referenced as + * follows: + * + * UPLO = 'U' or 'u' Only the upper triangular part of A + * is to be referenced. + * + * UPLO = 'L' or 'l' Only the lower triangular part of A + * is to be referenced. + * + * Unchanged on exit. + * + * N - INTEGER. + * On entry, N specifies the order of the matrix A. + * N must be at least zero. + * Unchanged on exit. + * + * ALPHA - DOUBLE PRECISION. + * On entry, ALPHA specifies the scalar alpha. + * Unchanged on exit. + * + * X - DOUBLE PRECISION array of dimension at least + * ( 1 + ( n - 1 )*abs( INCX ) ). + * Before entry, the incremented array X must contain the n + * element vector x. + * Unchanged on exit. + * + * INCX - INTEGER. + * On entry, INCX specifies the increment for the elements of + * X. INCX must not be zero. + * Unchanged on exit. + * + * Y - DOUBLE PRECISION array of dimension at least + * ( 1 + ( n - 1 )*abs( INCY ) ). + * Before entry, the incremented array Y must contain the n + * element vector y. + * Unchanged on exit. + * + * INCY - INTEGER. + * On entry, INCY specifies the increment for the elements of + * Y. INCY must not be zero. + * Unchanged on exit. + * + * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). + * Before entry with UPLO = 'U' or 'u', the leading n by n + * upper triangular part of the array A must contain the upper + * triangular part of the symmetric matrix and the strictly + * lower triangular part of A is not referenced. On exit, the + * upper triangular part of the array A is overwritten by the + * upper triangular part of the updated matrix. + * Before entry with UPLO = 'L' or 'l', the leading n by n + * lower triangular part of the array A must contain the lower + * triangular part of the symmetric matrix and the strictly + * upper triangular part of A is not referenced. On exit, the + * lower triangular part of the array A is overwritten by the + * lower triangular part of the updated matrix. + * + * LDA - INTEGER. + * On entry, LDA specifies the first dimension of A as declared + * in the calling (sub) program. LDA must be at least + * max( 1, n ). + * Unchanged on exit. + * + * + * Level 2 Blas routine. + * + * -- Written on 22-October-1986. + * Jack Dongarra, Argonne National Lab. + * Jeremy Du Croz, Nag Central Office. + * Sven Hammarling, Nag Central Office. + * Richard Hanson, Sandia National Labs. + * + * + * .. Parameters .. + DOUBLE PRECISION ZERO + PARAMETER ( ZERO = 0.0D+0 ) + * .. Local Scalars .. + DOUBLE PRECISION TEMP1, TEMP2 + INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY + * .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME + * .. External Subroutines .. + EXTERNAL XERBLA + * .. Intrinsic Functions .. + INTRINSIC MAX + * .. + * .. Executable Statements .. + * + * Test the input parameters. + * + INFO = 0 + IF ( .NOT.LSAME( UPLO, 'U' ).AND. + $ .NOT.LSAME( UPLO, 'L' ) )THEN + INFO = 1 + ELSE IF( N.LT.0 )THEN + INFO = 2 + ELSE IF( INCX.EQ.0 )THEN + INFO = 5 + ELSE IF( INCY.EQ.0 )THEN + INFO = 7 + ELSE IF( LDA.LT.MAX( 1, N ) )THEN + INFO = 9 + END IF + IF( INFO.NE.0 )THEN + CALL XERBLA( 'DSYR2 ', INFO ) + RETURN + END IF + * + * Quick return if possible. + * + IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) + $ RETURN + * + * Set up the start points in X and Y if the increments are not both + * unity. + * + IF( ( INCX.NE.1 ).OR.( INCY.NE.1 ) )THEN + IF( INCX.GT.0 )THEN + KX = 1 + ELSE + KX = 1 - ( N - 1 )*INCX + END IF + IF( INCY.GT.0 )THEN + KY = 1 + ELSE + KY = 1 - ( N - 1 )*INCY + END IF + JX = KX + JY = KY + END IF + * + * Start the operations. In this version the elements of A are + * accessed sequentially with one pass through the triangular part + * of A. + * + IF( LSAME( UPLO, 'U' ) )THEN + * + * Form A when A is stored in the upper triangle. + * + IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN + DO 20, J = 1, N + IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN + TEMP1 = ALPHA*Y( J ) + TEMP2 = ALPHA*X( J ) + DO 10, I = 1, J + A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2 + 10 CONTINUE + END IF + 20 CONTINUE + ELSE + DO 40, J = 1, N + IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN + TEMP1 = ALPHA*Y( JY ) + TEMP2 = ALPHA*X( JX ) + IX = KX + IY = KY + DO 30, I = 1, J + A( I, J ) = A( I, J ) + X( IX )*TEMP1 + $ + Y( IY )*TEMP2 + IX = IX + INCX + IY = IY + INCY + 30 CONTINUE + END IF + JX = JX + INCX + JY = JY + INCY + 40 CONTINUE + END IF + ELSE + * + * Form A when A is stored in the lower triangle. + * + IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN + DO 60, J = 1, N + IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN + TEMP1 = ALPHA*Y( J ) + TEMP2 = ALPHA*X( J ) + DO 50, I = J, N + A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2 + 50 CONTINUE + END IF + 60 CONTINUE + ELSE + DO 80, J = 1, N + IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN + TEMP1 = ALPHA*Y( JY ) + TEMP2 = ALPHA*X( JX ) + IX = JX + IY = JY + DO 70, I = J, N + A( I, J ) = A( I, J ) + X( IX )*TEMP1 + $ + Y( IY )*TEMP2 + IX = IX + INCX + IY = IY + INCY + 70 CONTINUE + END IF + JX = JX + INCX + JY = JY + INCY + 80 CONTINUE + END IF + END IF + * + RETURN + * + * End of DSYR2 . + * + END diff -cNr octave-2.0.7/libcruft/blas/dsyr2k.f octave-2.0.8/libcruft/blas/dsyr2k.f *** octave-2.0.7/libcruft/blas/dsyr2k.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/blas/dsyr2k.f Mon Jun 23 14:10:50 1997 *************** *** 0 **** --- 1,330 ---- + * + ************************************************************************ + * + SUBROUTINE DSYR2K( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, + $ BETA, C, LDC ) + * .. Scalar Arguments .. + CHARACTER*1 UPLO, TRANS + INTEGER N, K, LDA, LDB, LDC + DOUBLE PRECISION ALPHA, BETA + * .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) + * .. + * + * Purpose + * ======= + * + * DSYR2K performs one of the symmetric rank 2k operations + * + * C := alpha*A*B' + alpha*B*A' + beta*C, + * + * or + * + * C := alpha*A'*B + alpha*B'*A + beta*C, + * + * where alpha and beta are scalars, C is an n by n symmetric matrix + * and A and B are n by k matrices in the first case and k by n + * matrices in the second case. + * + * Parameters + * ========== + * + * UPLO - CHARACTER*1. + * On entry, UPLO specifies whether the upper or lower + * triangular part of the array C is to be referenced as + * follows: + * + * UPLO = 'U' or 'u' Only the upper triangular part of C + * is to be referenced. + * + * UPLO = 'L' or 'l' Only the lower triangular part of C + * is to be referenced. + * + * Unchanged on exit. + * + * TRANS - CHARACTER*1. + * On entry, TRANS specifies the operation to be performed as + * follows: + * + * TRANS = 'N' or 'n' C := alpha*A*B' + alpha*B*A' + + * beta*C. + * + * TRANS = 'T' or 't' C := alpha*A'*B + alpha*B'*A + + * beta*C. + * + * TRANS = 'C' or 'c' C := alpha*A'*B + alpha*B'*A + + * beta*C. + * + * Unchanged on exit. + * + * N - INTEGER. + * On entry, N specifies the order of the matrix C. N must be + * at least zero. + * Unchanged on exit. + * + * K - INTEGER. + * On entry with TRANS = 'N' or 'n', K specifies the number + * of columns of the matrices A and B, and on entry with + * TRANS = 'T' or 't' or 'C' or 'c', K specifies the number + * of rows of the matrices A and B. K must be at least zero. + * Unchanged on exit. + * + * ALPHA - DOUBLE PRECISION. + * On entry, ALPHA specifies the scalar alpha. + * Unchanged on exit. + * + * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is + * k when TRANS = 'N' or 'n', and is n otherwise. + * Before entry with TRANS = 'N' or 'n', the leading n by k + * part of the array A must contain the matrix A, otherwise + * the leading k by n part of the array A must contain the + * matrix A. + * Unchanged on exit. + * + * LDA - INTEGER. + * On entry, LDA specifies the first dimension of A as declared + * in the calling (sub) program. When TRANS = 'N' or 'n' + * then LDA must be at least max( 1, n ), otherwise LDA must + * be at least max( 1, k ). + * Unchanged on exit. + * + * B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is + * k when TRANS = 'N' or 'n', and is n otherwise. + * Before entry with TRANS = 'N' or 'n', the leading n by k + * part of the array B must contain the matrix B, otherwise + * the leading k by n part of the array B must contain the + * matrix B. + * Unchanged on exit. + * + * LDB - INTEGER. + * On entry, LDB specifies the first dimension of B as declared + * in the calling (sub) program. When TRANS = 'N' or 'n' + * then LDB must be at least max( 1, n ), otherwise LDB must + * be at least max( 1, k ). + * Unchanged on exit. + * + * BETA - DOUBLE PRECISION. + * On entry, BETA specifies the scalar beta. + * Unchanged on exit. + * + * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). + * Before entry with UPLO = 'U' or 'u', the leading n by n + * upper triangular part of the array C must contain the upper + * triangular part of the symmetric matrix and the strictly + * lower triangular part of C is not referenced. On exit, the + * upper triangular part of the array C is overwritten by the + * upper triangular part of the updated matrix. + * Before entry with UPLO = 'L' or 'l', the leading n by n + * lower triangular part of the array C must contain the lower + * triangular part of the symmetric matrix and the strictly + * upper triangular part of C is not referenced. On exit, the + * lower triangular part of the array C is overwritten by the + * lower triangular part of the updated matrix. + * + * LDC - INTEGER. + * On entry, LDC specifies the first dimension of C as declared + * in the calling (sub) program. LDC must be at least + * max( 1, n ). + * Unchanged on exit. + * + * + * Level 3 Blas routine. + * + * + * -- Written on 8-February-1989. + * Jack Dongarra, Argonne National Laboratory. + * Iain Duff, AERE Harwell. + * Jeremy Du Croz, Numerical Algorithms Group Ltd. + * Sven Hammarling, Numerical Algorithms Group Ltd. + * + * + * .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME + * .. External Subroutines .. + EXTERNAL XERBLA + * .. Intrinsic Functions .. + INTRINSIC MAX + * .. Local Scalars .. + LOGICAL UPPER + INTEGER I, INFO, J, L, NROWA + DOUBLE PRECISION TEMP1, TEMP2 + * .. Parameters .. + DOUBLE PRECISION ONE , ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) + * .. + * .. Executable Statements .. + * + * Test the input parameters. + * + IF( LSAME( TRANS, 'N' ) )THEN + NROWA = N + ELSE + NROWA = K + END IF + UPPER = LSAME( UPLO, 'U' ) + * + INFO = 0 + IF( ( .NOT.UPPER ).AND. + $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN + INFO = 1 + ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND. + $ ( .NOT.LSAME( TRANS, 'T' ) ).AND. + $ ( .NOT.LSAME( TRANS, 'C' ) ) )THEN + INFO = 2 + ELSE IF( N .LT.0 )THEN + INFO = 3 + ELSE IF( K .LT.0 )THEN + INFO = 4 + ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN + INFO = 7 + ELSE IF( LDB.LT.MAX( 1, NROWA ) )THEN + INFO = 9 + ELSE IF( LDC.LT.MAX( 1, N ) )THEN + INFO = 12 + END IF + IF( INFO.NE.0 )THEN + CALL XERBLA( 'DSYR2K', INFO ) + RETURN + END IF + * + * Quick return if possible. + * + IF( ( N.EQ.0 ).OR. + $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) + $ RETURN + * + * And when alpha.eq.zero. + * + IF( ALPHA.EQ.ZERO )THEN + IF( UPPER )THEN + IF( BETA.EQ.ZERO )THEN + DO 20, J = 1, N + DO 10, I = 1, J + C( I, J ) = ZERO + 10 CONTINUE + 20 CONTINUE + ELSE + DO 40, J = 1, N + DO 30, I = 1, J + C( I, J ) = BETA*C( I, J ) + 30 CONTINUE + 40 CONTINUE + END IF + ELSE + IF( BETA.EQ.ZERO )THEN + DO 60, J = 1, N + DO 50, I = J, N + C( I, J ) = ZERO + 50 CONTINUE + 60 CONTINUE + ELSE + DO 80, J = 1, N + DO 70, I = J, N + C( I, J ) = BETA*C( I, J ) + 70 CONTINUE + 80 CONTINUE + END IF + END IF + RETURN + END IF + * + * Start the operations. + * + IF( LSAME( TRANS, 'N' ) )THEN + * + * Form C := alpha*A*B' + alpha*B*A' + C. + * + IF( UPPER )THEN + DO 130, J = 1, N + IF( BETA.EQ.ZERO )THEN + DO 90, I = 1, J + C( I, J ) = ZERO + 90 CONTINUE + ELSE IF( BETA.NE.ONE )THEN + DO 100, I = 1, J + C( I, J ) = BETA*C( I, J ) + 100 CONTINUE + END IF + DO 120, L = 1, K + IF( ( A( J, L ).NE.ZERO ).OR. + $ ( B( J, L ).NE.ZERO ) )THEN + TEMP1 = ALPHA*B( J, L ) + TEMP2 = ALPHA*A( J, L ) + DO 110, I = 1, J + C( I, J ) = C( I, J ) + + $ A( I, L )*TEMP1 + B( I, L )*TEMP2 + 110 CONTINUE + END IF + 120 CONTINUE + 130 CONTINUE + ELSE + DO 180, J = 1, N + IF( BETA.EQ.ZERO )THEN + DO 140, I = J, N + C( I, J ) = ZERO + 140 CONTINUE + ELSE IF( BETA.NE.ONE )THEN + DO 150, I = J, N + C( I, J ) = BETA*C( I, J ) + 150 CONTINUE + END IF + DO 170, L = 1, K + IF( ( A( J, L ).NE.ZERO ).OR. + $ ( B( J, L ).NE.ZERO ) )THEN + TEMP1 = ALPHA*B( J, L ) + TEMP2 = ALPHA*A( J, L ) + DO 160, I = J, N + C( I, J ) = C( I, J ) + + $ A( I, L )*TEMP1 + B( I, L )*TEMP2 + 160 CONTINUE + END IF + 170 CONTINUE + 180 CONTINUE + END IF + ELSE + * + * Form C := alpha*A'*B + alpha*B'*A + C. + * + IF( UPPER )THEN + DO 210, J = 1, N + DO 200, I = 1, J + TEMP1 = ZERO + TEMP2 = ZERO + DO 190, L = 1, K + TEMP1 = TEMP1 + A( L, I )*B( L, J ) + TEMP2 = TEMP2 + B( L, I )*A( L, J ) + 190 CONTINUE + IF( BETA.EQ.ZERO )THEN + C( I, J ) = ALPHA*TEMP1 + ALPHA*TEMP2 + ELSE + C( I, J ) = BETA *C( I, J ) + + $ ALPHA*TEMP1 + ALPHA*TEMP2 + END IF + 200 CONTINUE + 210 CONTINUE + ELSE + DO 240, J = 1, N + DO 230, I = J, N + TEMP1 = ZERO + TEMP2 = ZERO + DO 220, L = 1, K + TEMP1 = TEMP1 + A( L, I )*B( L, J ) + TEMP2 = TEMP2 + B( L, I )*A( L, J ) + 220 CONTINUE + IF( BETA.EQ.ZERO )THEN + C( I, J ) = ALPHA*TEMP1 + ALPHA*TEMP2 + ELSE + C( I, J ) = BETA *C( I, J ) + + $ ALPHA*TEMP1 + ALPHA*TEMP2 + END IF + 230 CONTINUE + 240 CONTINUE + END IF + END IF + * + RETURN + * + * End of DSYR2K. + * + END diff -cNr octave-2.0.7/libcruft/blas/zhemv.f octave-2.0.8/libcruft/blas/zhemv.f *** octave-2.0.7/libcruft/blas/zhemv.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/blas/zhemv.f Mon Jun 23 14:11:00 1997 *************** *** 0 **** --- 1,269 ---- + * + ************************************************************************ + * + SUBROUTINE ZHEMV ( UPLO, N, ALPHA, A, LDA, X, INCX, + $ BETA, Y, INCY ) + * .. Scalar Arguments .. + COMPLEX*16 ALPHA, BETA + INTEGER INCX, INCY, LDA, N + CHARACTER*1 UPLO + * .. Array Arguments .. + COMPLEX*16 A( LDA, * ), X( * ), Y( * ) + * .. + * + * Purpose + * ======= + * + * ZHEMV performs the matrix-vector operation + * + * y := alpha*A*x + beta*y, + * + * where alpha and beta are scalars, x and y are n element vectors and + * A is an n by n hermitian matrix. + * + * Parameters + * ========== + * + * UPLO - CHARACTER*1. + * On entry, UPLO specifies whether the upper or lower + * triangular part of the array A is to be referenced as + * follows: + * + * UPLO = 'U' or 'u' Only the upper triangular part of A + * is to be referenced. + * + * UPLO = 'L' or 'l' Only the lower triangular part of A + * is to be referenced. + * + * Unchanged on exit. + * + * N - INTEGER. + * On entry, N specifies the order of the matrix A. + * N must be at least zero. + * Unchanged on exit. + * + * ALPHA - COMPLEX*16 . + * On entry, ALPHA specifies the scalar alpha. + * Unchanged on exit. + * + * A - COMPLEX*16 array of DIMENSION ( LDA, n ). + * Before entry with UPLO = 'U' or 'u', the leading n by n + * upper triangular part of the array A must contain the upper + * triangular part of the hermitian matrix and the strictly + * lower triangular part of A is not referenced. + * Before entry with UPLO = 'L' or 'l', the leading n by n + * lower triangular part of the array A must contain the lower + * triangular part of the hermitian matrix and the strictly + * upper triangular part of A is not referenced. + * Note that the imaginary parts of the diagonal elements need + * not be set and are assumed to be zero. + * Unchanged on exit. + * + * LDA - INTEGER. + * On entry, LDA specifies the first dimension of A as declared + * in the calling (sub) program. LDA must be at least + * max( 1, n ). + * Unchanged on exit. + * + * X - COMPLEX*16 array of dimension at least + * ( 1 + ( n - 1 )*abs( INCX ) ). + * Before entry, the incremented array X must contain the n + * element vector x. + * Unchanged on exit. + * + * INCX - INTEGER. + * On entry, INCX specifies the increment for the elements of + * X. INCX must not be zero. + * Unchanged on exit. + * + * BETA - COMPLEX*16 . + * On entry, BETA specifies the scalar beta. When BETA is + * supplied as zero then Y need not be set on input. + * Unchanged on exit. + * + * Y - COMPLEX*16 array of dimension at least + * ( 1 + ( n - 1 )*abs( INCY ) ). + * Before entry, the incremented array Y must contain the n + * element vector y. On exit, Y is overwritten by the updated + * vector y. + * + * INCY - INTEGER. + * On entry, INCY specifies the increment for the elements of + * Y. INCY must not be zero. + * Unchanged on exit. + * + * + * Level 2 Blas routine. + * + * -- Written on 22-October-1986. + * Jack Dongarra, Argonne National Lab. + * Jeremy Du Croz, Nag Central Office. + * Sven Hammarling, Nag Central Office. + * Richard Hanson, Sandia National Labs. + * + * + * .. Parameters .. + COMPLEX*16 ONE + PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) + COMPLEX*16 ZERO + PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) + * .. Local Scalars .. + COMPLEX*16 TEMP1, TEMP2 + INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY + * .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME + * .. External Subroutines .. + EXTERNAL XERBLA + * .. Intrinsic Functions .. + INTRINSIC DCONJG, MAX, DBLE + * .. + * .. Executable Statements .. + * + * Test the input parameters. + * + INFO = 0 + IF ( .NOT.LSAME( UPLO, 'U' ).AND. + $ .NOT.LSAME( UPLO, 'L' ) )THEN + INFO = 1 + ELSE IF( N.LT.0 )THEN + INFO = 2 + ELSE IF( LDA.LT.MAX( 1, N ) )THEN + INFO = 5 + ELSE IF( INCX.EQ.0 )THEN + INFO = 7 + ELSE IF( INCY.EQ.0 )THEN + INFO = 10 + END IF + IF( INFO.NE.0 )THEN + CALL XERBLA( 'ZHEMV ', INFO ) + RETURN + END IF + * + * Quick return if possible. + * + IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) + $ RETURN + * + * Set up the start points in X and Y. + * + IF( INCX.GT.0 )THEN + KX = 1 + ELSE + KX = 1 - ( N - 1 )*INCX + END IF + IF( INCY.GT.0 )THEN + KY = 1 + ELSE + KY = 1 - ( N - 1 )*INCY + END IF + * + * Start the operations. In this version the elements of A are + * accessed sequentially with one pass through the triangular part + * of A. + * + * First form y := beta*y. + * + IF( BETA.NE.ONE )THEN + IF( INCY.EQ.1 )THEN + IF( BETA.EQ.ZERO )THEN + DO 10, I = 1, N + Y( I ) = ZERO + 10 CONTINUE + ELSE + DO 20, I = 1, N + Y( I ) = BETA*Y( I ) + 20 CONTINUE + END IF + ELSE + IY = KY + IF( BETA.EQ.ZERO )THEN + DO 30, I = 1, N + Y( IY ) = ZERO + IY = IY + INCY + 30 CONTINUE + ELSE + DO 40, I = 1, N + Y( IY ) = BETA*Y( IY ) + IY = IY + INCY + 40 CONTINUE + END IF + END IF + END IF + IF( ALPHA.EQ.ZERO ) + $ RETURN + IF( LSAME( UPLO, 'U' ) )THEN + * + * Form y when A is stored in upper triangle. + * + IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN + DO 60, J = 1, N + TEMP1 = ALPHA*X( J ) + TEMP2 = ZERO + DO 50, I = 1, J - 1 + Y( I ) = Y( I ) + TEMP1*A( I, J ) + TEMP2 = TEMP2 + DCONJG( A( I, J ) )*X( I ) + 50 CONTINUE + Y( J ) = Y( J ) + TEMP1*DBLE( A( J, J ) ) + ALPHA*TEMP2 + 60 CONTINUE + ELSE + JX = KX + JY = KY + DO 80, J = 1, N + TEMP1 = ALPHA*X( JX ) + TEMP2 = ZERO + IX = KX + IY = KY + DO 70, I = 1, J - 1 + Y( IY ) = Y( IY ) + TEMP1*A( I, J ) + TEMP2 = TEMP2 + DCONJG( A( I, J ) )*X( IX ) + IX = IX + INCX + IY = IY + INCY + 70 CONTINUE + Y( JY ) = Y( JY ) + TEMP1*DBLE( A( J, J ) ) + ALPHA*TEMP2 + JX = JX + INCX + JY = JY + INCY + 80 CONTINUE + END IF + ELSE + * + * Form y when A is stored in lower triangle. + * + IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN + DO 100, J = 1, N + TEMP1 = ALPHA*X( J ) + TEMP2 = ZERO + Y( J ) = Y( J ) + TEMP1*DBLE( A( J, J ) ) + DO 90, I = J + 1, N + Y( I ) = Y( I ) + TEMP1*A( I, J ) + TEMP2 = TEMP2 + DCONJG( A( I, J ) )*X( I ) + 90 CONTINUE + Y( J ) = Y( J ) + ALPHA*TEMP2 + 100 CONTINUE + ELSE + JX = KX + JY = KY + DO 120, J = 1, N + TEMP1 = ALPHA*X( JX ) + TEMP2 = ZERO + Y( JY ) = Y( JY ) + TEMP1*DBLE( A( J, J ) ) + IX = JX + IY = JY + DO 110, I = J + 1, N + IX = IX + INCX + IY = IY + INCY + Y( IY ) = Y( IY ) + TEMP1*A( I, J ) + TEMP2 = TEMP2 + DCONJG( A( I, J ) )*X( IX ) + 110 CONTINUE + Y( JY ) = Y( JY ) + ALPHA*TEMP2 + JX = JX + INCX + JY = JY + INCY + 120 CONTINUE + END IF + END IF + * + RETURN + * + * End of ZHEMV . + * + END diff -cNr octave-2.0.7/libcruft/blas/zher2.f octave-2.0.8/libcruft/blas/zher2.f *** octave-2.0.7/libcruft/blas/zher2.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/blas/zher2.f Mon Jun 23 14:11:00 1997 *************** *** 0 **** --- 1,252 ---- + * + ************************************************************************ + * + SUBROUTINE ZHER2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA ) + * .. Scalar Arguments .. + COMPLEX*16 ALPHA + INTEGER INCX, INCY, LDA, N + CHARACTER*1 UPLO + * .. Array Arguments .. + COMPLEX*16 A( LDA, * ), X( * ), Y( * ) + * .. + * + * Purpose + * ======= + * + * ZHER2 performs the hermitian rank 2 operation + * + * A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A, + * + * where alpha is a scalar, x and y are n element vectors and A is an n + * by n hermitian matrix. + * + * Parameters + * ========== + * + * UPLO - CHARACTER*1. + * On entry, UPLO specifies whether the upper or lower + * triangular part of the array A is to be referenced as + * follows: + * + * UPLO = 'U' or 'u' Only the upper triangular part of A + * is to be referenced. + * + * UPLO = 'L' or 'l' Only the lower triangular part of A + * is to be referenced. + * + * Unchanged on exit. + * + * N - INTEGER. + * On entry, N specifies the order of the matrix A. + * N must be at least zero. + * Unchanged on exit. + * + * ALPHA - COMPLEX*16 . + * On entry, ALPHA specifies the scalar alpha. + * Unchanged on exit. + * + * X - COMPLEX*16 array of dimension at least + * ( 1 + ( n - 1 )*abs( INCX ) ). + * Before entry, the incremented array X must contain the n + * element vector x. + * Unchanged on exit. + * + * INCX - INTEGER. + * On entry, INCX specifies the increment for the elements of + * X. INCX must not be zero. + * Unchanged on exit. + * + * Y - COMPLEX*16 array of dimension at least + * ( 1 + ( n - 1 )*abs( INCY ) ). + * Before entry, the incremented array Y must contain the n + * element vector y. + * Unchanged on exit. + * + * INCY - INTEGER. + * On entry, INCY specifies the increment for the elements of + * Y. INCY must not be zero. + * Unchanged on exit. + * + * A - COMPLEX*16 array of DIMENSION ( LDA, n ). + * Before entry with UPLO = 'U' or 'u', the leading n by n + * upper triangular part of the array A must contain the upper + * triangular part of the hermitian matrix and the strictly + * lower triangular part of A is not referenced. On exit, the + * upper triangular part of the array A is overwritten by the + * upper triangular part of the updated matrix. + * Before entry with UPLO = 'L' or 'l', the leading n by n + * lower triangular part of the array A must contain the lower + * triangular part of the hermitian matrix and the strictly + * upper triangular part of A is not referenced. On exit, the + * lower triangular part of the array A is overwritten by the + * lower triangular part of the updated matrix. + * Note that the imaginary parts of the diagonal elements need + * not be set, they are assumed to be zero, and on exit they + * are set to zero. + * + * LDA - INTEGER. + * On entry, LDA specifies the first dimension of A as declared + * in the calling (sub) program. LDA must be at least + * max( 1, n ). + * Unchanged on exit. + * + * + * Level 2 Blas routine. + * + * -- Written on 22-October-1986. + * Jack Dongarra, Argonne National Lab. + * Jeremy Du Croz, Nag Central Office. + * Sven Hammarling, Nag Central Office. + * Richard Hanson, Sandia National Labs. + * + * + * .. Parameters .. + COMPLEX*16 ZERO + PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) + * .. Local Scalars .. + COMPLEX*16 TEMP1, TEMP2 + INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY + * .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME + * .. External Subroutines .. + EXTERNAL XERBLA + * .. Intrinsic Functions .. + INTRINSIC DCONJG, MAX, DBLE + * .. + * .. Executable Statements .. + * + * Test the input parameters. + * + INFO = 0 + IF ( .NOT.LSAME( UPLO, 'U' ).AND. + $ .NOT.LSAME( UPLO, 'L' ) )THEN + INFO = 1 + ELSE IF( N.LT.0 )THEN + INFO = 2 + ELSE IF( INCX.EQ.0 )THEN + INFO = 5 + ELSE IF( INCY.EQ.0 )THEN + INFO = 7 + ELSE IF( LDA.LT.MAX( 1, N ) )THEN + INFO = 9 + END IF + IF( INFO.NE.0 )THEN + CALL XERBLA( 'ZHER2 ', INFO ) + RETURN + END IF + * + * Quick return if possible. + * + IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) + $ RETURN + * + * Set up the start points in X and Y if the increments are not both + * unity. + * + IF( ( INCX.NE.1 ).OR.( INCY.NE.1 ) )THEN + IF( INCX.GT.0 )THEN + KX = 1 + ELSE + KX = 1 - ( N - 1 )*INCX + END IF + IF( INCY.GT.0 )THEN + KY = 1 + ELSE + KY = 1 - ( N - 1 )*INCY + END IF + JX = KX + JY = KY + END IF + * + * Start the operations. In this version the elements of A are + * accessed sequentially with one pass through the triangular part + * of A. + * + IF( LSAME( UPLO, 'U' ) )THEN + * + * Form A when A is stored in the upper triangle. + * + IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN + DO 20, J = 1, N + IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN + TEMP1 = ALPHA*DCONJG( Y( J ) ) + TEMP2 = DCONJG( ALPHA*X( J ) ) + DO 10, I = 1, J - 1 + A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2 + 10 CONTINUE + A( J, J ) = DBLE( A( J, J ) ) + + $ DBLE( X( J )*TEMP1 + Y( J )*TEMP2 ) + ELSE + A( J, J ) = DBLE( A( J, J ) ) + END IF + 20 CONTINUE + ELSE + DO 40, J = 1, N + IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN + TEMP1 = ALPHA*DCONJG( Y( JY ) ) + TEMP2 = DCONJG( ALPHA*X( JX ) ) + IX = KX + IY = KY + DO 30, I = 1, J - 1 + A( I, J ) = A( I, J ) + X( IX )*TEMP1 + $ + Y( IY )*TEMP2 + IX = IX + INCX + IY = IY + INCY + 30 CONTINUE + A( J, J ) = DBLE( A( J, J ) ) + + $ DBLE( X( JX )*TEMP1 + Y( JY )*TEMP2 ) + ELSE + A( J, J ) = DBLE( A( J, J ) ) + END IF + JX = JX + INCX + JY = JY + INCY + 40 CONTINUE + END IF + ELSE + * + * Form A when A is stored in the lower triangle. + * + IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN + DO 60, J = 1, N + IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN + TEMP1 = ALPHA*DCONJG( Y( J ) ) + TEMP2 = DCONJG( ALPHA*X( J ) ) + A( J, J ) = DBLE( A( J, J ) ) + + $ DBLE( X( J )*TEMP1 + Y( J )*TEMP2 ) + DO 50, I = J + 1, N + A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2 + 50 CONTINUE + ELSE + A( J, J ) = DBLE( A( J, J ) ) + END IF + 60 CONTINUE + ELSE + DO 80, J = 1, N + IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN + TEMP1 = ALPHA*DCONJG( Y( JY ) ) + TEMP2 = DCONJG( ALPHA*X( JX ) ) + A( J, J ) = DBLE( A( J, J ) ) + + $ DBLE( X( JX )*TEMP1 + Y( JY )*TEMP2 ) + IX = JX + IY = JY + DO 70, I = J + 1, N + IX = IX + INCX + IY = IY + INCY + A( I, J ) = A( I, J ) + X( IX )*TEMP1 + $ + Y( IY )*TEMP2 + 70 CONTINUE + ELSE + A( J, J ) = DBLE( A( J, J ) ) + END IF + JX = JX + INCX + JY = JY + INCY + 80 CONTINUE + END IF + END IF + * + RETURN + * + * End of ZHER2 . + * + END diff -cNr octave-2.0.7/libcruft/blas/zher2k.f octave-2.0.8/libcruft/blas/zher2k.f *** octave-2.0.7/libcruft/blas/zher2k.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/blas/zher2k.f Mon Jun 23 14:11:00 1997 *************** *** 0 **** --- 1,372 ---- + SUBROUTINE ZHER2K( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, BETA, + $ C, LDC ) + * .. Scalar Arguments .. + CHARACTER TRANS, UPLO + INTEGER K, LDA, LDB, LDC, N + DOUBLE PRECISION BETA + COMPLEX*16 ALPHA + * .. + * .. Array Arguments .. + COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ) + * .. + * + * Purpose + * ======= + * + * ZHER2K performs one of the hermitian rank 2k operations + * + * C := alpha*A*conjg( B' ) + conjg( alpha )*B*conjg( A' ) + beta*C, + * + * or + * + * C := alpha*conjg( A' )*B + conjg( alpha )*conjg( B' )*A + beta*C, + * + * where alpha and beta are scalars with beta real, C is an n by n + * hermitian matrix and A and B are n by k matrices in the first case + * and k by n matrices in the second case. + * + * Parameters + * ========== + * + * UPLO - CHARACTER*1. + * On entry, UPLO specifies whether the upper or lower + * triangular part of the array C is to be referenced as + * follows: + * + * UPLO = 'U' or 'u' Only the upper triangular part of C + * is to be referenced. + * + * UPLO = 'L' or 'l' Only the lower triangular part of C + * is to be referenced. + * + * Unchanged on exit. + * + * TRANS - CHARACTER*1. + * On entry, TRANS specifies the operation to be performed as + * follows: + * + * TRANS = 'N' or 'n' C := alpha*A*conjg( B' ) + + * conjg( alpha )*B*conjg( A' ) + + * beta*C. + * + * TRANS = 'C' or 'c' C := alpha*conjg( A' )*B + + * conjg( alpha )*conjg( B' )*A + + * beta*C. + * + * Unchanged on exit. + * + * N - INTEGER. + * On entry, N specifies the order of the matrix C. N must be + * at least zero. + * Unchanged on exit. + * + * K - INTEGER. + * On entry with TRANS = 'N' or 'n', K specifies the number + * of columns of the matrices A and B, and on entry with + * TRANS = 'C' or 'c', K specifies the number of rows of the + * matrices A and B. K must be at least zero. + * Unchanged on exit. + * + * ALPHA - COMPLEX*16 . + * On entry, ALPHA specifies the scalar alpha. + * Unchanged on exit. + * + * A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is + * k when TRANS = 'N' or 'n', and is n otherwise. + * Before entry with TRANS = 'N' or 'n', the leading n by k + * part of the array A must contain the matrix A, otherwise + * the leading k by n part of the array A must contain the + * matrix A. + * Unchanged on exit. + * + * LDA - INTEGER. + * On entry, LDA specifies the first dimension of A as declared + * in the calling (sub) program. When TRANS = 'N' or 'n' + * then LDA must be at least max( 1, n ), otherwise LDA must + * be at least max( 1, k ). + * Unchanged on exit. + * + * B - COMPLEX*16 array of DIMENSION ( LDB, kb ), where kb is + * k when TRANS = 'N' or 'n', and is n otherwise. + * Before entry with TRANS = 'N' or 'n', the leading n by k + * part of the array B must contain the matrix B, otherwise + * the leading k by n part of the array B must contain the + * matrix B. + * Unchanged on exit. + * + * LDB - INTEGER. + * On entry, LDB specifies the first dimension of B as declared + * in the calling (sub) program. When TRANS = 'N' or 'n' + * then LDB must be at least max( 1, n ), otherwise LDB must + * be at least max( 1, k ). + * Unchanged on exit. + * + * BETA - DOUBLE PRECISION . + * On entry, BETA specifies the scalar beta. + * Unchanged on exit. + * + * C - COMPLEX*16 array of DIMENSION ( LDC, n ). + * Before entry with UPLO = 'U' or 'u', the leading n by n + * upper triangular part of the array C must contain the upper + * triangular part of the hermitian matrix and the strictly + * lower triangular part of C is not referenced. On exit, the + * upper triangular part of the array C is overwritten by the + * upper triangular part of the updated matrix. + * Before entry with UPLO = 'L' or 'l', the leading n by n + * lower triangular part of the array C must contain the lower + * triangular part of the hermitian matrix and the strictly + * upper triangular part of C is not referenced. On exit, the + * lower triangular part of the array C is overwritten by the + * lower triangular part of the updated matrix. + * Note that the imaginary parts of the diagonal elements need + * not be set, they are assumed to be zero, and on exit they + * are set to zero. + * + * LDC - INTEGER. + * On entry, LDC specifies the first dimension of C as declared + * in the calling (sub) program. LDC must be at least + * max( 1, n ). + * Unchanged on exit. + * + * + * Level 3 Blas routine. + * + * -- Written on 8-February-1989. + * Jack Dongarra, Argonne National Laboratory. + * Iain Duff, AERE Harwell. + * Jeremy Du Croz, Numerical Algorithms Group Ltd. + * Sven Hammarling, Numerical Algorithms Group Ltd. + * + * -- Modified 8-Nov-93 to set C(J,J) to DBLE( C(J,J) ) when BETA = 1. + * Ed Anderson, Cray Research Inc. + * + * + * .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME + * .. + * .. External Subroutines .. + EXTERNAL XERBLA + * .. + * .. Intrinsic Functions .. + INTRINSIC DBLE, DCONJG, MAX + * .. + * .. Local Scalars .. + LOGICAL UPPER + INTEGER I, INFO, J, L, NROWA + COMPLEX*16 TEMP1, TEMP2 + * .. + * .. Parameters .. + DOUBLE PRECISION ONE + PARAMETER ( ONE = 1.0D+0 ) + COMPLEX*16 ZERO + PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) + * .. + * .. Executable Statements .. + * + * Test the input parameters. + * + IF( LSAME( TRANS, 'N' ) ) THEN + NROWA = N + ELSE + NROWA = K + END IF + UPPER = LSAME( UPLO, 'U' ) + * + INFO = 0 + IF( ( .NOT.UPPER ) .AND. ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN + INFO = 1 + ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ) .AND. + $ ( .NOT.LSAME( TRANS, 'C' ) ) ) THEN + INFO = 2 + ELSE IF( N.LT.0 ) THEN + INFO = 3 + ELSE IF( K.LT.0 ) THEN + INFO = 4 + ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN + INFO = 7 + ELSE IF( LDB.LT.MAX( 1, NROWA ) ) THEN + INFO = 9 + ELSE IF( LDC.LT.MAX( 1, N ) ) THEN + INFO = 12 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZHER2K', INFO ) + RETURN + END IF + * + * Quick return if possible. + * + IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND. + $ ( BETA.EQ.ONE ) ) )RETURN + * + * And when alpha.eq.zero. + * + IF( ALPHA.EQ.ZERO ) THEN + IF( UPPER ) THEN + IF( BETA.EQ.DBLE( ZERO ) ) THEN + DO 20 J = 1, N + DO 10 I = 1, J + C( I, J ) = ZERO + 10 CONTINUE + 20 CONTINUE + ELSE + DO 40 J = 1, N + DO 30 I = 1, J - 1 + C( I, J ) = BETA*C( I, J ) + 30 CONTINUE + C( J, J ) = BETA*DBLE( C( J, J ) ) + 40 CONTINUE + END IF + ELSE + IF( BETA.EQ.DBLE( ZERO ) ) THEN + DO 60 J = 1, N + DO 50 I = J, N + C( I, J ) = ZERO + 50 CONTINUE + 60 CONTINUE + ELSE + DO 80 J = 1, N + C( J, J ) = BETA*DBLE( C( J, J ) ) + DO 70 I = J + 1, N + C( I, J ) = BETA*C( I, J ) + 70 CONTINUE + 80 CONTINUE + END IF + END IF + RETURN + END IF + * + * Start the operations. + * + IF( LSAME( TRANS, 'N' ) ) THEN + * + * Form C := alpha*A*conjg( B' ) + conjg( alpha )*B*conjg( A' ) + + * C. + * + IF( UPPER ) THEN + DO 130 J = 1, N + IF( BETA.EQ.DBLE( ZERO ) ) THEN + DO 90 I = 1, J + C( I, J ) = ZERO + 90 CONTINUE + ELSE IF( BETA.NE.ONE ) THEN + DO 100 I = 1, J - 1 + C( I, J ) = BETA*C( I, J ) + 100 CONTINUE + C( J, J ) = BETA*DBLE( C( J, J ) ) + ELSE + C( J, J ) = DBLE( C( J, J ) ) + END IF + DO 120 L = 1, K + IF( ( A( J, L ).NE.ZERO ) .OR. ( B( J, L ).NE.ZERO ) ) + $ THEN + TEMP1 = ALPHA*DCONJG( B( J, L ) ) + TEMP2 = DCONJG( ALPHA*A( J, L ) ) + DO 110 I = 1, J - 1 + C( I, J ) = C( I, J ) + A( I, L )*TEMP1 + + $ B( I, L )*TEMP2 + 110 CONTINUE + C( J, J ) = DBLE( C( J, J ) ) + + $ DBLE( A( J, L )*TEMP1+B( J, L )*TEMP2 ) + END IF + 120 CONTINUE + 130 CONTINUE + ELSE + DO 180 J = 1, N + IF( BETA.EQ.DBLE( ZERO ) ) THEN + DO 140 I = J, N + C( I, J ) = ZERO + 140 CONTINUE + ELSE IF( BETA.NE.ONE ) THEN + DO 150 I = J + 1, N + C( I, J ) = BETA*C( I, J ) + 150 CONTINUE + C( J, J ) = BETA*DBLE( C( J, J ) ) + ELSE + C( J, J ) = DBLE( C( J, J ) ) + END IF + DO 170 L = 1, K + IF( ( A( J, L ).NE.ZERO ) .OR. ( B( J, L ).NE.ZERO ) ) + $ THEN + TEMP1 = ALPHA*DCONJG( B( J, L ) ) + TEMP2 = DCONJG( ALPHA*A( J, L ) ) + DO 160 I = J + 1, N + C( I, J ) = C( I, J ) + A( I, L )*TEMP1 + + $ B( I, L )*TEMP2 + 160 CONTINUE + C( J, J ) = DBLE( C( J, J ) ) + + $ DBLE( A( J, L )*TEMP1+B( J, L )*TEMP2 ) + END IF + 170 CONTINUE + 180 CONTINUE + END IF + ELSE + * + * Form C := alpha*conjg( A' )*B + conjg( alpha )*conjg( B' )*A + + * C. + * + IF( UPPER ) THEN + DO 210 J = 1, N + DO 200 I = 1, J + TEMP1 = ZERO + TEMP2 = ZERO + DO 190 L = 1, K + TEMP1 = TEMP1 + DCONJG( A( L, I ) )*B( L, J ) + TEMP2 = TEMP2 + DCONJG( B( L, I ) )*A( L, J ) + 190 CONTINUE + IF( I.EQ.J ) THEN + IF( BETA.EQ.DBLE( ZERO ) ) THEN + C( J, J ) = DBLE( ALPHA*TEMP1+DCONJG( ALPHA )* + $ TEMP2 ) + ELSE + C( J, J ) = BETA*DBLE( C( J, J ) ) + + $ DBLE( ALPHA*TEMP1+DCONJG( ALPHA )* + $ TEMP2 ) + END IF + ELSE + IF( BETA.EQ.DBLE( ZERO ) ) THEN + C( I, J ) = ALPHA*TEMP1 + DCONJG( ALPHA )*TEMP2 + ELSE + C( I, J ) = BETA*C( I, J ) + ALPHA*TEMP1 + + $ DCONJG( ALPHA )*TEMP2 + END IF + END IF + 200 CONTINUE + 210 CONTINUE + ELSE + DO 240 J = 1, N + DO 230 I = J, N + TEMP1 = ZERO + TEMP2 = ZERO + DO 220 L = 1, K + TEMP1 = TEMP1 + DCONJG( A( L, I ) )*B( L, J ) + TEMP2 = TEMP2 + DCONJG( B( L, I ) )*A( L, J ) + 220 CONTINUE + IF( I.EQ.J ) THEN + IF( BETA.EQ.DBLE( ZERO ) ) THEN + C( J, J ) = DBLE( ALPHA*TEMP1+DCONJG( ALPHA )* + $ TEMP2 ) + ELSE + C( J, J ) = BETA*DBLE( C( J, J ) ) + + $ DBLE( ALPHA*TEMP1+DCONJG( ALPHA )* + $ TEMP2 ) + END IF + ELSE + IF( BETA.EQ.DBLE( ZERO ) ) THEN + C( I, J ) = ALPHA*TEMP1 + DCONJG( ALPHA )*TEMP2 + ELSE + C( I, J ) = BETA*C( I, J ) + ALPHA*TEMP1 + + $ DCONJG( ALPHA )*TEMP2 + END IF + END IF + 230 CONTINUE + 240 CONTINUE + END IF + END IF + * + RETURN + * + * End of ZHER2K. + * + END diff -cNr octave-2.0.7/libcruft/lapack/dlae2.f octave-2.0.8/libcruft/lapack/dlae2.f *** octave-2.0.7/libcruft/lapack/dlae2.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/dlae2.f Mon Jun 23 14:00:35 1997 *************** *** 0 **** --- 1,124 ---- + SUBROUTINE DLAE2( A, B, C, RT1, RT2 ) + * + * -- LAPACK auxiliary routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * October 31, 1992 + * + * .. Scalar Arguments .. + DOUBLE PRECISION A, B, C, RT1, RT2 + * .. + * + * Purpose + * ======= + * + * DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix + * [ A B ] + * [ B C ]. + * On return, RT1 is the eigenvalue of larger absolute value, and RT2 + * is the eigenvalue of smaller absolute value. + * + * Arguments + * ========= + * + * A (input) DOUBLE PRECISION + * The (1,1) element of the 2-by-2 matrix. + * + * B (input) DOUBLE PRECISION + * The (1,2) and (2,1) elements of the 2-by-2 matrix. + * + * C (input) DOUBLE PRECISION + * The (2,2) element of the 2-by-2 matrix. + * + * RT1 (output) DOUBLE PRECISION + * The eigenvalue of larger absolute value. + * + * RT2 (output) DOUBLE PRECISION + * The eigenvalue of smaller absolute value. + * + * Further Details + * =============== + * + * RT1 is accurate to a few ulps barring over/underflow. + * + * RT2 may be inaccurate if there is massive cancellation in the + * determinant A*C-B*B; higher precision or correctly rounded or + * correctly truncated arithmetic would be needed to compute RT2 + * accurately in all cases. + * + * Overflow is possible only if RT1 is within a factor of 5 of overflow. + * Underflow is harmless if the input data is 0 or exceeds + * underflow_threshold / macheps. + * + * ===================================================================== + * + * .. Parameters .. + DOUBLE PRECISION ONE + PARAMETER ( ONE = 1.0D0 ) + DOUBLE PRECISION TWO + PARAMETER ( TWO = 2.0D0 ) + DOUBLE PRECISION ZERO + PARAMETER ( ZERO = 0.0D0 ) + DOUBLE PRECISION HALF + PARAMETER ( HALF = 0.5D0 ) + * .. + * .. Local Scalars .. + DOUBLE PRECISION AB, ACMN, ACMX, ADF, DF, RT, SM, TB + * .. + * .. Intrinsic Functions .. + INTRINSIC ABS, SQRT + * .. + * .. Executable Statements .. + * + * Compute the eigenvalues + * + SM = A + C + DF = A - C + ADF = ABS( DF ) + TB = B + B + AB = ABS( TB ) + IF( ABS( A ).GT.ABS( C ) ) THEN + ACMX = A + ACMN = C + ELSE + ACMX = C + ACMN = A + END IF + IF( ADF.GT.AB ) THEN + RT = ADF*SQRT( ONE+( AB / ADF )**2 ) + ELSE IF( ADF.LT.AB ) THEN + RT = AB*SQRT( ONE+( ADF / AB )**2 ) + ELSE + * + * Includes case AB=ADF=0 + * + RT = AB*SQRT( TWO ) + END IF + IF( SM.LT.ZERO ) THEN + RT1 = HALF*( SM-RT ) + * + * Order of execution important. + * To get fully accurate smaller eigenvalue, + * next line needs to be executed in higher precision. + * + RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B + ELSE IF( SM.GT.ZERO ) THEN + RT1 = HALF*( SM+RT ) + * + * Order of execution important. + * To get fully accurate smaller eigenvalue, + * next line needs to be executed in higher precision. + * + RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B + ELSE + * + * Includes case RT1 = RT2 = 0 + * + RT1 = HALF*RT + RT2 = -HALF*RT + END IF + RETURN + * + * End of DLAE2 + * + END diff -cNr octave-2.0.7/libcruft/lapack/dlaev2.f octave-2.0.8/libcruft/lapack/dlaev2.f *** octave-2.0.7/libcruft/lapack/dlaev2.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/dlaev2.f Mon Jun 23 14:00:35 1997 *************** *** 0 **** --- 1,170 ---- + SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 ) + * + * -- LAPACK auxiliary routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * October 31, 1992 + * + * .. Scalar Arguments .. + DOUBLE PRECISION A, B, C, CS1, RT1, RT2, SN1 + * .. + * + * Purpose + * ======= + * + * DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix + * [ A B ] + * [ B C ]. + * On return, RT1 is the eigenvalue of larger absolute value, RT2 is the + * eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right + * eigenvector for RT1, giving the decomposition + * + * [ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ] + * [-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ]. + * + * Arguments + * ========= + * + * A (input) DOUBLE PRECISION + * The (1,1) element of the 2-by-2 matrix. + * + * B (input) DOUBLE PRECISION + * The (1,2) element and the conjugate of the (2,1) element of + * the 2-by-2 matrix. + * + * C (input) DOUBLE PRECISION + * The (2,2) element of the 2-by-2 matrix. + * + * RT1 (output) DOUBLE PRECISION + * The eigenvalue of larger absolute value. + * + * RT2 (output) DOUBLE PRECISION + * The eigenvalue of smaller absolute value. + * + * CS1 (output) DOUBLE PRECISION + * SN1 (output) DOUBLE PRECISION + * The vector (CS1, SN1) is a unit right eigenvector for RT1. + * + * Further Details + * =============== + * + * RT1 is accurate to a few ulps barring over/underflow. + * + * RT2 may be inaccurate if there is massive cancellation in the + * determinant A*C-B*B; higher precision or correctly rounded or + * correctly truncated arithmetic would be needed to compute RT2 + * accurately in all cases. + * + * CS1 and SN1 are accurate to a few ulps barring over/underflow. + * + * Overflow is possible only if RT1 is within a factor of 5 of overflow. + * Underflow is harmless if the input data is 0 or exceeds + * underflow_threshold / macheps. + * + * ===================================================================== + * + * .. Parameters .. + DOUBLE PRECISION ONE + PARAMETER ( ONE = 1.0D0 ) + DOUBLE PRECISION TWO + PARAMETER ( TWO = 2.0D0 ) + DOUBLE PRECISION ZERO + PARAMETER ( ZERO = 0.0D0 ) + DOUBLE PRECISION HALF + PARAMETER ( HALF = 0.5D0 ) + * .. + * .. Local Scalars .. + INTEGER SGN1, SGN2 + DOUBLE PRECISION AB, ACMN, ACMX, ACS, ADF, CS, CT, DF, RT, SM, + $ TB, TN + * .. + * .. Intrinsic Functions .. + INTRINSIC ABS, SQRT + * .. + * .. Executable Statements .. + * + * Compute the eigenvalues + * + SM = A + C + DF = A - C + ADF = ABS( DF ) + TB = B + B + AB = ABS( TB ) + IF( ABS( A ).GT.ABS( C ) ) THEN + ACMX = A + ACMN = C + ELSE + ACMX = C + ACMN = A + END IF + IF( ADF.GT.AB ) THEN + RT = ADF*SQRT( ONE+( AB / ADF )**2 ) + ELSE IF( ADF.LT.AB ) THEN + RT = AB*SQRT( ONE+( ADF / AB )**2 ) + ELSE + * + * Includes case AB=ADF=0 + * + RT = AB*SQRT( TWO ) + END IF + IF( SM.LT.ZERO ) THEN + RT1 = HALF*( SM-RT ) + SGN1 = -1 + * + * Order of execution important. + * To get fully accurate smaller eigenvalue, + * next line needs to be executed in higher precision. + * + RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B + ELSE IF( SM.GT.ZERO ) THEN + RT1 = HALF*( SM+RT ) + SGN1 = 1 + * + * Order of execution important. + * To get fully accurate smaller eigenvalue, + * next line needs to be executed in higher precision. + * + RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B + ELSE + * + * Includes case RT1 = RT2 = 0 + * + RT1 = HALF*RT + RT2 = -HALF*RT + SGN1 = 1 + END IF + * + * Compute the eigenvector + * + IF( DF.GE.ZERO ) THEN + CS = DF + RT + SGN2 = 1 + ELSE + CS = DF - RT + SGN2 = -1 + END IF + ACS = ABS( CS ) + IF( ACS.GT.AB ) THEN + CT = -TB / CS + SN1 = ONE / SQRT( ONE+CT*CT ) + CS1 = CT*SN1 + ELSE + IF( AB.EQ.ZERO ) THEN + CS1 = ONE + SN1 = ZERO + ELSE + TN = -CS / TB + CS1 = ONE / SQRT( ONE+TN*TN ) + SN1 = TN*CS1 + END IF + END IF + IF( SGN1.EQ.SGN2 ) THEN + TN = CS1 + CS1 = -SN1 + SN1 = TN + END IF + RETURN + * + * End of DLAEV2 + * + END diff -cNr octave-2.0.7/libcruft/lapack/dlanst.f octave-2.0.8/libcruft/lapack/dlanst.f *** octave-2.0.7/libcruft/lapack/dlanst.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/dlanst.f Mon Jun 23 14:00:35 1997 *************** *** 0 **** --- 1,125 ---- + DOUBLE PRECISION FUNCTION DLANST( NORM, N, D, E ) + * + * -- LAPACK auxiliary routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * February 29, 1992 + * + * .. Scalar Arguments .. + CHARACTER NORM + INTEGER N + * .. + * .. Array Arguments .. + DOUBLE PRECISION D( * ), E( * ) + * .. + * + * Purpose + * ======= + * + * DLANST returns the value of the one norm, or the Frobenius norm, or + * the infinity norm, or the element of largest absolute value of a + * real symmetric tridiagonal matrix A. + * + * Description + * =========== + * + * DLANST returns the value + * + * DLANST = ( max(abs(A(i,j))), NORM = 'M' or 'm' + * ( + * ( norm1(A), NORM = '1', 'O' or 'o' + * ( + * ( normI(A), NORM = 'I' or 'i' + * ( + * ( normF(A), NORM = 'F', 'f', 'E' or 'e' + * + * where norm1 denotes the one norm of a matrix (maximum column sum), + * normI denotes the infinity norm of a matrix (maximum row sum) and + * normF denotes the Frobenius norm of a matrix (square root of sum of + * squares). Note that max(abs(A(i,j))) is not a matrix norm. + * + * Arguments + * ========= + * + * NORM (input) CHARACTER*1 + * Specifies the value to be returned in DLANST as described + * above. + * + * N (input) INTEGER + * The order of the matrix A. N >= 0. When N = 0, DLANST is + * set to zero. + * + * D (input) DOUBLE PRECISION array, dimension (N) + * The diagonal elements of A. + * + * E (input) DOUBLE PRECISION array, dimension (N-1) + * The (n-1) sub-diagonal or super-diagonal elements of A. + * + * ===================================================================== + * + * .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) + * .. + * .. Local Scalars .. + INTEGER I + DOUBLE PRECISION ANORM, SCALE, SUM + * .. + * .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME + * .. + * .. External Subroutines .. + EXTERNAL DLASSQ + * .. + * .. Intrinsic Functions .. + INTRINSIC ABS, MAX, SQRT + * .. + * .. Executable Statements .. + * + IF( N.LE.0 ) THEN + ANORM = ZERO + ELSE IF( LSAME( NORM, 'M' ) ) THEN + * + * Find max(abs(A(i,j))). + * + ANORM = ABS( D( N ) ) + DO 10 I = 1, N - 1 + ANORM = MAX( ANORM, ABS( D( I ) ) ) + ANORM = MAX( ANORM, ABS( E( I ) ) ) + 10 CONTINUE + ELSE IF( LSAME( NORM, 'O' ) .OR. NORM.EQ.'1' .OR. + $ LSAME( NORM, 'I' ) ) THEN + * + * Find norm1(A). + * + IF( N.EQ.1 ) THEN + ANORM = ABS( D( 1 ) ) + ELSE + ANORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ), + $ ABS( E( N-1 ) )+ABS( D( N ) ) ) + DO 20 I = 2, N - 1 + ANORM = MAX( ANORM, ABS( D( I ) )+ABS( E( I ) )+ + $ ABS( E( I-1 ) ) ) + 20 CONTINUE + END IF + ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN + * + * Find normF(A). + * + SCALE = ZERO + SUM = ONE + IF( N.GT.1 ) THEN + CALL DLASSQ( N-1, E, 1, SCALE, SUM ) + SUM = 2*SUM + END IF + CALL DLASSQ( N, D, 1, SCALE, SUM ) + ANORM = SCALE*SQRT( SUM ) + END IF + * + DLANST = ANORM + RETURN + * + * End of DLANST + * + END diff -cNr octave-2.0.7/libcruft/lapack/dlansy.f octave-2.0.8/libcruft/lapack/dlansy.f *** octave-2.0.7/libcruft/lapack/dlansy.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/dlansy.f Mon Jun 23 14:00:35 1997 *************** *** 0 **** --- 1,174 ---- + DOUBLE PRECISION FUNCTION DLANSY( NORM, UPLO, N, A, LDA, WORK ) + * + * -- LAPACK auxiliary routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * October 31, 1992 + * + * .. Scalar Arguments .. + CHARACTER NORM, UPLO + INTEGER LDA, N + * .. + * .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), WORK( * ) + * .. + * + * Purpose + * ======= + * + * DLANSY returns the value of the one norm, or the Frobenius norm, or + * the infinity norm, or the element of largest absolute value of a + * real symmetric matrix A. + * + * Description + * =========== + * + * DLANSY returns the value + * + * DLANSY = ( max(abs(A(i,j))), NORM = 'M' or 'm' + * ( + * ( norm1(A), NORM = '1', 'O' or 'o' + * ( + * ( normI(A), NORM = 'I' or 'i' + * ( + * ( normF(A), NORM = 'F', 'f', 'E' or 'e' + * + * where norm1 denotes the one norm of a matrix (maximum column sum), + * normI denotes the infinity norm of a matrix (maximum row sum) and + * normF denotes the Frobenius norm of a matrix (square root of sum of + * squares). Note that max(abs(A(i,j))) is not a matrix norm. + * + * Arguments + * ========= + * + * NORM (input) CHARACTER*1 + * Specifies the value to be returned in DLANSY as described + * above. + * + * UPLO (input) CHARACTER*1 + * Specifies whether the upper or lower triangular part of the + * symmetric matrix A is to be referenced. + * = 'U': Upper triangular part of A is referenced + * = 'L': Lower triangular part of A is referenced + * + * N (input) INTEGER + * The order of the matrix A. N >= 0. When N = 0, DLANSY is + * set to zero. + * + * A (input) DOUBLE PRECISION array, dimension (LDA,N) + * The symmetric matrix A. If UPLO = 'U', the leading n by n + * upper triangular part of A contains the upper triangular part + * of the matrix A, and the strictly lower triangular part of A + * is not referenced. If UPLO = 'L', the leading n by n lower + * triangular part of A contains the lower triangular part of + * the matrix A, and the strictly upper triangular part of A is + * not referenced. + * + * LDA (input) INTEGER + * The leading dimension of the array A. LDA >= max(N,1). + * + * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK), + * where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, + * WORK is not referenced. + * + * ===================================================================== + * + * .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) + * .. + * .. Local Scalars .. + INTEGER I, J + DOUBLE PRECISION ABSA, SCALE, SUM, VALUE + * .. + * .. External Subroutines .. + EXTERNAL DLASSQ + * .. + * .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME + * .. + * .. Intrinsic Functions .. + INTRINSIC ABS, MAX, SQRT + * .. + * .. Executable Statements .. + * + IF( N.EQ.0 ) THEN + VALUE = ZERO + ELSE IF( LSAME( NORM, 'M' ) ) THEN + * + * Find max(abs(A(i,j))). + * + VALUE = ZERO + IF( LSAME( UPLO, 'U' ) ) THEN + DO 20 J = 1, N + DO 10 I = 1, J + VALUE = MAX( VALUE, ABS( A( I, J ) ) ) + 10 CONTINUE + 20 CONTINUE + ELSE + DO 40 J = 1, N + DO 30 I = J, N + VALUE = MAX( VALUE, ABS( A( I, J ) ) ) + 30 CONTINUE + 40 CONTINUE + END IF + ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR. + $ ( NORM.EQ.'1' ) ) THEN + * + * Find normI(A) ( = norm1(A), since A is symmetric). + * + VALUE = ZERO + IF( LSAME( UPLO, 'U' ) ) THEN + DO 60 J = 1, N + SUM = ZERO + DO 50 I = 1, J - 1 + ABSA = ABS( A( I, J ) ) + SUM = SUM + ABSA + WORK( I ) = WORK( I ) + ABSA + 50 CONTINUE + WORK( J ) = SUM + ABS( A( J, J ) ) + 60 CONTINUE + DO 70 I = 1, N + VALUE = MAX( VALUE, WORK( I ) ) + 70 CONTINUE + ELSE + DO 80 I = 1, N + WORK( I ) = ZERO + 80 CONTINUE + DO 100 J = 1, N + SUM = WORK( J ) + ABS( A( J, J ) ) + DO 90 I = J + 1, N + ABSA = ABS( A( I, J ) ) + SUM = SUM + ABSA + WORK( I ) = WORK( I ) + ABSA + 90 CONTINUE + VALUE = MAX( VALUE, SUM ) + 100 CONTINUE + END IF + ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN + * + * Find normF(A). + * + SCALE = ZERO + SUM = ONE + IF( LSAME( UPLO, 'U' ) ) THEN + DO 110 J = 2, N + CALL DLASSQ( J-1, A( 1, J ), 1, SCALE, SUM ) + 110 CONTINUE + ELSE + DO 120 J = 1, N - 1 + CALL DLASSQ( N-J, A( J+1, J ), 1, SCALE, SUM ) + 120 CONTINUE + END IF + SUM = 2*SUM + CALL DLASSQ( N, A, LDA+1, SCALE, SUM ) + VALUE = SCALE*SQRT( SUM ) + END IF + * + DLANSY = VALUE + RETURN + * + * End of DLANSY + * + END diff -cNr octave-2.0.7/libcruft/lapack/dlatrd.f octave-2.0.8/libcruft/lapack/dlatrd.f *** octave-2.0.7/libcruft/lapack/dlatrd.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/dlatrd.f Mon Jun 23 14:00:08 1997 *************** *** 0 **** --- 1,259 ---- + SUBROUTINE DLATRD( UPLO, N, NB, A, LDA, E, TAU, W, LDW ) + * + * -- LAPACK auxiliary routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * October 31, 1992 + * + * .. Scalar Arguments .. + CHARACTER UPLO + INTEGER LDA, LDW, N, NB + * .. + * .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), E( * ), TAU( * ), W( LDW, * ) + * .. + * + * Purpose + * ======= + * + * DLATRD reduces NB rows and columns of a real symmetric matrix A to + * symmetric tridiagonal form by an orthogonal similarity + * transformation Q' * A * Q, and returns the matrices V and W which are + * needed to apply the transformation to the unreduced part of A. + * + * If UPLO = 'U', DLATRD reduces the last NB rows and columns of a + * matrix, of which the upper triangle is supplied; + * if UPLO = 'L', DLATRD reduces the first NB rows and columns of a + * matrix, of which the lower triangle is supplied. + * + * This is an auxiliary routine called by DSYTRD. + * + * Arguments + * ========= + * + * UPLO (input) CHARACTER + * Specifies whether the upper or lower triangular part of the + * symmetric matrix A is stored: + * = 'U': Upper triangular + * = 'L': Lower triangular + * + * N (input) INTEGER + * The order of the matrix A. + * + * NB (input) INTEGER + * The number of rows and columns to be reduced. + * + * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) + * On entry, the symmetric matrix A. If UPLO = 'U', the leading + * n-by-n upper triangular part of A contains the upper + * triangular part of the matrix A, and the strictly lower + * triangular part of A is not referenced. If UPLO = 'L', the + * leading n-by-n lower triangular part of A contains the lower + * triangular part of the matrix A, and the strictly upper + * triangular part of A is not referenced. + * On exit: + * if UPLO = 'U', the last NB columns have been reduced to + * tridiagonal form, with the diagonal elements overwriting + * the diagonal elements of A; the elements above the diagonal + * with the array TAU, represent the orthogonal matrix Q as a + * product of elementary reflectors; + * if UPLO = 'L', the first NB columns have been reduced to + * tridiagonal form, with the diagonal elements overwriting + * the diagonal elements of A; the elements below the diagonal + * with the array TAU, represent the orthogonal matrix Q as a + * product of elementary reflectors. + * See Further Details. + * + * LDA (input) INTEGER + * The leading dimension of the array A. LDA >= (1,N). + * + * E (output) DOUBLE PRECISION array, dimension (N-1) + * If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal + * elements of the last NB columns of the reduced matrix; + * if UPLO = 'L', E(1:nb) contains the subdiagonal elements of + * the first NB columns of the reduced matrix. + * + * TAU (output) DOUBLE PRECISION array, dimension (N-1) + * The scalar factors of the elementary reflectors, stored in + * TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'. + * See Further Details. + * + * W (output) DOUBLE PRECISION array, dimension (LDW,NB) + * The n-by-nb matrix W required to update the unreduced part + * of A. + * + * LDW (input) INTEGER + * The leading dimension of the array W. LDW >= max(1,N). + * + * Further Details + * =============== + * + * If UPLO = 'U', the matrix Q is represented as a product of elementary + * reflectors + * + * Q = H(n) H(n-1) . . . H(n-nb+1). + * + * Each H(i) has the form + * + * H(i) = I - tau * v * v' + * + * where tau is a real scalar, and v is a real vector with + * v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i), + * and tau in TAU(i-1). + * + * If UPLO = 'L', the matrix Q is represented as a product of elementary + * reflectors + * + * Q = H(1) H(2) . . . H(nb). + * + * Each H(i) has the form + * + * H(i) = I - tau * v * v' + * + * where tau is a real scalar, and v is a real vector with + * v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i), + * and tau in TAU(i). + * + * The elements of the vectors v together form the n-by-nb matrix V + * which is needed, with W, to apply the transformation to the unreduced + * part of the matrix, using a symmetric rank-2k update of the form: + * A := A - V*W' - W*V'. + * + * The contents of A on exit are illustrated by the following examples + * with n = 5 and nb = 2: + * + * if UPLO = 'U': if UPLO = 'L': + * + * ( a a a v4 v5 ) ( d ) + * ( a a v4 v5 ) ( 1 d ) + * ( a 1 v5 ) ( v1 1 a ) + * ( d 1 ) ( v1 v2 a a ) + * ( d ) ( v1 v2 a a a ) + * + * where d denotes a diagonal element of the reduced matrix, a denotes + * an element of the original matrix that is unchanged, and vi denotes + * an element of the vector defining H(i). + * + * ===================================================================== + * + * .. Parameters .. + DOUBLE PRECISION ZERO, ONE, HALF + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, HALF = 0.5D+0 ) + * .. + * .. Local Scalars .. + INTEGER I, IW + DOUBLE PRECISION ALPHA + * .. + * .. External Subroutines .. + EXTERNAL DAXPY, DGEMV, DLARFG, DSCAL, DSYMV + * .. + * .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DDOT + EXTERNAL LSAME, DDOT + * .. + * .. Intrinsic Functions .. + INTRINSIC MIN + * .. + * .. Executable Statements .. + * + * Quick return if possible + * + IF( N.LE.0 ) + $ RETURN + * + IF( LSAME( UPLO, 'U' ) ) THEN + * + * Reduce last NB columns of upper triangle + * + DO 10 I = N, N - NB + 1, -1 + IW = I - N + NB + IF( I.LT.N ) THEN + * + * Update A(1:i,i) + * + CALL DGEMV( 'No transpose', I, N-I, -ONE, A( 1, I+1 ), + $ LDA, W( I, IW+1 ), LDW, ONE, A( 1, I ), 1 ) + CALL DGEMV( 'No transpose', I, N-I, -ONE, W( 1, IW+1 ), + $ LDW, A( I, I+1 ), LDA, ONE, A( 1, I ), 1 ) + END IF + IF( I.GT.1 ) THEN + * + * Generate elementary reflector H(i) to annihilate + * A(1:i-2,i) + * + CALL DLARFG( I-1, A( I-1, I ), A( 1, I ), 1, TAU( I-1 ) ) + E( I-1 ) = A( I-1, I ) + A( I-1, I ) = ONE + * + * Compute W(1:i-1,i) + * + CALL DSYMV( 'Upper', I-1, ONE, A, LDA, A( 1, I ), 1, + $ ZERO, W( 1, IW ), 1 ) + IF( I.LT.N ) THEN + CALL DGEMV( 'Transpose', I-1, N-I, ONE, W( 1, IW+1 ), + $ LDW, A( 1, I ), 1, ZERO, W( I+1, IW ), 1 ) + CALL DGEMV( 'No transpose', I-1, N-I, -ONE, + $ A( 1, I+1 ), LDA, W( I+1, IW ), 1, ONE, + $ W( 1, IW ), 1 ) + CALL DGEMV( 'Transpose', I-1, N-I, ONE, A( 1, I+1 ), + $ LDA, A( 1, I ), 1, ZERO, W( I+1, IW ), 1 ) + CALL DGEMV( 'No transpose', I-1, N-I, -ONE, + $ W( 1, IW+1 ), LDW, W( I+1, IW ), 1, ONE, + $ W( 1, IW ), 1 ) + END IF + CALL DSCAL( I-1, TAU( I-1 ), W( 1, IW ), 1 ) + ALPHA = -HALF*TAU( I-1 )*DDOT( I-1, W( 1, IW ), 1, + $ A( 1, I ), 1 ) + CALL DAXPY( I-1, ALPHA, A( 1, I ), 1, W( 1, IW ), 1 ) + END IF + * + 10 CONTINUE + ELSE + * + * Reduce first NB columns of lower triangle + * + DO 20 I = 1, NB + * + * Update A(i:n,i) + * + CALL DGEMV( 'No transpose', N-I+1, I-1, -ONE, A( I, 1 ), + $ LDA, W( I, 1 ), LDW, ONE, A( I, I ), 1 ) + CALL DGEMV( 'No transpose', N-I+1, I-1, -ONE, W( I, 1 ), + $ LDW, A( I, 1 ), LDA, ONE, A( I, I ), 1 ) + IF( I.LT.N ) THEN + * + * Generate elementary reflector H(i) to annihilate + * A(i+2:n,i) + * + CALL DLARFG( N-I, A( I+1, I ), A( MIN( I+2, N ), I ), 1, + $ TAU( I ) ) + E( I ) = A( I+1, I ) + A( I+1, I ) = ONE + * + * Compute W(i+1:n,i) + * + CALL DSYMV( 'Lower', N-I, ONE, A( I+1, I+1 ), LDA, + $ A( I+1, I ), 1, ZERO, W( I+1, I ), 1 ) + CALL DGEMV( 'Transpose', N-I, I-1, ONE, W( I+1, 1 ), LDW, + $ A( I+1, I ), 1, ZERO, W( 1, I ), 1 ) + CALL DGEMV( 'No transpose', N-I, I-1, -ONE, A( I+1, 1 ), + $ LDA, W( 1, I ), 1, ONE, W( I+1, I ), 1 ) + CALL DGEMV( 'Transpose', N-I, I-1, ONE, A( I+1, 1 ), LDA, + $ A( I+1, I ), 1, ZERO, W( 1, I ), 1 ) + CALL DGEMV( 'No transpose', N-I, I-1, -ONE, W( I+1, 1 ), + $ LDW, W( 1, I ), 1, ONE, W( I+1, I ), 1 ) + CALL DSCAL( N-I, TAU( I ), W( I+1, I ), 1 ) + ALPHA = -HALF*TAU( I )*DDOT( N-I, W( I+1, I ), 1, + $ A( I+1, I ), 1 ) + CALL DAXPY( N-I, ALPHA, A( I+1, I ), 1, W( I+1, I ), 1 ) + END IF + * + 20 CONTINUE + END IF + * + RETURN + * + * End of DLATRD + * + END diff -cNr octave-2.0.7/libcruft/lapack/dorg2l.f octave-2.0.8/libcruft/lapack/dorg2l.f *** octave-2.0.7/libcruft/lapack/dorg2l.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/dorg2l.f Mon Jun 23 14:00:08 1997 *************** *** 0 **** --- 1,128 ---- + SUBROUTINE DORG2L( M, N, K, A, LDA, TAU, WORK, INFO ) + * + * -- LAPACK routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * February 29, 1992 + * + * .. Scalar Arguments .. + INTEGER INFO, K, LDA, M, N + * .. + * .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) + * .. + * + * Purpose + * ======= + * + * DORG2L generates an m by n real matrix Q with orthonormal columns, + * which is defined as the last n columns of a product of k elementary + * reflectors of order m + * + * Q = H(k) . . . H(2) H(1) + * + * as returned by DGEQLF. + * + * Arguments + * ========= + * + * M (input) INTEGER + * The number of rows of the matrix Q. M >= 0. + * + * N (input) INTEGER + * The number of columns of the matrix Q. M >= N >= 0. + * + * K (input) INTEGER + * The number of elementary reflectors whose product defines the + * matrix Q. N >= K >= 0. + * + * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) + * On entry, the (n-k+i)-th column must contain the vector which + * defines the elementary reflector H(i), for i = 1,2,...,k, as + * returned by DGEQLF in the last k columns of its array + * argument A. + * On exit, the m by n matrix Q. + * + * LDA (input) INTEGER + * The first dimension of the array A. LDA >= max(1,M). + * + * TAU (input) DOUBLE PRECISION array, dimension (K) + * TAU(i) must contain the scalar factor of the elementary + * reflector H(i), as returned by DGEQLF. + * + * WORK (workspace) DOUBLE PRECISION array, dimension (N) + * + * INFO (output) INTEGER + * = 0: successful exit + * < 0: if INFO = -i, the i-th argument has an illegal value + * + * ===================================================================== + * + * .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) + * .. + * .. Local Scalars .. + INTEGER I, II, J, L + * .. + * .. External Subroutines .. + EXTERNAL DLARF, DSCAL, XERBLA + * .. + * .. Intrinsic Functions .. + INTRINSIC MAX + * .. + * .. Executable Statements .. + * + * Test the input arguments + * + INFO = 0 + IF( M.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 .OR. N.GT.M ) THEN + INFO = -2 + ELSE IF( K.LT.0 .OR. K.GT.N ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -5 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DORG2L', -INFO ) + RETURN + END IF + * + * Quick return if possible + * + IF( N.LE.0 ) + $ RETURN + * + * Initialise columns 1:n-k to columns of the unit matrix + * + DO 20 J = 1, N - K + DO 10 L = 1, M + A( L, J ) = ZERO + 10 CONTINUE + A( M-N+J, J ) = ONE + 20 CONTINUE + * + DO 40 I = 1, K + II = N - K + I + * + * Apply H(i) to A(1:m-k+i,1:n-k+i) from the left + * + A( M-N+II, II ) = ONE + CALL DLARF( 'Left', M-N+II, II-1, A( 1, II ), 1, TAU( I ), A, + $ LDA, WORK ) + CALL DSCAL( M-N+II-1, -TAU( I ), A( 1, II ), 1 ) + A( M-N+II, II ) = ONE - TAU( I ) + * + * Set A(m-k+i+1:m,n-k+i) to zero + * + DO 30 L = M - N + II + 1, M + A( L, II ) = ZERO + 30 CONTINUE + 40 CONTINUE + RETURN + * + * End of DORG2L + * + END diff -cNr octave-2.0.7/libcruft/lapack/dorgql.f octave-2.0.8/libcruft/lapack/dorgql.f *** octave-2.0.7/libcruft/lapack/dorgql.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/dorgql.f Mon Jun 23 14:00:08 1997 *************** *** 0 **** --- 1,205 ---- + SUBROUTINE DORGQL( M, N, K, A, LDA, TAU, WORK, LWORK, INFO ) + * + * -- LAPACK routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * September 30, 1994 + * + * .. Scalar Arguments .. + INTEGER INFO, K, LDA, LWORK, M, N + * .. + * .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( LWORK ) + * .. + * + * Purpose + * ======= + * + * DORGQL generates an M-by-N real matrix Q with orthonormal columns, + * which is defined as the last N columns of a product of K elementary + * reflectors of order M + * + * Q = H(k) . . . H(2) H(1) + * + * as returned by DGEQLF. + * + * Arguments + * ========= + * + * M (input) INTEGER + * The number of rows of the matrix Q. M >= 0. + * + * N (input) INTEGER + * The number of columns of the matrix Q. M >= N >= 0. + * + * K (input) INTEGER + * The number of elementary reflectors whose product defines the + * matrix Q. N >= K >= 0. + * + * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) + * On entry, the (n-k+i)-th column must contain the vector which + * defines the elementary reflector H(i), for i = 1,2,...,k, as + * returned by DGEQLF in the last k columns of its array + * argument A. + * On exit, the M-by-N matrix Q. + * + * LDA (input) INTEGER + * The first dimension of the array A. LDA >= max(1,M). + * + * TAU (input) DOUBLE PRECISION array, dimension (K) + * TAU(i) must contain the scalar factor of the elementary + * reflector H(i), as returned by DGEQLF. + * + * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) + * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. + * + * LWORK (input) INTEGER + * The dimension of the array WORK. LWORK >= max(1,N). + * For optimum performance LWORK >= N*NB, where NB is the + * optimal blocksize. + * + * INFO (output) INTEGER + * = 0: successful exit + * < 0: if INFO = -i, the i-th argument has an illegal value + * + * ===================================================================== + * + * .. Parameters .. + DOUBLE PRECISION ZERO + PARAMETER ( ZERO = 0.0D+0 ) + * .. + * .. Local Scalars .. + INTEGER I, IB, IINFO, IWS, J, KK, L, LDWORK, NB, NBMIN, + $ NX + * .. + * .. External Subroutines .. + EXTERNAL DLARFB, DLARFT, DORG2L, XERBLA + * .. + * .. Intrinsic Functions .. + INTRINSIC MAX, MIN + * .. + * .. External Functions .. + INTEGER ILAENV + EXTERNAL ILAENV + * .. + * .. Executable Statements .. + * + * Test the input arguments + * + INFO = 0 + IF( M.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 .OR. N.GT.M ) THEN + INFO = -2 + ELSE IF( K.LT.0 .OR. K.GT.N ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -5 + ELSE IF( LWORK.LT.MAX( 1, N ) ) THEN + INFO = -8 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DORGQL', -INFO ) + RETURN + END IF + * + * Quick return if possible + * + IF( N.LE.0 ) THEN + WORK( 1 ) = 1 + RETURN + END IF + * + * Determine the block size. + * + NB = ILAENV( 1, 'DORGQL', ' ', M, N, K, -1 ) + NBMIN = 2 + NX = 0 + IWS = N + IF( NB.GT.1 .AND. NB.LT.K ) THEN + * + * Determine when to cross over from blocked to unblocked code. + * + NX = MAX( 0, ILAENV( 3, 'DORGQL', ' ', M, N, K, -1 ) ) + IF( NX.LT.K ) THEN + * + * Determine if workspace is large enough for blocked code. + * + LDWORK = N + IWS = LDWORK*NB + IF( LWORK.LT.IWS ) THEN + * + * Not enough workspace to use optimal NB: reduce NB and + * determine the minimum value of NB. + * + NB = LWORK / LDWORK + NBMIN = MAX( 2, ILAENV( 2, 'DORGQL', ' ', M, N, K, -1 ) ) + END IF + END IF + END IF + * + IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN + * + * Use blocked code after the first block. + * The last kk columns are handled by the block method. + * + KK = MIN( K, ( ( K-NX+NB-1 ) / NB )*NB ) + * + * Set A(m-kk+1:m,1:n-kk) to zero. + * + DO 20 J = 1, N - KK + DO 10 I = M - KK + 1, M + A( I, J ) = ZERO + 10 CONTINUE + 20 CONTINUE + ELSE + KK = 0 + END IF + * + * Use unblocked code for the first or only block. + * + CALL DORG2L( M-KK, N-KK, K-KK, A, LDA, TAU, WORK, IINFO ) + * + IF( KK.GT.0 ) THEN + * + * Use blocked code + * + DO 50 I = K - KK + 1, K, NB + IB = MIN( NB, K-I+1 ) + IF( N-K+I.GT.1 ) THEN + * + * Form the triangular factor of the block reflector + * H = H(i+ib-1) . . . H(i+1) H(i) + * + CALL DLARFT( 'Backward', 'Columnwise', M-K+I+IB-1, IB, + $ A( 1, N-K+I ), LDA, TAU( I ), WORK, LDWORK ) + * + * Apply H to A(1:m-k+i+ib-1,1:n-k+i-1) from the left + * + CALL DLARFB( 'Left', 'No transpose', 'Backward', + $ 'Columnwise', M-K+I+IB-1, N-K+I-1, IB, + $ A( 1, N-K+I ), LDA, WORK, LDWORK, A, LDA, + $ WORK( IB+1 ), LDWORK ) + END IF + * + * Apply H to rows 1:m-k+i+ib-1 of current block + * + CALL DORG2L( M-K+I+IB-1, IB, IB, A( 1, N-K+I ), LDA, + $ TAU( I ), WORK, IINFO ) + * + * Set rows m-k+i+ib:m of current block to zero + * + DO 40 J = N - K + I, N - K + I + IB - 1 + DO 30 L = M - K + I + IB, M + A( L, J ) = ZERO + 30 CONTINUE + 40 CONTINUE + 50 CONTINUE + END IF + * + WORK( 1 ) = IWS + RETURN + * + * End of DORGQL + * + END diff -cNr octave-2.0.7/libcruft/lapack/dorgtr.f octave-2.0.8/libcruft/lapack/dorgtr.f *** octave-2.0.7/libcruft/lapack/dorgtr.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/dorgtr.f Mon Jun 23 14:00:08 1997 *************** *** 0 **** --- 1,163 ---- + SUBROUTINE DORGTR( UPLO, N, A, LDA, TAU, WORK, LWORK, INFO ) + * + * -- LAPACK routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * September 30, 1994 + * + * .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, LWORK, N + * .. + * .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( LWORK ) + * .. + * + * Purpose + * ======= + * + * DORGTR generates a real orthogonal matrix Q which is defined as the + * product of n-1 elementary reflectors of order N, as returned by + * DSYTRD: + * + * if UPLO = 'U', Q = H(n-1) . . . H(2) H(1), + * + * if UPLO = 'L', Q = H(1) H(2) . . . H(n-1). + * + * Arguments + * ========= + * + * UPLO (input) CHARACTER*1 + * = 'U': Upper triangle of A contains elementary reflectors + * from DSYTRD; + * = 'L': Lower triangle of A contains elementary reflectors + * from DSYTRD. + * + * N (input) INTEGER + * The order of the matrix Q. N >= 0. + * + * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) + * On entry, the vectors which define the elementary reflectors, + * as returned by DSYTRD. + * On exit, the N-by-N orthogonal matrix Q. + * + * LDA (input) INTEGER + * The leading dimension of the array A. LDA >= max(1,N). + * + * TAU (input) DOUBLE PRECISION array, dimension (N-1) + * TAU(i) must contain the scalar factor of the elementary + * reflector H(i), as returned by DSYTRD. + * + * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) + * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. + * + * LWORK (input) INTEGER + * The dimension of the array WORK. LWORK >= max(1,N-1). + * For optimum performance LWORK >= (N-1)*NB, where NB is + * the optimal blocksize. + * + * INFO (output) INTEGER + * = 0: successful exit + * < 0: if INFO = -i, the i-th argument had an illegal value + * + * ===================================================================== + * + * .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) + * .. + * .. Local Scalars .. + LOGICAL UPPER + INTEGER I, IINFO, J + * .. + * .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME + * .. + * .. External Subroutines .. + EXTERNAL DORGQL, DORGQR, XERBLA + * .. + * .. Intrinsic Functions .. + INTRINSIC MAX + * .. + * .. Executable Statements .. + * + * Test the input arguments + * + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + ELSE IF( LWORK.LT.MAX( 1, N-1 ) ) THEN + INFO = -7 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DORGTR', -INFO ) + RETURN + END IF + * + * Quick return if possible + * + IF( N.EQ.0 ) THEN + WORK( 1 ) = 1 + RETURN + END IF + * + IF( UPPER ) THEN + * + * Q was determined by a call to DSYTRD with UPLO = 'U' + * + * Shift the vectors which define the elementary reflectors one + * column to the left, and set the last row and column of Q to + * those of the unit matrix + * + DO 20 J = 1, N - 1 + DO 10 I = 1, J - 1 + A( I, J ) = A( I, J+1 ) + 10 CONTINUE + A( N, J ) = ZERO + 20 CONTINUE + DO 30 I = 1, N - 1 + A( I, N ) = ZERO + 30 CONTINUE + A( N, N ) = ONE + * + * Generate Q(1:n-1,1:n-1) + * + CALL DORGQL( N-1, N-1, N-1, A, LDA, TAU, WORK, LWORK, IINFO ) + * + ELSE + * + * Q was determined by a call to DSYTRD with UPLO = 'L'. + * + * Shift the vectors which define the elementary reflectors one + * column to the right, and set the first row and column of Q to + * those of the unit matrix + * + DO 50 J = N, 2, -1 + A( 1, J ) = ZERO + DO 40 I = J + 1, N + A( I, J ) = A( I, J-1 ) + 40 CONTINUE + 50 CONTINUE + A( 1, 1 ) = ONE + DO 60 I = 2, N + A( I, 1 ) = ZERO + 60 CONTINUE + IF( N.GT.1 ) THEN + * + * Generate Q(2:n,2:n) + * + CALL DORGQR( N-1, N-1, N-1, A( 2, 2 ), LDA, TAU, WORK, + $ LWORK, IINFO ) + END IF + END IF + RETURN + * + * End of DORGTR + * + END diff -cNr octave-2.0.7/libcruft/lapack/dsteqr.f octave-2.0.8/libcruft/lapack/dsteqr.f *** octave-2.0.7/libcruft/lapack/dsteqr.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/dsteqr.f Mon Jun 23 14:00:08 1997 *************** *** 0 **** --- 1,501 ---- + SUBROUTINE DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO ) + * + * -- LAPACK routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * September 30, 1994 + * + * .. Scalar Arguments .. + CHARACTER COMPZ + INTEGER INFO, LDZ, N + * .. + * .. Array Arguments .. + DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * ) + * .. + * + * Purpose + * ======= + * + * DSTEQR computes all eigenvalues and, optionally, eigenvectors of a + * symmetric tridiagonal matrix using the implicit QL or QR method. + * The eigenvectors of a full or band symmetric matrix can also be found + * if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to + * tridiagonal form. + * + * Arguments + * ========= + * + * COMPZ (input) CHARACTER*1 + * = 'N': Compute eigenvalues only. + * = 'V': Compute eigenvalues and eigenvectors of the original + * symmetric matrix. On entry, Z must contain the + * orthogonal matrix used to reduce the original matrix + * to tridiagonal form. + * = 'I': Compute eigenvalues and eigenvectors of the + * tridiagonal matrix. Z is initialized to the identity + * matrix. + * + * N (input) INTEGER + * The order of the matrix. N >= 0. + * + * D (input/output) DOUBLE PRECISION array, dimension (N) + * On entry, the diagonal elements of the tridiagonal matrix. + * On exit, if INFO = 0, the eigenvalues in ascending order. + * + * E (input/output) DOUBLE PRECISION array, dimension (N-1) + * On entry, the (n-1) subdiagonal elements of the tridiagonal + * matrix. + * On exit, E has been destroyed. + * + * Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N) + * On entry, if COMPZ = 'V', then Z contains the orthogonal + * matrix used in the reduction to tridiagonal form. + * On exit, if INFO = 0, then if COMPZ = 'V', Z contains the + * orthonormal eigenvectors of the original symmetric matrix, + * and if COMPZ = 'I', Z contains the orthonormal eigenvectors + * of the symmetric tridiagonal matrix. + * If COMPZ = 'N', then Z is not referenced. + * + * LDZ (input) INTEGER + * The leading dimension of the array Z. LDZ >= 1, and if + * eigenvectors are desired, then LDZ >= max(1,N). + * + * WORK (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2)) + * If COMPZ = 'N', then WORK is not referenced. + * + * INFO (output) INTEGER + * = 0: successful exit + * < 0: if INFO = -i, the i-th argument had an illegal value + * > 0: the algorithm has failed to find all the eigenvalues in + * a total of 30*N iterations; if INFO = i, then i + * elements of E have not converged to zero; on exit, D + * and E contain the elements of a symmetric tridiagonal + * matrix which is orthogonally similar to the original + * matrix. + * + * ===================================================================== + * + * .. Parameters .. + DOUBLE PRECISION ZERO, ONE, TWO, THREE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, + $ THREE = 3.0D0 ) + INTEGER MAXIT + PARAMETER ( MAXIT = 30 ) + * .. + * .. Local Scalars .. + INTEGER I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND, + $ LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1, + $ NM1, NMAXIT + DOUBLE PRECISION ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2, + $ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST + * .. + * .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DLAMCH, DLANST, DLAPY2 + EXTERNAL LSAME, DLAMCH, DLANST, DLAPY2 + * .. + * .. External Subroutines .. + EXTERNAL DLAE2, DLAEV2, DLARTG, DLASCL, DLASET, DLASR, + $ DLASRT, DSWAP, XERBLA + * .. + * .. Intrinsic Functions .. + INTRINSIC ABS, MAX, SIGN, SQRT + * .. + * .. Executable Statements .. + * + * Test the input parameters. + * + INFO = 0 + * + IF( LSAME( COMPZ, 'N' ) ) THEN + ICOMPZ = 0 + ELSE IF( LSAME( COMPZ, 'V' ) ) THEN + ICOMPZ = 1 + ELSE IF( LSAME( COMPZ, 'I' ) ) THEN + ICOMPZ = 2 + ELSE + ICOMPZ = -1 + END IF + IF( ICOMPZ.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, + $ N ) ) ) THEN + INFO = -6 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DSTEQR', -INFO ) + RETURN + END IF + * + * Quick return if possible + * + IF( N.EQ.0 ) + $ RETURN + * + IF( N.EQ.1 ) THEN + IF( ICOMPZ.EQ.2 ) + $ Z( 1, 1 ) = ONE + RETURN + END IF + * + * Determine the unit roundoff and over/underflow thresholds. + * + EPS = DLAMCH( 'E' ) + EPS2 = EPS**2 + SAFMIN = DLAMCH( 'S' ) + SAFMAX = ONE / SAFMIN + SSFMAX = SQRT( SAFMAX ) / THREE + SSFMIN = SQRT( SAFMIN ) / EPS2 + * + * Compute the eigenvalues and eigenvectors of the tridiagonal + * matrix. + * + IF( ICOMPZ.EQ.2 ) + $ CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ ) + * + NMAXIT = N*MAXIT + JTOT = 0 + * + * Determine where the matrix splits and choose QL or QR iteration + * for each block, according to whether top or bottom diagonal + * element is smaller. + * + L1 = 1 + NM1 = N - 1 + * + 10 CONTINUE + IF( L1.GT.N ) + $ GO TO 160 + IF( L1.GT.1 ) + $ E( L1-1 ) = ZERO + IF( L1.LE.NM1 ) THEN + DO 20 M = L1, NM1 + TST = ABS( E( M ) ) + IF( TST.EQ.ZERO ) + $ GO TO 30 + IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+ + $ 1 ) ) ) )*EPS ) THEN + E( M ) = ZERO + GO TO 30 + END IF + 20 CONTINUE + END IF + M = N + * + 30 CONTINUE + L = L1 + LSV = L + LEND = M + LENDSV = LEND + L1 = M + 1 + IF( LEND.EQ.L ) + $ GO TO 10 + * + * Scale submatrix in rows and columns L to LEND + * + ANORM = DLANST( 'I', LEND-L+1, D( L ), E( L ) ) + ISCALE = 0 + IF( ANORM.EQ.ZERO ) + $ GO TO 10 + IF( ANORM.GT.SSFMAX ) THEN + ISCALE = 1 + CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N, + $ INFO ) + CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N, + $ INFO ) + ELSE IF( ANORM.LT.SSFMIN ) THEN + ISCALE = 2 + CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N, + $ INFO ) + CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N, + $ INFO ) + END IF + * + * Choose between QL and QR iteration + * + IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN + LEND = LSV + L = LENDSV + END IF + * + IF( LEND.GT.L ) THEN + * + * QL Iteration + * + * Look for small subdiagonal element. + * + 40 CONTINUE + IF( L.NE.LEND ) THEN + LENDM1 = LEND - 1 + DO 50 M = L, LENDM1 + TST = ABS( E( M ) )**2 + IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+ + $ SAFMIN )GO TO 60 + 50 CONTINUE + END IF + * + M = LEND + * + 60 CONTINUE + IF( M.LT.LEND ) + $ E( M ) = ZERO + P = D( L ) + IF( M.EQ.L ) + $ GO TO 80 + * + * If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 + * to compute its eigensystem. + * + IF( M.EQ.L+1 ) THEN + IF( ICOMPZ.GT.0 ) THEN + CALL DLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S ) + WORK( L ) = C + WORK( N-1+L ) = S + CALL DLASR( 'R', 'V', 'B', N, 2, WORK( L ), + $ WORK( N-1+L ), Z( 1, L ), LDZ ) + ELSE + CALL DLAE2( D( L ), E( L ), D( L+1 ), RT1, RT2 ) + END IF + D( L ) = RT1 + D( L+1 ) = RT2 + E( L ) = ZERO + L = L + 2 + IF( L.LE.LEND ) + $ GO TO 40 + GO TO 140 + END IF + * + IF( JTOT.EQ.NMAXIT ) + $ GO TO 140 + JTOT = JTOT + 1 + * + * Form shift. + * + G = ( D( L+1 )-P ) / ( TWO*E( L ) ) + R = DLAPY2( G, ONE ) + G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) ) + * + S = ONE + C = ONE + P = ZERO + * + * Inner loop + * + MM1 = M - 1 + DO 70 I = MM1, L, -1 + F = S*E( I ) + B = C*E( I ) + CALL DLARTG( G, F, C, S, R ) + IF( I.NE.M-1 ) + $ E( I+1 ) = R + G = D( I+1 ) - P + R = ( D( I )-G )*S + TWO*C*B + P = S*R + D( I+1 ) = G + P + G = C*R - B + * + * If eigenvectors are desired, then save rotations. + * + IF( ICOMPZ.GT.0 ) THEN + WORK( I ) = C + WORK( N-1+I ) = -S + END IF + * + 70 CONTINUE + * + * If eigenvectors are desired, then apply saved rotations. + * + IF( ICOMPZ.GT.0 ) THEN + MM = M - L + 1 + CALL DLASR( 'R', 'V', 'B', N, MM, WORK( L ), WORK( N-1+L ), + $ Z( 1, L ), LDZ ) + END IF + * + D( L ) = D( L ) - P + E( L ) = G + GO TO 40 + * + * Eigenvalue found. + * + 80 CONTINUE + D( L ) = P + * + L = L + 1 + IF( L.LE.LEND ) + $ GO TO 40 + GO TO 140 + * + ELSE + * + * QR Iteration + * + * Look for small superdiagonal element. + * + 90 CONTINUE + IF( L.NE.LEND ) THEN + LENDP1 = LEND + 1 + DO 100 M = L, LENDP1, -1 + TST = ABS( E( M-1 ) )**2 + IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+ + $ SAFMIN )GO TO 110 + 100 CONTINUE + END IF + * + M = LEND + * + 110 CONTINUE + IF( M.GT.LEND ) + $ E( M-1 ) = ZERO + P = D( L ) + IF( M.EQ.L ) + $ GO TO 130 + * + * If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 + * to compute its eigensystem. + * + IF( M.EQ.L-1 ) THEN + IF( ICOMPZ.GT.0 ) THEN + CALL DLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S ) + WORK( M ) = C + WORK( N-1+M ) = S + CALL DLASR( 'R', 'V', 'F', N, 2, WORK( M ), + $ WORK( N-1+M ), Z( 1, L-1 ), LDZ ) + ELSE + CALL DLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 ) + END IF + D( L-1 ) = RT1 + D( L ) = RT2 + E( L-1 ) = ZERO + L = L - 2 + IF( L.GE.LEND ) + $ GO TO 90 + GO TO 140 + END IF + * + IF( JTOT.EQ.NMAXIT ) + $ GO TO 140 + JTOT = JTOT + 1 + * + * Form shift. + * + G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) ) + R = DLAPY2( G, ONE ) + G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) ) + * + S = ONE + C = ONE + P = ZERO + * + * Inner loop + * + LM1 = L - 1 + DO 120 I = M, LM1 + F = S*E( I ) + B = C*E( I ) + CALL DLARTG( G, F, C, S, R ) + IF( I.NE.M ) + $ E( I-1 ) = R + G = D( I ) - P + R = ( D( I+1 )-G )*S + TWO*C*B + P = S*R + D( I ) = G + P + G = C*R - B + * + * If eigenvectors are desired, then save rotations. + * + IF( ICOMPZ.GT.0 ) THEN + WORK( I ) = C + WORK( N-1+I ) = S + END IF + * + 120 CONTINUE + * + * If eigenvectors are desired, then apply saved rotations. + * + IF( ICOMPZ.GT.0 ) THEN + MM = L - M + 1 + CALL DLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ), + $ Z( 1, M ), LDZ ) + END IF + * + D( L ) = D( L ) - P + E( LM1 ) = G + GO TO 90 + * + * Eigenvalue found. + * + 130 CONTINUE + D( L ) = P + * + L = L - 1 + IF( L.GE.LEND ) + $ GO TO 90 + GO TO 140 + * + END IF + * + * Undo scaling if necessary + * + 140 CONTINUE + IF( ISCALE.EQ.1 ) THEN + CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1, + $ D( LSV ), N, INFO ) + CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ), + $ N, INFO ) + ELSE IF( ISCALE.EQ.2 ) THEN + CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1, + $ D( LSV ), N, INFO ) + CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ), + $ N, INFO ) + END IF + * + * Check for no convergence to an eigenvalue after a total + * of N*MAXIT iterations. + * + IF( JTOT.LT.NMAXIT ) + $ GO TO 10 + DO 150 I = 1, N - 1 + IF( E( I ).NE.ZERO ) + $ INFO = INFO + 1 + 150 CONTINUE + GO TO 190 + * + * Order eigenvalues and eigenvectors. + * + 160 CONTINUE + IF( ICOMPZ.EQ.0 ) THEN + * + * Use Quick Sort + * + CALL DLASRT( 'I', N, D, INFO ) + * + ELSE + * + * Use Selection Sort to minimize swaps of eigenvectors + * + DO 180 II = 2, N + I = II - 1 + K = I + P = D( I ) + DO 170 J = II, N + IF( D( J ).LT.P ) THEN + K = J + P = D( J ) + END IF + 170 CONTINUE + IF( K.NE.I ) THEN + D( K ) = D( I ) + D( I ) = P + CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 ) + END IF + 180 CONTINUE + END IF + * + 190 CONTINUE + RETURN + * + * End of DSTEQR + * + END diff -cNr octave-2.0.7/libcruft/lapack/dsterf.f octave-2.0.8/libcruft/lapack/dsterf.f *** octave-2.0.7/libcruft/lapack/dsterf.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/dsterf.f Mon Jun 23 14:00:08 1997 *************** *** 0 **** --- 1,381 ---- + SUBROUTINE DSTERF( N, D, E, INFO ) + * + * -- LAPACK routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * September 30, 1994 + * + * .. Scalar Arguments .. + INTEGER INFO, N + * .. + * .. Array Arguments .. + DOUBLE PRECISION D( * ), E( * ) + * .. + * + * Purpose + * ======= + * + * DSTERF computes all eigenvalues of a symmetric tridiagonal matrix + * using the Pal-Walker-Kahan variant of the QL or QR algorithm. + * + * Arguments + * ========= + * + * N (input) INTEGER + * The order of the matrix. N >= 0. + * + * D (input/output) DOUBLE PRECISION array, dimension (N) + * On entry, the n diagonal elements of the tridiagonal matrix. + * On exit, if INFO = 0, the eigenvalues in ascending order. + * + * E (input/output) DOUBLE PRECISION array, dimension (N-1) + * On entry, the (n-1) subdiagonal elements of the tridiagonal + * matrix. + * On exit, E has been destroyed. + * + * INFO (output) INTEGER + * = 0: successful exit + * < 0: if INFO = -i, the i-th argument had an illegal value + * > 0: the algorithm failed to find all of the eigenvalues in + * a total of 30*N iterations; if INFO = i, then i + * elements of E have not converged to zero. + * + * ===================================================================== + * + * .. Parameters .. + DOUBLE PRECISION ZERO, ONE, TWO, THREE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, + $ THREE = 3.0D0 ) + INTEGER MAXIT + PARAMETER ( MAXIT = 30 ) + * .. + * .. Local Scalars .. + INTEGER I, ISCALE, JTOT, L, L1, LEND, LENDM1, LENDP1, + $ LENDSV, LM1, LSV, M, MM1, NM1, NMAXIT + DOUBLE PRECISION ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC, + $ OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN, + $ SIGMA, SSFMAX, SSFMIN, TST + * .. + * .. External Functions .. + DOUBLE PRECISION DLAMCH, DLANST, DLAPY2 + EXTERNAL DLAMCH, DLANST, DLAPY2 + * .. + * .. External Subroutines .. + EXTERNAL DLAE2, DLASCL, DLASRT, XERBLA + * .. + * .. Intrinsic Functions .. + INTRINSIC ABS, SIGN, SQRT + * .. + * .. Executable Statements .. + * + * Test the input parameters. + * + INFO = 0 + * + * Quick return if possible + * + IF( N.LT.0 ) THEN + INFO = -1 + CALL XERBLA( 'DSTERF', -INFO ) + RETURN + END IF + IF( N.LE.1 ) + $ RETURN + * + * Determine the unit roundoff for this environment. + * + EPS = DLAMCH( 'E' ) + EPS2 = EPS**2 + SAFMIN = DLAMCH( 'S' ) + SAFMAX = ONE / SAFMIN + SSFMAX = SQRT( SAFMAX ) / THREE + SSFMIN = SQRT( SAFMIN ) / EPS2 + * + * Compute the eigenvalues of the tridiagonal matrix. + * + NMAXIT = N*MAXIT + SIGMA = ZERO + JTOT = 0 + * + * Determine where the matrix splits and choose QL or QR iteration + * for each block, according to whether top or bottom diagonal + * element is smaller. + * + L1 = 1 + NM1 = N - 1 + * + 10 CONTINUE + IF( L1.GT.N ) + $ GO TO 170 + IF( L1.GT.1 ) + $ E( L1-1 ) = ZERO + IF( L1.LE.NM1 ) THEN + DO 20 M = L1, NM1 + TST = ABS( E( M ) ) + IF( TST.EQ.ZERO ) + $ GO TO 30 + IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+ + $ 1 ) ) ) )*EPS ) THEN + E( M ) = ZERO + GO TO 30 + END IF + 20 CONTINUE + END IF + M = N + * + 30 CONTINUE + L = L1 + LSV = L + LEND = M + LENDSV = LEND + L1 = M + 1 + IF( LEND.EQ.L ) + $ GO TO 10 + * + * Scale submatrix in rows and columns L to LEND + * + ANORM = DLANST( 'I', LEND-L+1, D( L ), E( L ) ) + ISCALE = 0 + IF( ANORM.GT.SSFMAX ) THEN + ISCALE = 1 + CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N, + $ INFO ) + CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N, + $ INFO ) + ELSE IF( ANORM.LT.SSFMIN ) THEN + ISCALE = 2 + CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N, + $ INFO ) + CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N, + $ INFO ) + END IF + * + DO 40 I = L, LEND - 1 + E( I ) = E( I )**2 + 40 CONTINUE + * + * Choose between QL and QR iteration + * + IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN + LEND = LSV + L = LENDSV + END IF + * + IF( LEND.GE.L ) THEN + * + * QL Iteration + * + * Look for small subdiagonal element. + * + 50 CONTINUE + IF( L.NE.LEND ) THEN + LENDM1 = LEND - 1 + DO 60 M = L, LENDM1 + TST = ABS( E( M ) ) + IF( TST.LE.EPS2*ABS( D( M )*D( M+1 ) ) ) + $ GO TO 70 + 60 CONTINUE + END IF + * + M = LEND + * + 70 CONTINUE + IF( M.LT.LEND ) + $ E( M ) = ZERO + P = D( L ) + IF( M.EQ.L ) + $ GO TO 90 + * + * If remaining matrix is 2 by 2, use DLAE2 to compute its + * eigenvalues. + * + IF( M.EQ.L+1 ) THEN + RTE = SQRT( E( L ) ) + CALL DLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 ) + D( L ) = RT1 + D( L+1 ) = RT2 + E( L ) = ZERO + L = L + 2 + IF( L.LE.LEND ) + $ GO TO 50 + GO TO 150 + END IF + * + IF( JTOT.EQ.NMAXIT ) + $ GO TO 150 + JTOT = JTOT + 1 + * + * Form shift. + * + RTE = SQRT( E( L ) ) + SIGMA = ( D( L+1 )-P ) / ( TWO*RTE ) + R = DLAPY2( SIGMA, ONE ) + SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) ) + * + C = ONE + S = ZERO + GAMMA = D( M ) - SIGMA + P = GAMMA*GAMMA + * + * Inner loop + * + MM1 = M - 1 + DO 80 I = MM1, L, -1 + BB = E( I ) + R = P + BB + IF( I.NE.M-1 ) + $ E( I+1 ) = S*R + OLDC = C + C = P / R + S = BB / R + OLDGAM = GAMMA + ALPHA = D( I ) + GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM + D( I+1 ) = OLDGAM + ( ALPHA-GAMMA ) + IF( C.NE.ZERO ) THEN + P = ( GAMMA*GAMMA ) / C + ELSE + P = OLDC*BB + END IF + 80 CONTINUE + * + E( L ) = S*P + D( L ) = SIGMA + GAMMA + GO TO 50 + * + * Eigenvalue found. + * + 90 CONTINUE + D( L ) = P + * + L = L + 1 + IF( L.LE.LEND ) + $ GO TO 50 + GO TO 150 + * + ELSE + * + * QR Iteration + * + * Look for small superdiagonal element. + * + 100 CONTINUE + IF( L.NE.LEND ) THEN + LENDP1 = LEND + 1 + DO 110 M = L, LENDP1, -1 + TST = ABS( E( M-1 ) ) + IF( TST.LE.EPS2*ABS( D( M )*D( M-1 ) ) ) + $ GO TO 120 + 110 CONTINUE + END IF + * + M = LEND + * + 120 CONTINUE + IF( M.GT.LEND ) + $ E( M-1 ) = ZERO + P = D( L ) + IF( M.EQ.L ) + $ GO TO 140 + * + * If remaining matrix is 2 by 2, use DLAE2 to compute its + * eigenvalues. + * + IF( M.EQ.L-1 ) THEN + RTE = SQRT( E( L-1 ) ) + CALL DLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 ) + D( L ) = RT1 + D( L-1 ) = RT2 + E( L-1 ) = ZERO + L = L - 2 + IF( L.GE.LEND ) + $ GO TO 100 + GO TO 150 + END IF + * + IF( JTOT.EQ.NMAXIT ) + $ GO TO 150 + JTOT = JTOT + 1 + * + * Form shift. + * + RTE = SQRT( E( L-1 ) ) + SIGMA = ( D( L-1 )-P ) / ( TWO*RTE ) + R = DLAPY2( SIGMA, ONE ) + SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) ) + * + C = ONE + S = ZERO + GAMMA = D( M ) - SIGMA + P = GAMMA*GAMMA + * + * Inner loop + * + LM1 = L - 1 + DO 130 I = M, LM1 + BB = E( I ) + R = P + BB + IF( I.NE.M ) + $ E( I-1 ) = S*R + OLDC = C + C = P / R + S = BB / R + OLDGAM = GAMMA + ALPHA = D( I+1 ) + GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM + D( I ) = OLDGAM + ( ALPHA-GAMMA ) + IF( C.NE.ZERO ) THEN + P = ( GAMMA*GAMMA ) / C + ELSE + P = OLDC*BB + END IF + 130 CONTINUE + * + E( LM1 ) = S*P + D( L ) = SIGMA + GAMMA + GO TO 100 + * + * Eigenvalue found. + * + 140 CONTINUE + D( L ) = P + * + L = L - 1 + IF( L.GE.LEND ) + $ GO TO 100 + GO TO 150 + * + END IF + * + * Undo scaling if necessary + * + 150 CONTINUE + IF( ISCALE.EQ.1 ) + $ CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1, + $ D( LSV ), N, INFO ) + IF( ISCALE.EQ.2 ) + $ CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1, + $ D( LSV ), N, INFO ) + * + * Check for no convergence to an eigenvalue after a total + * of N*MAXIT iterations. + * + IF( JTOT.EQ.NMAXIT ) THEN + DO 160 I = 1, N - 1 + IF( E( I ).NE.ZERO ) + $ INFO = INFO + 1 + 160 CONTINUE + RETURN + END IF + GO TO 10 + * + * Sort eigenvalues in increasing order. + * + 170 CONTINUE + CALL DLASRT( 'I', N, D, INFO ) + * + RETURN + * + * End of DSTERF + * + END diff -cNr octave-2.0.7/libcruft/lapack/dsyev.f octave-2.0.8/libcruft/lapack/dsyev.f *** octave-2.0.7/libcruft/lapack/dsyev.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/dsyev.f Mon Jun 23 14:00:08 1997 *************** *** 0 **** --- 1,198 ---- + SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO ) + * + * -- LAPACK driver routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * September 30, 1994 + * + * .. Scalar Arguments .. + CHARACTER JOBZ, UPLO + INTEGER INFO, LDA, LWORK, N + * .. + * .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * ) + * .. + * + * Purpose + * ======= + * + * DSYEV computes all eigenvalues and, optionally, eigenvectors of a + * real symmetric matrix A. + * + * Arguments + * ========= + * + * JOBZ (input) CHARACTER*1 + * = 'N': Compute eigenvalues only; + * = 'V': Compute eigenvalues and eigenvectors. + * + * UPLO (input) CHARACTER*1 + * = 'U': Upper triangle of A is stored; + * = 'L': Lower triangle of A is stored. + * + * N (input) INTEGER + * The order of the matrix A. N >= 0. + * + * A (input/output) DOUBLE PRECISION array, dimension (LDA, N) + * On entry, the symmetric matrix A. If UPLO = 'U', the + * leading N-by-N upper triangular part of A contains the + * upper triangular part of the matrix A. If UPLO = 'L', + * the leading N-by-N lower triangular part of A contains + * the lower triangular part of the matrix A. + * On exit, if JOBZ = 'V', then if INFO = 0, A contains the + * orthonormal eigenvectors of the matrix A. + * If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') + * or the upper triangle (if UPLO='U') of A, including the + * diagonal, is destroyed. + * + * LDA (input) INTEGER + * The leading dimension of the array A. LDA >= max(1,N). + * + * W (output) DOUBLE PRECISION array, dimension (N) + * If INFO = 0, the eigenvalues in ascending order. + * + * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) + * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. + * + * LWORK (input) INTEGER + * The length of the array WORK. LWORK >= max(1,3*N-1). + * For optimal efficiency, LWORK >= (NB+2)*N, + * where NB is the blocksize for DSYTRD returned by ILAENV. + * + * INFO (output) INTEGER + * = 0: successful exit + * < 0: if INFO = -i, the i-th argument had an illegal value + * > 0: if INFO = i, the algorithm failed to converge; i + * off-diagonal elements of an intermediate tridiagonal + * form did not converge to zero. + * + * ===================================================================== + * + * .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) + * .. + * .. Local Scalars .. + LOGICAL LOWER, WANTZ + INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE, + $ LLWORK, LOPT + DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, + $ SMLNUM + * .. + * .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DLAMCH, DLANSY + EXTERNAL LSAME, DLAMCH, DLANSY + * .. + * .. External Subroutines .. + EXTERNAL DLASCL, DORGTR, DSCAL, DSTEQR, DSTERF, DSYTRD, + $ XERBLA + * .. + * .. Intrinsic Functions .. + INTRINSIC MAX, SQRT + * .. + * .. Executable Statements .. + * + * Test the input parameters. + * + WANTZ = LSAME( JOBZ, 'V' ) + LOWER = LSAME( UPLO, 'L' ) + * + INFO = 0 + IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN + INFO = -1 + ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + ELSE IF( LWORK.LT.MAX( 1, 3*N-1 ) ) THEN + INFO = -8 + END IF + * + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DSYEV ', -INFO ) + RETURN + END IF + * + * Quick return if possible + * + IF( N.EQ.0 ) THEN + WORK( 1 ) = 1 + RETURN + END IF + * + IF( N.EQ.1 ) THEN + W( 1 ) = A( 1, 1 ) + WORK( 1 ) = 3 + IF( WANTZ ) + $ A( 1, 1 ) = ONE + RETURN + END IF + * + * Get machine constants. + * + SAFMIN = DLAMCH( 'Safe minimum' ) + EPS = DLAMCH( 'Precision' ) + SMLNUM = SAFMIN / EPS + BIGNUM = ONE / SMLNUM + RMIN = SQRT( SMLNUM ) + RMAX = SQRT( BIGNUM ) + * + * Scale matrix to allowable range, if necessary. + * + ANRM = DLANSY( 'M', UPLO, N, A, LDA, WORK ) + ISCALE = 0 + IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN + ISCALE = 1 + SIGMA = RMIN / ANRM + ELSE IF( ANRM.GT.RMAX ) THEN + ISCALE = 1 + SIGMA = RMAX / ANRM + END IF + IF( ISCALE.EQ.1 ) + $ CALL DLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO ) + * + * Call DSYTRD to reduce symmetric matrix to tridiagonal form. + * + INDE = 1 + INDTAU = INDE + N + INDWRK = INDTAU + N + LLWORK = LWORK - INDWRK + 1 + CALL DSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ), + $ WORK( INDWRK ), LLWORK, IINFO ) + LOPT = 2*N + WORK( INDWRK ) + * + * For eigenvalues only, call DSTERF. For eigenvectors, first call + * DORGTR to generate the orthogonal matrix, then call DSTEQR. + * + IF( .NOT.WANTZ ) THEN + CALL DSTERF( N, W, WORK( INDE ), INFO ) + ELSE + CALL DORGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ), + $ LLWORK, IINFO ) + CALL DSTEQR( JOBZ, N, W, WORK( INDE ), A, LDA, WORK( INDTAU ), + $ INFO ) + END IF + * + * If matrix was scaled, then rescale eigenvalues appropriately. + * + IF( ISCALE.EQ.1 ) THEN + IF( INFO.EQ.0 ) THEN + IMAX = N + ELSE + IMAX = INFO - 1 + END IF + CALL DSCAL( IMAX, ONE / SIGMA, W, 1 ) + END IF + * + * Set WORK(1) to optimal workspace size. + * + WORK( 1 ) = MAX( 3*N-1, LOPT ) + * + RETURN + * + * End of DSYEV + * + END diff -cNr octave-2.0.7/libcruft/lapack/dsytd2.f octave-2.0.8/libcruft/lapack/dsytd2.f *** octave-2.0.7/libcruft/lapack/dsytd2.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/dsytd2.f Mon Jun 23 14:00:08 1997 *************** *** 0 **** --- 1,249 ---- + SUBROUTINE DSYTD2( UPLO, N, A, LDA, D, E, TAU, INFO ) + * + * -- LAPACK routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * October 31, 1992 + * + * .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N + * .. + * .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), TAU( * ) + * .. + * + * Purpose + * ======= + * + * DSYTD2 reduces a real symmetric matrix A to symmetric tridiagonal + * form T by an orthogonal similarity transformation: Q' * A * Q = T. + * + * Arguments + * ========= + * + * UPLO (input) CHARACTER*1 + * Specifies whether the upper or lower triangular part of the + * symmetric matrix A is stored: + * = 'U': Upper triangular + * = 'L': Lower triangular + * + * N (input) INTEGER + * The order of the matrix A. N >= 0. + * + * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) + * On entry, the symmetric matrix A. If UPLO = 'U', the leading + * n-by-n upper triangular part of A contains the upper + * triangular part of the matrix A, and the strictly lower + * triangular part of A is not referenced. If UPLO = 'L', the + * leading n-by-n lower triangular part of A contains the lower + * triangular part of the matrix A, and the strictly upper + * triangular part of A is not referenced. + * On exit, if UPLO = 'U', the diagonal and first superdiagonal + * of A are overwritten by the corresponding elements of the + * tridiagonal matrix T, and the elements above the first + * superdiagonal, with the array TAU, represent the orthogonal + * matrix Q as a product of elementary reflectors; if UPLO + * = 'L', the diagonal and first subdiagonal of A are over- + * written by the corresponding elements of the tridiagonal + * matrix T, and the elements below the first subdiagonal, with + * the array TAU, represent the orthogonal matrix Q as a product + * of elementary reflectors. See Further Details. + * + * LDA (input) INTEGER + * The leading dimension of the array A. LDA >= max(1,N). + * + * D (output) DOUBLE PRECISION array, dimension (N) + * The diagonal elements of the tridiagonal matrix T: + * D(i) = A(i,i). + * + * E (output) DOUBLE PRECISION array, dimension (N-1) + * The off-diagonal elements of the tridiagonal matrix T: + * E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. + * + * TAU (output) DOUBLE PRECISION array, dimension (N-1) + * The scalar factors of the elementary reflectors (see Further + * Details). + * + * INFO (output) INTEGER + * = 0: successful exit + * < 0: if INFO = -i, the i-th argument had an illegal value. + * + * Further Details + * =============== + * + * If UPLO = 'U', the matrix Q is represented as a product of elementary + * reflectors + * + * Q = H(n-1) . . . H(2) H(1). + * + * Each H(i) has the form + * + * H(i) = I - tau * v * v' + * + * where tau is a real scalar, and v is a real vector with + * v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in + * A(1:i-1,i+1), and tau in TAU(i). + * + * If UPLO = 'L', the matrix Q is represented as a product of elementary + * reflectors + * + * Q = H(1) H(2) . . . H(n-1). + * + * Each H(i) has the form + * + * H(i) = I - tau * v * v' + * + * where tau is a real scalar, and v is a real vector with + * v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i), + * and tau in TAU(i). + * + * The contents of A on exit are illustrated by the following examples + * with n = 5: + * + * if UPLO = 'U': if UPLO = 'L': + * + * ( d e v2 v3 v4 ) ( d ) + * ( d e v3 v4 ) ( e d ) + * ( d e v4 ) ( v1 e d ) + * ( d e ) ( v1 v2 e d ) + * ( d ) ( v1 v2 v3 e d ) + * + * where d and e denote diagonal and off-diagonal elements of T, and vi + * denotes an element of the vector defining H(i). + * + * ===================================================================== + * + * .. Parameters .. + DOUBLE PRECISION ONE, ZERO, HALF + PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0, + $ HALF = 1.0D0 / 2.0D0 ) + * .. + * .. Local Scalars .. + LOGICAL UPPER + INTEGER I + DOUBLE PRECISION ALPHA, TAUI + * .. + * .. External Subroutines .. + EXTERNAL DAXPY, DLARFG, DSYMV, DSYR2, XERBLA + * .. + * .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DDOT + EXTERNAL LSAME, DDOT + * .. + * .. Intrinsic Functions .. + INTRINSIC MAX, MIN + * .. + * .. Executable Statements .. + * + * Test the input parameters + * + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DSYTD2', -INFO ) + RETURN + END IF + * + * Quick return if possible + * + IF( N.LE.0 ) + $ RETURN + * + IF( UPPER ) THEN + * + * Reduce the upper triangle of A + * + DO 10 I = N - 1, 1, -1 + * + * Generate elementary reflector H(i) = I - tau * v * v' + * to annihilate A(1:i-1,i+1) + * + CALL DLARFG( I, A( I, I+1 ), A( 1, I+1 ), 1, TAUI ) + E( I ) = A( I, I+1 ) + * + IF( TAUI.NE.ZERO ) THEN + * + * Apply H(i) from both sides to A(1:i,1:i) + * + A( I, I+1 ) = ONE + * + * Compute x := tau * A * v storing x in TAU(1:i) + * + CALL DSYMV( UPLO, I, TAUI, A, LDA, A( 1, I+1 ), 1, ZERO, + $ TAU, 1 ) + * + * Compute w := x - 1/2 * tau * (x'*v) * v + * + ALPHA = -HALF*TAUI*DDOT( I, TAU, 1, A( 1, I+1 ), 1 ) + CALL DAXPY( I, ALPHA, A( 1, I+1 ), 1, TAU, 1 ) + * + * Apply the transformation as a rank-2 update: + * A := A - v * w' - w * v' + * + CALL DSYR2( UPLO, I, -ONE, A( 1, I+1 ), 1, TAU, 1, A, + $ LDA ) + * + A( I, I+1 ) = E( I ) + END IF + D( I+1 ) = A( I+1, I+1 ) + TAU( I ) = TAUI + 10 CONTINUE + D( 1 ) = A( 1, 1 ) + ELSE + * + * Reduce the lower triangle of A + * + DO 20 I = 1, N - 1 + * + * Generate elementary reflector H(i) = I - tau * v * v' + * to annihilate A(i+2:n,i) + * + CALL DLARFG( N-I, A( I+1, I ), A( MIN( I+2, N ), I ), 1, + $ TAUI ) + E( I ) = A( I+1, I ) + * + IF( TAUI.NE.ZERO ) THEN + * + * Apply H(i) from both sides to A(i+1:n,i+1:n) + * + A( I+1, I ) = ONE + * + * Compute x := tau * A * v storing y in TAU(i:n-1) + * + CALL DSYMV( UPLO, N-I, TAUI, A( I+1, I+1 ), LDA, + $ A( I+1, I ), 1, ZERO, TAU( I ), 1 ) + * + * Compute w := x - 1/2 * tau * (x'*v) * v + * + ALPHA = -HALF*TAUI*DDOT( N-I, TAU( I ), 1, A( I+1, I ), + $ 1 ) + CALL DAXPY( N-I, ALPHA, A( I+1, I ), 1, TAU( I ), 1 ) + * + * Apply the transformation as a rank-2 update: + * A := A - v * w' - w * v' + * + CALL DSYR2( UPLO, N-I, -ONE, A( I+1, I ), 1, TAU( I ), 1, + $ A( I+1, I+1 ), LDA ) + * + A( I+1, I ) = E( I ) + END IF + D( I ) = A( I, I ) + TAU( I ) = TAUI + 20 CONTINUE + D( N ) = A( N, N ) + END IF + * + RETURN + * + * End of DSYTD2 + * + END diff -cNr octave-2.0.7/libcruft/lapack/dsytrd.f octave-2.0.8/libcruft/lapack/dsytrd.f *** octave-2.0.7/libcruft/lapack/dsytrd.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/dsytrd.f Mon Jun 23 14:00:08 1997 *************** *** 0 **** --- 1,279 ---- + SUBROUTINE DSYTRD( UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO ) + * + * -- LAPACK routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * September 30, 1994 + * + * .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, LWORK, N + * .. + * .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), TAU( * ), + $ WORK( * ) + * .. + * + * Purpose + * ======= + * + * DSYTRD reduces a real symmetric matrix A to real symmetric + * tridiagonal form T by an orthogonal similarity transformation: + * Q**T * A * Q = T. + * + * Arguments + * ========= + * + * UPLO (input) CHARACTER*1 + * = 'U': Upper triangle of A is stored; + * = 'L': Lower triangle of A is stored. + * + * N (input) INTEGER + * The order of the matrix A. N >= 0. + * + * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) + * On entry, the symmetric matrix A. If UPLO = 'U', the leading + * N-by-N upper triangular part of A contains the upper + * triangular part of the matrix A, and the strictly lower + * triangular part of A is not referenced. If UPLO = 'L', the + * leading N-by-N lower triangular part of A contains the lower + * triangular part of the matrix A, and the strictly upper + * triangular part of A is not referenced. + * On exit, if UPLO = 'U', the diagonal and first superdiagonal + * of A are overwritten by the corresponding elements of the + * tridiagonal matrix T, and the elements above the first + * superdiagonal, with the array TAU, represent the orthogonal + * matrix Q as a product of elementary reflectors; if UPLO + * = 'L', the diagonal and first subdiagonal of A are over- + * written by the corresponding elements of the tridiagonal + * matrix T, and the elements below the first subdiagonal, with + * the array TAU, represent the orthogonal matrix Q as a product + * of elementary reflectors. See Further Details. + * + * LDA (input) INTEGER + * The leading dimension of the array A. LDA >= max(1,N). + * + * D (output) DOUBLE PRECISION array, dimension (N) + * The diagonal elements of the tridiagonal matrix T: + * D(i) = A(i,i). + * + * E (output) DOUBLE PRECISION array, dimension (N-1) + * The off-diagonal elements of the tridiagonal matrix T: + * E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. + * + * TAU (output) DOUBLE PRECISION array, dimension (N-1) + * The scalar factors of the elementary reflectors (see Further + * Details). + * + * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) + * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. + * + * LWORK (input) INTEGER + * The dimension of the array WORK. LWORK >= 1. + * For optimum performance LWORK >= N*NB, where NB is the + * optimal blocksize. + * + * INFO (output) INTEGER + * = 0: successful exit + * < 0: if INFO = -i, the i-th argument had an illegal value + * + * Further Details + * =============== + * + * If UPLO = 'U', the matrix Q is represented as a product of elementary + * reflectors + * + * Q = H(n-1) . . . H(2) H(1). + * + * Each H(i) has the form + * + * H(i) = I - tau * v * v' + * + * where tau is a real scalar, and v is a real vector with + * v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in + * A(1:i-1,i+1), and tau in TAU(i). + * + * If UPLO = 'L', the matrix Q is represented as a product of elementary + * reflectors + * + * Q = H(1) H(2) . . . H(n-1). + * + * Each H(i) has the form + * + * H(i) = I - tau * v * v' + * + * where tau is a real scalar, and v is a real vector with + * v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i), + * and tau in TAU(i). + * + * The contents of A on exit are illustrated by the following examples + * with n = 5: + * + * if UPLO = 'U': if UPLO = 'L': + * + * ( d e v2 v3 v4 ) ( d ) + * ( d e v3 v4 ) ( e d ) + * ( d e v4 ) ( v1 e d ) + * ( d e ) ( v1 v2 e d ) + * ( d ) ( v1 v2 v3 e d ) + * + * where d and e denote diagonal and off-diagonal elements of T, and vi + * denotes an element of the vector defining H(i). + * + * ===================================================================== + * + * .. Parameters .. + DOUBLE PRECISION ONE + PARAMETER ( ONE = 1.0D0 ) + * .. + * .. Local Scalars .. + LOGICAL UPPER + INTEGER I, IINFO, IWS, J, KK, LDWORK, NB, NBMIN, NX + * .. + * .. External Subroutines .. + EXTERNAL DLATRD, DSYR2K, DSYTD2, XERBLA + * .. + * .. Intrinsic Functions .. + INTRINSIC MAX + * .. + * .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV + * .. + * .. Executable Statements .. + * + * Test the input parameters + * + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + ELSE IF( LWORK.LT.1 ) THEN + INFO = -9 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DSYTRD', -INFO ) + RETURN + END IF + * + * Quick return if possible + * + IF( N.EQ.0 ) THEN + WORK( 1 ) = 1 + RETURN + END IF + * + * Determine the block size. + * + NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 ) + NX = N + IWS = 1 + IF( NB.GT.1 .AND. NB.LT.N ) THEN + * + * Determine when to cross over from blocked to unblocked code + * (last block is always handled by unblocked code). + * + NX = MAX( NB, ILAENV( 3, 'DSYTRD', UPLO, N, -1, -1, -1 ) ) + IF( NX.LT.N ) THEN + * + * Determine if workspace is large enough for blocked code. + * + LDWORK = N + IWS = LDWORK*NB + IF( LWORK.LT.IWS ) THEN + * + * Not enough workspace to use optimal NB: determine the + * minimum value of NB, and reduce NB or force use of + * unblocked code by setting NX = N. + * + NB = MAX( LWORK / LDWORK, 1 ) + NBMIN = ILAENV( 2, 'DSYTRD', UPLO, N, -1, -1, -1 ) + IF( NB.LT.NBMIN ) + $ NX = N + END IF + ELSE + NX = N + END IF + ELSE + NB = 1 + END IF + * + IF( UPPER ) THEN + * + * Reduce the upper triangle of A. + * Columns 1:kk are handled by the unblocked method. + * + KK = N - ( ( N-NX+NB-1 ) / NB )*NB + DO 20 I = N - NB + 1, KK + 1, -NB + * + * Reduce columns i:i+nb-1 to tridiagonal form and form the + * matrix W which is needed to update the unreduced part of + * the matrix + * + CALL DLATRD( UPLO, I+NB-1, NB, A, LDA, E, TAU, WORK, + $ LDWORK ) + * + * Update the unreduced submatrix A(1:i-1,1:i-1), using an + * update of the form: A := A - V*W' - W*V' + * + CALL DSYR2K( UPLO, 'No transpose', I-1, NB, -ONE, A( 1, I ), + $ LDA, WORK, LDWORK, ONE, A, LDA ) + * + * Copy superdiagonal elements back into A, and diagonal + * elements into D + * + DO 10 J = I, I + NB - 1 + A( J-1, J ) = E( J-1 ) + D( J ) = A( J, J ) + 10 CONTINUE + 20 CONTINUE + * + * Use unblocked code to reduce the last or only block + * + CALL DSYTD2( UPLO, KK, A, LDA, D, E, TAU, IINFO ) + ELSE + * + * Reduce the lower triangle of A + * + DO 40 I = 1, N - NX, NB + * + * Reduce columns i:i+nb-1 to tridiagonal form and form the + * matrix W which is needed to update the unreduced part of + * the matrix + * + CALL DLATRD( UPLO, N-I+1, NB, A( I, I ), LDA, E( I ), + $ TAU( I ), WORK, LDWORK ) + * + * Update the unreduced submatrix A(i+ib:n,i+ib:n), using + * an update of the form: A := A - V*W' - W*V' + * + CALL DSYR2K( UPLO, 'No transpose', N-I-NB+1, NB, -ONE, + $ A( I+NB, I ), LDA, WORK( NB+1 ), LDWORK, ONE, + $ A( I+NB, I+NB ), LDA ) + * + * Copy subdiagonal elements back into A, and diagonal + * elements into D + * + DO 30 J = I, I + NB - 1 + A( J+1, J ) = E( J ) + D( J ) = A( J, J ) + 30 CONTINUE + 40 CONTINUE + * + * Use unblocked code to reduce the last or only block + * + CALL DSYTD2( UPLO, N-I+1, A( I, I ), LDA, D( I ), E( I ), + $ TAU( I ), IINFO ) + END IF + * + WORK( 1 ) = IWS + RETURN + * + * End of DSYTRD + * + END diff -cNr octave-2.0.7/libcruft/lapack/zheev.f octave-2.0.8/libcruft/lapack/zheev.f *** octave-2.0.7/libcruft/lapack/zheev.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/zheev.f Mon Jun 23 13:59:25 1997 *************** *** 0 **** --- 1,205 ---- + SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, + $ INFO ) + * + * -- LAPACK driver routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * September 30, 1994 + * + * .. Scalar Arguments .. + CHARACTER JOBZ, UPLO + INTEGER INFO, LDA, LWORK, N + * .. + * .. Array Arguments .. + DOUBLE PRECISION RWORK( * ), W( * ) + COMPLEX*16 A( LDA, * ), WORK( * ) + * .. + * + * Purpose + * ======= + * + * ZHEEV computes all eigenvalues and, optionally, eigenvectors of a + * complex Hermitian matrix A. + * + * Arguments + * ========= + * + * JOBZ (input) CHARACTER*1 + * = 'N': Compute eigenvalues only; + * = 'V': Compute eigenvalues and eigenvectors. + * + * UPLO (input) CHARACTER*1 + * = 'U': Upper triangle of A is stored; + * = 'L': Lower triangle of A is stored. + * + * N (input) INTEGER + * The order of the matrix A. N >= 0. + * + * A (input/output) COMPLEX*16 array, dimension (LDA, N) + * On entry, the Hermitian matrix A. If UPLO = 'U', the + * leading N-by-N upper triangular part of A contains the + * upper triangular part of the matrix A. If UPLO = 'L', + * the leading N-by-N lower triangular part of A contains + * the lower triangular part of the matrix A. + * On exit, if JOBZ = 'V', then if INFO = 0, A contains the + * orthonormal eigenvectors of the matrix A. + * If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') + * or the upper triangle (if UPLO='U') of A, including the + * diagonal, is destroyed. + * + * LDA (input) INTEGER + * The leading dimension of the array A. LDA >= max(1,N). + * + * W (output) DOUBLE PRECISION array, dimension (N) + * If INFO = 0, the eigenvalues in ascending order. + * + * WORK (workspace/output) COMPLEX*16 array, dimension (LWORK) + * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. + * + * LWORK (input) INTEGER + * The length of the array WORK. LWORK >= max(1,2*N-1). + * For optimal efficiency, LWORK >= (NB+1)*N, + * where NB is the blocksize for ZHETRD returned by ILAENV. + * + * RWORK (workspace) DOUBLE PRECISION array, dimension (max(1, 3*N-2)) + * + * INFO (output) INTEGER + * = 0: successful exit + * < 0: if INFO = -i, the i-th argument had an illegal value + * > 0: if INFO = i, the algorithm failed to converge; i + * off-diagonal elements of an intermediate tridiagonal + * form did not converge to zero. + * + * ===================================================================== + * + * .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) + COMPLEX*16 CONE + PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) ) + * .. + * .. Local Scalars .. + LOGICAL LOWER, WANTZ + INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE, + $ LLWORK, LOPT + DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, + $ SMLNUM + * .. + * .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DLAMCH, ZLANHE + EXTERNAL LSAME, DLAMCH, ZLANHE + * .. + * .. External Subroutines .. + EXTERNAL DSCAL, DSTERF, XERBLA, ZHETRD, ZLASCL, ZSTEQR, + $ ZUNGTR + * .. + * .. Intrinsic Functions .. + INTRINSIC MAX, SQRT + * .. + * .. Executable Statements .. + * + * Test the input parameters. + * + WANTZ = LSAME( JOBZ, 'V' ) + LOWER = LSAME( UPLO, 'L' ) + * + INFO = 0 + IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN + INFO = -1 + ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + ELSE IF( LWORK.LT.MAX( 1, 2*N-1 ) ) THEN + INFO = -8 + END IF + * + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZHEEV ', -INFO ) + RETURN + END IF + * + * Quick return if possible + * + IF( N.EQ.0 ) THEN + WORK( 1 ) = 1 + RETURN + END IF + * + IF( N.EQ.1 ) THEN + W( 1 ) = A( 1, 1 ) + WORK( 1 ) = 3 + IF( WANTZ ) + $ A( 1, 1 ) = CONE + RETURN + END IF + * + * Get machine constants. + * + SAFMIN = DLAMCH( 'Safe minimum' ) + EPS = DLAMCH( 'Precision' ) + SMLNUM = SAFMIN / EPS + BIGNUM = ONE / SMLNUM + RMIN = SQRT( SMLNUM ) + RMAX = SQRT( BIGNUM ) + * + * Scale matrix to allowable range, if necessary. + * + ANRM = ZLANHE( 'M', UPLO, N, A, LDA, RWORK ) + ISCALE = 0 + IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN + ISCALE = 1 + SIGMA = RMIN / ANRM + ELSE IF( ANRM.GT.RMAX ) THEN + ISCALE = 1 + SIGMA = RMAX / ANRM + END IF + IF( ISCALE.EQ.1 ) + $ CALL ZLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO ) + * + * Call ZHETRD to reduce Hermitian matrix to tridiagonal form. + * + INDE = 1 + INDTAU = 1 + INDWRK = INDTAU + N + LLWORK = LWORK - INDWRK + 1 + CALL ZHETRD( UPLO, N, A, LDA, W, RWORK( INDE ), WORK( INDTAU ), + $ WORK( INDWRK ), LLWORK, IINFO ) + LOPT = N + WORK( INDWRK ) + * + * For eigenvalues only, call DSTERF. For eigenvectors, first call + * ZUNGTR to generate the unitary matrix, then call ZSTEQR. + * + IF( .NOT.WANTZ ) THEN + CALL DSTERF( N, W, RWORK( INDE ), INFO ) + ELSE + CALL ZUNGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ), + $ LLWORK, IINFO ) + INDWRK = INDE + N + CALL ZSTEQR( JOBZ, N, W, RWORK( INDE ), A, LDA, + $ RWORK( INDWRK ), INFO ) + END IF + * + * If matrix was scaled, then rescale eigenvalues appropriately. + * + IF( ISCALE.EQ.1 ) THEN + IF( INFO.EQ.0 ) THEN + IMAX = N + ELSE + IMAX = INFO - 1 + END IF + CALL DSCAL( IMAX, ONE / SIGMA, W, 1 ) + END IF + * + * Set WORK(1) to optimal complex workspace size. + * + WORK( 1 ) = MAX( 2*N-1, LOPT ) + * + RETURN + * + * End of ZHEEV + * + END diff -cNr octave-2.0.7/libcruft/lapack/zhetd2.f octave-2.0.8/libcruft/lapack/zhetd2.f *** octave-2.0.7/libcruft/lapack/zhetd2.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/zhetd2.f Mon Jun 23 13:59:25 1997 *************** *** 0 **** --- 1,255 ---- + SUBROUTINE ZHETD2( UPLO, N, A, LDA, D, E, TAU, INFO ) + * + * -- LAPACK routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * September 30, 1994 + * + * .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N + * .. + * .. Array Arguments .. + DOUBLE PRECISION D( * ), E( * ) + COMPLEX*16 A( LDA, * ), TAU( * ) + * .. + * + * Purpose + * ======= + * + * ZHETD2 reduces a complex Hermitian matrix A to real symmetric + * tridiagonal form T by a unitary similarity transformation: + * Q' * A * Q = T. + * + * Arguments + * ========= + * + * UPLO (input) CHARACTER*1 + * Specifies whether the upper or lower triangular part of the + * Hermitian matrix A is stored: + * = 'U': Upper triangular + * = 'L': Lower triangular + * + * N (input) INTEGER + * The order of the matrix A. N >= 0. + * + * A (input/output) COMPLEX*16 array, dimension (LDA,N) + * On entry, the Hermitian matrix A. If UPLO = 'U', the leading + * n-by-n upper triangular part of A contains the upper + * triangular part of the matrix A, and the strictly lower + * triangular part of A is not referenced. If UPLO = 'L', the + * leading n-by-n lower triangular part of A contains the lower + * triangular part of the matrix A, and the strictly upper + * triangular part of A is not referenced. + * On exit, if UPLO = 'U', the diagonal and first superdiagonal + * of A are overwritten by the corresponding elements of the + * tridiagonal matrix T, and the elements above the first + * superdiagonal, with the array TAU, represent the unitary + * matrix Q as a product of elementary reflectors; if UPLO + * = 'L', the diagonal and first subdiagonal of A are over- + * written by the corresponding elements of the tridiagonal + * matrix T, and the elements below the first subdiagonal, with + * the array TAU, represent the unitary matrix Q as a product + * of elementary reflectors. See Further Details. + * + * LDA (input) INTEGER + * The leading dimension of the array A. LDA >= max(1,N). + * + * D (output) DOUBLE PRECISION array, dimension (N) + * The diagonal elements of the tridiagonal matrix T: + * D(i) = A(i,i). + * + * E (output) DOUBLE PRECISION array, dimension (N-1) + * The off-diagonal elements of the tridiagonal matrix T: + * E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. + * + * TAU (output) COMPLEX*16 array, dimension (N-1) + * The scalar factors of the elementary reflectors (see Further + * Details). + * + * INFO (output) INTEGER + * = 0: successful exit + * < 0: if INFO = -i, the i-th argument had an illegal value. + * + * Further Details + * =============== + * + * If UPLO = 'U', the matrix Q is represented as a product of elementary + * reflectors + * + * Q = H(n-1) . . . H(2) H(1). + * + * Each H(i) has the form + * + * H(i) = I - tau * v * v' + * + * where tau is a complex scalar, and v is a complex vector with + * v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in + * A(1:i-1,i+1), and tau in TAU(i). + * + * If UPLO = 'L', the matrix Q is represented as a product of elementary + * reflectors + * + * Q = H(1) H(2) . . . H(n-1). + * + * Each H(i) has the form + * + * H(i) = I - tau * v * v' + * + * where tau is a complex scalar, and v is a complex vector with + * v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i), + * and tau in TAU(i). + * + * The contents of A on exit are illustrated by the following examples + * with n = 5: + * + * if UPLO = 'U': if UPLO = 'L': + * + * ( d e v2 v3 v4 ) ( d ) + * ( d e v3 v4 ) ( e d ) + * ( d e v4 ) ( v1 e d ) + * ( d e ) ( v1 v2 e d ) + * ( d ) ( v1 v2 v3 e d ) + * + * where d and e denote diagonal and off-diagonal elements of T, and vi + * denotes an element of the vector defining H(i). + * + * ===================================================================== + * + * .. Parameters .. + COMPLEX*16 ONE, ZERO, HALF + PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ), + $ ZERO = ( 0.0D+0, 0.0D+0 ), + $ HALF = ( 0.5D+0, 0.0D+0 ) ) + * .. + * .. Local Scalars .. + LOGICAL UPPER + INTEGER I + COMPLEX*16 ALPHA, TAUI + * .. + * .. External Subroutines .. + EXTERNAL XERBLA, ZAXPY, ZHEMV, ZHER2, ZLARFG + * .. + * .. External Functions .. + LOGICAL LSAME + COMPLEX*16 ZDOTC + EXTERNAL LSAME, ZDOTC + * .. + * .. Intrinsic Functions .. + INTRINSIC DBLE, MAX, MIN + * .. + * .. Executable Statements .. + * + * Test the input parameters + * + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZHETD2', -INFO ) + RETURN + END IF + * + * Quick return if possible + * + IF( N.LE.0 ) + $ RETURN + * + IF( UPPER ) THEN + * + * Reduce the upper triangle of A + * + A( N, N ) = DBLE( A( N, N ) ) + DO 10 I = N - 1, 1, -1 + * + * Generate elementary reflector H(i) = I - tau * v * v' + * to annihilate A(1:i-1,i+1) + * + ALPHA = A( I, I+1 ) + CALL ZLARFG( I, ALPHA, A( 1, I+1 ), 1, TAUI ) + E( I ) = ALPHA + * + IF( TAUI.NE.ZERO ) THEN + * + * Apply H(i) from both sides to A(1:i,1:i) + * + A( I, I+1 ) = ONE + * + * Compute x := tau * A * v storing x in TAU(1:i) + * + CALL ZHEMV( UPLO, I, TAUI, A, LDA, A( 1, I+1 ), 1, ZERO, + $ TAU, 1 ) + * + * Compute w := x - 1/2 * tau * (x'*v) * v + * + ALPHA = -HALF*TAUI*ZDOTC( I, TAU, 1, A( 1, I+1 ), 1 ) + CALL ZAXPY( I, ALPHA, A( 1, I+1 ), 1, TAU, 1 ) + * + * Apply the transformation as a rank-2 update: + * A := A - v * w' - w * v' + * + CALL ZHER2( UPLO, I, -ONE, A( 1, I+1 ), 1, TAU, 1, A, + $ LDA ) + * + END IF + A( I, I+1 ) = E( I ) + D( I+1 ) = A( I+1, I+1 ) + TAU( I ) = TAUI + 10 CONTINUE + D( 1 ) = A( 1, 1 ) + ELSE + * + * Reduce the lower triangle of A + * + A( 1, 1 ) = DBLE( A( 1, 1 ) ) + DO 20 I = 1, N - 1 + * + * Generate elementary reflector H(i) = I - tau * v * v' + * to annihilate A(i+2:n,i) + * + ALPHA = A( I+1, I ) + CALL ZLARFG( N-I, ALPHA, A( MIN( I+2, N ), I ), 1, TAUI ) + E( I ) = ALPHA + * + IF( TAUI.NE.ZERO ) THEN + * + * Apply H(i) from both sides to A(i+1:n,i+1:n) + * + A( I+1, I ) = ONE + * + * Compute x := tau * A * v storing y in TAU(i:n-1) + * + CALL ZHEMV( UPLO, N-I, TAUI, A( I+1, I+1 ), LDA, + $ A( I+1, I ), 1, ZERO, TAU( I ), 1 ) + * + * Compute w := x - 1/2 * tau * (x'*v) * v + * + ALPHA = -HALF*TAUI*ZDOTC( N-I, TAU( I ), 1, A( I+1, I ), + $ 1 ) + CALL ZAXPY( N-I, ALPHA, A( I+1, I ), 1, TAU( I ), 1 ) + * + * Apply the transformation as a rank-2 update: + * A := A - v * w' - w * v' + * + CALL ZHER2( UPLO, N-I, -ONE, A( I+1, I ), 1, TAU( I ), 1, + $ A( I+1, I+1 ), LDA ) + * + END IF + A( I+1, I ) = E( I ) + D( I ) = A( I, I ) + TAU( I ) = TAUI + 20 CONTINUE + D( N ) = A( N, N ) + END IF + * + RETURN + * + * End of ZHETD2 + * + END diff -cNr octave-2.0.7/libcruft/lapack/zhetrd.f octave-2.0.8/libcruft/lapack/zhetrd.f *** octave-2.0.7/libcruft/lapack/zhetrd.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/zhetrd.f Mon Jun 23 13:59:25 1997 *************** *** 0 **** --- 1,281 ---- + SUBROUTINE ZHETRD( UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO ) + * + * -- LAPACK routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * September 30, 1994 + * + * .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, LWORK, N + * .. + * .. Array Arguments .. + DOUBLE PRECISION D( * ), E( * ) + COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * ) + * .. + * + * Purpose + * ======= + * + * ZHETRD reduces a complex Hermitian matrix A to real symmetric + * tridiagonal form T by a unitary similarity transformation: + * Q**H * A * Q = T. + * + * Arguments + * ========= + * + * UPLO (input) CHARACTER*1 + * = 'U': Upper triangle of A is stored; + * = 'L': Lower triangle of A is stored. + * + * N (input) INTEGER + * The order of the matrix A. N >= 0. + * + * A (input/output) COMPLEX*16 array, dimension (LDA,N) + * On entry, the Hermitian matrix A. If UPLO = 'U', the leading + * N-by-N upper triangular part of A contains the upper + * triangular part of the matrix A, and the strictly lower + * triangular part of A is not referenced. If UPLO = 'L', the + * leading N-by-N lower triangular part of A contains the lower + * triangular part of the matrix A, and the strictly upper + * triangular part of A is not referenced. + * On exit, if UPLO = 'U', the diagonal and first superdiagonal + * of A are overwritten by the corresponding elements of the + * tridiagonal matrix T, and the elements above the first + * superdiagonal, with the array TAU, represent the unitary + * matrix Q as a product of elementary reflectors; if UPLO + * = 'L', the diagonal and first subdiagonal of A are over- + * written by the corresponding elements of the tridiagonal + * matrix T, and the elements below the first subdiagonal, with + * the array TAU, represent the unitary matrix Q as a product + * of elementary reflectors. See Further Details. + * + * LDA (input) INTEGER + * The leading dimension of the array A. LDA >= max(1,N). + * + * D (output) DOUBLE PRECISION array, dimension (N) + * The diagonal elements of the tridiagonal matrix T: + * D(i) = A(i,i). + * + * E (output) DOUBLE PRECISION array, dimension (N-1) + * The off-diagonal elements of the tridiagonal matrix T: + * E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. + * + * TAU (output) COMPLEX*16 array, dimension (N-1) + * The scalar factors of the elementary reflectors (see Further + * Details). + * + * WORK (workspace/output) COMPLEX*16 array, dimension (LWORK) + * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. + * + * LWORK (input) INTEGER + * The dimension of the array WORK. LWORK >= 1. + * For optimum performance LWORK >= N*NB, where NB is the + * optimal blocksize. + * + * INFO (output) INTEGER + * = 0: successful exit + * < 0: if INFO = -i, the i-th argument had an illegal value + * + * Further Details + * =============== + * + * If UPLO = 'U', the matrix Q is represented as a product of elementary + * reflectors + * + * Q = H(n-1) . . . H(2) H(1). + * + * Each H(i) has the form + * + * H(i) = I - tau * v * v' + * + * where tau is a complex scalar, and v is a complex vector with + * v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in + * A(1:i-1,i+1), and tau in TAU(i). + * + * If UPLO = 'L', the matrix Q is represented as a product of elementary + * reflectors + * + * Q = H(1) H(2) . . . H(n-1). + * + * Each H(i) has the form + * + * H(i) = I - tau * v * v' + * + * where tau is a complex scalar, and v is a complex vector with + * v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i), + * and tau in TAU(i). + * + * The contents of A on exit are illustrated by the following examples + * with n = 5: + * + * if UPLO = 'U': if UPLO = 'L': + * + * ( d e v2 v3 v4 ) ( d ) + * ( d e v3 v4 ) ( e d ) + * ( d e v4 ) ( v1 e d ) + * ( d e ) ( v1 v2 e d ) + * ( d ) ( v1 v2 v3 e d ) + * + * where d and e denote diagonal and off-diagonal elements of T, and vi + * denotes an element of the vector defining H(i). + * + * ===================================================================== + * + * .. Parameters .. + DOUBLE PRECISION ONE + PARAMETER ( ONE = 1.0D+0 ) + COMPLEX*16 CONE + PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) + * .. + * .. Local Scalars .. + LOGICAL UPPER + INTEGER I, IINFO, IWS, J, KK, LDWORK, NB, NBMIN, NX + * .. + * .. External Subroutines .. + EXTERNAL XERBLA, ZHER2K, ZHETD2, ZLATRD + * .. + * .. Intrinsic Functions .. + INTRINSIC MAX + * .. + * .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV + * .. + * .. Executable Statements .. + * + * Test the input parameters + * + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + ELSE IF( LWORK.LT.1 ) THEN + INFO = -9 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZHETRD', -INFO ) + RETURN + END IF + * + * Quick return if possible + * + IF( N.EQ.0 ) THEN + WORK( 1 ) = 1 + RETURN + END IF + * + * Determine the block size. + * + NB = ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 ) + NX = N + IWS = 1 + IF( NB.GT.1 .AND. NB.LT.N ) THEN + * + * Determine when to cross over from blocked to unblocked code + * (last block is always handled by unblocked code). + * + NX = MAX( NB, ILAENV( 3, 'ZHETRD', UPLO, N, -1, -1, -1 ) ) + IF( NX.LT.N ) THEN + * + * Determine if workspace is large enough for blocked code. + * + LDWORK = N + IWS = LDWORK*NB + IF( LWORK.LT.IWS ) THEN + * + * Not enough workspace to use optimal NB: determine the + * minimum value of NB, and reduce NB or force use of + * unblocked code by setting NX = N. + * + NB = MAX( LWORK / LDWORK, 1 ) + NBMIN = ILAENV( 2, 'ZHETRD', UPLO, N, -1, -1, -1 ) + IF( NB.LT.NBMIN ) + $ NX = N + END IF + ELSE + NX = N + END IF + ELSE + NB = 1 + END IF + * + IF( UPPER ) THEN + * + * Reduce the upper triangle of A. + * Columns 1:kk are handled by the unblocked method. + * + KK = N - ( ( N-NX+NB-1 ) / NB )*NB + DO 20 I = N - NB + 1, KK + 1, -NB + * + * Reduce columns i:i+nb-1 to tridiagonal form and form the + * matrix W which is needed to update the unreduced part of + * the matrix + * + CALL ZLATRD( UPLO, I+NB-1, NB, A, LDA, E, TAU, WORK, + $ LDWORK ) + * + * Update the unreduced submatrix A(1:i-1,1:i-1), using an + * update of the form: A := A - V*W' - W*V' + * + CALL ZHER2K( UPLO, 'No transpose', I-1, NB, -CONE, + $ A( 1, I ), LDA, WORK, LDWORK, ONE, A, LDA ) + * + * Copy superdiagonal elements back into A, and diagonal + * elements into D + * + DO 10 J = I, I + NB - 1 + A( J-1, J ) = E( J-1 ) + D( J ) = A( J, J ) + 10 CONTINUE + 20 CONTINUE + * + * Use unblocked code to reduce the last or only block + * + CALL ZHETD2( UPLO, KK, A, LDA, D, E, TAU, IINFO ) + ELSE + * + * Reduce the lower triangle of A + * + DO 40 I = 1, N - NX, NB + * + * Reduce columns i:i+nb-1 to tridiagonal form and form the + * matrix W which is needed to update the unreduced part of + * the matrix + * + CALL ZLATRD( UPLO, N-I+1, NB, A( I, I ), LDA, E( I ), + $ TAU( I ), WORK, LDWORK ) + * + * Update the unreduced submatrix A(i+nb:n,i+nb:n), using + * an update of the form: A := A - V*W' - W*V' + * + CALL ZHER2K( UPLO, 'No transpose', N-I-NB+1, NB, -CONE, + $ A( I+NB, I ), LDA, WORK( NB+1 ), LDWORK, ONE, + $ A( I+NB, I+NB ), LDA ) + * + * Copy subdiagonal elements back into A, and diagonal + * elements into D + * + DO 30 J = I, I + NB - 1 + A( J+1, J ) = E( J ) + D( J ) = A( J, J ) + 30 CONTINUE + 40 CONTINUE + * + * Use unblocked code to reduce the last or only block + * + CALL ZHETD2( UPLO, N-I+1, A( I, I ), LDA, D( I ), E( I ), + $ TAU( I ), IINFO ) + END IF + * + WORK( 1 ) = IWS + RETURN + * + * End of ZHETRD + * + END diff -cNr octave-2.0.7/libcruft/lapack/zlanhe.f octave-2.0.8/libcruft/lapack/zlanhe.f *** octave-2.0.7/libcruft/lapack/zlanhe.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/zlanhe.f Mon Jun 23 13:59:25 1997 *************** *** 0 **** --- 1,188 ---- + DOUBLE PRECISION FUNCTION ZLANHE( NORM, UPLO, N, A, LDA, WORK ) + * + * -- LAPACK auxiliary routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * October 31, 1992 + * + * .. Scalar Arguments .. + CHARACTER NORM, UPLO + INTEGER LDA, N + * .. + * .. Array Arguments .. + DOUBLE PRECISION WORK( * ) + COMPLEX*16 A( LDA, * ) + * .. + * + * Purpose + * ======= + * + * ZLANHE returns the value of the one norm, or the Frobenius norm, or + * the infinity norm, or the element of largest absolute value of a + * complex hermitian matrix A. + * + * Description + * =========== + * + * ZLANHE returns the value + * + * ZLANHE = ( max(abs(A(i,j))), NORM = 'M' or 'm' + * ( + * ( norm1(A), NORM = '1', 'O' or 'o' + * ( + * ( normI(A), NORM = 'I' or 'i' + * ( + * ( normF(A), NORM = 'F', 'f', 'E' or 'e' + * + * where norm1 denotes the one norm of a matrix (maximum column sum), + * normI denotes the infinity norm of a matrix (maximum row sum) and + * normF denotes the Frobenius norm of a matrix (square root of sum of + * squares). Note that max(abs(A(i,j))) is not a matrix norm. + * + * Arguments + * ========= + * + * NORM (input) CHARACTER*1 + * Specifies the value to be returned in ZLANHE as described + * above. + * + * UPLO (input) CHARACTER*1 + * Specifies whether the upper or lower triangular part of the + * hermitian matrix A is to be referenced. + * = 'U': Upper triangular part of A is referenced + * = 'L': Lower triangular part of A is referenced + * + * N (input) INTEGER + * The order of the matrix A. N >= 0. When N = 0, ZLANHE is + * set to zero. + * + * A (input) COMPLEX*16 array, dimension (LDA,N) + * The hermitian matrix A. If UPLO = 'U', the leading n by n + * upper triangular part of A contains the upper triangular part + * of the matrix A, and the strictly lower triangular part of A + * is not referenced. If UPLO = 'L', the leading n by n lower + * triangular part of A contains the lower triangular part of + * the matrix A, and the strictly upper triangular part of A is + * not referenced. Note that the imaginary parts of the diagonal + * elements need not be set and are assumed to be zero. + * + * LDA (input) INTEGER + * The leading dimension of the array A. LDA >= max(N,1). + * + * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK), + * where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, + * WORK is not referenced. + * + * ===================================================================== + * + * .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) + * .. + * .. Local Scalars .. + INTEGER I, J + DOUBLE PRECISION ABSA, SCALE, SUM, VALUE + * .. + * .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME + * .. + * .. External Subroutines .. + EXTERNAL ZLASSQ + * .. + * .. Intrinsic Functions .. + INTRINSIC ABS, DBLE, MAX, SQRT + * .. + * .. Executable Statements .. + * + IF( N.EQ.0 ) THEN + VALUE = ZERO + ELSE IF( LSAME( NORM, 'M' ) ) THEN + * + * Find max(abs(A(i,j))). + * + VALUE = ZERO + IF( LSAME( UPLO, 'U' ) ) THEN + DO 20 J = 1, N + DO 10 I = 1, J - 1 + VALUE = MAX( VALUE, ABS( A( I, J ) ) ) + 10 CONTINUE + VALUE = MAX( VALUE, ABS( DBLE( A( J, J ) ) ) ) + 20 CONTINUE + ELSE + DO 40 J = 1, N + VALUE = MAX( VALUE, ABS( DBLE( A( J, J ) ) ) ) + DO 30 I = J + 1, N + VALUE = MAX( VALUE, ABS( A( I, J ) ) ) + 30 CONTINUE + 40 CONTINUE + END IF + ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR. + $ ( NORM.EQ.'1' ) ) THEN + * + * Find normI(A) ( = norm1(A), since A is hermitian). + * + VALUE = ZERO + IF( LSAME( UPLO, 'U' ) ) THEN + DO 60 J = 1, N + SUM = ZERO + DO 50 I = 1, J - 1 + ABSA = ABS( A( I, J ) ) + SUM = SUM + ABSA + WORK( I ) = WORK( I ) + ABSA + 50 CONTINUE + WORK( J ) = SUM + ABS( DBLE( A( J, J ) ) ) + 60 CONTINUE + DO 70 I = 1, N + VALUE = MAX( VALUE, WORK( I ) ) + 70 CONTINUE + ELSE + DO 80 I = 1, N + WORK( I ) = ZERO + 80 CONTINUE + DO 100 J = 1, N + SUM = WORK( J ) + ABS( DBLE( A( J, J ) ) ) + DO 90 I = J + 1, N + ABSA = ABS( A( I, J ) ) + SUM = SUM + ABSA + WORK( I ) = WORK( I ) + ABSA + 90 CONTINUE + VALUE = MAX( VALUE, SUM ) + 100 CONTINUE + END IF + ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN + * + * Find normF(A). + * + SCALE = ZERO + SUM = ONE + IF( LSAME( UPLO, 'U' ) ) THEN + DO 110 J = 2, N + CALL ZLASSQ( J-1, A( 1, J ), 1, SCALE, SUM ) + 110 CONTINUE + ELSE + DO 120 J = 1, N - 1 + CALL ZLASSQ( N-J, A( J+1, J ), 1, SCALE, SUM ) + 120 CONTINUE + END IF + SUM = 2*SUM + DO 130 I = 1, N + IF( DBLE( A( I, I ) ).NE.ZERO ) THEN + ABSA = ABS( DBLE( A( I, I ) ) ) + IF( SCALE.LT.ABSA ) THEN + SUM = ONE + SUM*( SCALE / ABSA )**2 + SCALE = ABSA + ELSE + SUM = SUM + ( ABSA / SCALE )**2 + END IF + END IF + 130 CONTINUE + VALUE = SCALE*SQRT( SUM ) + END IF + * + ZLANHE = VALUE + RETURN + * + * End of ZLANHE + * + END diff -cNr octave-2.0.7/libcruft/lapack/zlatrd.f octave-2.0.8/libcruft/lapack/zlatrd.f *** octave-2.0.7/libcruft/lapack/zlatrd.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/zlatrd.f Mon Jun 23 13:58:57 1997 *************** *** 0 **** --- 1,280 ---- + SUBROUTINE ZLATRD( UPLO, N, NB, A, LDA, E, TAU, W, LDW ) + * + * -- LAPACK auxiliary routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * September 30, 1994 + * + * .. Scalar Arguments .. + CHARACTER UPLO + INTEGER LDA, LDW, N, NB + * .. + * .. Array Arguments .. + DOUBLE PRECISION E( * ) + COMPLEX*16 A( LDA, * ), TAU( * ), W( LDW, * ) + * .. + * + * Purpose + * ======= + * + * ZLATRD reduces NB rows and columns of a complex Hermitian matrix A to + * Hermitian tridiagonal form by a unitary similarity + * transformation Q' * A * Q, and returns the matrices V and W which are + * needed to apply the transformation to the unreduced part of A. + * + * If UPLO = 'U', ZLATRD reduces the last NB rows and columns of a + * matrix, of which the upper triangle is supplied; + * if UPLO = 'L', ZLATRD reduces the first NB rows and columns of a + * matrix, of which the lower triangle is supplied. + * + * This is an auxiliary routine called by ZHETRD. + * + * Arguments + * ========= + * + * UPLO (input) CHARACTER + * Specifies whether the upper or lower triangular part of the + * Hermitian matrix A is stored: + * = 'U': Upper triangular + * = 'L': Lower triangular + * + * N (input) INTEGER + * The order of the matrix A. + * + * NB (input) INTEGER + * The number of rows and columns to be reduced. + * + * A (input/output) COMPLEX*16 array, dimension (LDA,N) + * On entry, the Hermitian matrix A. If UPLO = 'U', the leading + * n-by-n upper triangular part of A contains the upper + * triangular part of the matrix A, and the strictly lower + * triangular part of A is not referenced. If UPLO = 'L', the + * leading n-by-n lower triangular part of A contains the lower + * triangular part of the matrix A, and the strictly upper + * triangular part of A is not referenced. + * On exit: + * if UPLO = 'U', the last NB columns have been reduced to + * tridiagonal form, with the diagonal elements overwriting + * the diagonal elements of A; the elements above the diagonal + * with the array TAU, represent the unitary matrix Q as a + * product of elementary reflectors; + * if UPLO = 'L', the first NB columns have been reduced to + * tridiagonal form, with the diagonal elements overwriting + * the diagonal elements of A; the elements below the diagonal + * with the array TAU, represent the unitary matrix Q as a + * product of elementary reflectors. + * See Further Details. + * + * LDA (input) INTEGER + * The leading dimension of the array A. LDA >= max(1,N). + * + * E (output) DOUBLE PRECISION array, dimension (N-1) + * If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal + * elements of the last NB columns of the reduced matrix; + * if UPLO = 'L', E(1:nb) contains the subdiagonal elements of + * the first NB columns of the reduced matrix. + * + * TAU (output) COMPLEX*16 array, dimension (N-1) + * The scalar factors of the elementary reflectors, stored in + * TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'. + * See Further Details. + * + * W (output) COMPLEX*16 array, dimension (LDW,NB) + * The n-by-nb matrix W required to update the unreduced part + * of A. + * + * LDW (input) INTEGER + * The leading dimension of the array W. LDW >= max(1,N). + * + * Further Details + * =============== + * + * If UPLO = 'U', the matrix Q is represented as a product of elementary + * reflectors + * + * Q = H(n) H(n-1) . . . H(n-nb+1). + * + * Each H(i) has the form + * + * H(i) = I - tau * v * v' + * + * where tau is a complex scalar, and v is a complex vector with + * v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i), + * and tau in TAU(i-1). + * + * If UPLO = 'L', the matrix Q is represented as a product of elementary + * reflectors + * + * Q = H(1) H(2) . . . H(nb). + * + * Each H(i) has the form + * + * H(i) = I - tau * v * v' + * + * where tau is a complex scalar, and v is a complex vector with + * v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i), + * and tau in TAU(i). + * + * The elements of the vectors v together form the n-by-nb matrix V + * which is needed, with W, to apply the transformation to the unreduced + * part of the matrix, using a Hermitian rank-2k update of the form: + * A := A - V*W' - W*V'. + * + * The contents of A on exit are illustrated by the following examples + * with n = 5 and nb = 2: + * + * if UPLO = 'U': if UPLO = 'L': + * + * ( a a a v4 v5 ) ( d ) + * ( a a v4 v5 ) ( 1 d ) + * ( a 1 v5 ) ( v1 1 a ) + * ( d 1 ) ( v1 v2 a a ) + * ( d ) ( v1 v2 a a a ) + * + * where d denotes a diagonal element of the reduced matrix, a denotes + * an element of the original matrix that is unchanged, and vi denotes + * an element of the vector defining H(i). + * + * ===================================================================== + * + * .. Parameters .. + COMPLEX*16 ZERO, ONE, HALF + PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ), + $ ONE = ( 1.0D+0, 0.0D+0 ), + $ HALF = ( 0.5D+0, 0.0D+0 ) ) + * .. + * .. Local Scalars .. + INTEGER I, IW + COMPLEX*16 ALPHA + * .. + * .. External Subroutines .. + EXTERNAL ZAXPY, ZGEMV, ZHEMV, ZLACGV, ZLARFG, ZSCAL + * .. + * .. External Functions .. + LOGICAL LSAME + COMPLEX*16 ZDOTC + EXTERNAL LSAME, ZDOTC + * .. + * .. Intrinsic Functions .. + INTRINSIC DBLE, MIN + * .. + * .. Executable Statements .. + * + * Quick return if possible + * + IF( N.LE.0 ) + $ RETURN + * + IF( LSAME( UPLO, 'U' ) ) THEN + * + * Reduce last NB columns of upper triangle + * + DO 10 I = N, N - NB + 1, -1 + IW = I - N + NB + IF( I.LT.N ) THEN + * + * Update A(1:i,i) + * + A( I, I ) = DBLE( A( I, I ) ) + CALL ZLACGV( N-I, W( I, IW+1 ), LDW ) + CALL ZGEMV( 'No transpose', I, N-I, -ONE, A( 1, I+1 ), + $ LDA, W( I, IW+1 ), LDW, ONE, A( 1, I ), 1 ) + CALL ZLACGV( N-I, W( I, IW+1 ), LDW ) + CALL ZLACGV( N-I, A( I, I+1 ), LDA ) + CALL ZGEMV( 'No transpose', I, N-I, -ONE, W( 1, IW+1 ), + $ LDW, A( I, I+1 ), LDA, ONE, A( 1, I ), 1 ) + CALL ZLACGV( N-I, A( I, I+1 ), LDA ) + A( I, I ) = DBLE( A( I, I ) ) + END IF + IF( I.GT.1 ) THEN + * + * Generate elementary reflector H(i) to annihilate + * A(1:i-2,i) + * + ALPHA = A( I-1, I ) + CALL ZLARFG( I-1, ALPHA, A( 1, I ), 1, TAU( I-1 ) ) + E( I-1 ) = ALPHA + A( I-1, I ) = ONE + * + * Compute W(1:i-1,i) + * + CALL ZHEMV( 'Upper', I-1, ONE, A, LDA, A( 1, I ), 1, + $ ZERO, W( 1, IW ), 1 ) + IF( I.LT.N ) THEN + CALL ZGEMV( 'Conjugate transpose', I-1, N-I, ONE, + $ W( 1, IW+1 ), LDW, A( 1, I ), 1, ZERO, + $ W( I+1, IW ), 1 ) + CALL ZGEMV( 'No transpose', I-1, N-I, -ONE, + $ A( 1, I+1 ), LDA, W( I+1, IW ), 1, ONE, + $ W( 1, IW ), 1 ) + CALL ZGEMV( 'Conjugate transpose', I-1, N-I, ONE, + $ A( 1, I+1 ), LDA, A( 1, I ), 1, ZERO, + $ W( I+1, IW ), 1 ) + CALL ZGEMV( 'No transpose', I-1, N-I, -ONE, + $ W( 1, IW+1 ), LDW, W( I+1, IW ), 1, ONE, + $ W( 1, IW ), 1 ) + END IF + CALL ZSCAL( I-1, TAU( I-1 ), W( 1, IW ), 1 ) + ALPHA = -HALF*TAU( I-1 )*ZDOTC( I-1, W( 1, IW ), 1, + $ A( 1, I ), 1 ) + CALL ZAXPY( I-1, ALPHA, A( 1, I ), 1, W( 1, IW ), 1 ) + END IF + * + 10 CONTINUE + ELSE + * + * Reduce first NB columns of lower triangle + * + DO 20 I = 1, NB + * + * Update A(i:n,i) + * + A( I, I ) = DBLE( A( I, I ) ) + CALL ZLACGV( I-1, W( I, 1 ), LDW ) + CALL ZGEMV( 'No transpose', N-I+1, I-1, -ONE, A( I, 1 ), + $ LDA, W( I, 1 ), LDW, ONE, A( I, I ), 1 ) + CALL ZLACGV( I-1, W( I, 1 ), LDW ) + CALL ZLACGV( I-1, A( I, 1 ), LDA ) + CALL ZGEMV( 'No transpose', N-I+1, I-1, -ONE, W( I, 1 ), + $ LDW, A( I, 1 ), LDA, ONE, A( I, I ), 1 ) + CALL ZLACGV( I-1, A( I, 1 ), LDA ) + A( I, I ) = DBLE( A( I, I ) ) + IF( I.LT.N ) THEN + * + * Generate elementary reflector H(i) to annihilate + * A(i+2:n,i) + * + ALPHA = A( I+1, I ) + CALL ZLARFG( N-I, ALPHA, A( MIN( I+2, N ), I ), 1, + $ TAU( I ) ) + E( I ) = ALPHA + A( I+1, I ) = ONE + * + * Compute W(i+1:n,i) + * + CALL ZHEMV( 'Lower', N-I, ONE, A( I+1, I+1 ), LDA, + $ A( I+1, I ), 1, ZERO, W( I+1, I ), 1 ) + CALL ZGEMV( 'Conjugate transpose', N-I, I-1, ONE, + $ W( I+1, 1 ), LDW, A( I+1, I ), 1, ZERO, + $ W( 1, I ), 1 ) + CALL ZGEMV( 'No transpose', N-I, I-1, -ONE, A( I+1, 1 ), + $ LDA, W( 1, I ), 1, ONE, W( I+1, I ), 1 ) + CALL ZGEMV( 'Conjugate transpose', N-I, I-1, ONE, + $ A( I+1, 1 ), LDA, A( I+1, I ), 1, ZERO, + $ W( 1, I ), 1 ) + CALL ZGEMV( 'No transpose', N-I, I-1, -ONE, W( I+1, 1 ), + $ LDW, W( 1, I ), 1, ONE, W( I+1, I ), 1 ) + CALL ZSCAL( N-I, TAU( I ), W( I+1, I ), 1 ) + ALPHA = -HALF*TAU( I )*ZDOTC( N-I, W( I+1, I ), 1, + $ A( I+1, I ), 1 ) + CALL ZAXPY( N-I, ALPHA, A( I+1, I ), 1, W( I+1, I ), 1 ) + END IF + * + 20 CONTINUE + END IF + * + RETURN + * + * End of ZLATRD + * + END diff -cNr octave-2.0.7/libcruft/lapack/zsteqr.f octave-2.0.8/libcruft/lapack/zsteqr.f *** octave-2.0.7/libcruft/lapack/zsteqr.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/zsteqr.f Mon Jun 23 13:58:57 1997 *************** *** 0 **** --- 1,504 ---- + SUBROUTINE ZSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO ) + * + * -- LAPACK routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * September 30, 1994 + * + * .. Scalar Arguments .. + CHARACTER COMPZ + INTEGER INFO, LDZ, N + * .. + * .. Array Arguments .. + DOUBLE PRECISION D( * ), E( * ), WORK( * ) + COMPLEX*16 Z( LDZ, * ) + * .. + * + * Purpose + * ======= + * + * ZSTEQR computes all eigenvalues and, optionally, eigenvectors of a + * symmetric tridiagonal matrix using the implicit QL or QR method. + * The eigenvectors of a full or band complex Hermitian matrix can also + * be found if ZHETRD or ZHPTRD or ZHBTRD has been used to reduce this + * matrix to tridiagonal form. + * + * Arguments + * ========= + * + * COMPZ (input) CHARACTER*1 + * = 'N': Compute eigenvalues only. + * = 'V': Compute eigenvalues and eigenvectors of the original + * Hermitian matrix. On entry, Z must contain the + * unitary matrix used to reduce the original matrix + * to tridiagonal form. + * = 'I': Compute eigenvalues and eigenvectors of the + * tridiagonal matrix. Z is initialized to the identity + * matrix. + * + * N (input) INTEGER + * The order of the matrix. N >= 0. + * + * D (input/output) DOUBLE PRECISION array, dimension (N) + * On entry, the diagonal elements of the tridiagonal matrix. + * On exit, if INFO = 0, the eigenvalues in ascending order. + * + * E (input/output) DOUBLE PRECISION array, dimension (N-1) + * On entry, the (n-1) subdiagonal elements of the tridiagonal + * matrix. + * On exit, E has been destroyed. + * + * Z (input/output) COMPLEX*16 array, dimension (LDZ, N) + * On entry, if COMPZ = 'V', then Z contains the unitary + * matrix used in the reduction to tridiagonal form. + * On exit, if INFO = 0, then if COMPZ = 'V', Z contains the + * orthonormal eigenvectors of the original Hermitian matrix, + * and if COMPZ = 'I', Z contains the orthonormal eigenvectors + * of the symmetric tridiagonal matrix. + * If COMPZ = 'N', then Z is not referenced. + * + * LDZ (input) INTEGER + * The leading dimension of the array Z. LDZ >= 1, and if + * eigenvectors are desired, then LDZ >= max(1,N). + * + * WORK (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2)) + * If COMPZ = 'N', then WORK is not referenced. + * + * INFO (output) INTEGER + * = 0: successful exit + * < 0: if INFO = -i, the i-th argument had an illegal value + * > 0: the algorithm has failed to find all the eigenvalues in + * a total of 30*N iterations; if INFO = i, then i + * elements of E have not converged to zero; on exit, D + * and E contain the elements of a symmetric tridiagonal + * matrix which is unitarily similar to the original + * matrix. + * + * ===================================================================== + * + * .. Parameters .. + DOUBLE PRECISION ZERO, ONE, TWO, THREE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, + $ THREE = 3.0D0 ) + COMPLEX*16 CZERO, CONE + PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ), + $ CONE = ( 1.0D0, 0.0D0 ) ) + INTEGER MAXIT + PARAMETER ( MAXIT = 30 ) + * .. + * .. Local Scalars .. + INTEGER I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND, + $ LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1, + $ NM1, NMAXIT + DOUBLE PRECISION ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2, + $ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST + * .. + * .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DLAMCH, DLANST, DLAPY2 + EXTERNAL LSAME, DLAMCH, DLANST, DLAPY2 + * .. + * .. External Subroutines .. + EXTERNAL DLAE2, DLAEV2, DLARTG, DLASCL, DLASRT, XERBLA, + $ ZLASET, ZLASR, ZSWAP + * .. + * .. Intrinsic Functions .. + INTRINSIC ABS, MAX, SIGN, SQRT + * .. + * .. Executable Statements .. + * + * Test the input parameters. + * + INFO = 0 + * + IF( LSAME( COMPZ, 'N' ) ) THEN + ICOMPZ = 0 + ELSE IF( LSAME( COMPZ, 'V' ) ) THEN + ICOMPZ = 1 + ELSE IF( LSAME( COMPZ, 'I' ) ) THEN + ICOMPZ = 2 + ELSE + ICOMPZ = -1 + END IF + IF( ICOMPZ.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, + $ N ) ) ) THEN + INFO = -6 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZSTEQR', -INFO ) + RETURN + END IF + * + * Quick return if possible + * + IF( N.EQ.0 ) + $ RETURN + * + IF( N.EQ.1 ) THEN + IF( ICOMPZ.EQ.2 ) + $ Z( 1, 1 ) = CONE + RETURN + END IF + * + * Determine the unit roundoff and over/underflow thresholds. + * + EPS = DLAMCH( 'E' ) + EPS2 = EPS**2 + SAFMIN = DLAMCH( 'S' ) + SAFMAX = ONE / SAFMIN + SSFMAX = SQRT( SAFMAX ) / THREE + SSFMIN = SQRT( SAFMIN ) / EPS2 + * + * Compute the eigenvalues and eigenvectors of the tridiagonal + * matrix. + * + IF( ICOMPZ.EQ.2 ) + $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ ) + * + NMAXIT = N*MAXIT + JTOT = 0 + * + * Determine where the matrix splits and choose QL or QR iteration + * for each block, according to whether top or bottom diagonal + * element is smaller. + * + L1 = 1 + NM1 = N - 1 + * + 10 CONTINUE + IF( L1.GT.N ) + $ GO TO 160 + IF( L1.GT.1 ) + $ E( L1-1 ) = ZERO + IF( L1.LE.NM1 ) THEN + DO 20 M = L1, NM1 + TST = ABS( E( M ) ) + IF( TST.EQ.ZERO ) + $ GO TO 30 + IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+ + $ 1 ) ) ) )*EPS ) THEN + E( M ) = ZERO + GO TO 30 + END IF + 20 CONTINUE + END IF + M = N + * + 30 CONTINUE + L = L1 + LSV = L + LEND = M + LENDSV = LEND + L1 = M + 1 + IF( LEND.EQ.L ) + $ GO TO 10 + * + * Scale submatrix in rows and columns L to LEND + * + ANORM = DLANST( 'I', LEND-L+1, D( L ), E( L ) ) + ISCALE = 0 + IF( ANORM.EQ.ZERO ) + $ GO TO 10 + IF( ANORM.GT.SSFMAX ) THEN + ISCALE = 1 + CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N, + $ INFO ) + CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N, + $ INFO ) + ELSE IF( ANORM.LT.SSFMIN ) THEN + ISCALE = 2 + CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N, + $ INFO ) + CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N, + $ INFO ) + END IF + * + * Choose between QL and QR iteration + * + IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN + LEND = LSV + L = LENDSV + END IF + * + IF( LEND.GT.L ) THEN + * + * QL Iteration + * + * Look for small subdiagonal element. + * + 40 CONTINUE + IF( L.NE.LEND ) THEN + LENDM1 = LEND - 1 + DO 50 M = L, LENDM1 + TST = ABS( E( M ) )**2 + IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+ + $ SAFMIN )GO TO 60 + 50 CONTINUE + END IF + * + M = LEND + * + 60 CONTINUE + IF( M.LT.LEND ) + $ E( M ) = ZERO + P = D( L ) + IF( M.EQ.L ) + $ GO TO 80 + * + * If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 + * to compute its eigensystem. + * + IF( M.EQ.L+1 ) THEN + IF( ICOMPZ.GT.0 ) THEN + CALL DLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S ) + WORK( L ) = C + WORK( N-1+L ) = S + CALL ZLASR( 'R', 'V', 'B', N, 2, WORK( L ), + $ WORK( N-1+L ), Z( 1, L ), LDZ ) + ELSE + CALL DLAE2( D( L ), E( L ), D( L+1 ), RT1, RT2 ) + END IF + D( L ) = RT1 + D( L+1 ) = RT2 + E( L ) = ZERO + L = L + 2 + IF( L.LE.LEND ) + $ GO TO 40 + GO TO 140 + END IF + * + IF( JTOT.EQ.NMAXIT ) + $ GO TO 140 + JTOT = JTOT + 1 + * + * Form shift. + * + G = ( D( L+1 )-P ) / ( TWO*E( L ) ) + R = DLAPY2( G, ONE ) + G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) ) + * + S = ONE + C = ONE + P = ZERO + * + * Inner loop + * + MM1 = M - 1 + DO 70 I = MM1, L, -1 + F = S*E( I ) + B = C*E( I ) + CALL DLARTG( G, F, C, S, R ) + IF( I.NE.M-1 ) + $ E( I+1 ) = R + G = D( I+1 ) - P + R = ( D( I )-G )*S + TWO*C*B + P = S*R + D( I+1 ) = G + P + G = C*R - B + * + * If eigenvectors are desired, then save rotations. + * + IF( ICOMPZ.GT.0 ) THEN + WORK( I ) = C + WORK( N-1+I ) = -S + END IF + * + 70 CONTINUE + * + * If eigenvectors are desired, then apply saved rotations. + * + IF( ICOMPZ.GT.0 ) THEN + MM = M - L + 1 + CALL ZLASR( 'R', 'V', 'B', N, MM, WORK( L ), WORK( N-1+L ), + $ Z( 1, L ), LDZ ) + END IF + * + D( L ) = D( L ) - P + E( L ) = G + GO TO 40 + * + * Eigenvalue found. + * + 80 CONTINUE + D( L ) = P + * + L = L + 1 + IF( L.LE.LEND ) + $ GO TO 40 + GO TO 140 + * + ELSE + * + * QR Iteration + * + * Look for small superdiagonal element. + * + 90 CONTINUE + IF( L.NE.LEND ) THEN + LENDP1 = LEND + 1 + DO 100 M = L, LENDP1, -1 + TST = ABS( E( M-1 ) )**2 + IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+ + $ SAFMIN )GO TO 110 + 100 CONTINUE + END IF + * + M = LEND + * + 110 CONTINUE + IF( M.GT.LEND ) + $ E( M-1 ) = ZERO + P = D( L ) + IF( M.EQ.L ) + $ GO TO 130 + * + * If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 + * to compute its eigensystem. + * + IF( M.EQ.L-1 ) THEN + IF( ICOMPZ.GT.0 ) THEN + CALL DLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S ) + WORK( M ) = C + WORK( N-1+M ) = S + CALL ZLASR( 'R', 'V', 'F', N, 2, WORK( M ), + $ WORK( N-1+M ), Z( 1, L-1 ), LDZ ) + ELSE + CALL DLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 ) + END IF + D( L-1 ) = RT1 + D( L ) = RT2 + E( L-1 ) = ZERO + L = L - 2 + IF( L.GE.LEND ) + $ GO TO 90 + GO TO 140 + END IF + * + IF( JTOT.EQ.NMAXIT ) + $ GO TO 140 + JTOT = JTOT + 1 + * + * Form shift. + * + G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) ) + R = DLAPY2( G, ONE ) + G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) ) + * + S = ONE + C = ONE + P = ZERO + * + * Inner loop + * + LM1 = L - 1 + DO 120 I = M, LM1 + F = S*E( I ) + B = C*E( I ) + CALL DLARTG( G, F, C, S, R ) + IF( I.NE.M ) + $ E( I-1 ) = R + G = D( I ) - P + R = ( D( I+1 )-G )*S + TWO*C*B + P = S*R + D( I ) = G + P + G = C*R - B + * + * If eigenvectors are desired, then save rotations. + * + IF( ICOMPZ.GT.0 ) THEN + WORK( I ) = C + WORK( N-1+I ) = S + END IF + * + 120 CONTINUE + * + * If eigenvectors are desired, then apply saved rotations. + * + IF( ICOMPZ.GT.0 ) THEN + MM = L - M + 1 + CALL ZLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ), + $ Z( 1, M ), LDZ ) + END IF + * + D( L ) = D( L ) - P + E( LM1 ) = G + GO TO 90 + * + * Eigenvalue found. + * + 130 CONTINUE + D( L ) = P + * + L = L - 1 + IF( L.GE.LEND ) + $ GO TO 90 + GO TO 140 + * + END IF + * + * Undo scaling if necessary + * + 140 CONTINUE + IF( ISCALE.EQ.1 ) THEN + CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1, + $ D( LSV ), N, INFO ) + CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ), + $ N, INFO ) + ELSE IF( ISCALE.EQ.2 ) THEN + CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1, + $ D( LSV ), N, INFO ) + CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ), + $ N, INFO ) + END IF + * + * Check for no convergence to an eigenvalue after a total + * of N*MAXIT iterations. + * + IF( JTOT.EQ.NMAXIT ) THEN + DO 150 I = 1, N - 1 + IF( E( I ).NE.ZERO ) + $ INFO = INFO + 1 + 150 CONTINUE + RETURN + END IF + GO TO 10 + * + * Order eigenvalues and eigenvectors. + * + 160 CONTINUE + IF( ICOMPZ.EQ.0 ) THEN + * + * Use Quick Sort + * + CALL DLASRT( 'I', N, D, INFO ) + * + ELSE + * + * Use Selection Sort to minimize swaps of eigenvectors + * + DO 180 II = 2, N + I = II - 1 + K = I + P = D( I ) + DO 170 J = II, N + IF( D( J ).LT.P ) THEN + K = J + P = D( J ) + END IF + 170 CONTINUE + IF( K.NE.I ) THEN + D( K ) = D( I ) + D( I ) = P + CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 ) + END IF + 180 CONTINUE + END IF + RETURN + * + * End of ZSTEQR + * + END diff -cNr octave-2.0.7/libcruft/lapack/zung2l.f octave-2.0.8/libcruft/lapack/zung2l.f *** octave-2.0.7/libcruft/lapack/zung2l.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/zung2l.f Mon Jun 23 13:58:57 1997 *************** *** 0 **** --- 1,129 ---- + SUBROUTINE ZUNG2L( M, N, K, A, LDA, TAU, WORK, INFO ) + * + * -- LAPACK routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * September 30, 1994 + * + * .. Scalar Arguments .. + INTEGER INFO, K, LDA, M, N + * .. + * .. Array Arguments .. + COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * ) + * .. + * + * Purpose + * ======= + * + * ZUNG2L generates an m by n complex matrix Q with orthonormal columns, + * which is defined as the last n columns of a product of k elementary + * reflectors of order m + * + * Q = H(k) . . . H(2) H(1) + * + * as returned by ZGEQLF. + * + * Arguments + * ========= + * + * M (input) INTEGER + * The number of rows of the matrix Q. M >= 0. + * + * N (input) INTEGER + * The number of columns of the matrix Q. M >= N >= 0. + * + * K (input) INTEGER + * The number of elementary reflectors whose product defines the + * matrix Q. N >= K >= 0. + * + * A (input/output) COMPLEX*16 array, dimension (LDA,N) + * On entry, the (n-k+i)-th column must contain the vector which + * defines the elementary reflector H(i), for i = 1,2,...,k, as + * returned by ZGEQLF in the last k columns of its array + * argument A. + * On exit, the m-by-n matrix Q. + * + * LDA (input) INTEGER + * The first dimension of the array A. LDA >= max(1,M). + * + * TAU (input) COMPLEX*16 array, dimension (K) + * TAU(i) must contain the scalar factor of the elementary + * reflector H(i), as returned by ZGEQLF. + * + * WORK (workspace) COMPLEX*16 array, dimension (N) + * + * INFO (output) INTEGER + * = 0: successful exit + * < 0: if INFO = -i, the i-th argument has an illegal value + * + * ===================================================================== + * + * .. Parameters .. + COMPLEX*16 ONE, ZERO + PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ), + $ ZERO = ( 0.0D+0, 0.0D+0 ) ) + * .. + * .. Local Scalars .. + INTEGER I, II, J, L + * .. + * .. External Subroutines .. + EXTERNAL XERBLA, ZLARF, ZSCAL + * .. + * .. Intrinsic Functions .. + INTRINSIC MAX + * .. + * .. Executable Statements .. + * + * Test the input arguments + * + INFO = 0 + IF( M.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 .OR. N.GT.M ) THEN + INFO = -2 + ELSE IF( K.LT.0 .OR. K.GT.N ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -5 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZUNG2L', -INFO ) + RETURN + END IF + * + * Quick return if possible + * + IF( N.LE.0 ) + $ RETURN + * + * Initialise columns 1:n-k to columns of the unit matrix + * + DO 20 J = 1, N - K + DO 10 L = 1, M + A( L, J ) = ZERO + 10 CONTINUE + A( M-N+J, J ) = ONE + 20 CONTINUE + * + DO 40 I = 1, K + II = N - K + I + * + * Apply H(i) to A(1:m-k+i,1:n-k+i) from the left + * + A( M-N+II, II ) = ONE + CALL ZLARF( 'Left', M-N+II, II-1, A( 1, II ), 1, TAU( I ), A, + $ LDA, WORK ) + CALL ZSCAL( M-N+II-1, -TAU( I ), A( 1, II ), 1 ) + A( M-N+II, II ) = ONE - TAU( I ) + * + * Set A(m-k+i+1:m,n-k+i) to zero + * + DO 30 L = M - N + II + 1, M + A( L, II ) = ZERO + 30 CONTINUE + 40 CONTINUE + RETURN + * + * End of ZUNG2L + * + END diff -cNr octave-2.0.7/libcruft/lapack/zungql.f octave-2.0.8/libcruft/lapack/zungql.f *** octave-2.0.7/libcruft/lapack/zungql.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/zungql.f Mon Jun 23 13:58:57 1997 *************** *** 0 **** --- 1,205 ---- + SUBROUTINE ZUNGQL( M, N, K, A, LDA, TAU, WORK, LWORK, INFO ) + * + * -- LAPACK routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * September 30, 1994 + * + * .. Scalar Arguments .. + INTEGER INFO, K, LDA, LWORK, M, N + * .. + * .. Array Arguments .. + COMPLEX*16 A( LDA, * ), TAU( * ), WORK( LWORK ) + * .. + * + * Purpose + * ======= + * + * ZUNGQL generates an M-by-N complex matrix Q with orthonormal columns, + * which is defined as the last N columns of a product of K elementary + * reflectors of order M + * + * Q = H(k) . . . H(2) H(1) + * + * as returned by ZGEQLF. + * + * Arguments + * ========= + * + * M (input) INTEGER + * The number of rows of the matrix Q. M >= 0. + * + * N (input) INTEGER + * The number of columns of the matrix Q. M >= N >= 0. + * + * K (input) INTEGER + * The number of elementary reflectors whose product defines the + * matrix Q. N >= K >= 0. + * + * A (input/output) COMPLEX*16 array, dimension (LDA,N) + * On entry, the (n-k+i)-th column must contain the vector which + * defines the elementary reflector H(i), for i = 1,2,...,k, as + * returned by ZGEQLF in the last k columns of its array + * argument A. + * On exit, the M-by-N matrix Q. + * + * LDA (input) INTEGER + * The first dimension of the array A. LDA >= max(1,M). + * + * TAU (input) COMPLEX*16 array, dimension (K) + * TAU(i) must contain the scalar factor of the elementary + * reflector H(i), as returned by ZGEQLF. + * + * WORK (workspace/output) COMPLEX*16 array, dimension (LWORK) + * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. + * + * LWORK (input) INTEGER + * The dimension of the array WORK. LWORK >= max(1,N). + * For optimum performance LWORK >= N*NB, where NB is the + * optimal blocksize. + * + * INFO (output) INTEGER + * = 0: successful exit + * < 0: if INFO = -i, the i-th argument has an illegal value + * + * ===================================================================== + * + * .. Parameters .. + COMPLEX*16 ZERO + PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) + * .. + * .. Local Scalars .. + INTEGER I, IB, IINFO, IWS, J, KK, L, LDWORK, NB, NBMIN, + $ NX + * .. + * .. External Subroutines .. + EXTERNAL XERBLA, ZLARFB, ZLARFT, ZUNG2L + * .. + * .. Intrinsic Functions .. + INTRINSIC MAX, MIN + * .. + * .. External Functions .. + INTEGER ILAENV + EXTERNAL ILAENV + * .. + * .. Executable Statements .. + * + * Test the input arguments + * + INFO = 0 + IF( M.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 .OR. N.GT.M ) THEN + INFO = -2 + ELSE IF( K.LT.0 .OR. K.GT.N ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -5 + ELSE IF( LWORK.LT.MAX( 1, N ) ) THEN + INFO = -8 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZUNGQL', -INFO ) + RETURN + END IF + * + * Quick return if possible + * + IF( N.LE.0 ) THEN + WORK( 1 ) = 1 + RETURN + END IF + * + * Determine the block size. + * + NB = ILAENV( 1, 'ZUNGQL', ' ', M, N, K, -1 ) + NBMIN = 2 + NX = 0 + IWS = N + IF( NB.GT.1 .AND. NB.LT.K ) THEN + * + * Determine when to cross over from blocked to unblocked code. + * + NX = MAX( 0, ILAENV( 3, 'ZUNGQL', ' ', M, N, K, -1 ) ) + IF( NX.LT.K ) THEN + * + * Determine if workspace is large enough for blocked code. + * + LDWORK = N + IWS = LDWORK*NB + IF( LWORK.LT.IWS ) THEN + * + * Not enough workspace to use optimal NB: reduce NB and + * determine the minimum value of NB. + * + NB = LWORK / LDWORK + NBMIN = MAX( 2, ILAENV( 2, 'ZUNGQL', ' ', M, N, K, -1 ) ) + END IF + END IF + END IF + * + IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN + * + * Use blocked code after the first block. + * The last kk columns are handled by the block method. + * + KK = MIN( K, ( ( K-NX+NB-1 ) / NB )*NB ) + * + * Set A(m-kk+1:m,1:n-kk) to zero. + * + DO 20 J = 1, N - KK + DO 10 I = M - KK + 1, M + A( I, J ) = ZERO + 10 CONTINUE + 20 CONTINUE + ELSE + KK = 0 + END IF + * + * Use unblocked code for the first or only block. + * + CALL ZUNG2L( M-KK, N-KK, K-KK, A, LDA, TAU, WORK, IINFO ) + * + IF( KK.GT.0 ) THEN + * + * Use blocked code + * + DO 50 I = K - KK + 1, K, NB + IB = MIN( NB, K-I+1 ) + IF( N-K+I.GT.1 ) THEN + * + * Form the triangular factor of the block reflector + * H = H(i+ib-1) . . . H(i+1) H(i) + * + CALL ZLARFT( 'Backward', 'Columnwise', M-K+I+IB-1, IB, + $ A( 1, N-K+I ), LDA, TAU( I ), WORK, LDWORK ) + * + * Apply H to A(1:m-k+i+ib-1,1:n-k+i-1) from the left + * + CALL ZLARFB( 'Left', 'No transpose', 'Backward', + $ 'Columnwise', M-K+I+IB-1, N-K+I-1, IB, + $ A( 1, N-K+I ), LDA, WORK, LDWORK, A, LDA, + $ WORK( IB+1 ), LDWORK ) + END IF + * + * Apply H to rows 1:m-k+i+ib-1 of current block + * + CALL ZUNG2L( M-K+I+IB-1, IB, IB, A( 1, N-K+I ), LDA, + $ TAU( I ), WORK, IINFO ) + * + * Set rows m-k+i+ib:m of current block to zero + * + DO 40 J = N - K + I, N - K + I + IB - 1 + DO 30 L = M - K + I + IB, M + A( L, J ) = ZERO + 30 CONTINUE + 40 CONTINUE + 50 CONTINUE + END IF + * + WORK( 1 ) = IWS + RETURN + * + * End of ZUNGQL + * + END diff -cNr octave-2.0.7/libcruft/lapack/zungtr.f octave-2.0.8/libcruft/lapack/zungtr.f *** octave-2.0.7/libcruft/lapack/zungtr.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/lapack/zungtr.f Mon Jun 23 13:58:57 1997 *************** *** 0 **** --- 1,164 ---- + SUBROUTINE ZUNGTR( UPLO, N, A, LDA, TAU, WORK, LWORK, INFO ) + * + * -- LAPACK routine (version 2.0) -- + * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + * Courant Institute, Argonne National Lab, and Rice University + * September 30, 1994 + * + * .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, LWORK, N + * .. + * .. Array Arguments .. + COMPLEX*16 A( LDA, * ), TAU( * ), WORK( LWORK ) + * .. + * + * Purpose + * ======= + * + * ZUNGTR generates a complex unitary matrix Q which is defined as the + * product of n-1 elementary reflectors of order N, as returned by + * ZHETRD: + * + * if UPLO = 'U', Q = H(n-1) . . . H(2) H(1), + * + * if UPLO = 'L', Q = H(1) H(2) . . . H(n-1). + * + * Arguments + * ========= + * + * UPLO (input) CHARACTER*1 + * = 'U': Upper triangle of A contains elementary reflectors + * from ZHETRD; + * = 'L': Lower triangle of A contains elementary reflectors + * from ZHETRD. + * + * N (input) INTEGER + * The order of the matrix Q. N >= 0. + * + * A (input/output) COMPLEX*16 array, dimension (LDA,N) + * On entry, the vectors which define the elementary reflectors, + * as returned by ZHETRD. + * On exit, the N-by-N unitary matrix Q. + * + * LDA (input) INTEGER + * The leading dimension of the array A. LDA >= N. + * + * TAU (input) COMPLEX*16 array, dimension (N-1) + * TAU(i) must contain the scalar factor of the elementary + * reflector H(i), as returned by ZHETRD. + * + * WORK (workspace/output) COMPLEX*16 array, dimension (LWORK) + * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. + * + * LWORK (input) INTEGER + * The dimension of the array WORK. LWORK >= N-1. + * For optimum performance LWORK >= (N-1)*NB, where NB is + * the optimal blocksize. + * + * INFO (output) INTEGER + * = 0: successful exit + * < 0: if INFO = -i, the i-th argument had an illegal value + * + * ===================================================================== + * + * .. Parameters .. + COMPLEX*16 ZERO, ONE + PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ), + $ ONE = ( 1.0D+0, 0.0D+0 ) ) + * .. + * .. Local Scalars .. + LOGICAL UPPER + INTEGER I, IINFO, J + * .. + * .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME + * .. + * .. External Subroutines .. + EXTERNAL XERBLA, ZUNGQL, ZUNGQR + * .. + * .. Intrinsic Functions .. + INTRINSIC MAX + * .. + * .. Executable Statements .. + * + * Test the input arguments + * + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + ELSE IF( LWORK.LT.MAX( 1, N-1 ) ) THEN + INFO = -7 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZUNGTR', -INFO ) + RETURN + END IF + * + * Quick return if possible + * + IF( N.EQ.0 ) THEN + WORK( 1 ) = 1 + RETURN + END IF + * + IF( UPPER ) THEN + * + * Q was determined by a call to ZHETRD with UPLO = 'U' + * + * Shift the vectors which define the elementary reflectors one + * column to the left, and set the last row and column of Q to + * those of the unit matrix + * + DO 20 J = 1, N - 1 + DO 10 I = 1, J - 1 + A( I, J ) = A( I, J+1 ) + 10 CONTINUE + A( N, J ) = ZERO + 20 CONTINUE + DO 30 I = 1, N - 1 + A( I, N ) = ZERO + 30 CONTINUE + A( N, N ) = ONE + * + * Generate Q(1:n-1,1:n-1) + * + CALL ZUNGQL( N-1, N-1, N-1, A, LDA, TAU, WORK, LWORK, IINFO ) + * + ELSE + * + * Q was determined by a call to ZHETRD with UPLO = 'L'. + * + * Shift the vectors which define the elementary reflectors one + * column to the right, and set the first row and column of Q to + * those of the unit matrix + * + DO 50 J = N, 2, -1 + A( 1, J ) = ZERO + DO 40 I = J + 1, N + A( I, J ) = A( I, J-1 ) + 40 CONTINUE + 50 CONTINUE + A( 1, 1 ) = ONE + DO 60 I = 2, N + A( I, 1 ) = ZERO + 60 CONTINUE + IF( N.GT.1 ) THEN + * + * Generate Q(2:n,2:n) + * + CALL ZUNGQR( N-1, N-1, N-1, A( 2, 2 ), LDA, TAU, WORK, + $ LWORK, IINFO ) + END IF + END IF + RETURN + * + * End of ZUNGTR + * + END diff -cNr octave-2.0.7/libcruft/slatec-fn/xdgamma.f octave-2.0.8/libcruft/slatec-fn/xdgamma.f *** octave-2.0.7/libcruft/slatec-fn/xdgamma.f Wed Dec 31 18:00:00 1969 --- octave-2.0.8/libcruft/slatec-fn/xdgamma.f Fri Jun 6 21:52:26 1997 *************** *** 0 **** --- 1,5 ---- + subroutine xdgamma (x, result) + double precision x, result, dgamma + result = dgamma (x) + return + end diff -cNr octave-2.0.7/liboctave/Array2-idx.h octave-2.0.8/liboctave/Array2-idx.h *** octave-2.0.7/liboctave/Array2-idx.h Mon May 19 16:13:58 1997 --- octave-2.0.8/liboctave/Array2-idx.h Tue Jun 10 11:05:05 1997 *************** *** 118,124 **** if (idx.is_colon ()) { result_nr = nr * nc; ! result_nc = 1; } else if (idx.one_zero_only ()) { --- 118,124 ---- if (idx.is_colon ()) { result_nr = nr * nc; ! result_nc = result_nr ? 1 : 0; } else if (idx.one_zero_only ()) { diff -cNr octave-2.0.7/liboctave/ChangeLog octave-2.0.8/liboctave/ChangeLog *** octave-2.0.7/liboctave/ChangeLog Wed Jun 4 12:44:18 1997 --- octave-2.0.8/liboctave/ChangeLog Mon Jun 23 14:27:10 1997 *************** *** 1,3 **** --- 1,30 ---- + Mon Jun 23 14:26:56 1997 John W. Eaton + + * Version 2.0.8 released. + + Mon Jun 23 13:49:52 1997 John W. Eaton + + * EIG.cc (EIG::hermitian_init (const ComplexMatrix&)): New function. + (EIG::init (const ComplexMatrix&)): Call it if arg is hermitian. + (EIG::symmetric_init (const Matrix&)): New function. + (EIG::init (const Matrix&)): Call it if arg is symmetric. + (is_hermitian, is_symmetric): New static functions. + + Sun Jun 15 21:06:37 1997 John W. Eaton + + * lo-mappers.cc (acos (const Complex&)): Select branch that is + compatible with Matlab. + + Tue Jun 10 10:58:05 1997 John W. Eaton + + * Array2-idx.h: Correctly handle empty matrices indexed by a + single colon. + + Fri Jun 6 21:54:37 1997 John W. Eaton + + * lo-mappers.cc (xlgamma): Use F77_XFCN function to call dlgams. + (xgamma): Likewise, for calling xdgamma. + Wed Jun 4 12:43:30 1997 John W. Eaton * Version 2.0.7 released. diff -cNr octave-2.0.7/liboctave/EIG.cc octave-2.0.8/liboctave/EIG.cc *** octave-2.0.7/liboctave/EIG.cc Sat Mar 2 19:16:15 1996 --- octave-2.0.8/liboctave/EIG.cc Mon Jun 23 13:56:42 1997 *************** *** 1,6 **** /* ! Copyright (C) 1996 John W. Eaton This file is part of Octave. --- 1,6 ---- /* ! Copyright (C) 1996, 1997 John W. Eaton This file is part of Octave. *************** *** 29,37 **** #endif #include "EIG.h" #include "f77-fcn.h" #include "lo-error.h" - #include "mx-inlines.cc" extern "C" { --- 29,37 ---- #endif #include "EIG.h" + #include "dColVector.h" #include "f77-fcn.h" #include "lo-error.h" extern "C" { *************** *** 46,56 **** --- 46,105 ---- Complex*, const int&, Complex*, const int&, Complex*, const int&, double*, int&, long, long); + + int F77_FCN (dsyev, DSYEV) (const char*, const char*, const int&, + double*, const int&, double*, double*, + const int&, int&, long, long); + + int F77_FCN (zheev, ZHEEV) (const char*, const char*, const int&, + Complex*, const int&, double*, Complex*, + const int&, double*, int&, long, long); + } + + static bool + is_hermitian (const ComplexMatrix m) + { + int nr = m.rows (); + int nc = m.cols (); + + if (nr > 0 && nr == nc) + { + for (int i = 0; i < nr; i++) + for (int j = i; j < nc; j++) + if (m.elem (i, j) != conj (m.elem (j, i))) + return false; + + return true; + } + + return false; + } + + static bool + is_symmetric (const Matrix& m) + { + int nr = m.rows (); + int nc = m.cols (); + + if (nr > 0 && nr == nc) + { + for (int i = 0; i < nr; i++) + for (int j = i+1; j < nc; j++) + if (m.elem (i, j) != m.elem (j, i)) + return false; + + return true; + } + + return false; } int EIG::init (const Matrix& a) { + if (is_symmetric (a)) + return symmetric_init (a); + int n = a.rows (); if (n != a.cols ()) *************** *** 59,65 **** return -1; } ! int info; Matrix atmp = a; double *tmp_data = atmp.fortran_vec (); --- 108,114 ---- return -1; } ! int info = 0; Matrix atmp = a; double *tmp_data = atmp.fortran_vec (); *************** *** 73,80 **** Matrix vr (n, n); double *pvr = vr.fortran_vec (); ! int lwork = 8*n; Array work (lwork); double *pwork = work.fortran_vec (); --- 122,131 ---- Matrix vr (n, n); double *pvr = vr.fortran_vec (); ! // XXX FIXME XXX -- it might be possible to choose a better value of ! // lwork that would result in more efficient computations. + int lwork = 8*n; Array work (lwork); double *pwork = work.fortran_vec (); *************** *** 84,122 **** F77_XFCN (dgeev, DGEEV, ("N", "V", n, tmp_data, n, pwr, pwi, dummy, idummy, pvr, n, pwork, lwork, info, 1L, 1L)); ! if (f77_exception_encountered) (*current_liboctave_error_handler) ("unrecoverable error in dgeev"); else { ! lambda.resize (n); ! v.resize (n, n); ! ! for (int j = 0; j < n; j++) { ! if (wi.elem (j) == 0.0) ! { ! lambda.elem (j) = Complex (wr.elem (j)); ! for (int i = 0; i < n; i++) ! v.elem (i, j) = vr.elem (i, j); ! } ! else { ! if (j+1 >= n) { ! (*current_liboctave_error_handler) ("EIG: internal error"); ! return -1; } ! ! for (int i = 0; i < n; i++) { ! lambda.elem (j) = Complex (wr.elem (j), wi.elem (j)); ! lambda.elem (j+1) = Complex (wr.elem (j+1), wi.elem (j+1)); ! double real_part = vr.elem (i, j); ! double imag_part = vr.elem (i, j+1); ! v.elem (i, j) = Complex (real_part, imag_part); ! v.elem (i, j+1) = Complex (real_part, -imag_part); } - j++; } } } --- 135,179 ---- F77_XFCN (dgeev, DGEEV, ("N", "V", n, tmp_data, n, pwr, pwi, dummy, idummy, pvr, n, pwork, lwork, info, 1L, 1L)); ! if (f77_exception_encountered || info < 0) (*current_liboctave_error_handler) ("unrecoverable error in dgeev"); else { ! if (info > 0) ! (*current_liboctave_error_handler) ("dgeev failed to converge"); ! else { ! lambda.resize (n); ! v.resize (n, n); ! ! for (int j = 0; j < n; j++) { ! if (wi.elem (j) == 0.0) { ! lambda.elem (j) = Complex (wr.elem (j)); ! for (int i = 0; i < n; i++) ! v.elem (i, j) = vr.elem (i, j); } ! else { ! if (j+1 >= n) ! { ! (*current_liboctave_error_handler) ! ("EIG: internal error"); ! return -1; ! } ! ! for (int i = 0; i < n; i++) ! { ! lambda.elem(j) = Complex (wr.elem(j), wi.elem(j)); ! lambda.elem(j+1) = Complex (wr.elem(j+1), wi.elem(j+1)); ! double real_part = vr.elem (i, j); ! double imag_part = vr.elem (i, j+1); ! v.elem (i, j) = Complex (real_part, imag_part); ! v.elem (i, j+1) = Complex (real_part, -imag_part); ! } ! j++; } } } } *************** *** 125,131 **** } int ! EIG::init (const ComplexMatrix& a) { int n = a.rows (); --- 182,188 ---- } int ! EIG::symmetric_init (const Matrix& a) { int n = a.rows (); *************** *** 135,167 **** return -1; } ! int info; ! lambda.resize (n); ! v.resize (n, n); ! ! Complex *pw = lambda.fortran_vec (); ! Complex *pvr = v.fortran_vec (); ComplexMatrix atmp = a; Complex *tmp_data = atmp.fortran_vec (); ! int lwork = 8*n; Array work (lwork); Complex *pwork = work.fortran_vec (); ! Array rwork (4*n); double *prwork = rwork.fortran_vec (); Complex *dummy = 0; int idummy = 1; F77_XFCN (zgeev, ZGEEV, ("N", "V", n, tmp_data, n, pw, dummy, idummy, ! pvr, n, pwork, lwork, prwork, info, 1L, 1L)); ! if (f77_exception_encountered) (*current_liboctave_error_handler) ("unrecoverable error in zgeev"); return info; } --- 192,332 ---- return -1; } ! int info = 0; ! ! Matrix atmp = a; ! double *tmp_data = atmp.fortran_vec (); ! ! Array wr (n); ! double *pwr = wr.fortran_vec (); ! ! // XXX FIXME XXX -- it might be possible to choose a better value of ! // lwork that would result in more efficient computations. ! ! int lwork = 8*n; ! Array work (lwork); ! double *pwork = work.fortran_vec (); ! ! F77_XFCN (dsyev, DSYEV, ("V", "U", n, tmp_data, n, pwr, pwork, ! lwork, info, 1L, 1L)); ! ! if (f77_exception_encountered || info < 0) ! (*current_liboctave_error_handler) ("unrecoverable error in dsyev"); ! else ! { ! if (info > 0) ! (*current_liboctave_error_handler) ("dsyev failed to converge"); ! else ! { ! lambda.resize (n); ! ! for (int j = 0; j < n; j++) ! lambda.elem (j) = Complex (wr.elem (j)); ! ! v = atmp; ! } ! } ! ! return info; ! } ! ! int ! EIG::init (const ComplexMatrix& a) ! { ! if (is_hermitian (a)) ! return hermitian_init (a); ! ! int n = a.rows (); ! ! if (n != a.cols ()) ! { ! (*current_liboctave_error_handler) ("EIG requires square matrix"); ! return -1; ! } ! int info = 0; ComplexMatrix atmp = a; Complex *tmp_data = atmp.fortran_vec (); ! ComplexColumnVector w (n); ! Complex *pw = w.fortran_vec (); ! ! ComplexMatrix vtmp (n, n); ! Complex *pv = vtmp.fortran_vec (); ! ! // XXX FIXME XXX -- it might be possible to choose a better value of ! // lwork that would result in more efficient computations. + int lwork = 8*n; Array work (lwork); Complex *pwork = work.fortran_vec (); ! int lrwork = 2*n; ! Array rwork (lrwork); double *prwork = rwork.fortran_vec (); Complex *dummy = 0; int idummy = 1; F77_XFCN (zgeev, ZGEEV, ("N", "V", n, tmp_data, n, pw, dummy, idummy, ! pv, n, pwork, lwork, prwork, info, 1L, 1L)); ! if (f77_exception_encountered || info < 0) (*current_liboctave_error_handler) ("unrecoverable error in zgeev"); + else if (info > 0) + (*current_liboctave_error_handler) ("zgeev failed to converge"); + else + { + lambda = w; + v = vtmp; + } + + return info; + } + + int + EIG::hermitian_init (const ComplexMatrix& a) + { + int n = a.rows (); + + if (n != a.cols ()) + { + (*current_liboctave_error_handler) ("EIG requires square matrix"); + return -1; + } + + int info = 0; + + ComplexMatrix atmp = a; + Complex *tmp_data = atmp.fortran_vec (); + + ColumnVector w (n); + double *pw = w.fortran_vec (); + + // XXX FIXME XXX -- it might be possible to choose a better value of + // lwork that would result in more efficient computations. + + int lwork = 8*n; + Array work (lwork); + Complex *pwork = work.fortran_vec (); + + int lrwork = 3*n; + Array rwork (lrwork); + double *prwork = rwork.fortran_vec (); + + F77_XFCN (zheev, ZHEEV, ("V", "U", n, tmp_data, n, pw, pwork, + lwork, prwork, info, 1L, 1L)); + + if (f77_exception_encountered || info < 0) + (*current_liboctave_error_handler) ("unrecoverable error in zheev"); + else if (info > 0) + (*current_liboctave_error_handler) ("zheev failed to converge"); + else + { + lambda = w; + v = atmp; + } return info; } diff -cNr octave-2.0.7/liboctave/EIG.h octave-2.0.8/liboctave/EIG.h *** octave-2.0.7/liboctave/EIG.h Sat Mar 2 19:16:15 1996 --- octave-2.0.8/liboctave/EIG.h Mon Jun 23 13:49:21 1997 *************** *** 1,6 **** /* ! Copyright (C) 1996 John W. Eaton This file is part of Octave. --- 1,6 ---- /* ! Copyright (C) 1996, 1997 John W. Eaton This file is part of Octave. *************** *** 80,85 **** --- 80,88 ---- int init (const Matrix& a); int init (const ComplexMatrix& a); + + int symmetric_init (const Matrix& a); + int hermitian_init (const ComplexMatrix& a); }; #endif diff -cNr octave-2.0.7/liboctave/lo-mappers.cc octave-2.0.8/liboctave/lo-mappers.cc *** octave-2.0.7/liboctave/lo-mappers.cc Mon Jun 2 12:46:21 1997 --- octave-2.0.8/liboctave/lo-mappers.cc Fri Jun 20 10:40:14 1997 *************** *** 171,177 **** double xgamma (double x) { ! return F77_FCN (dgamma, DGAMMA) (x); } double --- 171,181 ---- double xgamma (double x) { ! double result; ! ! F77_XFCN (xdgamma, XDGAMMA, (x, result)); ! ! return result; } double *************** *** 192,198 **** double result; double sgngam; ! F77_FCN (dlgams, DLGAMS) (x, result, sgngam); return result; } --- 196,202 ---- double result; double sgngam; ! F77_XFCN (dlgams, DLGAMS, (x, result, sgngam)); return result; } *************** *** 231,238 **** acos (const Complex& x) { static Complex i (0, 1); ! Complex retval = -i * log (x + sqrt (x*x - 1.0)); ! return retval; } Complex --- 235,242 ---- acos (const Complex& x) { static Complex i (0, 1); ! ! return (real (x) * imag (x) < 0.0) ? i * acosh (x) : -i * acosh (x); } Complex diff -cNr octave-2.0.7/mkoctfile.in octave-2.0.8/mkoctfile.in *** octave-2.0.7/mkoctfile.in Mon May 19 17:20:32 1997 --- octave-2.0.8/mkoctfile.in Fri Jun 20 14:28:09 1997 *************** *** 3,56 **** # mkoctfile -- create a .oct file suitable for dynamic linking by # Octave. set -e ! if [ $# -eq 1 ]; then ! srcfile="$1" ! basnm=`echo $srcfile | sed 's,\.cc$,,'` ! objfile=$basnm.o ! octfile=$basnm.oct ! else ! echo "usage: mkoctfile file.cc" 1>&2 ! exit 1 fi ! # Configuration: these variables are filled in at configuration time. ! CPPFLAGS=%CPPFLAGS% ! INCFLAGS=%INCFLAGS% ! CXX=%CXX% ! CXX_VERSION=%CXX_VERSION% ! CXXFLAGS=%CXXFLAGS% ! CXXPICFLAG=%CXXPICFLAG% ! HOST_CXXFLAGS=%HOST_CXXFLAGS% ! NO_IMPLICIT_TEMPLATES=%NO_IMPLICIT_TEMPLATES% ! GCC_IEEE_FP_FLAG=%GCC_IEEE_FP_FLAG% ! ! LDFLAGS=%LDFLAGS% ! LIBFLAGS=%LIBFLAGS% ! RLD_FLAG=%RLD_FLAG% ! FLIBS=%FLIBS% ! LIBS=%LIBS% ! LEXLIB=%LEXLIB% ! CXXLIBS=%CXXLIBS% ! TERMLIBS=%TERMLIBS% ! LIBPLPLOT=%LIBPLPLOT% ! LIBDLFCN=%LIBDLFCN% ! ! # For now, leave -lglob out (glob/Makefile.in needs to be fixed to ! # install it. ! ! OCTAVE_LIBS="-loctinterp -loctave -ltinst -lcruft \ ! $LIBPLPLOT -lreadline -lkpathsea $LIBDLFCN" ! ! ALL_CXXFLAGS="$INCFLAGS $HOST_CXXFLAGS $NO_IMPLICIT_TEMPLATES \ ! $GCC_IEEE_FP_FLAG $CXXFLAGS" ! echo "making $objfile from $srcfile" ! $CXX -c $CPPFLAGS $CXXPICFLAG $ALL_CXXFLAGS $srcfile -o $objfile ! echo "making $octfile from $objfile" ! $SH_LD $SH_LDFLAGS -o $octfile $objfile --- 3,194 ---- # mkoctfile -- create a .oct file suitable for dynamic linking by # Octave. + # Exit immediately on any error. + set -e ! # Default values for these variables are filled in when Octave is ! # compiled. ! ! : ${CPPFLAGS=%CPPFLAGS%} ! : ${INCFLAGS=%INCFLAGS%} ! : ${F77=%F77%} ! : ${FFLAGS=%FFLAGS%} ! : ${FPICFLAG=%FPICFLAG%} ! : ${CC=%CC%} ! : ${CFLAGS=%CFLAGS%} ! : ${CPICFLAG=%CPICFLAG%} ! : ${CXX=%CXX%} ! : ${CXXFLAGS=%CXXFLAGS%} ! : ${CXXPICFLAG=%CXXPICFLAG%} ! : ${HOST_CXXFLAGS=%HOST_CXXFLAGS%} ! : ${NO_IMPLICIT_TEMPLATES=%NO_IMPLICIT_TEMPLATES%} ! : ${GCC_IEEE_FP_FLAG=%GCC_IEEE_FP_FLAG%} ! ! : ${SH_LD=%SH_LD%} ! : ${SH_LDFLAGS=%SH_LDFLAGS%} ! ! : ${ALL_FFLAGS="$FFLAGS"} ! ! : ${ALL_CFLAGS="$INCFLAGS $GCC_IEEE_FP_FLAG $CFLAGS"} ! ! : ${ALL_CXXFLAGS="$INCFLAGS $HOST_CXXFLAGS $NO_IMPLICIT_TEMPLATES $GCC_IEEE_FP_FLAG $CXXFLAGS"} ! ! # Local variables. ! ! usage_msg="usage: mkoctfile [options] file ..." ! ! cfiles= ! ccfiles= ! f77files= ! objfiles= ! octfiles= ! octfile= ! ldflags= ! dbg=: ! strip=false ! ! if [ $# -eq 0 ]; then ! echo $usage_msg ! exit 1; ! fi ! ! while [ $# -gt 0 ]; do ! file= ! case "$1" in ! *.c) ! file=$1 ! cfiles="$cfiles $file" ! ;; ! *.cc | *.C | *.cpp) ! file=$1 ! ccfiles="$ccfiles $file" ! ;; ! *.f | *.F) ! file=$1 ! f77files="$f77files $file" ! ;; ! *.o) ! file=$1 ! objfiles="$objfiles $file" ! ;; ! -d | --debug | -v | --verbose) ! dbg=echo ! ;; ! -h | -? | --help) ! echo $usage_msg ! cat << EOF ! ! Options: ! ! -h, -? --help Print this message. ! -lLIB Add library LIB to link command. ! -o FILE, --output FILE Output file name. Default extension is .oct. ! -s, --strip Strip output file. ! -v, --verbose Echo commands as they are executed. ! ! FILE Compile or link FILE. Recognized file types are: ! ! .c C source ! .cc C++ source ! .C C++ source ! .cpp C++ source ! .f Fortran source ! .F Fortran source ! .o object file ! ! EOF ! exit 0 ! ;; ! -l*) ! ldflags="$ldflags $1";; ! "") ! break ! ;; ! -o | --output) ! shift ! if [ $# -gt 0 ]; then ! octfile=`echo $1 | sed 's,\.[^.]*$,,'`.oct ! else ! echo "mkoctfile: output file name missing" ! fi ! ;; ! -s | --strip) ! strip=true ! ;; ! *) ! echo "mkoctfile: unrecognized argument $1" ! exit 1 ! ;; ! esac ! if [ -n "$file" ]; then ! if [ -z "$octfile" ]; then ! octfile=`echo $file | sed 's,\.[^.]*$,,'`.oct ! fi ! fi ! shift ! done ! ! # Compile Fortran, C, and C++ files. Add the name of each object file ! # that is produced to the overall list of object files. ! ! if [ -n "$f77files" ]; then ! for f in $f77files; do ! case $f in ! *.f) ! b=`echo $f | sed 's,\.f$,,'` ! ;; ! *.F) ! b=`echo $f | sed 's,\.F$,,'` ! ;; ! esac ! o=$b.o ! objfiles="$objfiles $o" ! $dbg $F77 -c $FPICFLAG $ALL_FFLAGS $f -o $o ! eval $F77 -c $FPICFLAG $ALL_FFLAGS $f -o $o ! done fi ! if [ -n "$cfiles" ]; then ! for f in $cfiles; do ! b=`echo $f | sed 's,\.c$,,'` ! o=$b.o ! objfiles="$objfiles $o" ! $dbg $CC -c $CPPFLAGS $CPICFLAG $ALL_CFLAGS $f -o $o ! eval $CC -c $CPPFLAGS $CPICFLAG $ALL_CFLAGS $f -o $o ! done ! fi ! ! if [ -n "$ccfiles" ]; then ! for f in $ccfiles; do ! case $f in ! *.cc) ! b=`echo $f | sed 's,\.cc$,,'` ! ;; ! *.C) ! b=`echo $f | sed 's,\.C$,,'` ! ;; ! *.cpp) ! b=`echo $f | sed 's,\.cpp$,,'` ! ;; ! esac ! o=$b.o ! objfiles="$objfiles $o" ! $dbg $CXX -c $CPPFLAGS $CXXPICFLAG $ALL_CXXFLAGS $f -o $o ! eval $CXX -c $CPPFLAGS $CXXPICFLAG $ALL_CXXFLAGS $f -o $o ! done ! fi ! # Link all the object files. ! $dbg $SH_LD $SH_LDFLAGS -o $octfile $objfiles ! eval $SH_LD $SH_LDFLAGS -o $octfile $objfiles ! # Maybe strip it. ! if $strip; then ! $dbg strip $octfile ! eval strip $octfile ! fi ! exit 0 diff -cNr octave-2.0.7/octMakefile.in octave-2.0.8/octMakefile.in *** octave-2.0.7/octMakefile.in Wed Jun 4 03:05:45 1997 --- octave-2.0.8/octMakefile.in Thu Jun 12 12:29:40 1997 *************** *** 31,37 **** MAKEINFO.PATCH ChangeLog ChangeLog.[0-9] # Complete directory trees to distribute. ! DISTDIRS = glob kpathsea make # plplot # Subdirectories in which to run `make all'. SUBDIRS = @INFO_DIR@ @PLPLOT_DIR@ @READLINE_DIR@ @DLFCN_DIR@ glob \ --- 31,37 ---- MAKEINFO.PATCH ChangeLog ChangeLog.[0-9] # Complete directory trees to distribute. ! DISTDIRS = glob kpathsea # plplot # Subdirectories in which to run `make all'. SUBDIRS = @INFO_DIR@ @PLPLOT_DIR@ @READLINE_DIR@ @DLFCN_DIR@ glob \ diff -cNr octave-2.0.7/scripts/ChangeLog octave-2.0.8/scripts/ChangeLog *** octave-2.0.7/scripts/ChangeLog Wed Jun 4 12:44:12 1997 --- octave-2.0.8/scripts/ChangeLog Mon Jun 23 14:27:08 1997 *************** *** 1,3 **** --- 1,11 ---- + Mon Jun 23 14:26:56 1997 John W. Eaton + + * Version 2.0.8 released. + + Wed Jun 18 10:24:00 1997 John W. Eaton + + * control/dlqr.m: Use ao, not a, to compute k. + Wed Jun 4 12:43:30 1997 John W. Eaton * Version 2.0.7 released. diff -cNr octave-2.0.7/scripts/control/dlqr.m octave-2.0.8/scripts/control/dlqr.m *** octave-2.0.7/scripts/control/dlqr.m Mon Jul 15 17:20:21 1996 --- octave-2.0.8/scripts/control/dlqr.m Wed Jun 18 10:24:36 1997 *************** *** 92,98 **** if (is_symmetric (q) && is_symmetric (r) ... && all (eig (q) >= 0) && all (eig (r) > 0)) p = dare (ao, b, qo, r); ! k = (r+b'*p*b)\b'*p*a + r\zz'; e = eig (a - b*k); else error ("dlqr: q (r) must be symmetric positive (semi) definite"); --- 92,98 ---- if (is_symmetric (q) && is_symmetric (r) ... && all (eig (q) >= 0) && all (eig (r) > 0)) p = dare (ao, b, qo, r); ! k = (r+b'*p*b)\b'*p*ao + r\zz'; e = eig (a - b*k); else error ("dlqr: q (r) must be symmetric positive (semi) definite"); diff -cNr octave-2.0.7/src/ChangeLog octave-2.0.8/src/ChangeLog *** octave-2.0.7/src/ChangeLog Wed Jun 4 12:44:14 1997 --- octave-2.0.8/src/ChangeLog Mon Jun 23 14:27:09 1997 *************** *** 1,3 **** --- 1,24 ---- + Mon Jun 23 14:26:56 1997 John W. Eaton + + * Version 2.0.8 released. + + Fri Jun 20 12:33:35 1997 John W. Eaton + + * toplev.cc (cmd_death_handler): New function. + (run_command_and_return_output): Insert pid of command in + octave_child_list along with pointer to cmd_death_handler so we + can get the exit status without having to block SIGCHLD. + (cleanup_iprocstream): Remove pid of command from octave_child_list. + + Sun Jun 15 16:11:13 1997 John W. Eaton + + * OPERATORS/op-cs-s.cc (ldiv): Doh, v1 is complex, v2 is real. + + Thu Jun 5 11:18:30 1997 John W. Eaton + + * dynamic-ld.cc (octave_shl_load_dynamic_loader::resolve_reference): + Call shl_findsym with type set to TYPE_UNDEFINED. + Wed Jun 4 12:43:30 1997 John W. Eaton * Version 2.0.7 released. diff -cNr octave-2.0.7/src/dynamic-ld.cc octave-2.0.8/src/dynamic-ld.cc *** octave-2.0.7/src/dynamic-ld.cc Fri May 30 16:04:02 1997 --- octave-2.0.8/src/dynamic-ld.cc Thu Jun 5 11:16:56 1997 *************** *** 124,130 **** if (handle) { ! int status = shl_findsym (&handle, nm, TYPE_PROCEDURE, &retval); if (status < 0) { --- 124,133 ---- if (handle) { ! // Don't use TYPE_PROCEDURE here. The man page says that future ! // versions of HP-UX may not support it. ! ! int status = shl_findsym (&handle, nm, TYPE_UNDEFINED, &retval); if (status < 0) { diff -cNr octave-2.0.7/src/op-cs-s.cc octave-2.0.8/src/op-cs-s.cc *** octave-2.0.7/src/op-cs-s.cc Mon Oct 14 23:34:28 1996 --- octave-2.0.8/src/op-cs-s.cc Sun Jun 15 20:54:31 1997 *************** *** 91,102 **** { CAST_BINOP_ARGS (const octave_complex&, const octave_scalar&); ! double d = v1.double_value (); if (d == 0.0) gripe_divide_by_zero (); ! return octave_value (v2.complex_value () / d); } static octave_value --- 91,102 ---- { CAST_BINOP_ARGS (const octave_complex&, const octave_scalar&); ! Complex d = v1.complex_value (); if (d == 0.0) gripe_divide_by_zero (); ! return octave_value (v2.double_value () / d); } static octave_value diff -cNr octave-2.0.7/src/sysdep.cc octave-2.0.8/src/sysdep.cc *** octave-2.0.7/src/sysdep.cc Mon Jan 27 13:20:52 1997 --- octave-2.0.8/src/sysdep.cc Mon Jun 23 09:40:10 1997 *************** *** 153,158 **** --- 153,160 ---- } #endif + // XXX FIXME XXX -- some systems define struct __exception. + #if defined (EXCEPTION_IN_MATH) extern "C" int diff -cNr octave-2.0.7/src/toplev.cc octave-2.0.8/src/toplev.cc *** octave-2.0.7/src/toplev.cc Tue May 27 09:59:25 1997 --- octave-2.0.8/src/toplev.cc Mon Jun 23 14:04:32 1997 *************** *** 654,663 **** // Execute a shell command. static void cleanup_iprocstream (void *p) { ! delete (iprocstream *) p; } static octave_value_list --- 654,675 ---- // Execute a shell command. + static int cmd_status = 0; + + static void + cmd_death_handler (pid_t, int status) + { + cmd_status = status; + } + static void cleanup_iprocstream (void *p) { ! iprocstream *cmd = (iprocstream *) p; ! ! octave_child_list::remove (cmd->pid ()); ! ! delete cmd; } static octave_value_list *************** *** 667,706 **** iprocstream *cmd = new iprocstream (cmd_str.c_str ()); ! add_unwind_protect (cleanup_iprocstream, cmd); ! int status = 127; ! ! if (cmd && *cmd) { ! ostrstream output_buf; ! char ch; ! while (cmd->get (ch)) ! output_buf.put (ch); ! status = cmd->close (); ! // The value in status is as returned by waitpid. If the ! // process exited normally, extract the actual exit status of ! // the command. Otherwise, return 127 as a failure code. ! if (WIFEXITED (status)) ! status = WEXITSTATUS (status); ! output_buf << ends; ! char *msg = output_buf.str (); ! retval(1) = (double) status; ! retval(0) = msg; ! delete [] msg; } else error ("unable to start subprocess for `%s'", cmd_str.c_str ()); - - run_unwind_protect (); return retval; } --- 679,730 ---- iprocstream *cmd = new iprocstream (cmd_str.c_str ()); ! cmd_status = -1; ! if (cmd) { ! octave_child_list::insert (cmd->pid (), cmd_death_handler); ! ! add_unwind_protect (cleanup_iprocstream, cmd); ! ! if (*cmd) ! { ! ostrstream output_buf; ! char ch; ! while (cmd->get (ch)) ! output_buf.put (ch); ! cmd->close (); ! // One way or another, cmd_death_handler should be called ! // when the process exits, and it will save the exit status ! // of the command in cmd_status. ! // The value in cmd_status is as returned by waitpid. If ! // the process exited normally, extract the actual exit ! // status of the command. Otherwise, return 127 as a ! // failure code. ! if (WIFEXITED (cmd_status)) ! cmd_status = WEXITSTATUS (cmd_status); ! else ! cmd_status = 127; ! output_buf << ends; ! char *msg = output_buf.str (); ! retval(1) = (double) cmd_status; ! retval(0) = msg; ! ! delete [] msg; ! } ! ! run_unwind_protect (); } else error ("unable to start subprocess for `%s'", cmd_str.c_str ()); return retval; } diff -cNr octave-2.0.7/src/version.h octave-2.0.8/src/version.h *** octave-2.0.7/src/version.h Mon Jun 2 13:58:42 1997 --- octave-2.0.8/src/version.h Thu Jun 12 12:32:18 1997 *************** *** 20,26 **** */ ! #define OCTAVE_VERSION "2.0.7" #define OCTAVE_COPYRIGHT \ "Copyright (C) 1996, 1997 John W. Eaton." --- 20,26 ---- */ ! #define OCTAVE_VERSION "2.0.8" #define OCTAVE_COPYRIGHT \ "Copyright (C) 1996, 1997 John W. Eaton." diff -cNr octave-2.0.7/test/ChangeLog octave-2.0.8/test/ChangeLog *** octave-2.0.7/test/ChangeLog Wed Jun 4 12:43:48 1997 --- octave-2.0.8/test/ChangeLog Mon Jun 23 14:27:05 1997 *************** *** 1,3 **** --- 1,7 ---- + Mon Jun 23 14:26:56 1997 John W. Eaton + + * Version 2.0.8 released. + Wed Jun 4 12:43:30 1997 John W. Eaton * Version 2.0.7 released. diff -cNr octave-2.0.7/test/octave.test/linalg/eig-1.m octave-2.0.8/test/octave.test/linalg/eig-1.m *** octave-2.0.7/test/octave.test/linalg/eig-1.m Mon Feb 24 15:56:44 1997 --- octave-2.0.8/test/octave.test/linalg/eig-1.m Mon Jun 23 14:18:02 1997 *************** *** 1 **** ! all (eig ([1, 2; 2, 1]) - [3; -1] < sqrt (eps)) --- 1 ---- ! all (eig ([1, 2; 2, 1]) - [-1; 3] < sqrt (eps)) diff -cNr octave-2.0.7/test/octave.test/linalg/eig-2.m octave-2.0.8/test/octave.test/linalg/eig-2.m *** octave-2.0.7/test/octave.test/linalg/eig-2.m Mon Feb 24 15:50:53 1997 --- octave-2.0.8/test/octave.test/linalg/eig-2.m Mon Jun 23 14:18:02 1997 *************** *** 1,4 **** [v, d] = eig ([1, 2; 2, 1]); x = 1 / sqrt (2); ! (abs (d - [3, 0; 0, -1] < sqrt (eps)) ! && abs (v - [x, -x; x, x] < sqrt (eps))) --- 1,4 ---- [v, d] = eig ([1, 2; 2, 1]); x = 1 / sqrt (2); ! (abs (d - [-1, 0; 0, 3] < sqrt (eps)) ! && abs (v - [-x, x; x, x] < sqrt (eps))) PATCH_EOF