From 0e9c6e15f1ca7817a4af86bf37bb260de8c84e19 Mon Sep 17 00:00:00 2001 From: Troy Dawson Date: Thu, 14 Apr 2022 13:16:42 -0700 Subject: [PATCH] SuperLU package is retired on c9s for CS-704 --- .gitignore | 6 - SuperLU-5.2.2-set_soname.patch | 11 - SuperLU-fedora-make.inc.in | 48 - SuperLU-fix_example_builds.patch | 54 - SuperLU.spec | 296 ---- dead.package | 1 + sources | 1 - superlu-removemc64.patch | 2662 ------------------------------ 8 files changed, 1 insertion(+), 3078 deletions(-) delete mode 100644 .gitignore delete mode 100644 SuperLU-5.2.2-set_soname.patch delete mode 100644 SuperLU-fedora-make.inc.in delete mode 100644 SuperLU-fix_example_builds.patch delete mode 100644 SuperLU.spec create mode 100644 dead.package delete mode 100644 sources delete mode 100644 superlu-removemc64.patch diff --git a/.gitignore b/.gitignore deleted file mode 100644 index adb7e63..0000000 --- a/.gitignore +++ /dev/null @@ -1,6 +0,0 @@ -/superlu_4.3.tar.gz -/superlu_5.1.1.tar.gz -/superlu_5.2.0.tar.gz -/superlu_5.2.1.tar.gz -/superlu-5.2.1.tar.gz -/superlu-5.2.2.tar.gz diff --git a/SuperLU-5.2.2-set_soname.patch b/SuperLU-5.2.2-set_soname.patch deleted file mode 100644 index dde86ae..0000000 --- a/SuperLU-5.2.2-set_soname.patch +++ /dev/null @@ -1,11 +0,0 @@ ---- a/SRC/CMakeLists.orig.txt 2018-04-25 15:17:39.786715000 +0200 -+++ b/SRC/CMakeLists.txt 2018-04-25 15:25:00.253317607 +0200 -@@ -242,7 +242,7 @@ - $ - ) - set_target_properties(superlu PROPERTIES -- VERSION ${PROJECT_VERSION} SOVERSION ${VERSION_MAJOR} -+ VERSION ${PROJECT_VERSION} SOVERSION ${VERSION_MAJOR}.${VERSION_MINOR} - ) - - # Define GNU standard installation directories diff --git a/SuperLU-fedora-make.inc.in b/SuperLU-fedora-make.inc.in deleted file mode 100644 index 6b998c0..0000000 --- a/SuperLU-fedora-make.inc.in +++ /dev/null @@ -1,48 +0,0 @@ -############################################################################ -# -# Program: SuperLU -# -# Module: make.inc -# -# Purpose: Top-level Definitions -# -# Creation date: May 10, 2015 -# -# Modified: -# -# -############################################################################ -# -# The name of the libraries to be created/linked to -# -SuperLUroot = ${CMAKE_SOURCE_DIR}/build -# -# SuperLUroot = ${CMAKE_INSTALL_PREFIX} -SUPERLULIB = $(SuperLUroot)/SRC/${PROJECT_NAME_LIB_EXPORT} - -#TMGLIB = libtmglib.a - -# BLASDEF = -DUSE_VENDOR_BLAS -BLASLIB = $(TPL_BLAS_LIBRARIES) -LIBS = $(SUPERLULIB) $(BLASLIB) - -# -# The archiver and the flag(s) to use when building archive (library) -# If your system has no ranlib, set RANLIB = echo. -# -ARCH = @CMAKE_AR@ -ARCHFLAGS = cr -RANLIB = @CMAKE_RANLIB@ - -CC = @CMAKE_C_COMPILER@ -CFLAGS = @CMAKE_C_FLAGS_RELEASE@ @CMAKE_C_FLAGS@ -NOOPTS = -FORTRAN = @CMAKE_Fortran_COMPILER@ - -LOADER = $(CC) -LOADOPTS = $(LDFLAGS) - -# -# The directory in which Matlab is installed -# -MATLAB = /Applications/MATLAB_R2015b.app diff --git a/SuperLU-fix_example_builds.patch b/SuperLU-fix_example_builds.patch deleted file mode 100644 index aa1c4b7..0000000 --- a/SuperLU-fix_example_builds.patch +++ /dev/null @@ -1,54 +0,0 @@ ---- a/EXAMPLE/Makefile.orig 2020-10-18 00:23:21.000000000 +0200 -+++ b/EXAMPLE/Makefile 2021-01-29 17:58:26.396762602 +0100 -@@ -101,10 +101,10 @@ - $(LOADER) $(LOADOPTS) $(SLINXEXM3) $(LIBS) -lm -o $@ - - sitersol: $(SITSOL) $(SUPERLULIB) -- $(LOADER) $(LOADOPTS) $(SITSOL) $(LIBS) -lm -o $@ -+ $(LOADER) $(LOADOPTS) $(SITSOL) $(LIBS) $(TPL_BLAS_LIBRARIES) -lm -o $@ - - sitersol1: $(SITSOL1) $(SUPERLULIB) -- $(LOADER) $(LOADOPTS) $(SITSOL1) $(LIBS) -lm -o $@ -+ $(LOADER) $(LOADOPTS) $(SITSOL1) $(LIBS) $(TPL_BLAS_LIBRARIES) -lm -o $@ - - dlinsol: $(DLINEXM) $(SUPERLULIB) - $(LOADER) $(LOADOPTS) $(DLINEXM) $(LIBS) -lm -o $@ -@@ -128,10 +128,10 @@ - $(LOADER) $(LOADOPTS) $(SUPERLUEXM) $(LIBS) -lm -o $@ - - ditersol: $(DITSOL) $(SUPERLULIB) -- $(LOADER) $(LOADOPTS) $(DITSOL) $(LIBS) -lm -o $@ -+ $(LOADER) $(LOADOPTS) $(DITSOL) $(LIBS) $(TPL_BLAS_LIBRARIES) -lm -o $@ - - ditersol1: $(DITSOL1) $(SUPERLULIB) -- $(LOADER) $(LOADOPTS) $(DITSOL1) $(LIBS) -lm -o $@ -+ $(LOADER) $(LOADOPTS) $(DITSOL1) $(LIBS) $(TPL_BLAS_LIBRARIES) -lm -o $@ - - clinsol: $(CLINEXM) $(SUPERLULIB) - $(LOADER) $(LOADOPTS) $(CLINEXM) $(LIBS) -lm -o $@ -@@ -152,10 +152,10 @@ - $(LOADER) $(LOADOPTS) $(CLINXEXM3) $(LIBS) -lm -o $@ - - citersol: $(CITSOL) $(SUPERLULIB) -- $(LOADER) $(LOADOPTS) $(CITSOL) $(LIBS) -lm -o $@ -+ $(LOADER) $(LOADOPTS) $(CITSOL) $(LIBS) $(TPL_BLAS_LIBRARIES) -lm -o $@ - - citersol1: $(CITSOL1) $(SUPERLULIB) -- $(LOADER) $(LOADOPTS) $(CITSOL1) $(LIBS) -lm -o $@ -+ $(LOADER) $(LOADOPTS) $(CITSOL1) $(LIBS) $(TPL_BLAS_LIBRARIES) -lm -o $@ - - zlinsol: $(ZLINEXM) $(SUPERLULIB) - $(LOADER) $(LOADOPTS) $(ZLINEXM) $(LIBS) -lm -o $@ -@@ -176,10 +176,10 @@ - $(LOADER) $(LOADOPTS) $(ZLINXEXM3) $(LIBS) -lm -o $@ - - zitersol: $(ZITSOL) $(SUPERLULIB) -- $(LOADER) $(LOADOPTS) $(ZITSOL) $(LIBS) -lm -o $@ -+ $(LOADER) $(LOADOPTS) $(ZITSOL) $(LIBS) $(TPL_BLAS_LIBRARIES) -lm -o $@ - - zitersol1: $(ZITSOL1) $(SUPERLULIB) -- $(LOADER) $(LOADOPTS) $(ZITSOL1) $(LIBS) -lm -o $@ -+ $(LOADER) $(LOADOPTS) $(ZITSOL1) $(LIBS) $(TPL_BLAS_LIBRARIES) -lm -o $@ - - .c.o: - $(CC) $(CFLAGS) $(CDEFS) $(INCLUDEDIR) -c $< $(VERBOSE) diff --git a/SuperLU.spec b/SuperLU.spec deleted file mode 100644 index fd5d7d2..0000000 --- a/SuperLU.spec +++ /dev/null @@ -1,296 +0,0 @@ -%global genname superlu -%global libver 5 - -Name: SuperLU -Version: 5.2.2 -Release: 2%{?dist} -Summary: Subroutines to solve sparse linear systems -License: BSD and GPLv2+ -URL: https://portal.nersc.gov/project/sparse/superlu/ -Source0: https://github.com/xiaoyeli/superlu/archive/v%{version}/%{genname}-%{version}.tar.gz - -# Use a pre-made configuration file for Make -Source1: %{name}-fedora-make.inc.in - -Patch0: %{genname}-removemc64.patch - -# Patch soname (5 -> 5.2) of shared library -Patch1: %{name}-%{version}-set_soname.patch - -# Fix ldflags of example files -Patch2: %{name}-fix_example_builds.patch - -%if 0%{?fedora} >= 33 -BuildRequires: pkgconfig(flexiblas) -%else -%ifarch %{openblas_arches} -BuildRequires: openblas-devel, openblas-srpm-macros -%else -BuildRequires: blas-devel -%endif -BuildRequires: atlas-devel -%if 0%{?epel} -BuildRequires: epel-rpm-macros -%endif -%endif -BuildRequires: make -BuildRequires: cmake3 -BuildRequires: gcc, gcc-gfortran -BuildRequires: csh - -%description -SuperLU contains a set of subroutines to solve a sparse linear system -A*X=B. It uses Gaussian elimination with partial pivoting (GEPP). -The columns of A may be preordered before factorization; the -preordering for sparsity is completely separate from the factorization. - -%package devel -Summary: Header files and libraries for SuperLU development -Requires: %{name}%{?_isa} = %{version}-%{release} - -%description devel -The %{name}-devel package contains the header files -and libraries for use with %{name} package. - -%package doc -Summary: Documentation and Examples for SuperLU -BuildArch: noarch -%description doc -The %{name}-doc package contains all the help HTML documentation. - -%prep -%autosetup -n %{genname}-%{version} -p1 - -rm -f make.inc -cp -pf %{SOURCE1} make.inc.in - -# Remove bundled BLAS -rm -rf CBLAS - -rm -fr SRC/mc64ad.f.bak -find . -type f | sed -e "/TESTING/d" | xargs chmod a-x - -# Remove the shippped executables from EXAMPLE -find EXAMPLE -type f | while read file -do - [ "$(file $file | awk '{print $2}')" = ELF ] && rm $file || : -done - -# Change optimization level -sed -e 's|-O0|-O2|g' -i SRC/CMakeLists.txt - -%build -mkdir -p build -%cmake3 -S . -B build \ - -Denable_internal_blaslib:BOOL=NO \ -%if 0%{?fedora} >= 33 - -DTPL_BLAS_LIBRARIES="`pkg-config --libs flexiblas`" \ -%else - -DTPL_BLAS_LIBRARIES=-lopenblas \ -%endif - -DCMAKE_BUILD_TYPE:STRING=Release \ - -DCMAKE_INSTALL_INCLUDEDIR:PATH=include/%{name} \ - -DCMAKE_INSTALL_LIBDIR:PATH=%{_lib} \ - -DCMAKE_SKIP_RPATH:BOOL=YES -DCMAKE_SKIP_INSTALL_RPATH:BOOL=YES -%make_build -C build - -# Compile Fortran example interface to use the C routines in SuperLU -%if 0%{?fedora} >= 33 -export TPL_BLAS_LIBRARIES="`pkg-config --libs flexiblas`" \ -%else -export TPL_BLAS_LIBRARIES=-lopenblas \ -%endif -export LDFLAGS="%{build_ldflags}" -export FFLAGS="%{build_fflags}" -export CFLAGS="%{build_cflags}" -make -C FORTRAN -make -C EXAMPLE - -%install -%make_install -C build - -%check -pushd build -ctest3 -V %{?_smp_mflags} -popd - -# Test Fortran example interface to use the C routines in SuperLU -pushd FORTRAN -export LD_LIBRARY_PATH=%{buildroot}%{_libdir} -./df77exm < ../EXAMPLE/g20.rua -./zf77exm < ../EXAMPLE/cg20.cua -./test_omp < ../EXAMPLE/g20.rua -make clean -popd -# Test sample programs -pushd EXAMPLE -export LD_LIBRARY_PATH=%{buildroot}%{_libdir} -./superlu -make clean -popd - -%files -%license License.txt -%{_libdir}/libsuperlu.so.%{libver}* - -%files devel -%{_includedir}/%{name}/ -%{_libdir}/libsuperlu.so -%{_libdir}/cmake/%{genname}/ -%{_libdir}/pkgconfig/%{genname}.pc - -%files doc -%license License.txt -%doc DOC - -%changelog -* Thu Apr 15 2021 Mohan Boddu - 5.2.2-2 -- Rebuilt for RHEL 9 BETA on Apr 15th 2021. Related: rhbz#1947937 - -* Thu Jan 28 2021 Antonio Trande - 5.2.2-1 -- Release 5.2.2 - -* Mon Jan 25 2021 Fedora Release Engineering - 5.2.1-15 -- Rebuilt for https://fedoraproject.org/wiki/Fedora_34_Mass_Rebuild - -* Sat Aug 01 2020 Fedora Release Engineering - 5.2.1-14 -- Second attempt - Rebuilt for - https://fedoraproject.org/wiki/Fedora_33_Mass_Rebuild - -* Mon Jul 27 2020 Fedora Release Engineering - 5.2.1-13 -- Rebuilt for https://fedoraproject.org/wiki/Fedora_33_Mass_Rebuild - -* Sat Jul 25 2020 Iñaki Úcar - 5.2.1-12 -- https://fedoraproject.org/wiki/Changes/FlexiBLAS_as_BLAS/LAPACK_manager - -* Tue Jul 21 2020 Merlin Mathesius - 5.2.1-11 -- Minor conditional fix for ELN -- Stick to cmake in-source building - -* Sun Apr 19 2020 Antonio Trande - 5.2.1-10 -- Doc sub-package provides its own license file - -* Sat Apr 18 2020 Antonio Trande - 5.2.1-9 -- Some minor fixes -- Compile/execute Fortran tests and examples -- Do not pack example's source code - -* Tue Jan 28 2020 Fedora Release Engineering - 5.2.1-8 -- Rebuilt for https://fedoraproject.org/wiki/Fedora_32_Mass_Rebuild - -* Wed Jul 24 2019 Fedora Release Engineering - 5.2.1-7 -- Rebuilt for https://fedoraproject.org/wiki/Fedora_31_Mass_Rebuild - -* Thu Jan 31 2019 Fedora Release Engineering - 5.2.1-6 -- Rebuilt for https://fedoraproject.org/wiki/Fedora_30_Mass_Rebuild - -* Fri Sep 14 2018 Antonio Trande - 5.2.1-5 -- Remove gcc-gfortran as required package - -* Thu Jul 12 2018 Fedora Release Engineering - 5.2.1-4 -- Rebuilt for https://fedoraproject.org/wiki/Fedora_29_Mass_Rebuild - -* Wed Apr 25 2018 Antonio Trande - 5.2.1-3 -- Patch soname (5 -> 5.2) of shared library - -* Sun Apr 15 2018 Antonio Trande - 5.2.1-2 -- Use cmake3 on rhel -- Use blas from system - -* 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 - -* Thu Feb 15 2018 Antonio Trande - 5.2.0-7 -- Use %%ldconfig_scriptlets - -* Wed Feb 07 2018 Fedora Release Engineering - 5.2.0-6 -- Rebuilt for https://fedoraproject.org/wiki/Fedora_28_Mass_Rebuild - -* Sat Oct 21 2017 Mukundan Ragavan - 5.2.0-5 -- Fix lib linking in EPEL (added conditionals; Thanks Antonio Trande)) -- Use license macro - -* Wed Aug 02 2017 Fedora Release Engineering - 5.2.0-4 -- Rebuilt for https://fedoraproject.org/wiki/Fedora_27_Binutils_Mass_Rebuild - -* Wed Jul 26 2017 Fedora Release Engineering - 5.2.0-3 -- Rebuilt for https://fedoraproject.org/wiki/Fedora_27_Mass_Rebuild - -* Fri Feb 10 2017 Fedora Release Engineering - 5.2.0-2 -- Rebuilt for https://fedoraproject.org/wiki/Fedora_26_Mass_Rebuild - -* Thu Apr 14 2016 Mukundan Ragavan - 5.2.0-1 -- Update to 5.2.0 -- spec file cleanup - -* Mon Mar 28 2016 Mukundan Ragavan - 5.1.1-2 -- Added -doc subpackage -- Added GPLv2 in the license field - -* Mon Mar 21 2016 Mukundan Ragavan - 5.1.1-1 -- Update to 5.1.1 -- Remove format security patch - not needed anymore -- Edit patches to be version specific -- Renamed patch4 to be consistent with others -- Minor spec file housekeeping - -* Wed Feb 03 2016 Fedora Release Engineering - 4.3-14 -- Rebuilt for https://fedoraproject.org/wiki/Fedora_24_Mass_Rebuild - -* Tue Jun 16 2015 Fedora Release Engineering - 4.3-13 -- Rebuilt for https://fedoraproject.org/wiki/Fedora_23_Mass_Rebuild - -* Sun Jan 25 2015 Mukundan Ragavan - 4.3-12 -- Fix spec file errors and remove backup files -- fixes 1084707 - -* Fri Aug 15 2014 Fedora Release Engineering - 4.3-11 -- Rebuilt for https://fedoraproject.org/wiki/Fedora_21_22_Mass_Rebuild - -* Sun Jul 13 2014 Mukundan Ragavan - 4.3-10 -- Removed non-free files, fixes bz#1114264 - -* Fri Jun 06 2014 Fedora Release Engineering - 4.3-9 -- Rebuilt for https://fedoraproject.org/wiki/Fedora_21_Mass_Rebuild - -* Mon Jan 06 2014 Björn Esser - 4.3-8 -- fixed FTBFS if "-Werror=format-security" flag is used (#1037343) -- devel-pkg must Requires: %%{name}%%{?_isa} -- apply proper LDFLAGS -- added needed bits for el5 -- reenable testsuite using Patch3 - -* Fri Oct 4 2013 Orion Poplawski - 4.3-7 -- Rebuild for atlas 3.10 -- Handle UnversionedDocDirs change - -* Fri Aug 02 2013 Fedora Release Engineering - 4.3-6 -- Rebuilt for https://fedoraproject.org/wiki/Fedora_20_Mass_Rebuild - -* Mon Mar 25 2013 Shakthi Kannan 4.3-5 -- Ship SuperLU examples - -* Wed Feb 13 2013 Fedora Release Engineering - 4.3-4 -- Rebuilt for https://fedoraproject.org/wiki/Fedora_19_Mass_Rebuild - -* Sat Aug 25 2012 Shakthi Kannan 4.3-3 -- Use README in main package and DOC in devel package -- chmod a-x on SRC/qselect.c -- Remove -latlas linking in prep section -- Added Patch comments -- Use name RPM macro in patch name - -* Wed Feb 01 2012 Shakthi Kannan 4.3-2 -- Use atlas library instead of blas. -- Use RPM_OPT_FLAGS and LIBS when building sources. -- Use macros as required for name and version. - -* Fri Jan 06 2012 Shakthi Kannan 4.3-1 -- First release. diff --git a/dead.package b/dead.package new file mode 100644 index 0000000..caf4d55 --- /dev/null +++ b/dead.package @@ -0,0 +1 @@ +SuperLU package is retired on c9s for CS-704 diff --git a/sources b/sources deleted file mode 100644 index d565943..0000000 --- a/sources +++ /dev/null @@ -1 +0,0 @@ -SHA512 (superlu-5.2.2.tar.gz) = 091928a3df3433b337ebdacdb28de341d6d29d655965de6ffd656a6de18cf11171555bfd3af73082af62b1cead6835b4c11e4ba524a32db7f7d28db47c9d490c diff --git a/superlu-removemc64.patch b/superlu-removemc64.patch deleted file mode 100644 index 58596b8..0000000 --- a/superlu-removemc64.patch +++ /dev/null @@ -1,2662 +0,0 @@ ---- 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(); -+} -+ -+void mc64ad_(int *a, int *b, int *c, int d[], int e[], double f[], -+ int *g, int h[], int *i, int j[], int *k, double l[], -+ int m[], int n[]) -+{ -+ fprintf(stderr, "SuperLU: MC64 functionality not available (it uses non-free code). Aborting.\n"); -+ abort(); -+}