diff --git a/.gitignore b/.gitignore index 848a7c3..03dfd16 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ /superlu_4.3.tar.gz /superlu_5.1.1.tar.gz /superlu_5.2.0.tar.gz +/superlu_5.2.1.tar.gz diff --git a/SuperLU-add-fpic.patch b/SuperLU-add-fpic.patch deleted file mode 100644 index 906e379..0000000 --- a/SuperLU-add-fpic.patch +++ /dev/null @@ -1,13 +0,0 @@ -diff -up SuperLU_4.3/MAKE_INC/make.linux.fpic SuperLU_4.3/MAKE_INC/make.linux ---- SuperLU_4.3/MAKE_INC/make.linux.fpic 2012-01-06 10:39:58.474356562 +0530 -+++ SuperLU_4.3/MAKE_INC/make.linux 2012-01-06 10:40:59.933356590 +0530 -@@ -46,7 +46,8 @@ ARCHFLAGS = cr - RANLIB = ranlib - - CC = gcc --CFLAGS = -O3 -+FPIC = -fPIC -+CFLAGS = -O3 $(FPIC) - NOOPTS = - FORTRAN = g77 - FFLAGS = -O2 diff --git a/SuperLU-build-shared-lib3.patch b/SuperLU-build-shared-lib3.patch deleted file mode 100644 index afaefe9..0000000 --- a/SuperLU-build-shared-lib3.patch +++ /dev/null @@ -1,39 +0,0 @@ -diff -up SuperLU_4.3/SRC/Makefile.fix SuperLU_4.3/SRC/Makefile ---- SuperLU_4.3/SRC/Makefile.fix 2012-02-01 17:15:17.711698876 +0530 -+++ SuperLU_4.3/SRC/Makefile 2012-02-01 17:17:07.794698927 +0530 -@@ -111,7 +111,7 @@ ZLUSRC = \ - ilu_zcolumn_dfs.o ilu_zpanel_dfs.o ilu_zcopy_to_ucol.o \ - ilu_zpivotL.o zdiagonal.o - --all: single double complex complex16 -+all: sharedlib - - single: $(SLUSRC) $(ALLAUX) $(LAAUX) $(SLASRC) $(SCLAUX) - $(ARCH) $(ARCHFLAGS) $(SUPERLULIB) \ -@@ -133,17 +133,22 @@ complex16: $(ZLUSRC) $(ALLAUX) $(LAAUX) - $(ZLUSRC) $(ALLAUX) $(LAAUX) $(ZLASRC) $(DZLAUX) - $(RANLIB) $(SUPERLULIB) - -+sharedlib: $(ALLAUX) $(LAAUX) $(SLASRC) $(DLASRC) $(CLASRC) $(ZLASRC) $(SCLAUX) $(DZLAUX) $(SLUSRC) $(DLUSRC) $(CLUSRC) $(ZLUSRC) -+ $(CC) $(CFLAGS) $(LIBS) -shared -Wl,-soname,libsuperlu.so.4.3 -o libsuperlu.so.4.3 \ -+ $(ALLAUX) $(LAAUX) $(SLASRC) $(DLASRC) $(CLASRC) $(ZLASRC) $(SCLAUX) \ -+ $(DZLAUX) $(SLUSRC) $(DLUSRC) $(CLUSRC) $(ZLUSRC) -+ ln -sf libsuperlu.so.4.3 libsuperlu.so - - ################################## - # Do not optimize these routines # - ################################## --slamch.o: slamch.c ; $(CC) -c $(NOOPTS) $(CDEFS) $< --dlamch.o: dlamch.c ; $(CC) -c $(NOOPTS) $(CDEFS) $< --superlu_timer.o: superlu_timer.c ; $(CC) -c $(NOOPTS) $< -+slamch.o: slamch.c ; $(CC) $(FPIC) $(LIBS) -c $(NOOPTS) $(CDEFS) $< -+dlamch.o: dlamch.c ; $(CC) $(FPIC) $(LIBS) -c $(NOOPTS) $(CDEFS) $< -+superlu_timer.o: superlu_timer.c ; $(CC) $(FPIC) $(LIBS) -c $(NOOPTS) $< - ################################## - - .c.o: -- $(CC) $(CFLAGS) $(CDEFS) $(BLASDEF) -c $< $(VERBOSE) -+ $(CC) $(CFLAGS) $(CDEFS) $(BLASDEF) $(LIBS) -c $< $(VERBOSE) - - .f.o: - $(FORTRAN) $(FFLAGS) -c $< diff --git a/SuperLU-fix-format-security.patch b/SuperLU-fix-format-security.patch deleted file mode 100644 index 5334a08..0000000 --- a/SuperLU-fix-format-security.patch +++ /dev/null @@ -1,13 +0,0 @@ -Index: SuperLU_4.3/SRC/util.c -=================================================================== ---- SuperLU_4.3.orig/SRC/util.c -+++ SuperLU_4.3/SRC/util.c -@@ -29,7 +29,7 @@ - - void superlu_abort_and_exit(char* msg) - { -- fprintf(stderr, msg); -+ fputs(stderr, msg); - exit (-1); - } - diff --git a/SuperLU-fix-testsuite.patch b/SuperLU-fix-testsuite.patch deleted file mode 100644 index e358875..0000000 --- a/SuperLU-fix-testsuite.patch +++ /dev/null @@ -1,52 +0,0 @@ -Index: SuperLU_4.3/TESTING/Makefile -=================================================================== ---- SuperLU_4.3.orig/TESTING/Makefile -+++ SuperLU_4.3/TESTING/Makefile -@@ -52,9 +52,9 @@ testmat: - - single: ./stest stest.out - --./stest: $(SLINTST) $(ALINTST) $(SUPERLULIB) $(TMGLIB) -+./stest: $(SLINTST) $(ALINTST) $(TMGLIB) - $(LOADER) $(LOADOPTS) $(SLINTST) $(ALINTST) \ -- $(TMGLIB) $(SUPERLULIB) $(BLASLIB) -lm -o $@ -+ $(TMGLIB) -Wl,-rpath,../SRC $(SUPERLULIB) $(BLASLIB) -lm -o $@ - - stest.out: stest stest.csh - @echo Testing SINGLE PRECISION linear equation routines -@@ -62,9 +62,9 @@ stest.out: stest stest.csh - - double: ./dtest dtest.out - --./dtest: $(DLINTST) $(ALINTST) $(SUPERLULIB) $(TMGLIB) -+./dtest: $(DLINTST) $(ALINTST) $(TMGLIB) - $(LOADER) $(LOADOPTS) $(DLINTST) $(ALINTST) \ -- $(TMGLIB) $(SUPERLULIB) $(BLASLIB) -lm -o $@ -+ $(TMGLIB) -Wl,-rpath,../SRC $(SUPERLULIB) $(BLASLIB) -lm -o $@ - - dtest.out: dtest dtest.csh - @echo Testing DOUBLE PRECISION linear equation routines -@@ -72,9 +72,9 @@ dtest.out: dtest dtest.csh - - complex: ./ctest ctest.out - --./ctest: $(CLINTST) $(ALINTST) $(SUPERLULIB) $(TMGLIB) -+./ctest: $(CLINTST) $(ALINTST) $(TMGLIB) - $(LOADER) $(LOADOPTS) $(CLINTST) $(ALINTST) \ -- $(TMGLIB) $(SUPERLULIB) $(BLASLIB) -lm -o $@ -+ $(TMGLIB) -Wl,-rpath,../SRC $(SUPERLULIB) $(BLASLIB) -lm -o $@ - - ctest.out: ctest ctest.csh - @echo Testing SINGLE COMPLEX linear equation routines -@@ -82,9 +82,9 @@ ctest.out: ctest ctest.csh - - complex16: ./ztest ztest.out - --./ztest: $(ZLINTST) $(ALINTST) $(SUPERLULIB) $(TMGLIB) -+./ztest: $(ZLINTST) $(ALINTST) $(TMGLIB) - $(LOADER) $(LOADOPTS) $(ZLINTST) $(ALINTST) \ -- $(TMGLIB) $(SUPERLULIB) $(BLASLIB) -lm -o $@ -+ $(TMGLIB) -Wl,-rpath,../SRC $(SUPERLULIB) $(BLASLIB) -lm -o $@ - - ztest.out: ztest ztest.csh - @echo Testing DOUBLE COMPLEX linear equation routines diff --git a/SuperLU.spec b/SuperLU.spec index 8785a00..9fa1efb 100644 --- a/SuperLU.spec +++ b/SuperLU.spec @@ -1,5 +1,5 @@ %global genname superlu -%global libver 5.1 +%global libver 5 ## The RPM macro for the linker flags does not exist on EPEL %if 0%{?rhel} < 7 @@ -7,24 +7,18 @@ %endif Name: SuperLU -Version: 5.2.0 -Release: 8%{?dist} +Version: 5.2.1 +Release: 1%{?dist} Summary: Subroutines to solve sparse linear systems License: BSD and GPLV2+ URL: http://crd-legacy.lbl.gov/~xiaoye/SuperLU/ Source0: http://crd-legacy.lbl.gov/~xiaoye/SuperLU/%{genname}_%{version}.tar.gz -# Build with -fPIC -Patch0: %{genname}-5x-add-fpic.patch -# Build shared library -Patch1: %{genname}-5x-build-shared-lib3.patch -# Fixes testsuite -Patch3: %{genname}-5x-fix-testsuite.patch -# remove non-free mc64 functionality -# patch obtained from the debian package -Patch4: %{genname}-removemc64.patch - -BuildRequires: atlas-devel, gcc +Patch0: %{genname}-cmake-includedir.patch +Patch1: %{genname}-removemc64.patch + +BuildRequires: atlas-devel, gcc, cmake +BuildRequires: gcc-gfortran BuildRequires: csh %description @@ -50,11 +44,7 @@ The %{name}-doc package contains all the help documentation along with C and FORTRAN examples. %prep -%setup -q -n %{name}_%{version} -%patch0 -p1 -%patch1 -p1 -%patch3 -p1 -%patch4 +%autosetup -n %{name}_%{version} -p1 rm -fr SRC/mc64ad.f.bak find . -type f | sed -e "/TESTING/d" | xargs chmod a-x @@ -63,46 +53,23 @@ find EXAMPLE -type f | while read file do [ "$(file $file | awk '{print $2}')" = ELF ] && rm $file || : done -cp -p MAKE_INC/make.linux make.inc -sed -i -e "s|-O3|$RPM_OPT_FLAGS|" \ - -e "s|\$(SUPERLULIB) ||" \ - -e "s|\$(HOME)/Dropbox/Codes/%{name}/%{name}|%{_builddir}/%{name}_%{version}|" \ - -e 's!lib/libsuperlu_5.1.a$!SRC/libsuperlu.so!' \ - -e 's!-shared!& %{__global_ldflags}!' \ -%if 0%{?rhel} && 0%{?rhel} < 7 - -e "s|-L/usr/lib -lblas|-L%{_libdir}/atlas -lf77blas|" \ -%else - -e "s|-L/usr/lib -lblas|-L%{_libdir}/atlas -lsatlas|" \ -%endif - make.inc %build -make %{?_smp_mflags} superlulib -make -C TESTING +%cmake . +%make_build %install -mkdir -p %{buildroot}%{_libdir} -mkdir -p %{buildroot}%{_includedir}/%{name} -install -p SRC/libsuperlu.so.%{libver} %{buildroot}%{_libdir} -install -p SRC/*.h %{buildroot}%{_includedir}/%{name} -chmod -x %{buildroot}%{_includedir}/%{name}/*.h -cp -Pp SRC/libsuperlu.so %{buildroot}%{_libdir} +%make_install %check -pushd TESTING -for _test in c d s z -do - chmod +x ${_test}test.csh - ./${_test}test.csh -done -popd +ctest -V %{?_smp_mflags} %ldconfig_scriptlets %files %doc README %license License.txt -%{_libdir}/libsuperlu.so.%{libver} +%{_libdir}/libsuperlu.so.%{libver}* %files devel %{_includedir}/%{name}/ @@ -112,6 +79,12 @@ popd %doc DOC EXAMPLE FORTRAN %changelog +* Fri Apr 13 2018 Rafael dos Santos - 5.2.1-1 +- Update to 5.2.1 +- Use cmake build system +- Drop obsolete patches +- Resolves #1547494 - build with standard Fedora flags + * Wed Feb 21 2018 Antonio Trande - 5.2.0-8 - Add gcc BR - Remove el5 bits diff --git a/sources b/sources index f4c2ffb..835c930 100644 --- a/sources +++ b/sources @@ -1 +1 @@ -ba3cfc93a93a8caed90bcf6f9b804ac7 superlu_5.2.0.tar.gz +SHA512 (superlu_5.2.1.tar.gz) = 30538b4c2809294b8f34646bce6445944f21a1dffaf3ec0a0f29a55d5261caa56e4279d7722bb95cc9d89450d36ded969617edc82ecce7d0f1dfb24040d80d07 diff --git a/superlu-5x-add-fpic.patch b/superlu-5x-add-fpic.patch deleted file mode 100644 index ab1e0cc..0000000 --- a/superlu-5x-add-fpic.patch +++ /dev/null @@ -1,12 +0,0 @@ ---- SuperLU_5.1.1/MAKE_INC/make.linux.orig 2016-03-26 08:43:03.038994570 -0400 -+++ SuperLU_5.1.1/MAKE_INC/make.linux 2016-03-26 08:44:09.141092611 -0400 -@@ -46,7 +46,8 @@ - RANLIB = ranlib - - CC = gcc --CFLAGS = -O3 -g -+FPIC = -fPIC -+CFLAGS = -O3 $(FPIC) - NOOPTS = - FORTRAN = gfortran #g77 - FFLAGS = -O2 -g -fopenmp diff --git a/superlu-5x-build-shared-lib3.patch b/superlu-5x-build-shared-lib3.patch deleted file mode 100644 index 88e7aa3..0000000 --- a/superlu-5x-build-shared-lib3.patch +++ /dev/null @@ -1,39 +0,0 @@ ---- SuperLU_5.1.1/SRC/Makefile.orig 2016-03-26 08:47:18.885503561 -0400 -+++ SuperLU_5.1.1/SRC/Makefile 2016-03-26 08:52:10.512497969 -0400 -@@ -98,7 +98,7 @@ - ilu_zcolumn_dfs.o ilu_zpanel_dfs.o ilu_zcopy_to_ucol.o \ - ilu_zpivotL.o zdiagonal.o zlacon2.o dzsum1.o izmax1.o - --all: single double complex complex16 -+all: sharedlib - - single: $(SLUSRC) $(ALLAUX) $(SCAUX) - $(ARCH) $(ARCHFLAGS) $(SUPERLULIB) $(SLUSRC) $(ALLAUX) $(SCAUX) -@@ -116,17 +116,23 @@ - $(ARCH) $(ARCHFLAGS) $(SUPERLULIB) $(ZLUSRC) $(ALLAUX) $(DZLAUX) - $(RANLIB) $(SUPERLULIB) - -+sharedlib: $(ALLAUX) $(LAAUX) $(SLASRC) $(DLASRC) $(CLASRC) $(ZLASRC) $(SCLAUX) $(DZLAUX) $(SLUSRC) $(DLUSRC) $(CLUSRC) $(ZLUSRC) -+ $(CC) $(CFLAGS) $(LIBS) -shared -Wl,-soname,libsuperlu.so.5.1 -o libsuperlu.so.5.1 \ -+ $(ALLAUX) $(LAAUX) $(SLASRC) $(DLASRC) $(CLASRC) $(ZLASRC) $(SCLAUX) \ -+ $(DZLAUX) $(SLUSRC) $(DLUSRC) $(CLUSRC) $(ZLUSRC) -+ ln -sf libsuperlu.so.5.1 libsuperlu.so -+ - - ################################## - # Do not optimize these routines # - ################################## --smach.o: smach.c ; $(CC) -c $(NOOPTS) $(CDEFS) $< --dmach.o: dmach.c ; $(CC) -c $(NOOPTS) $(CDEFS) $< --superlu_timer.o: superlu_timer.c ; $(CC) -c $(NOOPTS) $< -+slamch.o: slamch.c ; $(CC) $(FPIC) $(LIBS) -c $(NOOPTS) $(CDEFS) $< -+dlamch.o: dlamch.c ; $(CC) $(FPIC) $(LIBS) -c $(NOOPTS) $(CDEFS) $< -+superlu_timer.o: superlu_timer.c ; $(CC) $(FPIC) $(LIBS) -c $(NOOPTS) $< - ################################## - - .c.o: -- $(CC) $(CFLAGS) $(CDEFS) $(BLASDEF) -c $< $(VERBOSE) -+ $(CC) $(CFLAGS) $(CDEFS) $(BLASDEF) $(LIBS) -c $< $(VERBOSE) - - .f.o: - $(FORTRAN) $(FFLAGS) -c $< diff --git a/superlu-5x-fix-testsuite.patch b/superlu-5x-fix-testsuite.patch deleted file mode 100644 index 798690b..0000000 --- a/superlu-5x-fix-testsuite.patch +++ /dev/null @@ -1,50 +0,0 @@ ---- SuperLU_5.1.1/TESTING/Makefile.orig 2016-03-26 08:56:10.347188385 -0400 -+++ SuperLU_5.1.1/TESTING/Makefile 2016-03-26 08:58:39.763126530 -0400 -@@ -52,9 +52,9 @@ - - single: ./stest stest.out - --./stest: $(SLINTST) $(ALINTST) $(SUPERLULIB) $(TMGLIB) -+./stest: $(SLINTST) $(ALINTST) $(TMGLIB) - $(LOADER) $(LOADOPTS) $(SLINTST) $(ALINTST) \ -- $(TMGLIB) $(SUPERLULIB) $(BLASLIB) -lm -o $@ -+ $(TMGLIB) -Wl,-rpath,../SRC $(SUPERLULIB) $(BLASLIB) -lm -o $@ - - stest.out: stest stest.csh - @echo Testing SINGLE PRECISION linear equation routines -@@ -62,9 +62,9 @@ - - double: ./dtest dtest.out - --./dtest: $(DLINTST) $(ALINTST) $(SUPERLULIB) $(TMGLIB) -+./dtest: $(DLINTST) $(ALINTST) $(TMGLIB) - $(LOADER) $(LOADOPTS) $(DLINTST) $(ALINTST) \ -- $(TMGLIB) $(SUPERLULIB) $(BLASLIB) -lm -o $@ -+ $(TMGLIB) -Wl,-rpath,../SRC $(SUPERLULIB) $(BLASLIB) -lm -o $@ - - dtest.out: dtest dtest.csh - @echo Testing DOUBLE PRECISION linear equation routines -@@ -72,9 +72,9 @@ - - complex: ./ctest ctest.out - --./ctest: $(CLINTST) $(ALINTST) $(SUPERLULIB) $(TMGLIB) -+./ctest: $(CLINTST) $(ALINTST) $(TMGLIB) - $(LOADER) $(LOADOPTS) $(CLINTST) $(ALINTST) \ -- $(TMGLIB) $(SUPERLULIB) $(BLASLIB) -lm -o $@ -+ $(TMGLIB) -Wl,-rpath,../SRC $(SUPERLULIB) $(BLASLIB) -lm -o $@ - - ctest.out: ctest ctest.csh - @echo Testing SINGLE COMPLEX linear equation routines -@@ -82,9 +82,9 @@ - - complex16: ./ztest ztest.out - --./ztest: $(ZLINTST) $(ALINTST) $(SUPERLULIB) $(TMGLIB) -+./ztest: $(ZLINTST) $(ALINTST) $(TMGLIB) - $(LOADER) $(LOADOPTS) $(ZLINTST) $(ALINTST) \ -- $(TMGLIB) $(SUPERLULIB) $(BLASLIB) -lm -o $@ -+ $(TMGLIB) -Wl,-rpath,../SRC $(SUPERLULIB) $(BLASLIB) -lm -o $@ - - ztest.out: ztest ztest.csh - @echo Testing DOUBLE COMPLEX linear equation routines diff --git a/superlu-cmake-includedir.patch b/superlu-cmake-includedir.patch new file mode 100644 index 0000000..c955879 --- /dev/null +++ b/superlu-cmake-includedir.patch @@ -0,0 +1,9 @@ +--- SuperLU_5.2.1/SRC/CMakeLists.txt 2018-04-13 17:01:44.007761121 +0200 ++++ CMakeLists.txt 2018-04-13 17:01:30.573718261 +0200 +@@ -242,5 +242,5 @@ + + install(FILES ${headers} + # DESTINATION ${CMAKE_INSTALL_PREFIX}/include +- DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} ++ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/SuperLU + ) diff --git a/superlu-cmake-setup.patch b/superlu-cmake-setup.patch new file mode 100644 index 0000000..37066e0 --- /dev/null +++ b/superlu-cmake-setup.patch @@ -0,0 +1,14 @@ +--- /tmp/SuperLU_5.2.0/SRC/CMakeLists.txt 2016-04-09 03:53:17.000000000 +0200 ++++ SuperLU_5.2.0/SRC/CMakeLists.txt 2018-04-13 11:42:28.481583088 +0200 +@@ -233,9 +233,9 @@ + VERSION ${PROJECT_VERSION} SOVERSION ${VERSION_MAJOR} + ) + install(TARGETS superlu +- DESTINATION ${CMAKE_INSTALL_PREFIX}/lib ++ DESTINATION ${CMAKE_INSTALL_LIBDIR} + ) + + install(FILES ${headers} +- DESTINATION ${CMAKE_INSTALL_PREFIX}/include ++ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} + ) diff --git a/superlu-removemc64.patch b/superlu-removemc64.patch index 00d36af..58596b8 100644 --- a/superlu-removemc64.patch +++ b/superlu-removemc64.patch @@ -1,17 +1,2654 @@ -Description: Provide stubs for the two functions whose implementation is nonfree -Author: Sébastien Villemot -Forwarded: not-needed -Last-Update: 2013-11-27 ---- -This patch header follows DEP-3: http://dep.debian.net/deps/dep3/ ---- /dev/null -+++ b/SRC/mc64ad.c -@@ -0,0 +1,16 @@ +--- SuperLU_5.2.1/SRC/mc64ad.c 2016-05-22 17:58:44.000000000 +0200 ++++ mc64ad.c 2018-04-13 17:13:23.571981656 +0200 +@@ -1,2645 +1,16 @@ +-/* mc64ad.f -- translated by f2c (version 20100827). +- You must link the resulting object file with libf2c: +- on Microsoft Windows system, link with libf2c.lib; +- on Linux or Unix systems, link with .../path/to/libf2c.a -lm +- or, if you install libf2c.a in a standard place, with -lf2c -lm +- -- in that order, at the end of the command line, as in +- cc *.o -lf2c -lm +- Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., +#include +#include -+ + +- http://www.netlib.org/f2c/libf2c.zip +-*/ +- +-#include "slu_ddefs.h" +- +-#define abs(a) ((a) >= 0) ? (a) : -(a) +-#define min(a,b) ((a) < (b)) ? (a) : (b) +- +-#if 0 +-/* Table of constant values */ +-static int_t c__1 = 1; +-static int_t c__2 = 2; +-#endif +- +-/* CCCC COPYRIGHT (c) 1999 Council for the Central Laboratory of the */ +-/* CCCC Research Councils. All rights reserved. */ +-/* CCCC PACKAGE MC64A/AD */ +-/* CCCC AUTHORS Iain Duff (i.duff@rl.ac.uk) and Jacko Koster (jak@ii.uib.no) */ +-/* CCCC LAST UPDATE 20/09/99 */ +-/* CCCC */ +-/* *** Conditions on external use *** */ +- +-/* The user shall acknowledge the contribution of this */ +-/* package in any publication of material dependent upon the use of */ +-/* the package. The user shall use reasonable endeavours to notify */ +-/* the authors of the package of this publication. */ +- +-/* The user can modify this code but, at no time */ +-/* shall the right or title to all or any part of this package pass */ +-/* to the user. The user shall make available free of charge */ +-/* to the authors for any purpose all information relating to any */ +-/* alteration or addition made to this package for the purposes of */ +-/* extending the capabilities or enhancing the performance of this */ +-/* package. */ +- +-/* The user shall not pass this code directly to a third party without the */ +-/* express prior consent of the authors. Users wanting to licence their */ +-/* own copy of these routines should send email to hsl@stfc.ac.uk */ +- +-/* None of the comments from the Copyright notice up to and including this */ +-/* one shall be removed or altered in any way. */ +-/* ********************************************************************** */ +-/* Subroutine */ int_t mc64id_(int_t *icntl) +-{ +- int_t i__; +- +- +-/* *** Copyright (c) 1999 Council for the Central Laboratory of the */ +-/* Research Councils *** */ +-/* *** Although every effort has been made to ensure robustness and *** */ +-/* *** reliability of the subroutines in this MC64 suite, we *** */ +-/* *** disclaim any liability arising through the use or misuse of *** */ +-/* *** any of the subroutines. *** */ +-/* *** Any problems? Contact ... */ +-/* Iain Duff (I.Duff@rl.ac.uk) or Jacko Koster (jak@ii.uib.no) *** */ +- +-/* Purpose */ +-/* ======= */ +- +-/* The components of the array ICNTL control the action of MC64A/AD. */ +-/* Default values for these are set in this subroutine. */ +- +-/* Parameters */ +-/* ========== */ +- +- +-/* Local variables */ +- +-/* ICNTL(1) has default value 6. */ +-/* It is the output stream for error messages. If it */ +-/* is negative, these messages will be suppressed. */ +- +-/* ICNTL(2) has default value 6. */ +-/* It is the output stream for warning messages. */ +-/* If it is negative, these messages are suppressed. */ +- +-/* ICNTL(3) has default value -1. */ +-/* It is the output stream for monitoring printing. */ +-/* If it is negative, these messages are suppressed. */ +- +-/* ICNTL(4) has default value 0. */ +-/* If left at the defaut value, the incoming data is checked for */ +-/* out-of-range indices and duplicates. Setting ICNTL(4) to any */ +-/* other will avoid the checks but is likely to cause problems */ +-/* later if out-of-range indices or duplicates are present. */ +-/* The user should only set ICNTL(4) non-zero, if the data is */ +-/* known to avoid these problems. */ +- +-/* ICNTL(5) to ICNTL(10) are not used by MC64A/AD but are set to */ +-/* zero in this routine. */ +-/* Initialization of the ICNTL array. */ +- /* Parameter adjustments */ +- --icntl; +- +- /* Function Body */ +- icntl[1] = 6; +- icntl[2] = 6; +- icntl[3] = -1; +- for (i__ = 4; i__ <= 10; ++i__) { +- icntl[i__] = 0; +-/* L10: */ +- } +- return 0; +-} /* mc64id_ */ +- +-/* ********************************************************************** */ +-/* Subroutine */ int_t mc64ad_(int_t *job, int_t *n, int_t *ne, int_t * +- ip, int_t *irn, double *a, int_t *num, int_t *cperm, +- int_t *liw, int_t *iw, int_t *ldw, double *dw, int_t * +- icntl, int_t *info) +-{ +- /* System generated locals */ +- int_t i__1, i__2; +- double d__1, d__2; +- +- /* Builtin functions */ +- double log(double); +- +- /* Local variables */ +- int_t i__, j, k; +- double fact, rinf; +- +- extern /* Subroutine */ int_t mc21ad_(int_t *, int_t *, int_t *, +- int_t *, int_t *, int_t *, int_t *, int_t *), mc64bd_( +- int_t *, int_t *, int_t *, int_t *, double *, int_t +- *, int_t *, int_t *, int_t *, int_t *, int_t *, +- double *), mc64rd_(int_t *, int_t *, int_t *, int_t *, +- double *), mc64sd_(int_t *, int_t *, int_t *, int_t * +- , double *, int_t *, int_t *, int_t *, int_t *, +- int_t *, int_t *, int_t *, int_t *, int_t *), mc64wd_( +- int_t *, int_t *, int_t *, int_t *, double *, int_t +- *, int_t *, int_t *, int_t *, int_t *, int_t *, int_t +- *, double *, double *); +- +-/* *** Copyright (c) 1999 Council for the Central Laboratory of the */ +-/* Research Councils *** */ +-/* *** Although every effort has been made to ensure robustness and *** */ +-/* *** reliability of the subroutines in this MC64 suite, we *** */ +-/* *** disclaim any liability arising through the use or misuse of *** */ +-/* *** any of the subroutines. *** */ +-/* *** Any problems? Contact ... */ +-/* Iain Duff (I.Duff@rl.ac.uk) or Jacko Koster (jak@ii.uib.no) *** */ +- +-/* Purpose */ +-/* ======= */ +- +-/* This subroutine attempts to find a column permutation for an NxN */ +-/* sparse matrix A = {a_ij} that makes the permuted matrix have N */ +-/* entries on its diagonal. */ +-/* If the matrix is structurally nonsingular, the subroutine optionally */ +-/* returns a column permutation that maximizes the smallest element */ +-/* on the diagonal, maximizes the sum of the diagonal entries, or */ +-/* maximizes the product of the diagonal entries of the permuted matrix. */ +-/* For the latter option, the subroutine also finds scaling factors */ +-/* that may be used to scale the matrix so that the nonzero diagonal */ +-/* entries of the permuted matrix are one in absolute value and all the */ +-/* off-diagonal entries are less than or equal to one in absolute value. */ +-/* The natural logarithms of the scaling factors u(i), i=1..N, for the */ +-/* rows and v(j), j=1..N, for the columns are returned so that the */ +-/* scaled matrix B = {b_ij} has entries b_ij = a_ij * EXP(u_i + v_j). */ +- +-/* Parameters */ +-/* ========== */ +- +- +-/* JOB is an INT_T variable which must be set by the user to */ +-/* control the action. It is not altered by the subroutine. */ +-/* Possible values for JOB are: */ +-/* 1 Compute a column permutation of the matrix so that the */ +-/* permuted matrix has as many entries on its diagonal as possible. */ +-/* The values on the diagonal are of arbitrary size. HSL subroutine */ +-/* MC21A/AD is used for this. See [1]. */ +-/* 2 Compute a column permutation of the matrix so that the smallest */ +-/* value on the diagonal of the permuted matrix is maximized. */ +-/* See [3]. */ +-/* 3 Compute a column permutation of the matrix so that the smallest */ +-/* value on the diagonal of the permuted matrix is maximized. */ +-/* The algorithm differs from the one used for JOB = 2 and may */ +-/* have quite a different performance. See [2]. */ +-/* 4 Compute a column permutation of the matrix so that the sum */ +-/* of the diagonal entries of the permuted matrix is maximized. */ +-/* See [3]. */ +-/* 5 Compute a column permutation of the matrix so that the product */ +-/* of the diagonal entries of the permuted matrix is maximized */ +-/* and vectors to scale the matrix so that the nonzero diagonal */ +-/* entries of the permuted matrix are one in absolute value and */ +-/* all the off-diagonal entries are less than or equal to one in */ +-/* absolute value. See [3]. */ +-/* Restriction: 1 <= JOB <= 5. */ +- +-/* N is an INT_T variable which must be set by the user to the */ +-/* order of the matrix A. It is not altered by the subroutine. */ +-/* Restriction: N >= 1. */ +- +-/* NE is an INT_T variable which must be set by the user to the */ +-/* number of entries in the matrix. It is not altered by the */ +-/* subroutine. */ +-/* Restriction: NE >= 1. */ +- +-/* IP is an INT_T array of length N+1. */ +-/* IP(J), J=1..N, must be set by the user to the position in array IRN */ +-/* of the first row index of an entry in column J. IP(N+1) must be set */ +-/* to NE+1. It is not altered by the subroutine. */ +- +-/* IRN is an INT_T array of length NE. */ +-/* IRN(K), K=1..NE, must be set by the user to hold the row indices of */ +-/* the entries of the matrix. Those belonging to column J must be */ +-/* stored contiguously in the positions IP(J)..IP(J+1)-1. The ordering */ +-/* of the row indices within each column is unimportant. Repeated */ +-/* entries are not allowed. The array IRN is not altered by the */ +-/* subroutine. */ +- +-/* A is a REAL (DOUBLE PRECISION in the D-version) array of length NE. */ +-/* The user must set A(K), K=1..NE, to the numerical value of the */ +-/* entry that corresponds to IRN(K). */ +-/* It is not used by the subroutine when JOB = 1. */ +-/* It is not altered by the subroutine. */ +- +-/* NUM is an INT_T variable that need not be set by the user. */ +-/* On successful exit, NUM will be the number of entries on the */ +-/* diagonal of the permuted matrix. */ +-/* If NUM < N, the matrix is structurally singular. */ +- +-/* CPERM is an INT_T array of length N that need not be set by the */ +-/* user. On successful exit, CPERM contains the column permutation. */ +-/* Column CPERM(J) of the original matrix is column J in the permuted */ +-/* matrix, J=1..N. */ +- +-/* LIW is an INT_T variable that must be set by the user to */ +-/* the dimension of array IW. It is not altered by the subroutine. */ +-/* Restriction: */ +-/* JOB = 1 : LIW >= 5N */ +-/* JOB = 2 : LIW >= 4N */ +-/* JOB = 3 : LIW >= 10N + NE */ +-/* JOB = 4 : LIW >= 5N */ +-/* JOB = 5 : LIW >= 5N */ +- +-/* IW is an INT_T array of length LIW that is used for workspace. */ +- +-/* LDW is an INT_T variable that must be set by the user to the */ +-/* dimension of array DW. It is not altered by the subroutine. */ +-/* Restriction: */ +-/* JOB = 1 : LDW is not used */ +-/* JOB = 2 : LDW >= N */ +-/* JOB = 3 : LDW >= NE */ +-/* JOB = 4 : LDW >= 2N + NE */ +-/* JOB = 5 : LDW >= 3N + NE */ +- +-/* DW is a REAL (DOUBLE PRECISION in the D-version) array of length LDW */ +-/* that is used for workspace. If JOB = 5, on return, */ +-/* DW(i) contains u_i, i=1..N, and DW(N+j) contains v_j, j=1..N. */ +- +-/* ICNTL is an INT_T array of length 10. Its components control the */ +-/* output of MC64A/AD and must be set by the user before calling */ +-/* MC64A/AD. They are not altered by the subroutine. */ +- +-/* ICNTL(1) must be set to specify the output stream for */ +-/* error messages. If ICNTL(1) < 0, messages are suppressed. */ +-/* The default value set by MC46I/ID is 6. */ +- +-/* ICNTL(2) must be set by the user to specify the output stream for */ +-/* warning messages. If ICNTL(2) < 0, messages are suppressed. */ +-/* The default value set by MC46I/ID is 6. */ +- +-/* ICNTL(3) must be set by the user to specify the output stream for */ +-/* diagnostic messages. If ICNTL(3) < 0, messages are suppressed. */ +-/* The default value set by MC46I/ID is -1. */ +- +-/* ICNTL(4) must be set by the user to a value other than 0 to avoid */ +-/* checking of the input data. */ +-/* The default value set by MC46I/ID is 0. */ +- +-/* INFO is an INT_T array of length 10 which need not be set by the */ +-/* user. INFO(1) is set non-negative to indicate success. A negative */ +-/* value is returned if an error occurred, a positive value if a */ +-/* warning occurred. INFO(2) holds further information on the error. */ +-/* On exit from the subroutine, INFO(1) will take one of the */ +-/* following values: */ +-/* 0 : successful entry (for structurally nonsingular matrix). */ +-/* +1 : successful entry (for structurally singular matrix). */ +-/* +2 : the returned scaling factors are large and may cause */ +-/* overflow when used to scale the matrix. */ +-/* (For JOB = 5 entry only.) */ +-/* -1 : JOB < 1 or JOB > 5. Value of JOB held in INFO(2). */ +-/* -2 : N < 1. Value of N held in INFO(2). */ +-/* -3 : NE < 1. Value of NE held in INFO(2). */ +-/* -4 : the defined length LIW violates the restriction on LIW. */ +-/* Value of LIW required given by INFO(2). */ +-/* -5 : the defined length LDW violates the restriction on LDW. */ +-/* Value of LDW required given by INFO(2). */ +-/* -6 : entries are found whose row indices are out of range. INFO(2) */ +-/* contains the index of a column in which such an entry is found. */ +-/* -7 : repeated entries are found. INFO(2) contains the index of a */ +-/* column in which such entries are found. */ +-/* INFO(3) to INFO(10) are not currently used and are set to zero by */ +-/* the routine. */ +- +-/* References: */ +-/* [1] I. S. Duff, (1981), */ +-/* "Algorithm 575. Permutations for a zero-free diagonal", */ +-/* ACM Trans. Math. Software 7(3), 387-390. */ +-/* [2] I. S. Duff and J. Koster, (1998), */ +-/* "The design and use of algorithms for permuting large */ +-/* entries to the diagonal of sparse matrices", */ +-/* SIAM J. Matrix Anal. Appl., vol. 20, no. 4, pp. 889-901. */ +-/* [3] I. S. Duff and J. Koster, (1999), */ +-/* "On algorithms for permuting large entries to the diagonal */ +-/* of sparse matrices", */ +-/* Technical Report RAL-TR-1999-030, RAL, Oxfordshire, England. */ +-/* Local variables and parameters */ +-/* External routines and functions */ +-/* EXTERNAL FD05AD */ +-/* DOUBLE PRECISION FD05AD */ +-/* Intrinsic functions */ +-/* Set RINF to largest positive real number (infinity) */ +-/* XSL RINF = FD05AD(5) */ +- /* Parameter adjustments */ +- --cperm; +- --ip; +- --a; +- --irn; +- --iw; +- --dw; +- --icntl; +- --info; +- +- /* Function Body */ +- rinf = dmach("Overflow"); +-/* Check value of JOB */ +- if (*job < 1 || *job > 5) { +- info[1] = -1; +- info[2] = *job; +- if (icntl[1] >= 0) { +- printf(" ****** Error in MC64A/AD. INFO(1) = %2d" +- " because JOB = %d\n", info[1], *job); +- } +- goto L99; +- } +-/* Check value of N */ +- if (*n < 1) { +- info[1] = -2; +- info[2] = *n; +- if (icntl[1] >= 0) { +- printf(" ****** Error in MC64A/AD. INFO(1) = %2d" +- " because N = %d\n", info[1], *job); +- } +- goto L99; +- } +-/* Check value of NE */ +- if (*ne < 1) { +- info[1] = -3; +- info[2] = *ne; +- if (icntl[1] >= 0) { +- printf(" ****** Error in MC64A/AD. INFO(1) = %2d" +- " because NE = %d\n", info[1], *job); +- } +- goto L99; +- } +-/* Check LIW */ +- if (*job == 1) { +- k = *n * 5; +- } +- if (*job == 2) { +- k = *n << 2; +- } +- if (*job == 3) { +- k = *n * 10 + *ne; +- } +- if (*job == 4) { +- k = *n * 5; +- } +- if (*job == 5) { +- k = *n * 5; +- } +- if (*liw < k) { +- info[1] = -4; +- info[2] = k; +- if (icntl[1] >= 0) { +- printf(" ****** Error in MC64A/AD. INFO(1) = %2d" +- " LIW too small, must be at least %8d\n", info[1], k); +- } +- goto L99; +- } +-/* Check LDW */ +-/* If JOB = 1, do not check */ +- if (*job > 1) { +- if (*job == 2) { +- k = *n; +- } +- if (*job == 3) { +- k = *ne; +- } +- if (*job == 4) { +- k = (*n << 1) + *ne; +- } +- if (*job == 5) { +- k = *n * 3 + *ne; +- } +- if (*ldw < k) { +- info[1] = -5; +- info[2] = k; +- if (icntl[1] >= 0) { +- printf(" ****** Error in MC64A/AD. INFO(1) = %2d" +- " LDW too small, must be at least %8d\n", info[1], k); +- } +- goto L99; +- } +- } +- if (icntl[4] == 0) { +-/* Check row indices. Use IW(1:N) as workspace */ +- i__1 = *n; +- for (i__ = 1; i__ <= i__1; ++i__) { +- iw[i__] = 0; +-/* L3: */ +- } +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +- i__2 = ip[j + 1] - 1; +- for (k = ip[j]; k <= i__2; ++k) { +- i__ = irn[k]; +-/* Check for row indices that are out of range */ +- if (i__ < 1 || i__ > *n) { +- info[1] = -6; +- info[2] = j; +- if (icntl[1] >= 0) { +- printf(" ****** Error in MC64A/AD. INFO(1) = %2d Column %8d" +- " contains an entry with invalid row index %8d\n", +- info[1], j, i__); +- } +- goto L99; +- } +-/* Check for repeated row indices within a column */ +- if (iw[i__] == j) { +- info[1] = -7; +- info[2] = j; +- if (icntl[1] >= 0) { +- printf(" ****** Error in MC64A/AD. INFO(1) = %2d" +- " Column %8d" +- " contains two or more entries with row index %8d\n", +- info[1], j, i__); +- } +- goto L99; +- } else { +- iw[i__] = j; +- } +-/* L4: */ +- } +-/* L6: */ +- } +- } +-/* Print diagnostics on input */ +- if (icntl[3] >= 0) { +- printf(" ****** Input parameters for MC64A/AD: JOB = %8d," +- " N = %d, NE = %8d\n", *job, *n, *ne); +- printf(" IP(1:N+1) = "); +- for (j=1; j<=(*n+1); ++j) { +- printf("%8d", ip[j]); +- if (j%8 == 0) printf("\n"); +- } +- printf("\n IRN(1:NE) = "); +- for (j=1; j<=(*ne); ++j) { +- printf("%8d", irn[j]); +- if (j%8 == 0) printf("\n"); +- } +- printf("\n"); +- +- if (*job > 1) { +- printf(" A(1:NE) = "); +- for (j=1; j<=(*ne); ++j) { +- printf("%f14.4", a[j]); +- if (j%4 == 0) printf("\n"); +- } +- printf("\n"); +- } +- } +-/* Set components of INFO to zero */ +- for (i__ = 1; i__ <= 10; ++i__) { +- info[i__] = 0; +-/* L8: */ +- } +-/* Compute maximum matching with MC21A/AD */ +- if (*job == 1) { +-/* Put length of column J in IW(J) */ +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +- iw[j] = ip[j + 1] - ip[j]; +-/* L10: */ +- } +-/* IW(N+1:5N) is workspace */ +-#if 0 +- mc21ad_(n, &irn[1], ne, &ip[1], &iw[1], &cperm[1], num, &iw[*n+1]); +-#else +- printf(" ****** Warning from MC64A/AD. Need to link mc21ad.\n"); +-#endif +- goto L90; +- } +-/* Compute bottleneck matching */ +- if (*job == 2) { +-/* IW(1:5N), DW(1:N) are workspaces */ +- mc64bd_(n, ne, &ip[1], &irn[1], &a[1], &cperm[1], num, &iw[1], &iw[*n +- + 1], &iw[(*n << 1) + 1], &iw[*n * 3 + 1], &dw[1]); +- goto L90; +- } +-/* Compute bottleneck matching */ +- if (*job == 3) { +-/* Copy IRN(K) into IW(K), ABS(A(K)) into DW(K), K=1..NE */ +- i__1 = *ne; +- for (k = 1; k <= i__1; ++k) { +- iw[k] = irn[k]; +- dw[k] = (d__1 = a[k], abs(d__1)); +-/* L20: */ +- } +-/* Sort entries in each column by decreasing value. */ +- mc64rd_(n, ne, &ip[1], &iw[1], &dw[1]); +-/* IW(NE+1:NE+10N) is workspace */ +- mc64sd_(n, ne, &ip[1], &iw[1], &dw[1], &cperm[1], num, &iw[*ne + 1], & +- iw[*ne + *n + 1], &iw[*ne + (*n << 1) + 1], &iw[*ne + *n * 3 +- + 1], &iw[*ne + (*n << 2) + 1], &iw[*ne + *n * 5 + 1], &iw[* +- ne + *n * 6 + 1]); +- goto L90; +- } +- if (*job == 4) { +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +- fact = 0.; +- i__2 = ip[j + 1] - 1; +- for (k = ip[j]; k <= i__2; ++k) { +- if ((d__1 = a[k], abs(d__1)) > fact) { +- fact = (d__2 = a[k], abs(d__2)); +- } +-/* L30: */ +- } +- i__2 = ip[j + 1] - 1; +- for (k = ip[j]; k <= i__2; ++k) { +- dw[(*n << 1) + k] = fact - (d__1 = a[k], abs(d__1)); +-/* L40: */ +- } +-/* L50: */ +- } +-/* B = DW(2N+1:2N+NE); IW(1:5N) and DW(1:2N) are workspaces */ +- mc64wd_(n, ne, &ip[1], &irn[1], &dw[(*n << 1) + 1], &cperm[1], num, & +- iw[1], &iw[*n + 1], &iw[(*n << 1) + 1], &iw[*n * 3 + 1], &iw[( +- *n << 2) + 1], &dw[1], &dw[*n + 1]); +- goto L90; +- } +- if (*job == 5) { +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +- fact = 0.; +- i__2 = ip[j + 1] - 1; +- for (k = ip[j]; k <= i__2; ++k) { +- dw[*n * 3 + k] = (d__1 = a[k], abs(d__1)); +- if (dw[*n * 3 + k] > fact) { +- fact = dw[*n * 3 + k]; +- } +-/* L60: */ +- } +- dw[(*n << 1) + j] = fact; +- if (fact != 0.) { +- fact = log(fact); +- } else { +- fact = rinf / *n; +- } +- i__2 = ip[j + 1] - 1; +- for (k = ip[j]; k <= i__2; ++k) { +- if (dw[*n * 3 + k] != 0.) { +- dw[*n * 3 + k] = fact - log(dw[*n * 3 + k]); +- } else { +- dw[*n * 3 + k] = rinf / *n; +- } +-/* L70: */ +- } +-/* L75: */ +- } +-/* B = DW(3N+1:3N+NE); IW(1:5N) and DW(1:2N) are workspaces */ +- mc64wd_(n, ne, &ip[1], &irn[1], &dw[*n * 3 + 1], &cperm[1], num, &iw[ +- 1], &iw[*n + 1], &iw[(*n << 1) + 1], &iw[*n * 3 + 1], &iw[(*n +- << 2) + 1], &dw[1], &dw[*n + 1]); +- if (*num == *n) { +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +- if (dw[(*n << 1) + j] != 0.) { +- dw[*n + j] -= log(dw[(*n << 1) + j]); +- } else { +- dw[*n + j] = 0.; +- } +-/* L80: */ +- } +- } +-/* Check size of scaling factors */ +- fact = log(rinf) * .5f; +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +- if (dw[j] < fact && dw[*n + j] < fact) { +- goto L86; +- } +- info[1] = 2; +- goto L90; +-L86: +- ; +- } +-/* GO TO 90 */ +- } +-L90: +- if (info[1] == 0 && *num < *n) { +-/* Matrix is structurally singular, return with warning */ +- info[1] = 1; +- if (icntl[2] >= 0) { +- printf(" ****** Warning from MC64A/AD. INFO(1) = %2d" +- " The matrix is structurally singular.\n", info[1]); +- } +- } +- if (info[1] == 2) { +-/* Scaling factors are large, return with warning */ +- if (icntl[2] >= 0) { +- printf(" ****** Warning from MC64A/AD. INFO(1) = %2d\n" +- " Some scaling factors may be too large.\n", info[1]); +- } +- } +-/* Print diagnostics on output */ +- if (icntl[3] >= 0) { +- printf(" ****** Output parameters for MC64A/AD: INFO(1:2) = %8d%8d\n", +- info[1], info[2]); +- printf(" NUM = %8d", *num); +- printf(" CPERM(1:N) = "); +- for (j=1; j<=*n; ++j) { +- printf("%8d", cperm[j]); +- if (j%8 == 0) printf("\n"); +- } +- if (*job == 5) { +- printf("\n DW(1:N) = "); +- for (j=1; j<=*n; ++j) { +- printf("%11.3f", dw[j]); +- if (j%5 == 0) printf("\n"); +- } +- printf("\n DW(N+1:2N) = "); +- for (j=1; j<=*n; ++j) { +- printf("%11.3f", dw[*n+j]); +- if (j%5 == 0) printf("\n"); +- } +- printf("\n"); +- } +- } +-/* Return from subroutine. */ +-L99: +- return 0; +-} /* mc64ad_ */ +- +-/* ********************************************************************** */ +-/* Subroutine */ int_t mc64bd_(int_t *n, int_t *ne, int_t *ip, int_t * +- irn, double *a, int_t *iperm, int_t *num, int_t *jperm, +- int_t *pr, int_t *q, int_t *l, double *d__) +-{ +- /* System generated locals */ +- int_t i__1, i__2, i__3; +- double d__1, d__2, d__3; +- int_t c__1 = 1; +- +- /* Local variables */ +- int_t i__, j, k; +- double a0; +- int_t i0, q0; +- double ai, di; +- int_t ii, jj, kk; +- double bv; +- int_t up; +- double dq0; +- int_t kk1, kk2; +- double csp; +- int_t isp, jsp, low; +- double dnew; +- int_t jord, qlen, idum, jdum; +- double rinf; +- extern /* Subroutine */ int_t mc64dd_(int_t *, int_t *, int_t *, +- double *, int_t *, int_t *), mc64ed_(int_t *, int_t *, +- int_t *, double *, int_t *, int_t *), mc64fd_(int_t * +- , int_t *, int_t *, int_t *, double *, int_t *, int_t *); +- +- +-/* *** Copyright (c) 1999 Council for the Central Laboratory of the */ +-/* Research Councils *** */ +-/* *** Although every effort has been made to ensure robustness and *** */ +-/* *** reliability of the subroutines in this MC64 suite, we *** */ +-/* *** disclaim any liability arising through the use or misuse of *** */ +-/* *** any of the subroutines. *** */ +-/* *** Any problems? Contact ... */ +-/* Iain Duff (I.Duff@rl.ac.uk) or Jacko Koster (jak@ii.uib.no) *** */ +- +-/* N, NE, IP, IRN are described in MC64A/AD. */ +-/* A is a REAL (DOUBLE PRECISION in the D-version) array of length */ +-/* NE. A(K), K=1..NE, must be set to the value of the entry */ +-/* that corresponds to IRN(K). It is not altered. */ +-/* IPERM is an INT_T array of length N. On exit, it contains the */ +-/* matching: IPERM(I) = 0 or row I is matched to column IPERM(I). */ +-/* NUM is INT_T variable. On exit, it contains the cardinality of the */ +-/* matching stored in IPERM. */ +-/* IW is an INT_T work array of length 4N. */ +-/* DW is a REAL (DOUBLE PRECISION in D-version) work array of length N. */ +-/* Local variables */ +-/* Local parameters */ +-/* Intrinsic functions */ +-/* External subroutines and/or functions */ +-/* EXTERNAL FD05AD,MC64DD,MC64ED,MC64FD, DMACH */ +-/* DOUBLE PRECISION FD05AD, DMACH */ +-/* Set RINF to largest positive real number */ +-/* XSL RINF = FD05AD(5) */ +- /* Parameter adjustments */ +- --d__; +- --l; +- --q; +- --pr; +- --jperm; +- --iperm; +- --ip; +- --a; +- --irn; +- +- /* Function Body */ +- rinf = dmach("Overflow"); +-/* Initialization */ +- *num = 0; +- bv = rinf; +- i__1 = *n; +- for (k = 1; k <= i__1; ++k) { +- iperm[k] = 0; +- jperm[k] = 0; +- pr[k] = ip[k]; +- d__[k] = 0.; +-/* L10: */ +- } +-/* Scan columns of matrix; */ +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +- a0 = -1.; +- i__2 = ip[j + 1] - 1; +- for (k = ip[j]; k <= i__2; ++k) { +- i__ = irn[k]; +- ai = (d__1 = a[k], abs(d__1)); +- if (ai > d__[i__]) { +- d__[i__] = ai; +- } +- if (jperm[j] != 0) { +- goto L30; +- } +- if (ai >= bv) { +- a0 = bv; +- if (iperm[i__] != 0) { +- goto L30; +- } +- jperm[j] = i__; +- iperm[i__] = j; +- ++(*num); +- } else { +- if (ai <= a0) { +- goto L30; +- } +- a0 = ai; +- i0 = i__; +- } +-L30: +- ; +- } +- if (a0 != -1. && a0 < bv) { +- bv = a0; +- if (iperm[i0] != 0) { +- goto L20; +- } +- iperm[i0] = j; +- jperm[j] = i0; +- ++(*num); +- } +-L20: +- ; +- } +-/* Update BV with smallest of all the largest maximum absolute values */ +-/* of the rows. */ +- i__1 = *n; +- for (i__ = 1; i__ <= i__1; ++i__) { +-/* Computing MIN */ +- d__1 = bv, d__2 = d__[i__]; +- bv = min(d__1,d__2); +-/* L25: */ +- } +- if (*num == *n) { +- goto L1000; +- } +-/* Rescan unassigned columns; improve initial assignment */ +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +- if (jperm[j] != 0) { +- goto L95; +- } +- i__2 = ip[j + 1] - 1; +- for (k = ip[j]; k <= i__2; ++k) { +- i__ = irn[k]; +- ai = (d__1 = a[k], abs(d__1)); +- if (ai < bv) { +- goto L50; +- } +- if (iperm[i__] == 0) { +- goto L90; +- } +- jj = iperm[i__]; +- kk1 = pr[jj]; +- kk2 = ip[jj + 1] - 1; +- if (kk1 > kk2) { +- goto L50; +- } +- i__3 = kk2; +- for (kk = kk1; kk <= i__3; ++kk) { +- ii = irn[kk]; +- if (iperm[ii] != 0) { +- goto L70; +- } +- if ((d__1 = a[kk], abs(d__1)) >= bv) { +- goto L80; +- } +-L70: +- ; +- } +- pr[jj] = kk2 + 1; +-L50: +- ; +- } +- goto L95; +-L80: +- jperm[jj] = ii; +- iperm[ii] = jj; +- pr[jj] = kk + 1; +-L90: +- ++(*num); +- jperm[j] = i__; +- iperm[i__] = j; +- pr[j] = k + 1; +-L95: +- ; +- } +- if (*num == *n) { +- goto L1000; +- } +-/* Prepare for main loop */ +- i__1 = *n; +- for (i__ = 1; i__ <= i__1; ++i__) { +- d__[i__] = -1.; +- l[i__] = 0; +-/* L99: */ +- } +-/* Main loop ... each pass round this loop is similar to Dijkstra's */ +-/* algorithm for solving the single source shortest path problem */ +- i__1 = *n; +- for (jord = 1; jord <= i__1; ++jord) { +- if (jperm[jord] != 0) { +- goto L100; +- } +- qlen = 0; +- low = *n + 1; +- up = *n + 1; +-/* CSP is cost of shortest path to any unassigned row */ +-/* ISP is matrix position of unassigned row element in shortest path */ +-/* JSP is column index of unassigned row element in shortest path */ +- csp = -1.; +-/* Build shortest path tree starting from unassigned column JORD */ +- j = jord; +- pr[j] = -1; +-/* Scan column J */ +- i__2 = ip[j + 1] - 1; +- for (k = ip[j]; k <= i__2; ++k) { +- i__ = irn[k]; +- dnew = (d__1 = a[k], abs(d__1)); +- if (csp >= dnew) { +- goto L115; +- } +- if (iperm[i__] == 0) { +-/* Row I is unassigned; update shortest path info */ +- csp = dnew; +- isp = i__; +- jsp = j; +- if (csp >= bv) { +- goto L160; +- } +- } else { +- d__[i__] = dnew; +- if (dnew >= bv) { +-/* Add row I to Q2 */ +- --low; +- q[low] = i__; +- } else { +-/* Add row I to Q, and push it */ +- ++qlen; +- l[i__] = qlen; +- mc64dd_(&i__, n, &q[1], &d__[1], &l[1], &c__1); +- } +- jj = iperm[i__]; +- pr[jj] = j; +- } +-L115: +- ; +- } +- i__2 = *num; +- for (jdum = 1; jdum <= i__2; ++jdum) { +-/* If Q2 is empty, extract new rows from Q */ +- if (low == up) { +- if (qlen == 0) { +- goto L160; +- } +- i__ = q[1]; +- if (csp >= d__[i__]) { +- goto L160; +- } +- bv = d__[i__]; +- i__3 = *n; +- for (idum = 1; idum <= i__3; ++idum) { +- mc64ed_(&qlen, n, &q[1], &d__[1], &l[1], &c__1); +- l[i__] = 0; +- --low; +- q[low] = i__; +- if (qlen == 0) { +- goto L153; +- } +- i__ = q[1]; +- if (d__[i__] != bv) { +- goto L153; +- } +-/* L152: */ +- } +-/* End of dummy loop; this point is never reached */ +- } +-/* Move row Q0 */ +-L153: +- --up; +- q0 = q[up]; +- dq0 = d__[q0]; +- l[q0] = up; +-/* Scan column that matches with row Q0 */ +- j = iperm[q0]; +- i__3 = ip[j + 1] - 1; +- for (k = ip[j]; k <= i__3; ++k) { +- i__ = irn[k]; +-/* Update D(I) */ +- if (l[i__] >= up) { +- goto L155; +- } +-/* Computing MIN */ +- d__2 = dq0, d__3 = (d__1 = a[k], abs(d__1)); +- dnew = min(d__2,d__3); +- if (csp >= dnew) { +- goto L155; +- } +- if (iperm[i__] == 0) { +-/* Row I is unassigned; update shortest path info */ +- csp = dnew; +- isp = i__; +- jsp = j; +- if (csp >= bv) { +- goto L160; +- } +- } else { +- di = d__[i__]; +- if (di >= bv || di >= dnew) { +- goto L155; +- } +- d__[i__] = dnew; +- if (dnew >= bv) { +-/* Delete row I from Q (if necessary); add row I to Q2 */ +- if (di != -1.) { +- mc64fd_(&l[i__], &qlen, n, &q[1], &d__[1], &l[1], +- &c__1); +- } +- l[i__] = 0; +- --low; +- q[low] = i__; +- } else { +-/* Add row I to Q (if necessary); push row I up Q */ +- if (di == -1.) { +- ++qlen; +- l[i__] = qlen; +- } +- mc64dd_(&i__, n, &q[1], &d__[1], &l[1], &c__1); +- } +-/* Update tree */ +- jj = iperm[i__]; +- pr[jj] = j; +- } +-L155: +- ; +- } +-/* L150: */ +- } +-/* If CSP = MINONE, no augmenting path is found */ +-L160: +- if (csp == -1.) { +- goto L190; +- } +-/* Update bottleneck value */ +- bv = min(bv,csp); +-/* Find augmenting path by tracing backward in PR; update IPERM,JPERM */ +- ++(*num); +- i__ = isp; +- j = jsp; +- i__2 = *num + 1; +- for (jdum = 1; jdum <= i__2; ++jdum) { +- i0 = jperm[j]; +- jperm[j] = i__; +- iperm[i__] = j; +- j = pr[j]; +- if (j == -1) { +- goto L190; +- } +- i__ = i0; +-/* L170: */ +- } +-/* End of dummy loop; this point is never reached */ +-L190: +- i__2 = *n; +- for (kk = up; kk <= i__2; ++kk) { +- i__ = q[kk]; +- d__[i__] = -1.; +- l[i__] = 0; +-/* L191: */ +- } +- i__2 = up - 1; +- for (kk = low; kk <= i__2; ++kk) { +- i__ = q[kk]; +- d__[i__] = -1.; +-/* L192: */ +- } +- i__2 = qlen; +- for (kk = 1; kk <= i__2; ++kk) { +- i__ = q[kk]; +- d__[i__] = -1.; +- l[i__] = 0; +-/* L193: */ +- } +-L100: +- ; +- } +-/* End of main loop */ +-/* BV is bottleneck value of final matching */ +- if (*num == *n) { +- goto L1000; +- } +-/* Matrix is structurally singular, complete IPERM. */ +-/* JPERM, PR are work arrays */ +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +- jperm[j] = 0; +-/* L300: */ +- } +- k = 0; +- i__1 = *n; +- for (i__ = 1; i__ <= i__1; ++i__) { +- if (iperm[i__] == 0) { +- ++k; +- pr[k] = i__; +- } else { +- j = iperm[i__]; +- jperm[j] = i__; +- } +-/* L310: */ +- } +- k = 0; +- i__1 = *n; +- for (i__ = 1; i__ <= i__1; ++i__) { +- if (jperm[i__] != 0) { +- goto L320; +- } +- ++k; +- jdum = pr[k]; +- iperm[jdum] = i__; +-L320: +- ; +- } +-L1000: +- return 0; +-} /* mc64bd_ */ +- +-/* ********************************************************************** */ +-/* Subroutine */ int_t mc64dd_(int_t *i__, int_t *n, int_t *q, double +- *d__, int_t *l, int_t *iway) +-{ +- /* System generated locals */ +- int_t i__1; +- int_t c__1 = 1; +- +- /* Local variables */ +- double di; +- int_t qk, pos, idum, posk; +- +- +-/* *** Copyright (c) 1999 Council for the Central Laboratory of the */ +-/* Research Councils *** */ +-/* *** Although every effort has been made to ensure robustness and *** */ +-/* *** reliability of the subroutines in this MC64 suite, we *** */ +-/* *** disclaim any liability arising through the use or misuse of *** */ +-/* *** any of the subroutines. *** */ +-/* *** Any problems? Contact ... */ +-/* Iain Duff (I.Duff@rl.ac.uk) or Jacko Koster (jak@ii.uib.no) *** */ +- +-/* Variables N,Q,D,L are described in MC64B/BD */ +-/* IF IWAY is equal to 1, then */ +-/* node I is pushed from its current position upwards */ +-/* IF IWAY is not equal to 1, then */ +-/* node I is pushed from its current position downwards */ +-/* Local variables and parameters */ +- /* Parameter adjustments */ +- --l; +- --d__; +- --q; +- +- /* Function Body */ +- di = d__[*i__]; +- pos = l[*i__]; +-/* POS is index of current position of I in the tree */ +- if (*iway == 1) { +- i__1 = *n; +- for (idum = 1; idum <= i__1; ++idum) { +- if (pos <= 1) { +- goto L20; +- } +- posk = pos / 2; +- qk = q[posk]; +- if (di <= d__[qk]) { +- goto L20; +- } +- q[pos] = qk; +- l[qk] = pos; +- pos = posk; +-/* L10: */ +- } +-/* End of dummy loop; this point is never reached */ +- } else { +- i__1 = *n; +- for (idum = 1; idum <= i__1; ++idum) { +- if (pos <= 1) { +- goto L20; +- } +- posk = pos / 2; +- qk = q[posk]; +- if (di >= d__[qk]) { +- goto L20; +- } +- q[pos] = qk; +- l[qk] = pos; +- pos = posk; +-/* L15: */ +- } +-/* End of dummy loop; this point is never reached */ +- } +-/* End of dummy if; this point is never reached */ +-L20: +- q[pos] = *i__; +- l[*i__] = pos; +- return 0; +-} /* mc64dd_ */ +- +-/* ********************************************************************** */ +-/* Subroutine */ int_t mc64ed_(int_t *qlen, int_t *n, int_t *q, +- double *d__, int_t *l, int_t *iway) +-{ +- /* System generated locals */ +- int_t i__1; +- +- /* Local variables */ +- int_t i__; +- double di, dk, dr; +- int_t pos, idum, posk; +- +- +-/* *** Copyright (c) 1999 Council for the Central Laboratory of the */ +-/* Research Councils *** */ +-/* *** Although every effort has been made to ensure robustness and *** */ +-/* *** reliability of the subroutines in this MC64 suite, we *** */ +-/* *** disclaim any liability arising through the use or misuse of *** */ +-/* *** any of the subroutines. *** */ +-/* *** Any problems? Contact ... */ +-/* Iain Duff (I.Duff@rl.ac.uk) or Jacko Koster (jak@ii.uib.no) *** */ +- +-/* Variables QLEN,N,Q,D,L are described in MC64B/BD (IWAY = 1) or */ +-/* MC64W/WD (IWAY = 2) */ +-/* The root node is deleted from the binary heap. */ +-/* Local variables and parameters */ +-/* Move last element to begin of Q */ +- /* Parameter adjustments */ +- --l; +- --d__; +- --q; +- +- /* Function Body */ +- i__ = q[*qlen]; +- di = d__[i__]; +- --(*qlen); +- pos = 1; +- if (*iway == 1) { +- i__1 = *n; +- for (idum = 1; idum <= i__1; ++idum) { +- posk = pos << 1; +- if (posk > *qlen) { +- goto L20; +- } +- dk = d__[q[posk]]; +- if (posk < *qlen) { +- dr = d__[q[posk + 1]]; +- if (dk < dr) { +- ++posk; +- dk = dr; +- } +- } +- if (di >= dk) { +- goto L20; +- } +-/* Exchange old last element with larger priority child */ +- q[pos] = q[posk]; +- l[q[pos]] = pos; +- pos = posk; +-/* L10: */ +- } +-/* End of dummy loop; this point is never reached */ +- } else { +- i__1 = *n; +- for (idum = 1; idum <= i__1; ++idum) { +- posk = pos << 1; +- if (posk > *qlen) { +- goto L20; +- } +- dk = d__[q[posk]]; +- if (posk < *qlen) { +- dr = d__[q[posk + 1]]; +- if (dk > dr) { +- ++posk; +- dk = dr; +- } +- } +- if (di <= dk) { +- goto L20; +- } +-/* Exchange old last element with smaller child */ +- q[pos] = q[posk]; +- l[q[pos]] = pos; +- pos = posk; +-/* L15: */ +- } +-/* End of dummy loop; this point is never reached */ +- } +-/* End of dummy if; this point is never reached */ +-L20: +- q[pos] = i__; +- l[i__] = pos; +- return 0; +-} /* mc64ed_ */ +- +-/* ********************************************************************** */ +-/* Subroutine */ int_t mc64fd_(int_t *pos0, int_t *qlen, int_t *n, +- int_t *q, double *d__, int_t *l, int_t *iway) +-{ +- /* System generated locals */ +- int_t i__1; +- +- /* Local variables */ +- int_t i__; +- double di, dk, dr; +- int_t qk, pos, idum, posk; +- +- +-/* *** Copyright (c) 1999 Council for the Central Laboratory of the */ +-/* Research Councils *** */ +-/* *** Although every effort has been made to ensure robustness and *** */ +-/* *** reliability of the subroutines in this MC64 suite, we *** */ +-/* *** disclaim any liability arising through the use or misuse of *** */ +-/* *** any of the subroutines. *** */ +-/* *** Any problems? Contact ... */ +-/* Iain Duff (I.Duff@rl.ac.uk) or Jacko Koster (jak@ii.uib.no) *** */ +- +-/* Variables QLEN,N,Q,D,L are described in MC64B/BD (IWAY = 1) or */ +-/* MC64WD (IWAY = 2). */ +-/* Move last element in the heap */ +-/* Quick return, if possible */ +- /* Parameter adjustments */ +- --l; +- --d__; +- --q; +- +- /* Function Body */ +- if (*qlen == *pos0) { +- --(*qlen); +- return 0; +- } +-/* Move last element from queue Q to position POS0 */ +-/* POS is current position of node I in the tree */ +- i__ = q[*qlen]; +- di = d__[i__]; +- --(*qlen); +- pos = *pos0; +- if (*iway == 1) { +- i__1 = *n; +- for (idum = 1; idum <= i__1; ++idum) { +- if (pos <= 1) { +- goto L20; +- } +- posk = pos / 2; +- qk = q[posk]; +- if (di <= d__[qk]) { +- goto L20; +- } +- q[pos] = qk; +- l[qk] = pos; +- pos = posk; +-/* L10: */ +- } +-/* End of dummy loop; this point is never reached */ +-L20: +- q[pos] = i__; +- l[i__] = pos; +- i__1 = *n; +- for (idum = 1; idum <= i__1; ++idum) { +- posk = pos << 1; +- if (posk > *qlen) { +- goto L40; +- } +- dk = d__[q[posk]]; +- if (posk < *qlen) { +- dr = d__[q[posk + 1]]; +- if (dk < dr) { +- ++posk; +- dk = dr; +- } +- } +- if (di >= dk) { +- goto L40; +- } +- qk = q[posk]; +- q[pos] = qk; +- l[qk] = pos; +- pos = posk; +-/* L30: */ +- } +-/* End of dummy loop; this point is never reached */ +- } else { +- i__1 = *n; +- for (idum = 1; idum <= i__1; ++idum) { +- if (pos <= 1) { +- goto L34; +- } +- posk = pos / 2; +- qk = q[posk]; +- if (di >= d__[qk]) { +- goto L34; +- } +- q[pos] = qk; +- l[qk] = pos; +- pos = posk; +-/* L32: */ +- } +-/* End of dummy loop; this point is never reached */ +-L34: +- q[pos] = i__; +- l[i__] = pos; +- i__1 = *n; +- for (idum = 1; idum <= i__1; ++idum) { +- posk = pos << 1; +- if (posk > *qlen) { +- goto L40; +- } +- dk = d__[q[posk]]; +- if (posk < *qlen) { +- dr = d__[q[posk + 1]]; +- if (dk > dr) { +- ++posk; +- dk = dr; +- } +- } +- if (di <= dk) { +- goto L40; +- } +- qk = q[posk]; +- q[pos] = qk; +- l[qk] = pos; +- pos = posk; +-/* L36: */ +- } +-/* End of dummy loop; this point is never reached */ +- } +-/* End of dummy if; this point is never reached */ +-L40: +- q[pos] = i__; +- l[i__] = pos; +- return 0; +-} /* mc64fd_ */ +- +-/* ********************************************************************** */ +-/* Subroutine */ int_t mc64rd_(int_t *n, int_t *ne, int_t *ip, int_t * +- irn, double *a) +-{ +- /* System generated locals */ +- int_t i__1, i__2, i__3; +- +- /* Local variables */ +- int_t j, k, r__, s; +- double ha; +- int_t hi, td, mid, len, ipj; +- double key; +- int_t last, todo[50], first; +- +- +-/* *** Copyright (c) 1999 Council for the Central Laboratory of the */ +-/* Research Councils *** */ +-/* *** Although every effort has been made to ensure robustness and *** */ +-/* *** reliability of the subroutines in this MC64 suite, we *** */ +-/* *** disclaim any liability arising through the use or misuse of *** */ +-/* *** any of the subroutines. *** */ +-/* *** Any problems? Contact ... */ +-/* Iain Duff (I.Duff@rl.ac.uk) or Jacko Koster (jak@ii.uib.no) *** */ +- +-/* This subroutine sorts the entries in each column of the */ +-/* sparse matrix (defined by N,NE,IP,IRN,A) by decreasing */ +-/* numerical value. */ +-/* Local constants */ +-/* Local variables */ +-/* Local arrays */ +- /* Parameter adjustments */ +- --ip; +- --a; +- --irn; +- +- /* Function Body */ +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +- len = ip[j + 1] - ip[j]; +- if (len <= 1) { +- goto L100; +- } +- ipj = ip[j]; +-/* Sort array roughly with partial quicksort */ +- if (len < 15) { +- goto L400; +- } +- todo[0] = ipj; +- todo[1] = ipj + len; +- td = 2; +-L500: +- first = todo[td - 2]; +- last = todo[td - 1]; +-/* KEY is the smallest of two values present in interval [FIRST,LAST) */ +- key = a[(first + last) / 2]; +- i__2 = last - 1; +- for (k = first; k <= i__2; ++k) { +- ha = a[k]; +- if (ha == key) { +- goto L475; +- } +- if (ha > key) { +- goto L470; +- } +- key = ha; +- goto L470; +-L475: +- ; +- } +-/* Only one value found in interval, so it is already sorted */ +- td += -2; +- goto L425; +-/* Reorder interval [FIRST,LAST) such that entries before MID are gt KEY */ +-L470: +- mid = first; +- i__2 = last - 1; +- for (k = first; k <= i__2; ++k) { +- if (a[k] <= key) { +- goto L450; +- } +- ha = a[mid]; +- a[mid] = a[k]; +- a[k] = ha; +- hi = irn[mid]; +- irn[mid] = irn[k]; +- irn[k] = hi; +- ++mid; +-L450: +- ; +- } +-/* Both subintervals [FIRST,MID), [MID,LAST) are nonempty */ +-/* Stack the longest of the two subintervals first */ +- if (mid - first >= last - mid) { +- todo[td + 1] = last; +- todo[td] = mid; +- todo[td - 1] = mid; +-/* TODO(TD-1) = FIRST */ +- } else { +- todo[td + 1] = mid; +- todo[td] = first; +- todo[td - 1] = last; +- todo[td - 2] = mid; +- } +- td += 2; +-L425: +- if (td == 0) { +- goto L400; +- } +-/* There is still work to be done */ +- if (todo[td - 1] - todo[td - 2] >= 15) { +- goto L500; +- } +-/* Next interval is already short enough for straightforward insertion */ +- td += -2; +- goto L425; +-/* Complete sorting with straightforward insertion */ +-L400: +- i__2 = ipj + len - 1; +- for (r__ = ipj + 1; r__ <= i__2; ++r__) { +- if (a[r__ - 1] < a[r__]) { +- ha = a[r__]; +- hi = irn[r__]; +- a[r__] = a[r__ - 1]; +- irn[r__] = irn[r__ - 1]; +- i__3 = ipj + 1; +- for (s = r__ - 1; s >= i__3; --s) { +- if (a[s - 1] < ha) { +- a[s] = a[s - 1]; +- irn[s] = irn[s - 1]; +- } else { +- a[s] = ha; +- irn[s] = hi; +- goto L200; +- } +-/* L300: */ +- } +- a[ipj] = ha; +- irn[ipj] = hi; +- } +-L200: +- ; +- } +-L100: +- ; +- } +- return 0; +-} /* mc64rd_ */ +- +-/* ********************************************************************** */ +-/* Subroutine */ int_t mc64sd_(int_t *n, int_t *ne, int_t *ip, int_t * +- irn, double *a, int_t *iperm, int_t *numx, int_t *w, +- int_t *len, int_t *lenl, int_t *lenh, int_t *fc, int_t *iw, +- int_t *iw4) +void mc64id_(int *a) -+{ + { +- /* System generated locals */ +- int_t i__1, i__2, i__3, i__4; +- +- /* Local variables */ +- int_t i__, j, k, l, ii, mod, cnt, num; +- double bval, bmin, bmax, rinf; +- int_t nval, wlen, idum1, idum2, idum3; +- extern /* Subroutine */ int_t mc64qd_(int_t *, int_t *, int_t *, +- int_t *, int_t *, double *, int_t *, double *), +- mc64ud_(int_t *, int_t *, int_t *, int_t *, int_t *, +- int_t *, int_t *, int_t *, int_t *, int_t *, int_t *, +- int_t *, int_t *, int_t *, int_t *); +- +-/* *** Copyright (c) 1999 Council for the Central Laboratory of the */ +-/* Research Councils *** */ +-/* *** Although every effort has been made to ensure robustness and *** */ +-/* *** reliability of the subroutines in this MC64 suite, we *** */ +-/* *** disclaim any liability arising through the use or misuse of *** */ +-/* *** any of the subroutines. *** */ +-/* *** Any problems? Contact ... */ +-/* Iain Duff (I.Duff@rl.ac.uk) or Jacko Koster (jak@ii.uib.no) *** */ +- +-/* N, NE, IP, IRN, are described in MC64A/AD. */ +-/* A is a REAL (DOUBLE PRECISION in the D-version) array of length NE. */ +-/* A(K), K=1..NE, must be set to the value of the entry that */ +-/* corresponds to IRN(k). The entries in each column must be */ +-/* non-negative and ordered by decreasing value. */ +-/* IPERM is an INT_T array of length N. On exit, it contains the */ +-/* bottleneck matching: IPERM(I) - 0 or row I is matched to column */ +-/* IPERM(I). */ +-/* NUMX is an INT_T variable. On exit, it contains the cardinality */ +-/* of the matching stored in IPERM. */ +-/* IW is an INT_T work array of length 10N. */ +-/* FC is an int_t array of length N that contains the list of */ +-/* unmatched columns. */ +-/* LEN(J), LENL(J), LENH(J) are int_t arrays of length N that point */ +-/* to entries in matrix column J. */ +-/* In the matrix defined by the column parts IP(J)+LENL(J) we know */ +-/* a matching does not exist; in the matrix defined by the column */ +-/* parts IP(J)+LENH(J) we know one exists. */ +-/* LEN(J) lies between LENL(J) and LENH(J) and determines the matrix */ +-/* that is tested for a maximum matching. */ +-/* W is an int_t array of length N and contains the indices of the */ +-/* columns for which LENL ne LENH. */ +-/* WLEN is number of indices stored in array W. */ +-/* IW is int_t work array of length N. */ +-/* IW4 is int_t work array of length 4N used by MC64U/UD. */ +-/* EXTERNAL FD05AD,MC64QD,MC64UD */ +-/* DOUBLE PRECISION FD05AD */ +-/* BMIN and BMAX are such that a maximum matching exists for the input */ +-/* matrix in which all entries smaller than BMIN are dropped. */ +-/* For BMAX, a maximum matching does not exist. */ +-/* BVAL is a value between BMIN and BMAX. */ +-/* CNT is the number of calls made to MC64U/UD so far. */ +-/* NUM is the cardinality of last matching found. */ +-/* Set RINF to largest positive real number */ +-/* XSL RINF = FD05AD(5) */ +- /* Parameter adjustments */ +- --iw4; +- --iw; +- --fc; +- --lenh; +- --lenl; +- --len; +- --w; +- --iperm; +- --ip; +- --a; +- --irn; +- +- /* Function Body */ +- rinf = dmach("Overflow"); +-/* Compute a first maximum matching from scratch on whole matrix. */ +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +- fc[j] = j; +- iw[j] = 0; +- len[j] = ip[j + 1] - ip[j]; +-/* L20: */ +- } +-/* The first call to MC64U/UD */ +- cnt = 1; +- mod = 1; +- *numx = 0; +- mc64ud_(&cnt, &mod, n, &irn[1], ne, &ip[1], &len[1], &fc[1], &iw[1], numx, +- n, &iw4[1], &iw4[*n + 1], &iw4[(*n << 1) + 1], &iw4[*n * 3 + 1]); +-/* IW contains a maximum matching of length NUMX. */ +- num = *numx; +- if (num != *n) { +-/* Matrix is structurally singular */ +- bmax = rinf; +- } else { +-/* Matrix is structurally nonsingular, NUM=NUMX=N; */ +-/* Set BMAX just above the smallest of all the maximum absolute */ +-/* values of the columns */ +- bmax = rinf; +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +- bval = 0.f; +- i__2 = ip[j + 1] - 1; +- for (k = ip[j]; k <= i__2; ++k) { +- if (a[k] > bval) { +- bval = a[k]; +- } +-/* L25: */ +- } +- if (bval < bmax) { +- bmax = bval; +- } +-/* L30: */ +- } +- bmax *= 1.001f; +- } +-/* Initialize BVAL,BMIN */ +- bval = 0.f; +- bmin = 0.f; +-/* Initialize LENL,LEN,LENH,W,WLEN according to BMAX. */ +-/* Set LEN(J), LENH(J) just after last entry in column J. */ +-/* Set LENL(J) just after last entry in column J with value ge BMAX. */ +- wlen = 0; +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +- l = ip[j + 1] - ip[j]; +- lenh[j] = l; +- len[j] = l; +- i__2 = ip[j + 1] - 1; +- for (k = ip[j]; k <= i__2; ++k) { +- if (a[k] < bmax) { +- goto L46; +- } +-/* L45: */ +- } +-/* Column J is empty or all entries are ge BMAX */ +- k = ip[j + 1]; +-L46: +- lenl[j] = k - ip[j]; +-/* Add J to W if LENL(J) ne LENH(J) */ +- if (lenl[j] == l) { +- goto L48; +- } +- ++wlen; +- w[wlen] = j; +-L48: +- ; +- } +-/* Main loop */ +- i__1 = *ne; +- for (idum1 = 1; idum1 <= i__1; ++idum1) { +- if (num == *numx) { +-/* We have a maximum matching in IW; store IW in IPERM */ +- i__2 = *n; +- for (i__ = 1; i__ <= i__2; ++i__) { +- iperm[i__] = iw[i__]; +-/* L50: */ +- } +-/* Keep going round this loop until matching IW is no longer maximum. */ +- i__2 = *ne; +- for (idum2 = 1; idum2 <= i__2; ++idum2) { +- bmin = bval; +- if (bmax == bmin) { +- goto L99; +- } +-/* Find splitting value BVAL */ +- mc64qd_(&ip[1], &lenl[1], &len[1], &w[1], &wlen, &a[1], &nval, +- &bval); +- if (nval <= 1) { +- goto L99; +- } +-/* Set LEN such that all matrix entries with value lt BVAL are */ +-/* discarded. Store old LEN in LENH. Do this for all columns W(K). */ +-/* Each step, either K is incremented or WLEN is decremented. */ +- k = 1; +- i__3 = *n; +- for (idum3 = 1; idum3 <= i__3; ++idum3) { +- if (k > wlen) { +- goto L71; +- } +- j = w[k]; +- i__4 = ip[j] + lenl[j]; +- for (ii = ip[j] + len[j] - 1; ii >= i__4; --ii) { +- if (a[ii] >= bval) { +- goto L60; +- } +- i__ = irn[ii]; +- if (iw[i__] != j) { +- goto L55; +- } +-/* Remove entry from matching */ +- iw[i__] = 0; +- --num; +- fc[*n - num] = j; +-L55: +- ; +- } +-L60: +- lenh[j] = len[j]; +-/* IP(J)+LEN(J)-1 is last entry in column ge BVAL */ +- len[j] = ii - ip[j] + 1; +-/* If LENH(J) = LENL(J), remove J from W */ +- if (lenl[j] == lenh[j]) { +- w[k] = w[wlen]; +- --wlen; +- } else { +- ++k; +- } +-/* L70: */ +- } +-L71: +- if (num < *numx) { +- goto L81; +- } +-/* L80: */ +- } +-/* End of dummy loop; this point is never reached */ +-/* Set mode for next call to MC64U/UD */ +-L81: +- mod = 1; +- } else { +-/* We do not have a maximum matching in IW. */ +- bmax = bval; +-/* BMIN is the bottleneck value of a maximum matching; */ +-/* for BMAX the matching is not maximum, so BMAX>BMIN */ +-/* IF (BMAX .EQ. BMIN) GO TO 99 */ +-/* Find splitting value BVAL */ +- mc64qd_(&ip[1], &len[1], &lenh[1], &w[1], &wlen, &a[1], &nval, & +- bval); +- if (nval == 0 || bval == bmin) { +- goto L99; +- } +-/* Set LEN such that all matrix entries with value ge BVAL are */ +-/* inside matrix. Store old LEN in LENL. Do this for all columns W(K). */ +-/* Each step, either K is incremented or WLEN is decremented. */ +- k = 1; +- i__2 = *n; +- for (idum3 = 1; idum3 <= i__2; ++idum3) { +- if (k > wlen) { +- goto L88; +- } +- j = w[k]; +- i__3 = ip[j] + lenh[j] - 1; +- for (ii = ip[j] + len[j]; ii <= i__3; ++ii) { +- if (a[ii] < bval) { +- goto L86; +- } +-/* L85: */ +- } +-L86: +- lenl[j] = len[j]; +- len[j] = ii - ip[j]; +- if (lenl[j] == lenh[j]) { +- w[k] = w[wlen]; +- --wlen; +- } else { +- ++k; +- } +-/* L87: */ +- } +-/* End of dummy loop; this point is never reached */ +-/* Set mode for next call to MC64U/UD */ +-L88: +- mod = 0; +- } +- ++cnt; +- mc64ud_(&cnt, &mod, n, &irn[1], ne, &ip[1], &len[1], &fc[1], &iw[1], & +- num, numx, &iw4[1], &iw4[*n + 1], &iw4[(*n << 1) + 1], &iw4[* +- n * 3 + 1]); +-/* IW contains maximum matching of length NUM */ +-/* L90: */ +- } +-/* End of dummy loop; this point is never reached */ +-/* BMIN is bottleneck value of final matching */ +-L99: +- if (*numx == *n) { +- goto L1000; +- } +-/* The matrix is structurally singular, complete IPERM */ +-/* W, IW are work arrays */ +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +- w[j] = 0; +-/* L300: */ +- } +- k = 0; +- i__1 = *n; +- for (i__ = 1; i__ <= i__1; ++i__) { +- if (iperm[i__] == 0) { +- ++k; +- iw[k] = i__; +- } else { +- j = iperm[i__]; +- w[j] = i__; +- } +-/* L310: */ +- } +- k = 0; +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +- if (w[j] != 0) { +- goto L320; +- } +- ++k; +- idum1 = iw[k]; +- iperm[idum1] = j; +-L320: +- ; +- } +-L1000: +- return 0; +-} /* mc64sd_ */ +- +-/* ********************************************************************** */ +-/* Subroutine */ int_t mc64qd_(int_t *ip, int_t *lenl, int_t *lenh, +- int_t *w, int_t *wlen, double *a, int_t *nval, double * +- val) +-{ +- /* System generated locals */ +- int_t i__1, i__2, i__3; +- +- /* Local variables */ +- int_t j, k, s; +- double ha; +- int_t ii, pos; +- double split[10]; +- +- +-/* *** Copyright (c) 1999 Council for the Central Laboratory of the */ +-/* Research Councils *** */ +-/* *** Although every effort has been made to ensure robustness and *** */ +-/* *** reliability of the subroutines in this MC64 suite, we *** */ +-/* *** disclaim any liability arising through the use or misuse of *** */ +-/* *** any of the subroutines. *** */ +-/* *** Any problems? Contact ... */ +-/* Iain Duff (I.Duff@rl.ac.uk) or Jacko Koster (jak@ii.uib.no) *** */ +- +-/* This routine searches for at most XX different numerical values */ +-/* in the columns W(1:WLEN). XX>=2. */ +-/* Each column J is scanned between IP(J)+LENL(J) and IP(J)+LENH(J)-1 */ +-/* until XX values are found or all columns have been considered. */ +-/* On output, NVAL is the number of different values that is found */ +-/* and SPLIT(1:NVAL) contains the values in decreasing order. */ +-/* If NVAL > 0, the routine returns VAL = SPLIT((NVAL+1)/2). */ +- +-/* Scan columns in W(1:WLEN). For each encountered value, if value not */ +-/* already present in SPLIT(1:NVAL), insert value such that SPLIT */ +-/* remains sorted by decreasing value. */ +-/* The sorting is done by straightforward insertion; therefore the use */ +-/* of this routine should be avoided for large XX (XX < 20). */ +- /* Parameter adjustments */ +- --a; +- --w; +- --lenh; +- --lenl; +- --ip; +- +- /* Function Body */ +- *nval = 0; +- i__1 = *wlen; +- for (k = 1; k <= i__1; ++k) { +- j = w[k]; +- i__2 = ip[j] + lenh[j] - 1; +- for (ii = ip[j] + lenl[j]; ii <= i__2; ++ii) { +- ha = a[ii]; +- if (*nval == 0) { +- split[0] = ha; +- *nval = 1; +- } else { +-/* Check presence of HA in SPLIT */ +- for (s = *nval; s >= 1; --s) { +- if (split[s - 1] == ha) { +- goto L15; +- } +- if (split[s - 1] > ha) { +- pos = s + 1; +- goto L21; +- } +-/* L20: */ +- } +- pos = 1; +-/* The insertion */ +-L21: +- i__3 = pos; +- for (s = *nval; s >= i__3; --s) { +- split[s] = split[s - 1]; +-/* L22: */ +- } +- split[pos - 1] = ha; +- ++(*nval); +- } +-/* Exit loop if XX values are found */ +- if (*nval == 10) { +- goto L11; +- } +-L15: +- ; +- } +-/* L10: */ +- } +-/* Determine VAL */ +-L11: +- if (*nval > 0) { +- *val = split[(*nval + 1) / 2 - 1]; +- } +- return 0; +-} /* mc64qd_ */ +- +-/* ********************************************************************** */ +-/* Subroutine */ int_t mc64ud_(int_t *id, int_t *mod, int_t *n, int_t * +- irn, int_t *lirn, int_t *ip, int_t *lenc, int_t *fc, int_t * +- iperm, int_t *num, int_t *numx, int_t *pr, int_t *arp, +- int_t *cv, int_t *out) +-{ +- /* System generated locals */ +- int_t i__1, i__2, i__3, i__4; +- +- /* Local variables */ +- int_t i__, j, k, j1, ii, kk, id0, id1, in1, in2, nfc, num0, num1, num2, +- jord, last; +- +- +-/* *** Copyright (c) 1999 Council for the Central Laboratory of the */ +-/* Research Councils *** */ +-/* *** Although every effort has been made to ensure robustness and *** */ +-/* *** reliability of the subroutines in this MC64 suite, we *** */ +-/* *** disclaim any liability arising through the use or misuse of *** */ +-/* *** any of the subroutines. *** */ +-/* *** Any problems? Contact ... */ +-/* Iain Duff (I.Duff@rl.ac.uk) or Jacko Koster (jak@ii.uib.no) *** */ +- +-/* PR(J) is the previous column to J in the depth first search. */ +-/* Array PR is used as workspace in the sorting algorithm. */ +-/* Elements (I,IPERM(I)) I=1,..,N are entries at the end of the */ +-/* algorithm unless N assignments have not been made in which case */ +-/* N-NUM pairs (I,IPERM(I)) will not be entries in the matrix. */ +-/* CV(I) is the most recent loop number (ID+JORD) at which row I */ +-/* was visited. */ +-/* ARP(J) is the number of entries in column J which have been scanned */ +-/* when looking for a cheap assignment. */ +-/* OUT(J) is one less than the number of entries in column J which have */ +-/* not been scanned during one pass through the main loop. */ +-/* NUMX is maximum possible size of matching. */ +- /* Parameter adjustments */ +- --out; +- --cv; +- --arp; +- --pr; +- --iperm; +- --fc; +- --lenc; +- --ip; +- --irn; +- +- /* Function Body */ +- if (*id == 1) { +-/* The first call to MC64U/UD. */ +-/* Initialize CV and ARP; parameters MOD, NUMX are not accessed */ +- i__1 = *n; +- for (i__ = 1; i__ <= i__1; ++i__) { +- cv[i__] = 0; +- arp[i__] = 0; +-/* L5: */ +- } +- num1 = *n; +- num2 = *n; +- } else { +-/* Not the first call to MC64U/UD. */ +-/* Re-initialize ARP if entries were deleted since last call to MC64U/UD */ +- if (*mod == 1) { +- i__1 = *n; +- for (i__ = 1; i__ <= i__1; ++i__) { +- arp[i__] = 0; +-/* L8: */ +- } +- } +- num1 = *numx; +- num2 = *n - *numx; +- } +- num0 = *num; +-/* NUM0 is size of input matching */ +-/* NUM1 is maximum possible size of matching */ +-/* NUM2 is maximum allowed number of unassigned rows/columns */ +-/* NUM is size of current matching */ +-/* Quick return if possible */ +-/* IF (NUM.EQ.N) GO TO 199 */ +-/* NFC is number of rows/columns that could not be assigned */ +- nfc = 0; +-/* Integers ID0+1 to ID0+N are unique numbers for call ID to MC64U/UD, */ +-/* so 1st call uses 1..N, 2nd call uses N+1..2N, etc */ +- id0 = (*id - 1) * *n; +-/* Main loop. Each pass round this loop either results in a new */ +-/* assignment or gives a column with no assignment */ +- i__1 = *n; +- for (jord = num0 + 1; jord <= i__1; ++jord) { +-/* Each pass uses unique number ID1 */ +- id1 = id0 + jord; +-/* J is unmatched column */ +- j = fc[jord - num0]; +- pr[j] = -1; +- i__2 = jord; +- for (k = 1; k <= i__2; ++k) { +-/* Look for a cheap assignment */ +- if (arp[j] >= lenc[j]) { +- goto L30; +- } +- in1 = ip[j] + arp[j]; +- in2 = ip[j] + lenc[j] - 1; +- i__3 = in2; +- for (ii = in1; ii <= i__3; ++ii) { +- i__ = irn[ii]; +- if (iperm[i__] == 0) { +- goto L80; +- } +-/* L20: */ +- } +-/* No cheap assignment in row */ +- arp[j] = lenc[j]; +-/* Begin looking for assignment chain starting with row J */ +-L30: +- out[j] = lenc[j] - 1; +-/* Inner loop. Extends chain by one or backtracks */ +- i__3 = jord; +- for (kk = 1; kk <= i__3; ++kk) { +- in1 = out[j]; +- if (in1 < 0) { +- goto L50; +- } +- in2 = ip[j] + lenc[j] - 1; +- in1 = in2 - in1; +-/* Forward scan */ +- i__4 = in2; +- for (ii = in1; ii <= i__4; ++ii) { +- i__ = irn[ii]; +- if (cv[i__] == id1) { +- goto L40; +- } +-/* Column J has not yet been accessed during this pass */ +- j1 = j; +- j = iperm[i__]; +- cv[i__] = id1; +- pr[j] = j1; +- out[j1] = in2 - ii - 1; +- goto L70; +-L40: +- ; +- } +-/* Backtracking step. */ +-L50: +- j1 = pr[j]; +- if (j1 == -1) { +-/* No augmenting path exists for column J. */ +- ++nfc; +- fc[nfc] = j; +- if (nfc > num2) { +-/* A matching of maximum size NUM1 is not possible */ +- last = jord; +- goto L101; +- } +- goto L100; +- } +- j = j1; +-/* L60: */ +- } +-/* End of dummy loop; this point is never reached */ +-L70: +- ; +- } +-/* End of dummy loop; this point is never reached */ +-/* New assignment is made. */ +-L80: +- iperm[i__] = j; +- arp[j] = ii - ip[j] + 1; +- ++(*num); +- i__2 = jord; +- for (k = 1; k <= i__2; ++k) { +- j = pr[j]; +- if (j == -1) { +- goto L95; +- } +- ii = ip[j] + lenc[j] - out[j] - 2; +- i__ = irn[ii]; +- iperm[i__] = j; +-/* L90: */ +- } +-/* End of dummy loop; this point is never reached */ +-L95: +- if (*num == num1) { +-/* A matching of maximum size NUM1 is found */ +- last = jord; +- goto L101; +- } +- +-L100: +- ; +- } +-/* All unassigned columns have been considered */ +- last = *n; +-/* Now, a transversal is computed or is not possible. */ +-/* Complete FC before returning. */ +-L101: +- i__1 = *n; +- for (jord = last + 1; jord <= i__1; ++jord) { +- ++nfc; +- fc[nfc] = fc[jord - num0]; +-/* L110: */ +- } +-/* 199 RETURN */ +- return 0; +-} /* mc64ud_ */ +- +-/* ********************************************************************** */ +-/* Subroutine */ int_t mc64wd_(int_t *n, int_t *ne, int_t *ip, int_t * +- irn, double *a, int_t *iperm, int_t *num, int_t *jperm, +- int_t *out, int_t *pr, int_t *q, int_t *l, double *u, +- double *d__) +-{ +- /* System generated locals */ +- int_t i__1, i__2, i__3, c__2 = 2; +- +- /* Local variables */ +- int_t i__, j, k, i0, k0, k1, k2, q0; +- double di; +- int_t ii, jj, kk; +- double vj; +- int_t up; +- double dq0; +- int_t kk1, kk2; +- double csp; +- int_t isp, jsp, low; +- double dmin__, dnew; +- int_t jord, qlen, jdum; +- double rinf; +- extern /* Subroutine */ int_t mc64dd_(int_t *, int_t *, int_t *, +- double *, int_t *, int_t *), mc64ed_(int_t *, int_t *, +- int_t *, double *, int_t *, int_t *), mc64fd_(int_t * +- , int_t *, int_t *, int_t *, double *, int_t *, +- int_t *); +- +- +-/* *** Copyright (c) 1999 Council for the Central Laboratory of the */ +-/* Research Councils *** */ +-/* *** Although every effort has been made to ensure robustness and *** */ +-/* *** reliability of the subroutines in this MC64 suite, we *** */ +-/* *** disclaim any liability arising through the use or misuse of *** */ +-/* *** any of the subroutines. *** */ +-/* *** Any problems? Contact ... */ +-/* Iain Duff (I.Duff@rl.ac.uk) or Jacko Koster (jak@ii.uib.no) *** */ +- +-/* N, NE, IP, IRN are described in MC64A/AD. */ +-/* A is a REAL (DOUBLE PRECISION in the D-version) array of length NE. */ +-/* A(K), K=1..NE, must be set to the value of the entry that */ +-/* corresponds to IRN(K). It is not altered. */ +-/* All values A(K) must be non-negative. */ +-/* IPERM is an INT_T array of length N. On exit, it contains the */ +-/* weighted matching: IPERM(I) = 0 or row I is matched to column */ +-/* IPERM(I). */ +-/* NUM is an INT_T variable. On exit, it contains the cardinality of */ +-/* the matching stored in IPERM. */ +-/* IW is an INT_T work array of length 5N. */ +-/* DW is a REAL (DOUBLE PRECISION in the D-version) array of length 2N. */ +-/* On exit, U = D(1:N) contains the dual row variable and */ +-/* V = D(N+1:2N) contains the dual column variable. If the matrix */ +-/* is structurally nonsingular (NUM = N), the following holds: */ +-/* U(I)+V(J) <= A(I,J) if IPERM(I) |= J */ +-/* U(I)+V(J) = A(I,J) if IPERM(I) = J */ +-/* U(I) = 0 if IPERM(I) = 0 */ +-/* V(J) = 0 if there is no I for which IPERM(I) = J */ +-/* Local variables */ +-/* Local parameters */ +-/* External subroutines and/or functions */ +-/* EXTERNAL FD05AD,MC64DD,MC64ED,MC64FD */ +-/* DOUBLE PRECISION FD05AD */ +-/* Set RINF to largest positive real number */ +-/* XSL RINF = FD05AD(5) */ +- /* Parameter adjustments */ +- --d__; +- --u; +- --l; +- --q; +- --pr; +- --out; +- --jperm; +- --iperm; +- --ip; +- --a; +- --irn; +- +- /* Function Body */ +- rinf = dmach("Overflow"); +-/* Initialization */ +- *num = 0; +- i__1 = *n; +- for (k = 1; k <= i__1; ++k) { +- u[k] = rinf; +- d__[k] = 0.; +- iperm[k] = 0; +- jperm[k] = 0; +- pr[k] = ip[k]; +- l[k] = 0; +-/* L10: */ +- } +-/* Initialize U(I) */ +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +- i__2 = ip[j + 1] - 1; +- for (k = ip[j]; k <= i__2; ++k) { +- i__ = irn[k]; +- if (a[k] > u[i__]) { +- goto L20; +- } +- u[i__] = a[k]; +- iperm[i__] = j; +- l[i__] = k; +-L20: +- ; +- } +-/* L30: */ +- } +- i__1 = *n; +- for (i__ = 1; i__ <= i__1; ++i__) { +- j = iperm[i__]; +- if (j == 0) { +- goto L40; +- } +-/* Row I is not empty */ +- iperm[i__] = 0; +- if (jperm[j] != 0) { +- goto L40; +- } +-/* Assignment of column J to row I */ +- ++(*num); +- iperm[i__] = j; +- jperm[j] = l[i__]; +-L40: +- ; +- } +- if (*num == *n) { +- goto L1000; +- } +-/* Scan unassigned columns; improve assignment */ +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +-/* JPERM(J) ne 0 iff column J is already assigned */ +- if (jperm[j] != 0) { +- goto L95; +- } +- k1 = ip[j]; +- k2 = ip[j + 1] - 1; +-/* Continue only if column J is not empty */ +- if (k1 > k2) { +- goto L95; +- } +- vj = rinf; +- i__2 = k2; +- for (k = k1; k <= i__2; ++k) { +- i__ = irn[k]; +- di = a[k] - u[i__]; +- if (di > vj) { +- goto L50; +- } +- if (di < vj || di == rinf) { +- goto L55; +- } +- if (iperm[i__] != 0 || iperm[i0] == 0) { +- goto L50; +- } +-L55: +- vj = di; +- i0 = i__; +- k0 = k; +-L50: +- ; +- } +- d__[j] = vj; +- k = k0; +- i__ = i0; +- if (iperm[i__] == 0) { +- goto L90; +- } +- i__2 = k2; +- for (k = k0; k <= i__2; ++k) { +- i__ = irn[k]; +- if (a[k] - u[i__] > vj) { +- goto L60; +- } +- jj = iperm[i__]; +-/* Scan remaining part of assigned column JJ */ +- kk1 = pr[jj]; +- kk2 = ip[jj + 1] - 1; +- if (kk1 > kk2) { +- goto L60; +- } +- i__3 = kk2; +- for (kk = kk1; kk <= i__3; ++kk) { +- ii = irn[kk]; +- if (iperm[ii] > 0) { +- goto L70; +- } +- if (a[kk] - u[ii] <= d__[jj]) { +- goto L80; +- } +-L70: +- ; +- } +- pr[jj] = kk2 + 1; +-L60: +- ; +- } +- goto L95; +-L80: +- jperm[jj] = kk; +- iperm[ii] = jj; +- pr[jj] = kk + 1; +-L90: +- ++(*num); +- jperm[j] = k; +- iperm[i__] = j; +- pr[j] = k + 1; +-L95: +- ; +- } +- if (*num == *n) { +- goto L1000; +- } +-/* Prepare for main loop */ +- i__1 = *n; +- for (i__ = 1; i__ <= i__1; ++i__) { +- d__[i__] = rinf; +- l[i__] = 0; +-/* L99: */ +- } +-/* Main loop ... each pass round this loop is similar to Dijkstra's */ +-/* algorithm for solving the single source shortest path problem */ +- i__1 = *n; +- for (jord = 1; jord <= i__1; ++jord) { +- if (jperm[jord] != 0) { +- goto L100; +- } +-/* JORD is next unmatched column */ +-/* DMIN is the length of shortest path in the tree */ +- dmin__ = rinf; +- qlen = 0; +- low = *n + 1; +- up = *n + 1; +-/* CSP is the cost of the shortest augmenting path to unassigned row */ +-/* IRN(ISP). The corresponding column index is JSP. */ +- csp = rinf; +-/* Build shortest path tree starting from unassigned column (root) JORD */ +- j = jord; +- pr[j] = -1; +-/* Scan column J */ +- i__2 = ip[j + 1] - 1; +- for (k = ip[j]; k <= i__2; ++k) { +- i__ = irn[k]; +- dnew = a[k] - u[i__]; +- if (dnew >= csp) { +- goto L115; +- } +- if (iperm[i__] == 0) { +- csp = dnew; +- isp = k; +- jsp = j; +- } else { +- if (dnew < dmin__) { +- dmin__ = dnew; +- } +- d__[i__] = dnew; +- ++qlen; +- q[qlen] = k; +- } +-L115: +- ; +- } +-/* Initialize heap Q and Q2 with rows held in Q(1:QLEN) */ +- q0 = qlen; +- qlen = 0; +- i__2 = q0; +- for (kk = 1; kk <= i__2; ++kk) { +- k = q[kk]; +- i__ = irn[k]; +- if (csp <= d__[i__]) { +- d__[i__] = rinf; +- goto L120; +- } +- if (d__[i__] <= dmin__) { +- --low; +- q[low] = i__; +- l[i__] = low; +- } else { +- ++qlen; +- l[i__] = qlen; +- mc64dd_(&i__, n, &q[1], &d__[1], &l[1], &c__2); +- } +-/* Update tree */ +- jj = iperm[i__]; +- out[jj] = k; +- pr[jj] = j; +-L120: +- ; +- } +- i__2 = *num; +- for (jdum = 1; jdum <= i__2; ++jdum) { +-/* If Q2 is empty, extract rows from Q */ +- if (low == up) { +- if (qlen == 0) { +- goto L160; +- } +- i__ = q[1]; +- if (d__[i__] >= csp) { +- goto L160; +- } +- dmin__ = d__[i__]; +-L152: +- mc64ed_(&qlen, n, &q[1], &d__[1], &l[1], &c__2); +- --low; +- q[low] = i__; +- l[i__] = low; +- if (qlen == 0) { +- goto L153; +- } +- i__ = q[1]; +- if (d__[i__] > dmin__) { +- goto L153; +- } +- goto L152; +- } +-/* Q0 is row whose distance D(Q0) to the root is smallest */ +-L153: +- q0 = q[up - 1]; +- dq0 = d__[q0]; +-/* Exit loop if path to Q0 is longer than the shortest augmenting path */ +- if (dq0 >= csp) { +- goto L160; +- } +- --up; +-/* Scan column that matches with row Q0 */ +- j = iperm[q0]; +- vj = dq0 - a[jperm[j]] + u[q0]; +- i__3 = ip[j + 1] - 1; +- for (k = ip[j]; k <= i__3; ++k) { +- i__ = irn[k]; +- if (l[i__] >= up) { +- goto L155; +- } +-/* DNEW is new cost */ +- dnew = vj + a[k] - u[i__]; +-/* Do not update D(I) if DNEW ge cost of shortest path */ +- if (dnew >= csp) { +- goto L155; +- } +- if (iperm[i__] == 0) { +-/* Row I is unmatched; update shortest path info */ +- csp = dnew; +- isp = k; +- jsp = j; +- } else { +-/* Row I is matched; do not update D(I) if DNEW is larger */ +- di = d__[i__]; +- if (di <= dnew) { +- goto L155; +- } +- if (l[i__] >= low) { +- goto L155; +- } +- d__[i__] = dnew; +- if (dnew <= dmin__) { +- if (l[i__] != 0) { +- mc64fd_(&l[i__], &qlen, n, &q[1], &d__[1], &l[1], +- &c__2); +- } +- --low; +- q[low] = i__; +- l[i__] = low; +- } else { +- if (l[i__] == 0) { +- ++qlen; +- l[i__] = qlen; +- } +- mc64dd_(&i__, n, &q[1], &d__[1], &l[1], &c__2); +- } +-/* Update tree */ +- jj = iperm[i__]; +- out[jj] = k; +- pr[jj] = j; +- } +-L155: +- ; +- } +-/* L150: */ +- } +-/* If CSP = RINF, no augmenting path is found */ +-L160: +- if (csp == rinf) { +- goto L190; +- } +-/* Find augmenting path by tracing backward in PR; update IPERM,JPERM */ +- ++(*num); +- i__ = irn[isp]; +- iperm[i__] = jsp; +- jperm[jsp] = isp; +- j = jsp; +- i__2 = *num; +- for (jdum = 1; jdum <= i__2; ++jdum) { +- jj = pr[j]; +- if (jj == -1) { +- goto L180; +- } +- k = out[j]; +- i__ = irn[k]; +- iperm[i__] = jj; +- jperm[jj] = k; +- j = jj; +-/* L170: */ +- } +-/* End of dummy loop; this point is never reached */ +-/* Update U for rows in Q(UP:N) */ +-L180: +- i__2 = *n; +- for (kk = up; kk <= i__2; ++kk) { +- i__ = q[kk]; +- u[i__] = u[i__] + d__[i__] - csp; +-/* L185: */ +- } +-L190: +- i__2 = *n; +- for (kk = low; kk <= i__2; ++kk) { +- i__ = q[kk]; +- d__[i__] = rinf; +- l[i__] = 0; +-/* L191: */ +- } +- i__2 = qlen; +- for (kk = 1; kk <= i__2; ++kk) { +- i__ = q[kk]; +- d__[i__] = rinf; +- l[i__] = 0; +-/* L193: */ +- } +-L100: +- ; +- } +-/* End of main loop */ +-/* Set dual column variable in D(1:N) */ +-L1000: +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +- k = jperm[j]; +- if (k != 0) { +- d__[j] = a[k] - u[irn[k]]; +- } else { +- d__[j] = 0.; +- } +- if (iperm[j] == 0) { +- u[j] = 0.; +- } +-/* L200: */ +- } +- if (*num == *n) { +- goto L1100; +- } +-/* The matrix is structurally singular, complete IPERM. */ +-/* JPERM, OUT are work arrays */ +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +- jperm[j] = 0; +-/* L300: */ +- } +- k = 0; +- i__1 = *n; +- for (i__ = 1; i__ <= i__1; ++i__) { +- if (iperm[i__] == 0) { +- ++k; +- out[k] = i__; +- } else { +- j = iperm[i__]; +- jperm[j] = i__; +- } +-/* L310: */ +- } +- k = 0; +- i__1 = *n; +- for (j = 1; j <= i__1; ++j) { +- if (jperm[j] != 0) { +- goto L320; +- } +- ++k; +- jdum = out[k]; +- iperm[jdum] = j; +-L320: +- ; +- } +-L1100: +- return 0; +-} /* mc64wd_ */ +- +- + fprintf(stderr, "SuperLU: MC64 functionality not available (it uses non-free code). Aborting.\n"); + abort(); +}