/* 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.,

		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)

/* Table of constant values */

static int_t c__1 = 1;

static int_t c__2 = 2;

/* 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@aeat.co.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 = dlamch_("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", info[1],

		   " because JOB = %d\n", *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", info[1],

		   " because N = %d\n", *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", info[1],

		   " because NE = %d\n", *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", info[1],

		   " LIW too small, must be at least %8d\n", 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", info[1],

		       " LDW too small, must be at least %8d\n", 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",

			       info[1], " Column %8d", j, 

			       " contains an entry with invalid row index %8d\n", 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", info[1],

			       "        Column %8d", j,

			       " contains two or more entries with row index %8d\n", 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", info[1],

		   " The matrix is structurally singular.\n");

	}

    }

    if (info[1] == 2) {

/* Scaling factors are large, return with warning */

	if (icntl[2] >= 0) {

	    printf(" ****** Warning from MC64A/AD. INFO(1) = %2d\n", info[1],

		   "        Some scaling factors may be too large.\n");

	}

    }

/* 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        = ", *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;

    /* 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, DLAMCH */

/*      DOUBLE PRECISION FD05AD, DLAMCH */

/* 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 = dlamch_("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;

	    }