/* 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) #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) { /* 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_ */