PS_Method_CCS_test/libhsl/ma27d.c

/* ma27d.f -- translated by f2c (version 20200916).
   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 "f2c.h"

/* Table of constant values */

static integer c__1 = 1;
static integer c__2 = 2;

/* COPYRIGHT (c) 1982 AEA Technology */
/* ######DATE 20 September 2001 */
/*  September 2001: threadsafe version of MA27 */
/*  19/3/03. Array ICNTL in MA27GD made assumed size. */
/* Subroutine */ int ma27id_(integer *icntl, doublereal *cntl)
{
/* Stream number for error messages */
    /* Parameter adjustments */
    --cntl;
    --icntl;

    /* Function Body */
    icntl[1] = 6;
/* Stream number for diagnostic messages */
    icntl[2] = 6;
/* Control the level of diagnostic printing. */
/*   0 no printing */
/*   1 printing of scalar parameters and first parts of arrays. */
/*   2 printing of scalar parameters and whole of arrays. */
    icntl[3] = 0;
/* The largest integer such that all integers I in the range */
/* -ICNTL(4).LE.I.LE.ICNTL(4) can be handled by the shortest integer */
/* type in use. */
    icntl[4] = 2139062143;
/* Minimum number of eliminations in a step that is automatically */
/* accepted. if two adjacent steps can be combined and each has less */
/* eliminations then they are combined. */
    icntl[5] = 1;
/* Control whether direct or indirect access is used by MA27C/CD. */
/* Indirect access is employed in forward and back substitution */
/* respectively if the size of a block is less than */
/* ICNTL(IFRLVL+MIN(10,NPIV)) and ICNTL(IFRLVL+10+MIN(10,NPIV)) */
/* respectively, where NPIV is the number of pivots in the block. */
    icntl[6] = 32639;
    icntl[7] = 32639;
    icntl[8] = 32639;
    icntl[9] = 32639;
    icntl[10] = 14;
    icntl[11] = 9;
    icntl[12] = 8;
    icntl[13] = 8;
    icntl[14] = 9;
    icntl[15] = 10;
    icntl[16] = 32639;
    icntl[17] = 32639;
    icntl[18] = 32639;
    icntl[19] = 32689;
    icntl[20] = 24;
    icntl[21] = 11;
    icntl[22] = 9;
    icntl[23] = 8;
    icntl[24] = 9;
    icntl[25] = 10;
/* Not used */
    icntl[26] = 0;
    icntl[27] = 0;
    icntl[28] = 0;
    icntl[29] = 0;
    icntl[30] = 0;
/* Control threshold pivoting. */
    cntl[1] = .1;
/* If a column of the reduced matrix has relative density greater than */
/* CNTL(2), it is forced into the root. All such columns are taken to */
/* have sparsity pattern equal to their merged patterns, so the fill */
/* and operation counts may be overestimated. */
    cntl[2] = 1.;
/* An entry with absolute value less than CNTL(3) is never accepted as */
/* a 1x1 pivot or as the off-diagonal of a 2x2 pivot. */
    cntl[3] = 0.;
/* Not used */
    cntl[4] = 0.f;
    cntl[5] = 0.f;
    return 0;
} /* ma27id_ */

/* Subroutine */ int ma27ad_(integer *n, integer *nz, integer *irn, integer *
	icn, integer *iw, integer *liw, integer *ikeep, integer *iw1, integer 
	*nsteps, integer *iflag, integer *icntl, doublereal *cntl, integer *
	info, doublereal *ops)
{
    /* Format strings */
    static char fmt_10[] = "(/,/,\002 ENTERING MA27AD WITH      N     NZ    "
	    "  LIW  IFLAG\002,/,21x,i7,i7,i9,i7)";
    static char fmt_20[] = "(\002 MATRIX NON-ZEROS\002,/,4(i9,i6),/,(i9,i6,i"
	    "9,i6,i9,i6,i9,i6))";
    static char fmt_30[] = "(\002 IKEEP(.,1)=\002,10i6,/,(12x,10i6))";
    static char fmt_80[] = "(\002 **** ERROR RETURN FROM MA27AD **** INFO("
	    "1)=\002,i3)";
    static char fmt_90[] = "(\002 VALUE OF N OUT OF RANGE ... =\002,i10)";
    static char fmt_110[] = "(\002 VALUE OF NZ OUT OF RANGE .. =\002,i10)";
    static char fmt_150[] = "(\002 LIW TOO SMALL, MUST BE INCREASED FROM\002"
	    ",i10,\002 TO AT LEAST\002,i10)";
    static char fmt_170[] = "(/,\002 LEAVING MA27AD WITH NSTEPS  INFO(1)    "
	    "OPS IERROR\002,\002 NRLTOT NIRTOT\002,/,20x,2i7,f7.0,3i7,/,20x"
	    ",\002 NRLNEC NIRNEC NRLADU NIRADU  NCMPA\002,/,20x,6i7)";
    static char fmt_180[] = "(\002 IKEEP(.,2)=\002,10i6,/,(12x,10i6))";
    static char fmt_190[] = "(\002 IKEEP(.,3)=\002,10i6,/,(12x,10i6))";

    /* System generated locals */
    integer ikeep_dim1, ikeep_offset, iw1_dim1, iw1_offset, i__1;

    /* Builtin functions */
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);

    /* Local variables */
    static integer i__, k, l1, l2, iwfr, lliw;
    extern /* Subroutine */ int ma27gd_(integer *, integer *, integer *, 
	    integer *, integer *, integer *, integer *, integer *, integer *, 
	    integer *, integer *, integer *), ma27hd_(integer *, integer *, 
	    integer *, integer *, integer *, integer *, integer *, integer *, 
	    integer *, integer *, integer *, integer *, doublereal *), 
	    ma27jd_(integer *, integer *, integer *, integer *, integer *, 
	    integer *, integer *, integer *, integer *, integer *, integer *, 
	    integer *, integer *), ma27kd_(integer *, integer *, integer *, 
	    integer *, integer *, integer *, integer *, integer *, integer *, 
	    integer *), ma27ld_(integer *, integer *, integer *, integer *, 
	    integer *, integer *, integer *, integer *, integer *), ma27md_(
	    integer *, integer *, integer *, integer *, integer *, integer *, 
	    integer *, integer *, integer *, integer *, integer *, integer *, 
	    integer *, doublereal *);

    /* Fortran I/O blocks */
    static cilist io___2 = { 0, 0, 0, fmt_10, 0 };
    static cilist io___4 = { 0, 0, 0, fmt_20, 0 };
    static cilist io___5 = { 0, 0, 0, fmt_30, 0 };
    static cilist io___10 = { 0, 0, 0, fmt_80, 0 };
    static cilist io___11 = { 0, 0, 0, fmt_90, 0 };
    static cilist io___12 = { 0, 0, 0, fmt_80, 0 };
    static cilist io___13 = { 0, 0, 0, fmt_110, 0 };
    static cilist io___14 = { 0, 0, 0, fmt_80, 0 };
    static cilist io___15 = { 0, 0, 0, fmt_150, 0 };
    static cilist io___16 = { 0, 0, 0, fmt_170, 0 };
    static cilist io___17 = { 0, 0, 0, fmt_30, 0 };
    static cilist io___18 = { 0, 0, 0, fmt_180, 0 };
    static cilist io___19 = { 0, 0, 0, fmt_190, 0 };


/* THIS SUBROUTINE COMPUTES A MINIMUM DEGREE ORDERING OR ACCEPTS A GIVEN */
/*     ORDERING. IT COMPUTES ASSOCIATED ASSEMBLY AND ELIMINATION */
/*     INFORMATION FOR MA27B/BD. */
/* N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED. */
/* NZ MUST BE SET TO THE NUMBER OF NON-ZEROS INPUT. IT IS NOT */
/*     ALTERED. */
/* IRN(I),I=1,2,...,NZ MUST BE SET TO THE ROW NUMBERS OF THE */
/*     NON-ZEROS. IT IS NOT ALTERED UNLESS IT IS EQUIVALENCED */
/*     TO IW (SEE DESCRIPTION OF IW). */
/* ICN(I),I=1,2,...,NZ MUST BE SET TO THE COLUMN NUMBERS OF THE */
/*     NON-ZEROS. IT IS NOT ALTERED UNLESS IT IS EQUIVALENCED */
/*     TO IW (SEE DESCRIPTION OF IW). */
/* IW NEED NOT BE SET ON INPUT. IT IS USED FOR WORKSPACE. */
/*     IRN(1) MAY BE EQUIVALENCED TO IW(1) AND ICN(1) MAY BE */
/*     EQUIVALENCED TO IW(K), WHERE K.GT.NZ. */
/* LIW MUST BE SET TO THE LENGTH OF IW. IT MUST BE AT LEAST 2*NZ+3*N */
/*     FOR THE IFLAG=0 ENTRY AND AT LEAST NZ+3*N FOR THE IFLAG=1 */
/*     ENTRY. IT IS NOT ALTERED. */
/* IKEEP NEED NOT BE SET UNLESS AN ORDERING IS GIVEN, IN WHICH CASE */
/*     IKEEP(I,1) MUST BE SET TO THE POSITION OF VARIABLE I IN THE */
/*     ORDER. ON OUTPUT IKEEP CONTAINS INFORMATION NEEDED BY MA27B/BD. */
/*     IKEEP(I,1) HOLDS THE POSITION OF VARIABLE I IN THE PIVOT ORDER. */
/*     IKEEP(I,2), IKEEP(I,3) HOLD THE NUMBER OF ELIMINATIONS, ASSEMBLIES */
/*     AT MAJOR STEP I, I=1,2,...,NSTEPS. NOTE THAT WHEN AN ORDER IS */
/*     GIVEN IT MAY BE REPLACED BY ANOTHER ORDER THAT GIVES IDENTICAL */
/*     NUMERICAL RESULTS. */
/* IW1 IS USED FOR WORKSPACE. */
/* NSTEPS NEED NOT BE SET. ON OUTPUT IT CONTAINS THE NUMBER OF MAJOR */
/*     STEPS NEEDED FOR A DEFINITE MATRIX AND MUST BE PASSED UNCHANGED */
/*     TO MA27B/BD. */
/* IFLAG MUST SET TO ZERO IF THE USER WANTS THE PIVOT ORDER CHOSEN */
/*     AUTOMATICALLY AND TO ONE IF HE WANTS TO SPECIFY IT IN IKEEP. */
/* ICNTL is an INTEGER array of length 30 containing user options */
/*     with integer values (defaults set in MA27I/ID) */
/*   ICNTL(1) (LP) MUST BE SET TO THE STREAM NUMBER FOR ERROR MESSAGES. */
/*     ERROR MESSAGES ARE SUPPRESSED IF ICNTL(1) IS NOT POSITIVE. */
/*     IT IS NOT ALTERED. */
/*   ICNTL(2) (MP) MUST BE SET TO THE STREAM NUMBER FOR DIAGNOSTIC */
/*     MESSAGES.  DIAGNOSTIC MESSAGES ARE SUPPRESSED IF ICNTL(2) IS NOT */
/*     POSITIVE.  IT IS NOT ALTERED. */
/*   ICNTL(3) (LDIAG) CONTROLS THE LEVEL OF DIAGNOSTIC PRINTING. */
/*     0 NO PRINTING */
/*     1 PRINTING OF SCALAR PARAMETERS AND FIRST PARTS OF ARRAYS. */
/*     2 PRINTING OF SCALAR PARAMETERS AND WHOLE OF ARRAYS. */
/*   ICNTL(4) (IOVFLO) IS THE LARGEST INTEGER SUCH THAT ALL INTEGERS */
/*     I IN THE RANGE -IOVFLO.LE.I.LE.IOVFLO CAN BE HANDLED BY THE */
/*     SHORTEST INTEGER TYPE IN USE. */
/*   ICNT(5) (NEMIN) MUST BE SET TO THE MINIMUM NUMBER OF ELIMINATIONS */
/*     IN A STEP THAT IS AUTOMATICALLY ACCEPTED. IF TWO ADJACENT STEPS */
/*     CAN BE COMBINED AND EACH HAS LESS ELIMINATIONS THEN THEY ARE */
/*     COMBINED. */
/*   ICNTL(IFRLVL+I) I=1,20, (IFRLVL) MUST BE SET TO CONTROL WHETHER */
/*     DIRECT OR INDIRECT ACCESS IS USED BY MA27C/CD.  INDIRECT ACCESS */
/*     IS EMPLOYED IN FORWARD AND BACK SUBSTITUTION RESPECTIVELY IF THE */
/*     SIZE OF A BLOCK IS LESS THAN ICNTL(IFRLVL+(MIN(10,NPIV)) AND */
/*     ICNTL(IFRLVL+10+MIN(10,NPIV)) RESPECTIVELY, WHERE NPIV IS THE */
/*     NUMBER OF PIVOTS IN THE BLOCK. */
/*   ICNTL(I) I=26,30 are not used. */
/* CNTL is an DOUBLE PRECISION array of length 5 containing user options */
/*     with real values (defaults set in MA27I/ID) */
/*   CNTL(1) (U) IS USED TO CONTROL THRESHOLD PIVOTING. IT IS NOT */
/*     ALTERED. */
/*   CNTL(2) (FRATIO) has default value 1.0.  If a column of the */
/*      reduced matrix has relative density greater than CNTL(2), it */
/*      is forced into the root. All such columns are taken to have */
/*      sparsity pattern equal to their merged patterns, so the fill */
/*      and operation counts may be overestimated. */
/*   CNTL(3) (PIVTOL) has default value 0.0. An entry with absolute */
/*      value less than CNTL(3) is never accepted as a 1x1 pivot or */
/*      as the off-diagonal of a 2x2 pivot. */
/*   CNTL(I) I=4,5 are not used. */
/* INFO is an INTEGER array of length 20 which is used to return */
/*     information to the user. */
/*   INFO(1) (IFLAG) is an error return code, zero for success, greater */
/*      than zero for a warning and less than zero for errors, see */
/*      INFO(2). */
/*   INFO(2) (IERROR) HOLDS ADDITIONAL INFORMATION IN THE EVENT OF ERRORS. */
/*     IF INFO(1)=-3 INFO(2) HOLDS A LENGTH THAT MAY SUFFICE FOR IW. */
/*     IF INFO(1)=-4 INFO(2) HOLDS A LENGTH THAT MAY SUFFICE FOR A. */
/*     IF INFO(1)=-5 INFO(2) IS SET TO THE PIVOT STEP AT WHICH SINGULARITY */
/*                 WAS DETECTED. */
/*     IF INFO(1)=-6 INFO(2) IS SET TO THE PIVOT STEP AT WHICH A CHANGE OF */
/*                 PIVOT SIGN WAS FOUND. */
/*     IF INFO(1)= 1 INFO(2) HOLDS THE NUMBER OF FAULTY ENTRIES. */
/*     IF INFO(1)= 2 INFO(2) IS SET TO THE NUMBER OF SIGNS CHANGES IN */
/*                 THE PIVOTS. */
/*     IF INFO(1)=3 INFO(2) IS SET TO THE RANK OF THE MATRIX. */
/*   INFO(3) and INFO(4) (NRLTOT and NIRTOT) REAL AND INTEGER STRORAGE */
/*     RESPECTIVELY REQUIRED FOR THE FACTORIZATION IF NO COMPRESSES ARE */
/*     ALLOWED. */
/*   INFO(5) and INFO(6) (NRLNEC and NIRNEC) REAL AND INTEGER STORAGE */
/*     RESPECTIVELY REQUIRED FOR THE FACTORIZATION IF COMPRESSES ARE */
/*     ALLOWED AND THE MATRIX IS DEFINITE. */
/*   INFO(7) and INFO(8) (NRLADU and NIRADU) REAL AND INTEGER STORAGE */
/*     RESPECTIVELY FOR THE MATRIX FACTORS AS CALCULATED BY MA27A/AD */
/*     FOR THE DEFINITE CASE. */
/*   INFO(9) and INFO(10) (NRLBDU and NIRBDU) REAL AND INTEGER STORAGE */
/*     RESPECTIVELY FOR THE MATRIX FACTORS AS FOUND  BY MA27B/BD. */
/*   INFO(11) (NCMPA) ACCUMULATES THE NUMBER OF TIMES THE ARRAY IW IS */
/*     COMPRESSED BY MA27A/AD. */
/*   INFO(12) and INFO(13) (NCMPBR and NCMPBI) ACCUMULATE THE NUMBER */
/*     OF COMPRESSES OF THE REALS AND INTEGERS PERFORMED BY MA27B/BD. */
/*   INFO(14) (NTWO) IS USED BY MA27B/BD TO RECORD THE NUMBER OF 2*2 */
/*     PIVOTS USED. */
/*   INFO(15) (NEIG) IS USED BY ME27B/BD TO RECORD THE NUMBER OF */
/*     NEGATIVE EIGENVALUES OF A. */
/*   INFO(16) to INFO(20) are not used. */
/* OPS ACCUMULATES THE NO. OF MULTIPLY/ADD PAIRS NEEDED TO CREATE THE */
/*     TRIANGULAR FACTORIZATION, IN THE DEFINITE CASE. */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. External Subroutines .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */
    /* Parameter adjustments */
    iw1_dim1 = *n;
    iw1_offset = 1 + iw1_dim1;
    iw1 -= iw1_offset;
    ikeep_dim1 = *n;
    ikeep_offset = 1 + ikeep_dim1;
    ikeep -= ikeep_offset;
    --irn;
    --icn;
    --iw;
    --icntl;
    --cntl;
    --info;

    /* Function Body */
    for (i__ = 1; i__ <= 15; ++i__) {
	info[i__] = 0;
/* L5: */
    }
    if (icntl[3] <= 0 || icntl[2] <= 0) {
	goto L40;
    }
/* PRINT INPUT VARIABLES. */
    io___2.ciunit = icntl[2];
    s_wsfe(&io___2);
    do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&(*nz), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&(*liw), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&(*iflag), (ftnlen)sizeof(integer));
    e_wsfe();
    *nsteps = 0;
    k = min(8,*nz);
    if (icntl[3] > 1) {
	k = *nz;
    }
    if (k > 0) {
	io___4.ciunit = icntl[2];
	s_wsfe(&io___4);
	i__1 = k;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    do_fio(&c__1, (char *)&irn[i__], (ftnlen)sizeof(integer));
	    do_fio(&c__1, (char *)&icn[i__], (ftnlen)sizeof(integer));
	}
	e_wsfe();
    }
    k = min(10,*n);
    if (icntl[3] > 1) {
	k = *n;
    }
    if (*iflag == 1 && k > 0) {
	io___5.ciunit = icntl[2];
	s_wsfe(&io___5);
	i__1 = k;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    do_fio(&c__1, (char *)&ikeep[i__ + ikeep_dim1], (ftnlen)sizeof(
		    integer));
	}
	e_wsfe();
    }
L40:
    if (*n < 1 || *n > icntl[4]) {
	goto L70;
    }
    if (*nz < 0) {
	goto L100;
    }
    lliw = *liw - (*n << 1);
    l1 = lliw + 1;
    l2 = l1 + *n;
    if (*iflag == 1) {
	goto L50;
    }
    if (*liw < (*nz << 1) + *n * 3 + 1) {
	goto L130;
    }
/* SORT */
    ma27gd_(n, nz, &irn[1], &icn[1], &iw[1], &lliw, &iw1[iw1_offset], &iw1[(
	    iw1_dim1 << 1) + 1], &iw[l1], &iwfr, &icntl[1], &info[1]);
/* ANALYZE USING MINIMUM DEGREE ORDERING */
    ma27hd_(n, &iw1[iw1_offset], &iw[1], &lliw, &iwfr, &iw[l1], &iw[l2], &
	    ikeep[(ikeep_dim1 << 1) + 1], &ikeep[ikeep_dim1 * 3 + 1], &ikeep[
	    ikeep_offset], &icntl[4], &info[11], &cntl[2]);
    goto L60;
/* SORT USING GIVEN ORDER */
L50:
    if (*liw < *nz + *n * 3 + 1) {
	goto L120;
    }
    ma27jd_(n, nz, &irn[1], &icn[1], &ikeep[ikeep_offset], &iw[1], &lliw, &
	    iw1[iw1_offset], &iw1[(iw1_dim1 << 1) + 1], &iw[l1], &iwfr, &
	    icntl[1], &info[1]);
/* ANALYZE USING GIVEN ORDER */
    ma27kd_(n, &iw1[iw1_offset], &iw[1], &lliw, &iwfr, &ikeep[ikeep_offset], &
	    ikeep[(ikeep_dim1 << 1) + 1], &iw[l1], &iw[l2], &info[11]);
/* PERFORM DEPTH-FIRST SEARCH OF ASSEMBLY TREE */
L60:
    ma27ld_(n, &iw1[iw1_offset], &iw[l1], &ikeep[ikeep_offset], &ikeep[(
	    ikeep_dim1 << 1) + 1], &ikeep[ikeep_dim1 * 3 + 1], &iw[l2], 
	    nsteps, &icntl[5]);
/* EVALUATE STORAGE AND OPERATION COUNT REQUIRED BY MA27B/BD IN THE */
/*     DEFINITE CASE. */
/* SET IW(1) SO THAT ARRAYS IW AND IRN CAN BE TESTED FOR EQUIVALENCE */
/*     IN MA27M/MD. */
    if (*nz >= 1) {
	iw[1] = irn[1] + 1;
    }
    ma27md_(n, nz, &irn[1], &icn[1], &ikeep[ikeep_offset], &ikeep[ikeep_dim1 *
	     3 + 1], &ikeep[(ikeep_dim1 << 1) + 1], &iw[l2], nsteps, &iw1[
	    iw1_offset], &iw1[(iw1_dim1 << 1) + 1], &iw[1], &info[1], ops);
    goto L160;
L70:
    info[1] = -1;
    if (icntl[1] > 0) {
	io___10.ciunit = icntl[1];
	s_wsfe(&io___10);
	do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer));
	e_wsfe();
    }
    if (icntl[1] > 0) {
	io___11.ciunit = icntl[1];
	s_wsfe(&io___11);
	do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
	e_wsfe();
    }
    goto L160;
L100:
    info[1] = -2;
    if (icntl[1] > 0) {
	io___12.ciunit = icntl[1];
	s_wsfe(&io___12);
	do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer));
	e_wsfe();
    }
    if (icntl[1] > 0) {
	io___13.ciunit = icntl[1];
	s_wsfe(&io___13);
	do_fio(&c__1, (char *)&(*nz), (ftnlen)sizeof(integer));
	e_wsfe();
    }
    goto L160;
L120:
    info[2] = *nz + *n * 3 + 1;
    goto L140;
L130:
    info[2] = (*nz << 1) + *n * 3 + 1;
L140:
    info[1] = -3;
    if (icntl[1] > 0) {
	io___14.ciunit = icntl[1];
	s_wsfe(&io___14);
	do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer));
	e_wsfe();
    }
    if (icntl[1] > 0) {
	io___15.ciunit = icntl[1];
	s_wsfe(&io___15);
	do_fio(&c__1, (char *)&(*liw), (ftnlen)sizeof(integer));
	do_fio(&c__1, (char *)&info[2], (ftnlen)sizeof(integer));
	e_wsfe();
    }
L160:
    if (icntl[3] <= 0 || icntl[2] <= 0 || info[1] < 0) {
	goto L200;
    }
/* PRINT PARAMETER VALUES ON EXIT. */
    io___16.ciunit = icntl[2];
    s_wsfe(&io___16);
    do_fio(&c__1, (char *)&(*nsteps), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&(*ops), (ftnlen)sizeof(doublereal));
    do_fio(&c__1, (char *)&info[2], (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&info[3], (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&info[4], (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&info[5], (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&info[6], (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&info[7], (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&info[8], (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&info[11], (ftnlen)sizeof(integer));
    e_wsfe();
    k = min(9,*n);
    if (icntl[3] > 1) {
	k = *n;
    }
    if (k > 0) {
	io___17.ciunit = icntl[2];
	s_wsfe(&io___17);
	i__1 = k;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    do_fio(&c__1, (char *)&ikeep[i__ + ikeep_dim1], (ftnlen)sizeof(
		    integer));
	}
	e_wsfe();
    }
    k = min(k,*nsteps);
    if (k > 0) {
	io___18.ciunit = icntl[2];
	s_wsfe(&io___18);
	i__1 = k;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    do_fio(&c__1, (char *)&ikeep[i__ + (ikeep_dim1 << 1)], (ftnlen)
		    sizeof(integer));
	}
	e_wsfe();
    }
    if (k > 0) {
	io___19.ciunit = icntl[2];
	s_wsfe(&io___19);
	i__1 = k;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    do_fio(&c__1, (char *)&ikeep[i__ + ikeep_dim1 * 3], (ftnlen)
		    sizeof(integer));
	}
	e_wsfe();
    }
L200:
    return 0;
} /* ma27ad_ */

/* Subroutine */ int ma27bd_(integer *n, integer *nz, integer *irn, integer *
	icn, doublereal *a, integer *la, integer *iw, integer *liw, integer *
	ikeep, integer *nsteps, integer *maxfrt, integer *iw1, integer *icntl,
	 doublereal *cntl, integer *info)
{
    /* Format strings */
    static char fmt_10[] = "(/,/,\002 ENTERING MA27BD WITH      N     NZ    "
	    "   LA      LIW\002,\002 NSTEPS      U\002,/,21x,i7,i7,i9,i9,i7,1"
	    "pd10.2)";
    static char fmt_20[] = "(\002 MATRIX NON-ZEROS\002,/,1x,2(1p,d16.3,2i6),"
	    "/,(1x,1p,d16.3,2i6,1p,d16.3,2i6))";
    static char fmt_30[] = "(\002 IKEEP(.,1)=\002,10i6,/,(12x,10i6))";
    static char fmt_40[] = "(\002 IKEEP(.,2)=\002,10i6,/,(12x,10i6))";
    static char fmt_50[] = "(\002 IKEEP(.,3)=\002,10i6,/,(12x,10i6))";
    static char fmt_65[] = "(\002 *** WARNING MESSAGE FROM SUBROUTINE MA27B"
	    "D\002,\002  *** INFO(1) =\002,i2,/,5x,\002MATRIX IS SINGULAR. RA"
	    "NK=\002,i5)";
    static char fmt_80[] = "(\002 **** ERROR RETURN FROM MA27BD **** INFO("
	    "1)=\002,i3)";
    static char fmt_90[] = "(\002 VALUE OF N OUT OF RANGE ... =\002,i10)";
    static char fmt_110[] = "(\002 VALUE OF NZ OUT OF RANGE .. =\002,i10)";
    static char fmt_140[] = "(\002 LIW TOO SMALL, MUST BE INCREASED FROM\002"
	    ",i10,\002 TO\002,\002 AT LEAST\002,i10)";
    static char fmt_170[] = "(\002 LA TOO SMALL, MUST BE INCREASED FROM \002"
	    ",i10,\002 TO\002,\002 AT LEAST\002,i10)";
    static char fmt_190[] = "(\002 ZERO PIVOT AT STAGE\002,i10,\002 WHEN INP"
	    "UT MATRIX DECLARED DEFINITE\002)";
    static char fmt_210[] = "(\002 CHANGE IN SIGN OF PIVOT ENCOUNTERED\002"
	    ",\002 WHEN FACTORING ALLEGEDLY DEFINITE MATRIX\002)";
    static char fmt_230[] = "(/,\002 LEAVING MA27BD WITH\002,/,10x,\002  MAX"
	    "FRT  INFO(1) NRLBDU NIRBDU NCMPBR\002,\002 NCMPBI   NTWO IERRO"
	    "R\002,/,11x,8i7)";
    static char fmt_250[] = "(\002 BLOCK PIVOT =\002,i8,\002 NROWS =\002,i8"
	    ",\002 NCOLS =\002,i8)";
    static char fmt_260[] = "(\002 COLUMN INDICES =\002,10i6,/,(17x,10i6))";
    static char fmt_270[] = "(\002 REAL ENTRIES .. EACH ROW STARTS ON A NEW "
	    "LINE\002)";
    static char fmt_280[] = "(1p,5d13.3)";

    /* System generated locals */
    integer ikeep_dim1, ikeep_offset, i__1, i__2, i__3;

    /* Builtin functions */
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);

    /* Local variables */
    static integer i__, k, j1, j2, jj, kz, nz1, len, iblk, kblk, ipos;
    extern /* Subroutine */ int ma27nd_(integer *, integer *, integer *, 
	    doublereal *, integer *, integer *, integer *, integer *, integer 
	    *, integer *, integer *, integer *, integer *), ma27od_(integer *,
	     integer *, doublereal *, integer *, integer *, integer *, 
	    integer *, integer *, integer *, integer *, integer *, integer *, 
	    integer *, doublereal *, integer *);
    static integer iapos, ncols, irows, nrows;

    /* Fortran I/O blocks */
    static cilist io___20 = { 0, 0, 0, fmt_10, 0 };
    static cilist io___22 = { 0, 0, 0, fmt_20, 0 };
    static cilist io___24 = { 0, 0, 0, fmt_30, 0 };
    static cilist io___26 = { 0, 0, 0, fmt_40, 0 };
    static cilist io___27 = { 0, 0, 0, fmt_50, 0 };
    static cilist io___29 = { 0, 0, 0, fmt_65, 0 };
    static cilist io___30 = { 0, 0, 0, fmt_80, 0 };
    static cilist io___31 = { 0, 0, 0, fmt_90, 0 };
    static cilist io___32 = { 0, 0, 0, fmt_80, 0 };
    static cilist io___33 = { 0, 0, 0, fmt_110, 0 };
    static cilist io___34 = { 0, 0, 0, fmt_80, 0 };
    static cilist io___35 = { 0, 0, 0, fmt_140, 0 };
    static cilist io___36 = { 0, 0, 0, fmt_80, 0 };
    static cilist io___37 = { 0, 0, 0, fmt_170, 0 };
    static cilist io___38 = { 0, 0, 0, fmt_80, 0 };
    static cilist io___39 = { 0, 0, 0, "(A)", 0 };
    static cilist io___40 = { 0, 0, 0, fmt_80, 0 };
    static cilist io___41 = { 0, 0, 0, fmt_190, 0 };
    static cilist io___42 = { 0, 0, 0, fmt_80, 0 };
    static cilist io___43 = { 0, 0, 0, fmt_210, 0 };
    static cilist io___44 = { 0, 0, 0, fmt_230, 0 };
    static cilist io___52 = { 0, 0, 0, fmt_250, 0 };
    static cilist io___54 = { 0, 0, 0, fmt_260, 0 };
    static cilist io___56 = { 0, 0, 0, fmt_270, 0 };
    static cilist io___59 = { 0, 0, 0, fmt_280, 0 };


/* THIS SUBROUTINE COMPUTES THE FACTORISATION OF THE MATRIX INPUT IN */
/*     A,IRN,ICN USING INFORMATION (IN IKEEP) FROM MA27A/AD. */
/* N MUST BE SET TO THE ORDER OF THE MATRIX.  IT IS NOT ALTERED. */
/* NZ MUST BE SET TO THE NUMBER OF NON-ZEROS INPUT.  IT IS NOT */
/*     ALTERED. */
/* IRN,ICN,A.  ENTRY K (K=1,NZ) OF IRN,ICN MUST BE SET TO THE ROW */
/*     AND COLUMN INDEX RESPECTIVELY OF THE NON-ZERO IN A. */
/*     IRN AND ICN ARE UNALTERED BY MA27B/BD. */
/*     ON EXIT, ENTRIES 1 TO NRLBDU OF A HOLD REAL INFORMATION */
/*     ON THE FACTORS AND SHOULD BE PASSED UNCHANGED TO MA27C/CD. */
/* LA LENGTH OF ARRAY A.  AN INDICATION OF A SUITABLE VALUE, */
/*     SUFFICIENT FOR FACTORIZATION OF A DEFINITE MATRIX, WILL */
/*     HAVE BEEN PROVIDED IN NRLNEC AND NRLTOT BY MA27A/AD. */
/*     IT IS NOT ALTERED BY MA27B/BD. */
/* IW NEED NOT BE SET ON ENTRY.  USED AS A WORK ARRAY BY MA27B/BD. */
/*     ON EXIT, ENTRIES 1 TO NIRBDU HOLD INTEGER INFORMATION ON THE */
/*     FACTORS AND SHOULD BE PASSED UNCHANGED TO MA27C/CD. */
/* LIW LENGTH OF ARRAY IW.  AN INDICATION OF A SUITABLE VALUE WILL */
/*     HAVE BEEN PROVIDED IN NIRNEC AND NIRTOT BY MA27A/AD. */
/*     IT IS NOT ALTERED BY MA27B/BD. */
/* IKEEP MUST BE UNCHANGED SINCE THE CALL TO MA27A/AD. */
/*     IT IS NOT ALTERED BY MA27B/BD. */
/* NSTEPS MUST BE UNCHANGED SINCE THE CALL TO MA27A/AD. */
/*     IT IS NOT ALTERED BY MA27B/BD. */
/* MAXFRT NEED NOT BE SET AND MUST BE PASSED UNCHANGED TO MA27C/CD. */
/*     IT IS THE MAXIMUM SIZE OF THE FRONTAL MATRIX GENERATED DURING */
/*     THE DECOMPOSITION. */
/* IW1 USED AS WORKSPACE BY MA27B/BD. */
/* ICNTL is an INTEGER array of length 30, see MA27A/AD. */
/* CNTL is a DOUBLE PRECISION array of length 5, see MA27A/AD. */
/* INFO is an INTEGER array of length 20, see MA27A/AD. */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. External Subroutines .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */
    /* Parameter adjustments */
    --iw1;
    ikeep_dim1 = *n;
    ikeep_offset = 1 + ikeep_dim1;
    ikeep -= ikeep_offset;
    --irn;
    --icn;
    --a;
    --iw;
    --icntl;
    --cntl;
    --info;

    /* Function Body */
    info[1] = 0;
    if (icntl[3] <= 0 || icntl[2] <= 0) {
	goto L60;
    }
/* PRINT INPUT PARAMETERS. */
    io___20.ciunit = icntl[2];
    s_wsfe(&io___20);
    do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&(*nz), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&(*la), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&(*liw), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&(*nsteps), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&cntl[1], (ftnlen)sizeof(doublereal));
    e_wsfe();
    kz = min(6,*nz);
    if (icntl[3] > 1) {
	kz = *nz;
    }
    if (*nz > 0) {
	io___22.ciunit = icntl[2];
	s_wsfe(&io___22);
	i__1 = kz;
	for (k = 1; k <= i__1; ++k) {
	    do_fio(&c__1, (char *)&a[k], (ftnlen)sizeof(doublereal));
	    do_fio(&c__1, (char *)&irn[k], (ftnlen)sizeof(integer));
	    do_fio(&c__1, (char *)&icn[k], (ftnlen)sizeof(integer));
	}
	e_wsfe();
    }
    k = min(9,*n);
    if (icntl[3] > 1) {
	k = *n;
    }
    if (k > 0) {
	io___24.ciunit = icntl[2];
	s_wsfe(&io___24);
	i__1 = k;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    do_fio(&c__1, (char *)&ikeep[i__ + ikeep_dim1], (ftnlen)sizeof(
		    integer));
	}
	e_wsfe();
    }
    k = min(k,*nsteps);
    if (k > 0) {
	io___26.ciunit = icntl[2];
	s_wsfe(&io___26);
	i__1 = k;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    do_fio(&c__1, (char *)&ikeep[i__ + (ikeep_dim1 << 1)], (ftnlen)
		    sizeof(integer));
	}
	e_wsfe();
    }
    if (k > 0) {
	io___27.ciunit = icntl[2];
	s_wsfe(&io___27);
	i__1 = k;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    do_fio(&c__1, (char *)&ikeep[i__ + ikeep_dim1 * 3], (ftnlen)
		    sizeof(integer));
	}
	e_wsfe();
    }
L60:
    if (*n < 1 || *n > icntl[4]) {
	goto L70;
    }
    if (*nz < 0) {
	goto L100;
    }
    if (*liw < *nz) {
	goto L120;
    }
    if (*la < *nz + *n) {
	goto L150;
    }
    if (*nsteps < 1 || *nsteps > *n) {
	goto L175;
    }
/* SORT */
    ma27nd_(n, nz, &nz1, &a[1], la, &irn[1], &icn[1], &iw[1], liw, &ikeep[
	    ikeep_offset], &iw1[1], &icntl[1], &info[1]);
    if (info[1] == -3) {
	goto L130;
    }
    if (info[1] == -4) {
	goto L160;
    }
/* FACTORIZE */
    ma27od_(n, &nz1, &a[1], la, &iw[1], liw, &ikeep[ikeep_offset], &ikeep[
	    ikeep_dim1 * 3 + 1], nsteps, maxfrt, &ikeep[(ikeep_dim1 << 1) + 1]
	    , &iw1[1], &icntl[1], &cntl[1], &info[1]);
    if (info[1] == -3) {
	goto L130;
    }
    if (info[1] == -4) {
	goto L160;
    }
    if (info[1] == -5) {
	goto L180;
    }
    if (info[1] == -6) {
	goto L200;
    }
/* **** WARNING MESSAGE **** */
    if (info[1] == 3 && icntl[2] > 0) {
	io___29.ciunit = icntl[2];
	s_wsfe(&io___29);
	do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer));
	do_fio(&c__1, (char *)&info[2], (ftnlen)sizeof(integer));
	e_wsfe();
    }
    goto L220;
/* **** ERROR RETURNS **** */
L70:
    info[1] = -1;
    if (icntl[1] > 0) {
	io___30.ciunit = icntl[1];
	s_wsfe(&io___30);
	do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer));
	e_wsfe();
    }
    if (icntl[1] > 0) {
	io___31.ciunit = icntl[1];
	s_wsfe(&io___31);
	do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
	e_wsfe();
    }
    goto L220;
L100:
    info[1] = -2;
    if (icntl[1] > 0) {
	io___32.ciunit = icntl[1];
	s_wsfe(&io___32);
	do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer));
	e_wsfe();
    }
    if (icntl[1] > 0) {
	io___33.ciunit = icntl[1];
	s_wsfe(&io___33);
	do_fio(&c__1, (char *)&(*nz), (ftnlen)sizeof(integer));
	e_wsfe();
    }
    goto L220;
L120:
    info[1] = -3;
    info[2] = *nz;
L130:
    if (icntl[1] > 0) {
	io___34.ciunit = icntl[1];
	s_wsfe(&io___34);
	do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer));
	e_wsfe();
    }
    if (icntl[1] > 0) {
	io___35.ciunit = icntl[1];
	s_wsfe(&io___35);
	do_fio(&c__1, (char *)&(*liw), (ftnlen)sizeof(integer));
	do_fio(&c__1, (char *)&info[2], (ftnlen)sizeof(integer));
	e_wsfe();
    }
    goto L220;
L150:
    info[1] = -4;
    info[2] = *nz + *n;
L160:
    if (icntl[1] > 0) {
	io___36.ciunit = icntl[1];
	s_wsfe(&io___36);
	do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer));
	e_wsfe();
    }
    if (icntl[1] > 0) {
	io___37.ciunit = icntl[1];
	s_wsfe(&io___37);
	do_fio(&c__1, (char *)&(*la), (ftnlen)sizeof(integer));
	do_fio(&c__1, (char *)&info[2], (ftnlen)sizeof(integer));
	e_wsfe();
    }
    goto L220;
L175:
    info[1] = -7;
    if (icntl[1] > 0) {
	io___38.ciunit = icntl[1];
	s_wsfe(&io___38);
	do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer));
	e_wsfe();
    }
    if (icntl[1] > 0) {
	io___39.ciunit = icntl[1];
	s_wsfe(&io___39);
	do_fio(&c__1, " NSTEPS is out of range", (ftnlen)23);
	e_wsfe();
    }
    goto L220;
L180:
    if (icntl[1] > 0) {
	io___40.ciunit = icntl[1];
	s_wsfe(&io___40);
	do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer));
	e_wsfe();
    }
    if (icntl[1] > 0) {
	io___41.ciunit = icntl[1];
	s_wsfe(&io___41);
	do_fio(&c__1, (char *)&info[2], (ftnlen)sizeof(integer));
	e_wsfe();
    }
    goto L220;
L200:
    if (icntl[1] > 0) {
	io___42.ciunit = icntl[1];
	s_wsfe(&io___42);
	do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer));
	e_wsfe();
    }
    if (icntl[1] > 0) {
	io___43.ciunit = icntl[1];
	s_wsfe(&io___43);
	e_wsfe();
    }
L220:
    if (icntl[3] <= 0 || icntl[2] <= 0 || info[1] < 0) {
	goto L310;
    }
/* PRINT OUTPUT PARAMETERS. */
    io___44.ciunit = icntl[2];
    s_wsfe(&io___44);
    do_fio(&c__1, (char *)&(*maxfrt), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&info[9], (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&info[10], (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&info[12], (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&info[13], (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&info[14], (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&info[2], (ftnlen)sizeof(integer));
    e_wsfe();
    if (info[1] < 0) {
	goto L310;
    }
/* PRINT OUT MATRIX FACTORS FROM MA27B/BD. */
    kblk = abs(iw[1]);
    if (kblk == 0) {
	goto L310;
    }
    if (icntl[3] == 1) {
	kblk = 1;
    }
    ipos = 2;
    iapos = 1;
    i__1 = kblk;
    for (iblk = 1; iblk <= i__1; ++iblk) {
	ncols = iw[ipos];
	nrows = iw[ipos + 1];
	j1 = ipos + 2;
	if (ncols > 0) {
	    goto L240;
	}
	ncols = -ncols;
	nrows = 1;
	--j1;
L240:
	io___52.ciunit = icntl[2];
	s_wsfe(&io___52);
	do_fio(&c__1, (char *)&iblk, (ftnlen)sizeof(integer));
	do_fio(&c__1, (char *)&nrows, (ftnlen)sizeof(integer));
	do_fio(&c__1, (char *)&ncols, (ftnlen)sizeof(integer));
	e_wsfe();
	j2 = j1 + ncols - 1;
	ipos = j2 + 1;
	io___54.ciunit = icntl[2];
	s_wsfe(&io___54);
	i__2 = j2;
	for (jj = j1; jj <= i__2; ++jj) {
	    do_fio(&c__1, (char *)&iw[jj], (ftnlen)sizeof(integer));
	}
	e_wsfe();
	io___56.ciunit = icntl[2];
	s_wsfe(&io___56);
	e_wsfe();
	len = ncols;
	i__2 = nrows;
	for (irows = 1; irows <= i__2; ++irows) {
	    j1 = iapos;
	    j2 = iapos + len - 1;
	    io___59.ciunit = icntl[2];
	    s_wsfe(&io___59);
	    i__3 = j2;
	    for (jj = j1; jj <= i__3; ++jj) {
		do_fio(&c__1, (char *)&a[jj], (ftnlen)sizeof(doublereal));
	    }
	    e_wsfe();
	    --len;
	    iapos = j2 + 1;
/* L290: */
	}
/* L300: */
    }
L310:
    return 0;
} /* ma27bd_ */

/* Subroutine */ int ma27cd_(integer *n, doublereal *a, integer *la, integer *
	iw, integer *liw, doublereal *w, integer *maxfrt, doublereal *rhs, 
	integer *iw1, integer *nsteps, integer *icntl, integer *info)
{
    /* Format strings */
    static char fmt_10[] = "(/,/,\002 ENTERING MA27CD WITH      N     LA    "
	    "LIW MAXFRT\002,\002  NSTEPS\002,/,21x,5i7)";
    static char fmt_30[] = "(\002 BLOCK PIVOT =\002,i8,\002 NROWS =\002,i8"
	    ",\002 NCOLS =\002,i8)";
    static char fmt_40[] = "(\002 COLUMN INDICES =\002,10i6,/,(17x,10i6))";
    static char fmt_50[] = "(\002 REAL ENTRIES .. EACH ROW STARTS ON A NEW L"
	    "INE\002)";
    static char fmt_60[] = "(1p,5d13.3)";
    static char fmt_100[] = "(\002 RHS\002,1p,5d13.3,/,(4x,1p,5d13.3))";
    static char fmt_160[] = "(/,/,\002 LEAVING MA27CD WITH\002)";

    /* System generated locals */
    integer i__1, i__2, i__3;

    /* Builtin functions */
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);

    /* Local variables */
    static integer i__, k, j1, j2, jj, len, iblk, kblk, nblk, ipos;
    extern /* Subroutine */ int ma27qd_(integer *, doublereal *, integer *, 
	    integer *, integer *, doublereal *, integer *, doublereal *, 
	    integer *, integer *, integer *, integer *), ma27rd_(integer *, 
	    doublereal *, integer *, integer *, integer *, doublereal *, 
	    integer *, doublereal *, integer *, integer *, integer *, integer 
	    *);
    static integer iapos, ncols, latop, irows, nrows;

    /* Fortran I/O blocks */
    static cilist io___60 = { 0, 0, 0, fmt_10, 0 };
    static cilist io___68 = { 0, 0, 0, fmt_30, 0 };
    static cilist io___70 = { 0, 0, 0, fmt_40, 0 };
    static cilist io___72 = { 0, 0, 0, fmt_50, 0 };
    static cilist io___75 = { 0, 0, 0, fmt_60, 0 };
    static cilist io___77 = { 0, 0, 0, fmt_100, 0 };
    static cilist io___81 = { 0, 0, 0, fmt_160, 0 };
    static cilist io___82 = { 0, 0, 0, fmt_100, 0 };


/* THIS SUBROUTINE USES THE FACTORISATION OF THE MATRIX IN A,IW TO */
/*     SOLVE A SYSTEM OF EQUATIONS. */
/* N MUST BE SET TO THE ORDER OF THE MATRIX.  IT IS NOT ALTERED. */
/* A,IW HOLD INFORMATION ON THE FACTORS AND MUST BE UNCHANGED SINCE */
/*     THE CALL TO MA27B/BD. THEY ARE NOT ALTERED BY MA27C/CDD. */
/* LA,LIW MUST BE SET TO THE LENGTHS OF A,IW RESPECTIVELY.  THEY */
/*     ARE NOT ALTERED. */
/* W USED AS A WORK ARRAY. */
/* MAXFRT IS THE LENGTH OF W AND MUST BE PASSED UNCHANGED FROM THE */
/*     CALL TO MA27B/BD.  IT IS NOT ALTERED. */
/* RHS MUST BE SET TO THE RIGHT HAND SIDE FOR THE EQUATIONS BEING */
/*     SOLVED.  ON EXIT, THIS ARRAY WILL HOLD THE SOLUTION. */
/* IW1 USED AS A WORK ARRAY. */
/* NSTEPS IS THE LENGTH OF IW1 WHICH MUST BE AT LEAST THE ABSOLUTE */
/*     VALUE OF IW(1).  IT IS NOT ALTERED. */
/* ICNTL is an INTEGER array of length 30, see MA27A/AD. */
/* INFO is an INTEGER array of length 20, see MA27A/AD. */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. External Subroutines .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */
    /* Parameter adjustments */
    --rhs;
    --a;
    --iw;
    --w;
    --iw1;
    --icntl;
    --info;

    /* Function Body */
    info[1] = 0;
    if (icntl[3] <= 0 || icntl[2] <= 0) {
	goto L110;
    }
/* PRINT INPUT PARAMETERS. */
    io___60.ciunit = icntl[2];
    s_wsfe(&io___60);
    do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&(*la), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&(*liw), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&(*maxfrt), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&(*nsteps), (ftnlen)sizeof(integer));
    e_wsfe();
/* PRINT OUT MATRIX FACTORS FROM MA27B/BD. */
    kblk = abs(iw[1]);
    if (kblk == 0) {
	goto L90;
    }
    if (icntl[3] == 1) {
	kblk = 1;
    }
    ipos = 2;
    iapos = 1;
    i__1 = kblk;
    for (iblk = 1; iblk <= i__1; ++iblk) {
	ncols = iw[ipos];
	nrows = iw[ipos + 1];
	j1 = ipos + 2;
	if (ncols > 0) {
	    goto L20;
	}
	ncols = -ncols;
	nrows = 1;
	--j1;
L20:
	io___68.ciunit = icntl[2];
	s_wsfe(&io___68);
	do_fio(&c__1, (char *)&iblk, (ftnlen)sizeof(integer));
	do_fio(&c__1, (char *)&nrows, (ftnlen)sizeof(integer));
	do_fio(&c__1, (char *)&ncols, (ftnlen)sizeof(integer));
	e_wsfe();
	j2 = j1 + ncols - 1;
	ipos = j2 + 1;
	io___70.ciunit = icntl[2];
	s_wsfe(&io___70);
	i__2 = j2;
	for (jj = j1; jj <= i__2; ++jj) {
	    do_fio(&c__1, (char *)&iw[jj], (ftnlen)sizeof(integer));
	}
	e_wsfe();
	io___72.ciunit = icntl[2];
	s_wsfe(&io___72);
	e_wsfe();
	len = ncols;
	i__2 = nrows;
	for (irows = 1; irows <= i__2; ++irows) {
	    j1 = iapos;
	    j2 = iapos + len - 1;
	    io___75.ciunit = icntl[2];
	    s_wsfe(&io___75);
	    i__3 = j2;
	    for (jj = j1; jj <= i__3; ++jj) {
		do_fio(&c__1, (char *)&a[jj], (ftnlen)sizeof(doublereal));
	    }
	    e_wsfe();
	    --len;
	    iapos = j2 + 1;
/* L70: */
	}
/* L80: */
    }
L90:
    k = min(10,*n);
    if (icntl[3] > 1) {
	k = *n;
    }
    if (*n > 0) {
	io___77.ciunit = icntl[2];
	s_wsfe(&io___77);
	i__1 = k;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    do_fio(&c__1, (char *)&rhs[i__], (ftnlen)sizeof(doublereal));
	}
	e_wsfe();
    }
L110:
    if (iw[1] < 0) {
	goto L130;
    }
    nblk = iw[1];
    if (nblk > 0) {
	goto L140;
    }
/* WE HAVE A ZERO MATRIX */
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	rhs[i__] = 0.;
/* L120: */
    }
    goto L150;
L130:
    nblk = -iw[1];
/* FORWARD SUBSTITUTION */
L140:
    i__1 = *liw - 1;
    ma27qd_(n, &a[1], la, &iw[2], &i__1, &w[1], maxfrt, &rhs[1], &iw1[1], &
	    nblk, &latop, &icntl[1]);
/* BACK SUBSTITUTION. */
    i__1 = *liw - 1;
    ma27rd_(n, &a[1], la, &iw[2], &i__1, &w[1], maxfrt, &rhs[1], &iw1[1], &
	    nblk, &latop, &icntl[1]);
L150:
    if (icntl[3] <= 0 || icntl[2] <= 0) {
	goto L170;
    }
/* PRINT OUTPUT PARAMETERS. */
    io___81.ciunit = icntl[2];
    s_wsfe(&io___81);
    e_wsfe();
    if (*n > 0) {
	io___82.ciunit = icntl[2];
	s_wsfe(&io___82);
	i__1 = k;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    do_fio(&c__1, (char *)&rhs[i__], (ftnlen)sizeof(doublereal));
	}
	e_wsfe();
    }
L170:
    return 0;
} /* ma27cd_ */

/* Subroutine */ int ma27gd_(integer *n, integer *nz, integer *irn, integer *
	icn, integer *iw, integer *lw, integer *ipe, integer *iq, integer *
	flag__, integer *iwfr, integer *icntl, integer *info)
{
    /* Format strings */
    static char fmt_60[] = "(\002 *** WARNING MESSAGE FROM SUBROUTINE MA27A"
	    "D\002,\002  *** INFO(1) =\002,i2)";
    static char fmt_70[] = "(i6,\002TH NON-ZERO (IN ROW\002,i6,\002 AND COLU"
	    "MN\002,i6,\002) IGNORED\002)";

    /* System generated locals */
    integer i__1, i__2;

    /* Builtin functions */
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);

    /* Local variables */
    static integer i__, j, k, l, k1, k2, n1, id, jn, lr, last, ndup;

    /* Fortran I/O blocks */
    static cilist io___87 = { 0, 0, 0, fmt_60, 0 };
    static cilist io___88 = { 0, 0, 0, fmt_70, 0 };



/* SORT PRIOR TO CALLING ANALYSIS ROUTINE MA27H/HD. */

/* GIVEN THE POSITIONS OF THE OFF-DIAGONAL NON-ZEROS OF A SYMMETRIC */
/*     MATRIX, CONSTRUCT THE SPARSITY PATTERN OF THE OFF-DIAGONAL */
/*     PART OF THE WHOLE MATRIX (UPPER AND LOWER TRIANGULAR PARTS). */
/* EITHER ONE OF A PAIR (I,J),(J,I) MAY BE USED TO REPRESENT */
/*     THE PAIR. DIAGONAL ELEMENTS AND DUPLICATES ARE IGNORED. */

/* N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED. */
/* NZ MUST BE SET TO THE NUMBER OF NON-ZEROS INPUT. IT IS NOT */
/*     ALTERED. */
/* IRN(I),I=1,2,...,NZ MUST BE SET TO THE ROW NUMBERS OF THE */
/*     NON-ZEROS ON INPUT. IT IS NOT ALTERED UNLESS IT IS EQUIVALENCED */
/*     TO IW (SEE DESCRIPTION OF IW). */
/* ICN(I),I=1,2,...,NZ MUST BE SET TO THE COLUMN NUMBERS OF THE */
/*     NON-ZEROS ON INPUT. IT IS NOT ALTERED UNLESS IT IS EQUIVALENCED */
/*     TO IW (SEE DESCRIPTION OF IW). */
/* IW NEED NOT BE SET ON INPUT. ON OUTPUT IT CONTAINS LISTS OF */
/*     COLUMN INDICES, EACH LIST BEING HEADED BY ITS LENGTH. */
/*     IRN(1) MAY BE EQUIVALENCED TO IW(1) AND ICN(1) MAY BE */
/*     EQUIVALENCED TO IW(K), WHERE K.GT.NZ. */
/* LW MUST BE SET TO THE LENGTH OF IW. IT MUST BE AT LEAST 2*NZ+N. */
/*     IT IS NOT ALTERED. */
/* IPE NEED NOT BE SET ON INPUT. ON OUTPUT IPE(I) POINTS TO THE START OF */
/*     THE ENTRY IN IW FOR ROW I, OR IS ZERO IF THERE IS NO ENTRY. */
/* IQ NEED NOT BE SET.  ON OUTPUT IQ(I),I=1,N CONTAINS THE NUMBER OF */
/*     OFF-DIAGONAL N0N-ZEROS IN ROW I INCLUDING DUPLICATES. */
/* FLAG IS USED FOR WORKSPACE TO HOLD FLAGS TO PERMIT DUPLICATE ENTRIES */
/*     TO BE IDENTIFIED QUICKLY. */
/* IWFR NEED NOT BE SET ON INPUT. ON OUTPUT IT POINTS TO THE FIRST */
/*     UNUSED LOCATION IN IW. */
/* ICNTL is an INTEGER array of length 30, see MA27A/AD. */
/* INFO is an INTEGER array of length 20, see MA27A/AD. */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. Executable Statements .. */

/* INITIALIZE INFO(2) AND COUNT IN IPE THE */
/*     NUMBERS OF NON-ZEROS IN THE ROWS AND MOVE ROW AND COLUMN */
/*     NUMBERS INTO IW. */
    /* Parameter adjustments */
    --flag__;
    --iq;
    --ipe;
    --irn;
    --icn;
    --iw;
    --icntl;
    --info;

    /* Function Body */
    info[2] = 0;
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	ipe[i__] = 0;
/* L10: */
    }
    lr = *nz;
    if (*nz == 0) {
	goto L120;
    }
    i__1 = *nz;
    for (k = 1; k <= i__1; ++k) {
	i__ = irn[k];
	j = icn[k];
	if (i__ < j) {
	    if (i__ >= 1 && j <= *n) {
		goto L90;
	    }
	} else if (i__ > j) {
	    if (j >= 1 && i__ <= *n) {
		goto L90;
	    }
	} else {
	    if (i__ >= 1 && i__ <= *n) {
		goto L80;
	    }
	}
	++info[2];
	info[1] = 1;
	if (info[2] <= 1 && icntl[2] > 0) {
	    io___87.ciunit = icntl[2];
	    s_wsfe(&io___87);
	    do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer));
	    e_wsfe();
	}
	if (info[2] <= 10 && icntl[2] > 0) {
	    io___88.ciunit = icntl[2];
	    s_wsfe(&io___88);
	    do_fio(&c__1, (char *)&k, (ftnlen)sizeof(integer));
	    do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(integer));
	    do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
	    e_wsfe();
	}
L80:
	i__ = 0;
	j = 0;
	goto L100;
L90:
	++ipe[i__];
	++ipe[j];
L100:
	iw[k] = j;
	++lr;
	iw[lr] = i__;
/* L110: */
    }

/* ACCUMULATE ROW COUNTS TO GET POINTERS TO ROW STARTS IN BOTH IPE AND IQ */
/*     AND INITIALIZE FLAG */
L120:
    iq[1] = 1;
    n1 = *n - 1;
    if (n1 <= 0) {
	goto L140;
    }
    i__1 = n1;
    for (i__ = 1; i__ <= i__1; ++i__) {
	flag__[i__] = 0;
	if (ipe[i__] == 0) {
	    ipe[i__] = -1;
	}
	iq[i__ + 1] = ipe[i__] + iq[i__] + 1;
	ipe[i__] = iq[i__];
/* L130: */
    }
L140:
    last = ipe[*n] + iq[*n];
    flag__[*n] = 0;
    if (lr >= last) {
	goto L160;
    }
    k1 = lr + 1;
    i__1 = last;
    for (k = k1; k <= i__1; ++k) {
	iw[k] = 0;
/* L150: */
    }
L160:
    ipe[*n] = iq[*n];
    *iwfr = last + 1;
    if (*nz == 0) {
	goto L230;
    }

/* RUN THROUGH PUTTING THE MATRIX ELEMENTS IN THE RIGHT PLACE */
/*     BUT WITH SIGNS INVERTED. IQ IS USED FOR HOLDING RUNNING POINTERS */
/*     AND IS LEFT HOLDING POINTERS TO ROW ENDS. */
    i__1 = *nz;
    for (k = 1; k <= i__1; ++k) {
	j = iw[k];
	if (j <= 0) {
	    goto L220;
	}
	l = k;
	iw[k] = 0;
	i__2 = *nz;
	for (id = 1; id <= i__2; ++id) {
	    if (l > *nz) {
		goto L170;
	    }
	    l += *nz;
	    goto L180;
L170:
	    l -= *nz;
L180:
	    i__ = iw[l];
	    iw[l] = 0;
	    if (i__ < j) {
		goto L190;
	    }
	    l = iq[j] + 1;
	    iq[j] = l;
	    jn = iw[l];
	    iw[l] = -i__;
	    goto L200;
L190:
	    l = iq[i__] + 1;
	    iq[i__] = l;
	    jn = iw[l];
	    iw[l] = -j;
L200:
	    j = jn;
	    if (j <= 0) {
		goto L220;
	    }
/* L210: */
	}
L220:
	;
    }

/* RUN THROUGH RESTORING SIGNS, REMOVING DUPLICATES AND SETTING THE */
/*     MATE OF EACH NON-ZERO. */
/* NDUP COUNTS THE NUMBER OF DUPLICATE ELEMENTS. */
L230:
    ndup = 0;
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	k1 = ipe[i__] + 1;
	k2 = iq[i__];
	if (k1 <= k2) {
	    goto L240;
	}
/* ROW IS EMPTY. SET POINTER TO ZERO. */
	ipe[i__] = 0;
	iq[i__] = 0;
	goto L280;
/* ON ENTRY TO THIS LOOP FLAG(J).LT.I FOR J=1,2,...,N. DURING THE LOOP */
/*     FLAG(J) IS SET TO I IF A NON-ZERO IN COLUMN J IS FOUND. THIS */
/*     PERMITS DUPLICATES TO BE RECOGNIZED QUICKLY. */
L240:
	i__2 = k2;
	for (k = k1; k <= i__2; ++k) {
	    j = -iw[k];
	    if (j <= 0) {
		goto L270;
	    }
	    l = iq[j] + 1;
	    iq[j] = l;
	    iw[l] = i__;
	    iw[k] = j;
	    if (flag__[j] != i__) {
		goto L250;
	    }
	    ++ndup;
	    iw[l] = 0;
	    iw[k] = 0;
L250:
	    flag__[j] = i__;
/* L260: */
	}
L270:
	iq[i__] -= ipe[i__];
	if (ndup == 0) {
	    iw[k1 - 1] = iq[i__];
	}
L280:
	;
    }
    if (ndup == 0) {
	goto L310;
    }

/* COMPRESS IW TO REMOVE DUMMY ENTRIES CAUSED BY DUPLICATES. */
    *iwfr = 1;
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	k1 = ipe[i__] + 1;
	if (k1 == 1) {
	    goto L300;
	}
	k2 = iq[i__] + ipe[i__];
	l = *iwfr;
	ipe[i__] = *iwfr;
	++(*iwfr);
	i__2 = k2;
	for (k = k1; k <= i__2; ++k) {
	    if (iw[k] == 0) {
		goto L290;
	    }
	    iw[*iwfr] = iw[k];
	    ++(*iwfr);
L290:
	    ;
	}
	iw[l] = *iwfr - l - 1;
L300:
	;
    }
L310:
    return 0;
} /* ma27gd_ */

/* Subroutine */ int ma27hd_(integer *n, integer *ipe, integer *iw, integer *
	lw, integer *iwfr, integer *nv, integer *nxt, integer *lst, integer *
	ipd, integer *flag__, integer *iovflo, integer *ncmpa, doublereal *
	fratio)
{
    /* System generated locals */
    integer i__1, i__2, i__3, i__4, i__5;
    doublereal d__1;

    /* Builtin functions */
    integer i_dnnt(doublereal *);

    /* Local variables */
    static integer i__, k, l, k1, k2, id, ie, ke, md, me, ip, jp, kp, is, js, 
	    ks, ln, ls, ml, ms, np, ns, jp1, jp2, kp0, kp1, kp2, np0, idl, 
	    idn, len, nel, nflg, lwfr, root;
    extern /* Subroutine */ int ma27ud_(integer *, integer *, integer *, 
	    integer *, integer *, integer *);
    static integer limit, nvpiv, nvroot;


/* ANALYSIS SUBROUTINE */

/* GIVEN REPRESENTATION OF THE WHOLE MATRIX (EXCLUDING DIAGONAL) */
/*     PERFORM MINIMUM DEGREE ORDERING, CONSTRUCTING TREE POINTERS. */
/*     IT WORKS WITH SUPERVARIABLES WHICH ARE COLLECTIONS OF ONE OR MORE */
/*     VARIABLES, STARTING WITH SUPERVARIABLE I CONTAINING VARIABLE I FOR */
/*     I = 1,2,...,N. ALL VARIABLES IN A SUPERVARIABLE ARE ELIMINATED */
/*     TOGETHER. EACH SUPERVARIABLE HAS AS NUMERICAL NAME THAT OF ONE */
/*     OF ITS VARIABLES (ITS PRINCIPAL VARIABLE). */

/* N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED. */
/* IPE(I) MUST BE SET TO POINT TO THE POSITION IN IW OF THE */
/*     START OF ROW I OR HAVE THE VALUE ZERO IF ROW I HAS NO OFF- */
/*     DIAGONAL NON-ZEROS. DURING EXECUTION IT IS USED AS FOLLOWS. IF */
/*     SUPERVARIABLE I IS ABSORBED INTO SUPERVARIABLE J THEN IPE(I)=-J. */
/*     IF SUPERVARIABLE I IS ELIMINATED THEN IPE(I) EITHER POINTS TO THE */
/*     LIST OF SUPERVARIABLES FOR CREATED ELEMENT I OR IS ZERO IF */
/*     THE CREATED ELEMENT IS NULL. IF ELEMENT I */
/*     IS ABSORBED INTO ELEMENT J THEN IPE(I)=-J. */
/* IW MUST BE SET ON ENTRY TO HOLD LISTS OF VARIABLES BY */
/*     ROWS, EACH LIST BEING HEADED BY ITS LENGTH. */
/*     DURING EXECUTION THESE LISTS ARE REVISED AND HOLD */
/*     LISTS OF ELEMENTS AND SUPERVARIABLES. THE ELEMENTS */
/*     ALWAYS HEAD THE LISTS. WHEN A SUPERVARIABLE */
/*     IS ELIMINATED ITS LIST IS REPLACED BY A LIST OF SUPERVARIABLES */
/*     IN THE NEW ELEMENT. */
/* LW MUST BE SET TO THE LENGTH OF IW. IT IS NOT ALTERED. */
/* IWFR MUST BE SET TO THE POSITION IN IW OF THE FIRST FREE VARIABLE. */
/*     IT IS REVISED DURING EXECUTION AND CONTINUES TO HAVE THIS MEANING. */
/* NV(JS) NEED NOT BE SET. DURING EXECUTION IT IS ZERO IF */
/*     JS IS NOT A PRINCIPAL VARIABLE AND IF IT IS IT HOLDS */
/*     THE NUMBER OF VARIABLES IN SUPERVARIABLE JS. FOR ELIMINATED */
/*     VARIABLES IT IS SET TO THE DEGREE AT THE TIME OF ELIMINATION. */
/* NXT(JS) NEED NOT BE SET. DURING EXECUTION IT IS THE NEXT */
/*     SUPERVARIABLE HAVING THE SAME DEGREE AS JS, OR ZERO */
/*     IF IT IS LAST IN ITS LIST. */
/* LST(JS) NEED NOT BE SET. DURING EXECUTION IT IS THE */
/*     LAST SUPERVARIABLE HAVING THE SAME DEGREE AS JS OR */
/*     -(ITS DEGREE) IF IT IS FIRST IN ITS LIST. */
/* IPD(ID) NEED NOT BE SET. DURING EXECUTION IT */
/*     IS THE FIRST SUPERVARIABLE WITH DEGREE ID OR ZERO */
/*     IF THERE ARE NONE. */
/* FLAG IS USED AS WORKSPACE FOR ELEMENT AND SUPERVARIABLE FLAGS. */
/*     WHILE THE CODE IS FINDING THE DEGREE OF SUPERVARIABLE IS */
/*     FLAG HAS THE FOLLOWING VALUES. */
/*     A) FOR THE CURRENT PIVOT/NEW ELEMENT ME */
/*           FLAG(ME)=-1 */
/*     B) FOR VARIABLES JS */
/*           FLAG(JS)=-1 IF JS IS NOT A PRINCIPAL VARIABLE */
/*           FLAG(JS)=0 IF JS IS A SUPERVARIABLE IN THE NEW ELEMENT */
/*           FLAG(JS)=NFLG IF JS IS A SUPERVARIABLE NOT IN THE NEW */
/*                 ELEMENT THAT HAS BEEN COUNTED IN THE DEGREE */
/*                 CALCULATION */
/*           FLAG(JS).GT.NFLG IF JS IS A SUPERVARIABLE NOT IN THE NEW */
/*                 ELEMENT THAT HAS NOT BEEN COUNTED IN THE DEGREE */
/*                 CALCULATION */
/*     C) FOR ELEMENTS IE */
/*           FLAG(IE)=-1 IF ELEMENT IE HAS BEEN MERGED INTO ANOTHER */
/*           FLAG(IE)=-NFLG IF ELEMENT IE HAS BEEN USED IN THE DEGREE */
/*                 CALCULATION FOR IS. */
/*           FLAG(IE).LT.-NFLG IF ELEMENT IE HAS NOT BEEN USED IN THE */
/*                 DEGREE CALCULATION FOR IS */
/* IOVFLO see ICNTL(4) in MA27A/AD. */
/* NCMPA see INFO(11) in MA27A/AD. */
/* FRATIO see CNTL(2) in MA27A/AD. */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */
/*     .. Local Scalars .. */
/* LIMIT  Limit on number of variables for putting node in root. */
/* NVROOT Number of variables in the root node */
/* ROOT   Index of the root node (N+1 if none chosen yet). */
/*     .. */
/*     .. External Subroutines .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/* If a column of the reduced matrix has relative density greater than */
/* CNTL(2), it is forced into the root. All such columns are taken to */
/* have sparsity pattern equal to their merged patterns, so the fill */
/* and operation counts may be overestimated. */

/* IS,JS,KS,LS,MS,NS ARE USED TO REFER TO SUPERVARIABLES. */
/* IE,JE,KE ARE USED TO REFER TO ELEMENTS. */
/* IP,JP,KP,K,NP ARE USED TO POINT TO LISTS OF ELEMENTS. */
/*     OR SUPERVARIABLES. */
/* ID IS USED FOR THE DEGREE OF A SUPERVARIABLE. */
/* MD IS USED FOR THE CURRENT MINIMUM DEGREE. */
/* IDN IS USED FOR THE NO. OF VARIABLES IN A NEWLY CREATED ELEMENT */
/* NEL IS USED TO HOLD THE NO. OF VARIABLES THAT HAVE BEEN */
/*     ELIMINATED. */
/* ME=MS IS THE NAME OF THE SUPERVARIABLE ELIMINATED AND */
/*     OF THE ELEMENT CREATED IN THE MAIN LOOP. */
/* NFLG IS USED FOR THE CURRENT FLAG VALUE IN ARRAY FLAG. IT STARTS */
/*     WITH THE VALUE IOVFLO AND IS REDUCED BY 1 EACH TIME IT IS USED */
/*     UNTIL IT HAS THE VALUE 2 WHEN IT IS RESET TO THE VALUE IOVFLO. */

/*     .. Executable Statements .. */
/* INITIALIZATIONS */
    /* Parameter adjustments */
    --flag__;
    --ipd;
    --lst;
    --nxt;
    --nv;
    --ipe;
    --iw;

    /* Function Body */
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	ipd[i__] = 0;
	nv[i__] = 1;
	flag__[i__] = *iovflo;
/* L10: */
    }
    md = 1;
    *ncmpa = 0;
    nflg = *iovflo;
    nel = 0;
    root = *n + 1;
    nvroot = 0;

/* LINK TOGETHER VARIABLES HAVING SAME DEGREE */
    i__1 = *n;
    for (is = 1; is <= i__1; ++is) {
	k = ipe[is];
	if (k <= 0) {
	    goto L20;
	}
	id = iw[k] + 1;
	ns = ipd[id];
	if (ns > 0) {
	    lst[ns] = is;
	}
	nxt[is] = ns;
	ipd[id] = is;
	lst[is] = -id;
	goto L30;
/* WE HAVE A VARIABLE THAT CAN BE ELIMINATED AT ONCE BECAUSE THERE IS */
/*     NO OFF-DIAGONAL NON-ZERO IN ITS ROW. */
L20:
	++nel;
	flag__[is] = -1;
	nxt[is] = 0;
	lst[is] = 0;
L30:
	;
    }

/* START OF MAIN LOOP */

    i__1 = *n;
    for (ml = 1; ml <= i__1; ++ml) {
/* LEAVE LOOP IF ALL VARIABLES HAVE BEEN ELIMINATED. */
	if (nel + nvroot + 1 >= *n) {
	    goto L350;
	}

/* FIND NEXT SUPERVARIABLE FOR ELIMINATION. */
	i__2 = *n;
	for (id = md; id <= i__2; ++id) {
	    ms = ipd[id];
	    if (ms > 0) {
		goto L50;
	    }
/* L40: */
	}
L50:
	md = id;
/* NVPIV HOLDS THE NUMBER OF VARIABLES IN THE PIVOT. */
	nvpiv = nv[ms];

/* REMOVE CHOSEN VARIABLE FROM LINKED LIST */
	ns = nxt[ms];
	nxt[ms] = 0;
	lst[ms] = 0;
	if (ns > 0) {
	    lst[ns] = -id;
	}
	ipd[id] = ns;
	me = ms;
	nel += nvpiv;
/* IDN HOLDS THE DEGREE OF THE NEW ELEMENT. */
	idn = 0;

/* RUN THROUGH THE LIST OF THE PIVOTAL SUPERVARIABLE, SETTING TREE */
/*     POINTERS AND CONSTRUCTING NEW LIST OF SUPERVARIABLES. */
/* KP IS A POINTER TO THE CURRENT POSITION IN THE OLD LIST. */
	kp = ipe[me];
	flag__[ms] = -1;
/* IP POINTS TO THE START OF THE NEW LIST. */
	ip = *iwfr;
/* LEN HOLDS THE LENGTH OF THE LIST ASSOCIATED WITH THE PIVOT. */
	len = iw[kp];
	i__2 = len;
	for (kp1 = 1; kp1 <= i__2; ++kp1) {
	    ++kp;
	    ke = iw[kp];
/* JUMP IF KE IS AN ELEMENT THAT HAS NOT BEEN MERGED INTO ANOTHER. */
	    if (flag__[ke] <= -2) {
		goto L60;
	    }
/* JUMP IF KE IS AN ELEMENT THAT HAS BEEN MERGED INTO ANOTHER OR IS */
/*     A SUPERVARIABLE THAT HAS BEEN ELIMINATED. */
	    if (flag__[ke] <= 0) {
		if (ipe[ke] != -root) {
		    goto L140;
		}
/* KE has been merged into the root */
		ke = root;
		if (flag__[ke] <= 0) {
		    goto L140;
		}
	    }
/* WE HAVE A SUPERVARIABLE. PREPARE TO SEARCH REST OF LIST. */
	    jp = kp - 1;
	    ln = len - kp1 + 1;
	    ie = ms;
	    goto L70;
/* SEARCH VARIABLE LIST OF ELEMENT KE, USING JP AS A POINTER TO IT. */
L60:
	    ie = ke;
	    jp = ipe[ie];
	    ln = iw[jp];

/* SEARCH FOR DIFFERENT SUPERVARIABLES AND ADD THEM TO THE NEW LIST, */
/*     COMPRESSING WHEN NECESSARY. THIS LOOP IS EXECUTED ONCE FOR */
/*     EACH ELEMENT IN THE LIST AND ONCE FOR ALL THE SUPERVARIABLES */
/*     IN THE LIST. */
L70:
	    i__3 = ln;
	    for (jp1 = 1; jp1 <= i__3; ++jp1) {
		++jp;
		is = iw[jp];
/* JUMP IF IS IS NOT A PRINCIPAL VARIABLE OR HAS ALREADY BEEN COUNTED. */
		if (flag__[is] <= 0) {
		    if (ipe[is] == -root) {
/* IS has been merged into the root */
			is = root;
			iw[jp] = root;
			if (flag__[is] <= 0) {
			    goto L130;
			}
		    } else {
			goto L130;
		    }
		}
		flag__[is] = 0;
		if (*iwfr < *lw) {
		    goto L100;
		}
/* PREPARE FOR COMPRESSING IW BY ADJUSTING POINTERS AND */
/*     LENGTHS SO THAT THE LISTS BEING SEARCHED IN THE INNER AND OUTER */
/*     LOOPS CONTAIN ONLY THE REMAINING ENTRIES. */
		ipe[ms] = kp;
		iw[kp] = len - kp1;
		ipe[ie] = jp;
		iw[jp] = ln - jp1;
/* COMPRESS IW */
		i__4 = ip - 1;
		ma27ud_(n, &ipe[1], &iw[1], &i__4, &lwfr, ncmpa);
/* COPY NEW LIST FORWARD */
		jp2 = *iwfr - 1;
		*iwfr = lwfr;
		if (ip > jp2) {
		    goto L90;
		}
		i__4 = jp2;
		for (jp = ip; jp <= i__4; ++jp) {
		    iw[*iwfr] = iw[jp];
		    ++(*iwfr);
/* L80: */
		}
/* ADJUST POINTERS FOR THE NEW LIST AND THE LISTS BEING SEARCHED. */
L90:
		ip = lwfr;
		jp = ipe[ie];
		kp = ipe[me];
/* STORE IS IN NEW LIST. */
L100:
		iw[*iwfr] = is;
		idn += nv[is];
		++(*iwfr);
/* REMOVE IS FROM DEGREE LINKED LIST */
		ls = lst[is];
		lst[is] = 0;
		ns = nxt[is];
		nxt[is] = 0;
		if (ns > 0) {
		    lst[ns] = ls;
		}
		if (ls < 0) {
		    ls = -ls;
		    ipd[ls] = ns;
		} else if (ls > 0) {
		    nxt[ls] = ns;
		}
L130:
		;
	    }
/* JUMP IF WE HAVE JUST BEEN SEARCHING THE VARIABLES AT THE END OF */
/*     THE LIST OF THE PIVOT. */
	    if (ie == ms) {
		goto L150;
	    }
/* SET TREE POINTER AND FLAG TO INDICATE ELEMENT IE IS ABSORBED INTO */
/*     NEW ELEMENT ME. */
	    ipe[ie] = -me;
	    flag__[ie] = -1;
L140:
	    ;
	}
/* STORE THE DEGREE OF THE PIVOT. */
L150:
	nv[ms] = idn + nvpiv;
/* JUMP IF NEW ELEMENT IS NULL. */
	if (*iwfr == ip) {
	    goto L330;
	}
	k1 = ip;
	k2 = *iwfr - 1;

/* RUN THROUGH NEW LIST OF SUPERVARIABLES REVISING EACH ASSOCIATED LIST, */
/*     RECALCULATING DEGREES AND REMOVING DUPLICATES. */
	d__1 = *fratio * (*n - nel);
	limit = i_dnnt(&d__1);
	i__2 = k2;
	for (k = k1; k <= i__2; ++k) {
	    is = iw[k];
	    if (is == root) {
		goto L310;
	    }
	    if (nflg > 2) {
		goto L170;
	    }
/* RESET FLAG VALUES TO +/-IOVFLO. */
	    i__3 = *n;
	    for (i__ = 1; i__ <= i__3; ++i__) {
		if (flag__[i__] > 0) {
		    flag__[i__] = *iovflo;
		}
		if (flag__[i__] <= -2) {
		    flag__[i__] = -(*iovflo);
		}
/* L160: */
	    }
	    nflg = *iovflo;
/* REDUCE NFLG BY ONE TO CATER FOR THIS SUPERVARIABLE. */
L170:
	    --nflg;
/* BEGIN WITH THE DEGREE OF THE NEW ELEMENT. ITS VARIABLES MUST ALWAYS */
/*     BE COUNTED DURING THE DEGREE CALCULATION AND THEY ARE ALREADY */
/*     FLAGGED WITH THE VALUE 0. */
	    id = idn;
/* RUN THROUGH THE LIST ASSOCIATED WITH SUPERVARIABLE IS */
	    kp1 = ipe[is] + 1;
/* NP POINTS TO THE NEXT ENTRY IN THE REVISED LIST. */
	    np = kp1;
	    kp2 = iw[kp1 - 1] + kp1 - 1;
	    i__3 = kp2;
	    for (kp = kp1; kp <= i__3; ++kp) {
		ke = iw[kp];
/* TEST WHETHER KE IS AN ELEMENT, A REDUNDANT ENTRY OR A SUPERVARIABLE. */
		if (flag__[ke] == -1) {
		    if (ipe[ke] != -root) {
			goto L220;
		    }
/* KE has been merged into the root */
		    ke = root;
		    iw[kp] = root;
		    if (flag__[ke] == -1) {
			goto L220;
		    }
		}
		if (flag__[ke] >= 0) {
		    goto L230;
		}
/* SEARCH LIST OF ELEMENT KE, REVISING THE DEGREE WHEN NEW VARIABLES */
/*     FOUND. */
		jp1 = ipe[ke] + 1;
		jp2 = iw[jp1 - 1] + jp1 - 1;
		idl = id;
		i__4 = jp2;
		for (jp = jp1; jp <= i__4; ++jp) {
		    js = iw[jp];
/* JUMP IF JS HAS ALREADY BEEN COUNTED. */
		    if (flag__[js] <= nflg) {
			goto L190;
		    }
		    id += nv[js];
		    flag__[js] = nflg;
L190:
		    ;
		}
/* JUMP IF ONE OR MORE NEW SUPERVARIABLES WERE FOUND. */
		if (id > idl) {
		    goto L210;
		}
/* CHECK WHETHER EVERY VARIABLE OF ELEMENT KE IS IN NEW ELEMENT ME. */
		i__4 = jp2;
		for (jp = jp1; jp <= i__4; ++jp) {
		    js = iw[jp];
		    if (flag__[js] != 0) {
			goto L210;
		    }
/* L200: */
		}
/* SET TREE POINTER AND FLAG TO INDICATE THAT ELEMENT KE IS ABSORBED */
/*     INTO NEW ELEMENT ME. */
		ipe[ke] = -me;
		flag__[ke] = -1;
		goto L220;
/* STORE ELEMENT KE IN THE REVISED LIST FOR SUPERVARIABLE IS AND FLAG IT. */
L210:
		iw[np] = ke;
		flag__[ke] = -nflg;
		++np;
L220:
		;
	    }
	    np0 = np;
	    goto L250;
/* TREAT THE REST OF THE LIST ASSOCIATED WITH SUPERVARIABLE IS. IT */
/*     CONSISTS ENTIRELY OF SUPERVARIABLES. */
L230:
	    kp0 = kp;
	    np0 = np;
	    i__3 = kp2;
	    for (kp = kp0; kp <= i__3; ++kp) {
		ks = iw[kp];
		if (flag__[ks] <= nflg) {
		    if (ipe[ks] == -root) {
			ks = root;
			iw[kp] = root;
			if (flag__[ks] <= nflg) {
			    goto L240;
			}
		    } else {
			goto L240;
		    }
		}
/* ADD TO DEGREE, FLAG SUPERVARIABLE KS AND ADD IT TO NEW LIST. */
		id += nv[ks];
		flag__[ks] = nflg;
		iw[np] = ks;
		++np;
L240:
		;
	    }
/* MOVE FIRST SUPERVARIABLE TO END OF LIST, MOVE FIRST ELEMENT TO END */
/*     OF ELEMENT PART OF LIST AND ADD NEW ELEMENT TO FRONT OF LIST. */
L250:
	    if (id >= limit) {
		goto L295;
	    }
	    iw[np] = iw[np0];
	    iw[np0] = iw[kp1];
	    iw[kp1] = me;
/* STORE THE NEW LENGTH OF THE LIST. */
	    iw[kp1 - 1] = np - kp1 + 1;

/* CHECK WHETHER ROW IS IS IDENTICAL TO ANOTHER BY LOOKING IN LINKED */
/*     LIST OF SUPERVARIABLES WITH DEGREE ID AT THOSE WHOSE LISTS HAVE */
/*     FIRST ENTRY ME. NOTE THAT THOSE CONTAINING ME COME FIRST SO THE */
/*     SEARCH CAN BE TERMINATED WHEN A LIST NOT STARTING WITH ME IS */
/*     FOUND. */
	    js = ipd[id];
	    i__3 = *n;
	    for (l = 1; l <= i__3; ++l) {
		if (js <= 0) {
		    goto L300;
		}
		kp1 = ipe[js] + 1;
		if (iw[kp1] != me) {
		    goto L300;
		}
/* JS HAS SAME DEGREE AND IS ACTIVE. CHECK IF IDENTICAL TO IS. */
		kp2 = kp1 - 1 + iw[kp1 - 1];
		i__4 = kp2;
		for (kp = kp1; kp <= i__4; ++kp) {
		    ie = iw[kp];
/* JUMP IF IE IS A SUPERVARIABLE OR AN ELEMENT NOT IN THE LIST OF IS. */
		    if ((i__5 = flag__[ie], abs(i__5)) > nflg) {
			goto L270;
		    }
/* L260: */
		}
		goto L290;
L270:
		js = nxt[js];
/* L280: */
	    }
/* SUPERVARIABLE AMALGAMATION. ROW IS IS IDENTICAL TO ROW JS. */
/* REGARD ALL VARIABLES IN THE TWO SUPERVARIABLES AS BEING IN IS. SET */
/*     TREE POINTER, FLAG AND NV ENTRIES. */
L290:
	    ipe[js] = -is;
	    nv[is] += nv[js];
	    nv[js] = 0;
	    flag__[js] = -1;
/* REPLACE JS BY IS IN LINKED LIST. */
	    ns = nxt[js];
	    ls = lst[js];
	    if (ns > 0) {
		lst[ns] = is;
	    }
	    if (ls > 0) {
		nxt[ls] = is;
	    }
	    lst[is] = ls;
	    nxt[is] = ns;
	    lst[js] = 0;
	    nxt[js] = 0;
	    if (ipd[id] == js) {
		ipd[id] = is;
	    }
	    goto L310;
/* Treat IS as full. Merge it into the root node. */
L295:
	    if (nvroot == 0) {
		root = is;
		ipe[is] = 0;
	    } else {
		iw[k] = root;
		ipe[is] = -root;
		nv[root] += nv[is];
		nv[is] = 0;
		flag__[is] = -1;
	    }
	    nvroot = nv[root];
	    goto L310;
/* INSERT IS INTO LINKED LIST OF SUPERVARIABLES OF SAME DEGREE. */
L300:
	    ns = ipd[id];
	    if (ns > 0) {
		lst[ns] = is;
	    }
	    nxt[is] = ns;
	    ipd[id] = is;
	    lst[is] = -id;
	    md = min(md,id);
L310:
	    ;
	}

/* RESET FLAGS FOR SUPERVARIABLES IN NEWLY CREATED ELEMENT AND */
/*     REMOVE THOSE ABSORBED INTO OTHERS. */
	i__2 = k2;
	for (k = k1; k <= i__2; ++k) {
	    is = iw[k];
	    if (nv[is] == 0) {
		goto L320;
	    }
	    flag__[is] = nflg;
	    iw[ip] = is;
	    ++ip;
L320:
	    ;
	}
	*iwfr = k1;
	flag__[me] = -nflg;
/* MOVE FIRST ENTRY TO END TO MAKE ROOM FOR LENGTH. */
	iw[ip] = iw[k1];
	iw[k1] = ip - k1;
/* SET POINTER FOR NEW ELEMENT AND RESET IWFR. */
	ipe[me] = k1;
	*iwfr = ip + 1;
	goto L335;
L330:
	ipe[me] = 0;

L335:
/* L340: */
	;
    }

/* Absorb any remaining variables into the root */
L350:
    i__1 = *n;
    for (is = 1; is <= i__1; ++is) {
	if (nxt[is] != 0 || lst[is] != 0) {
	    if (nvroot == 0) {
		root = is;
		ipe[is] = 0;
	    } else {
		ipe[is] = -root;
	    }
	    nvroot += nv[is];
	    nv[is] = 0;
	}
/* L360: */
    }
/* Link any remaining elements to the root */
    i__1 = *n;
    for (ie = 1; ie <= i__1; ++ie) {
	if (ipe[ie] > 0) {
	    ipe[ie] = -root;
	}
/* L370: */
    }
    if (nvroot > 0) {
	nv[root] = nvroot;
    }
    return 0;
} /* ma27hd_ */

/* Subroutine */ int ma27ud_(integer *n, integer *ipe, integer *iw, integer *
	lw, integer *iwfr, integer *ncmpa)
{
    /* System generated locals */
    integer i__1, i__2;

    /* Local variables */
    static integer i__, k, k1, k2, ir, lwfr;

/* COMPRESS LISTS HELD BY MA27H/HD AND MA27K/KD IN IW AND ADJUST POINTERS */
/*     IN IPE TO CORRESPOND. */
/* N IS THE MATRIX ORDER. IT IS NOT ALTERED. */
/* IPE(I) POINTS TO THE POSITION IN IW OF THE START OF LIST I OR IS */
/*     ZERO IF THERE IS NO LIST I. ON EXIT IT POINTS TO THE NEW POSITION. */
/* IW HOLDS THE LISTS, EACH HEADED BY ITS LENGTH. ON OUTPUT THE SAME */
/*     LISTS ARE HELD, BUT THEY ARE NOW COMPRESSED TOGETHER. */
/* LW HOLDS THE LENGTH OF IW. IT IS NOT ALTERED. */
/* IWFR NEED NOT BE SET ON ENTRY. ON EXIT IT POINTS TO THE FIRST FREE */
/*     LOCATION IN IW. */
/*     ON RETURN IT IS SET TO THE FIRST FREE LOCATION IN IW. */
/* NCMPA see INFO(11) in MA27A/AD. */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. Executable Statements .. */
    /* Parameter adjustments */
    --ipe;
    --iw;

    /* Function Body */
    ++(*ncmpa);
/* PREPARE FOR COMPRESSING BY STORING THE LENGTHS OF THE */
/*     LISTS IN IPE AND SETTING THE FIRST ENTRY OF EACH LIST TO */
/*     -(LIST NUMBER). */
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	k1 = ipe[i__];
	if (k1 <= 0) {
	    goto L10;
	}
	ipe[i__] = iw[k1];
	iw[k1] = -i__;
L10:
	;
    }

/* COMPRESS */
/* IWFR POINTS JUST BEYOND THE END OF THE COMPRESSED FILE. */
/* LWFR POINTS JUST BEYOND THE END OF THE UNCOMPRESSED FILE. */
    *iwfr = 1;
    lwfr = *iwfr;
    i__1 = *n;
    for (ir = 1; ir <= i__1; ++ir) {
	if (lwfr > *lw) {
	    goto L70;
	}
/* SEARCH FOR THE NEXT NEGATIVE ENTRY. */
	i__2 = *lw;
	for (k = lwfr; k <= i__2; ++k) {
	    if (iw[k] < 0) {
		goto L30;
	    }
/* L20: */
	}
	goto L70;
/* PICK UP ENTRY NUMBER, STORE LENGTH IN NEW POSITION, SET NEW POINTER */
/*     AND PREPARE TO COPY LIST. */
L30:
	i__ = -iw[k];
	iw[*iwfr] = ipe[i__];
	ipe[i__] = *iwfr;
	k1 = k + 1;
	k2 = k + iw[*iwfr];
	++(*iwfr);
	if (k1 > k2) {
	    goto L50;
	}
/* COPY LIST TO NEW POSITION. */
	i__2 = k2;
	for (k = k1; k <= i__2; ++k) {
	    iw[*iwfr] = iw[k];
	    ++(*iwfr);
/* L40: */
	}
L50:
	lwfr = k2 + 1;
/* L60: */
    }
L70:
    return 0;
} /* ma27ud_ */

/* Subroutine */ int ma27jd_(integer *n, integer *nz, integer *irn, integer *
	icn, integer *perm, integer *iw, integer *lw, integer *ipe, integer *
	iq, integer *flag__, integer *iwfr, integer *icntl, integer *info)
{
    /* Format strings */
    static char fmt_60[] = "(\002 *** WARNING MESSAGE FROM SUBROUTINE MA27A"
	    "D\002,\002  *** INFO(1) =\002,i2)";
    static char fmt_70[] = "(i6,\002TH NON-ZERO (IN ROW\002,i6,\002 AND COLU"
	    "MN\002,i6,\002) IGNORED\002)";

    /* System generated locals */
    integer i__1, i__2;

    /* Builtin functions */
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);

    /* Local variables */
    static integer i__, j, k, l, k1, k2, id, in, len, lbig, jdummy;

    /* Fortran I/O blocks */
    static cilist io___144 = { 0, 0, 0, fmt_60, 0 };
    static cilist io___145 = { 0, 0, 0, fmt_70, 0 };



/* SORT PRIOR TO CALLING ANALYSIS ROUTINE MA27K/KD. */

/* GIVEN THE POSITIONS OF THE OFF-DIAGONAL NON-ZEROS OF A SYMMETRIC */
/*     MATRIX AND A PERMUTATION, CONSTRUCT THE SPARSITY PATTERN */
/*     OF THE STRICTLY UPPER TRIANGULAR PART OF THE PERMUTED MATRIX. */
/*     EITHER ONE OF A PAIR (I,J),(J,I) MAY BE USED TO REPRESENT */
/*     THE PAIR. DIAGONAL ELEMENTS ARE IGNORED. NO CHECK IS MADE */
/*     FOR DUPLICATE ELEMENTS UNLESS ANY ROW HAS MORE THAN ICNTL(4) */
/*     NON-ZEROS, IN WHICH CASE DUPLICATES ARE REMOVED. */

/* N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED. */
/* NZ MUST BE SET TO THE NUMBER OF NON-ZEROS INPUT. IT IS NOT */
/*     ALTERED. */
/* IRN(I),I=1,2,...,NZ MUST BE SET TO THE ROW INDICES OF THE */
/*     NON-ZEROS ON INPUT. IT IS NOT ALTERED UNLESS EQUIVALENCED WITH IW. */
/*     IRN(1) MAY BE EQUIVALENCED WITH IW(1). */
/* ICN(I),I=1,2,...,NZ MUST BE SET TO THE COLUMN INDICES OF THE */
/*     NON-ZEROS ON INPUT. IT IS NOT ALTERED UNLESS EQUIVALENCED */
/*     WITH IW.ICN(1) MAY BE EQUIVELENCED WITH IW(K),K.GT.NZ. */
/* PERM(I) MUST BE SET TO HOLD THE POSITION OF VARIABLE I IN THE */
/*     PERMUTED ORDER. IT IS NOT ALTERED. */
/* IW NEED NOT BE SET ON INPUT. ON OUTPUT IT CONTAINS LISTS OF */
/*     COLUMN NUMBERS, EACH LIST BEING HEADED BY ITS LENGTH. */
/* LW MUST BE SET TO THE LENGTH OF IW. IT MUST BE AT LEAST */
/*     MAX(NZ,N+(NO. OF OFF-DIAGONAL NON-ZEROS)). IT IS NOT ALTERED. */
/* IPE NEED NOT BE SET ON INPUT. ON OUTPUT IPE(I) POINTS TO THE START OF */
/*     THE ENTRY IN IW FOR ROW I, OR IS ZERO IF THERE IS NO ENTRY. */
/* IQ NEED NOT BE SET. ON OUTPUT IQ(I) CONTAINS THE NUMBER OF */
/*     OFF-DIAGONAL NON-ZEROS IN ROW I, INCLUDING DUPLICATES. */
/* FLAG IS USED FOR WORKSPACE TO HOLD FLAGS TO PERMIT DUPLICATE */
/*     ENTRIES TO BE IDENTIFIED QUICKLY. */
/* IWFR NEED NOT BE SET ON INPUT. ON OUTPUT IT POINTS TO THE FIRST */
/*     UNUSED LOCATION IN IW. */
/* ICNTL is an INTEGER array of length 30, see MA27A/AD. */
/* INFO is an INTEGER array of length 20, see MA27A/AD. */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */

/* INITIALIZE INFO(1), INFO(2) AND IQ */
    /* Parameter adjustments */
    --flag__;
    --iq;
    --ipe;
    --perm;
    --irn;
    --icn;
    --iw;
    --icntl;
    --info;

    /* Function Body */
    info[1] = 0;
    info[2] = 0;
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	iq[i__] = 0;
/* L10: */
    }

/* COUNT THE NUMBERS OF NON-ZEROS IN THE ROWS, PRINT WARNINGS ABOUT */
/*     OUT-OF-RANGE INDICES AND TRANSFER GENUINE ROW NUMBERS */
/*     (NEGATED) INTO IW. */
    if (*nz == 0) {
	goto L110;
    }
    i__1 = *nz;
    for (k = 1; k <= i__1; ++k) {
	i__ = irn[k];
	j = icn[k];
	iw[k] = -i__;
	if (i__ < j) {
	    if (i__ >= 1 && j <= *n) {
		goto L80;
	    }
	} else if (i__ > j) {
	    if (j >= 1 && i__ <= *n) {
		goto L80;
	    }
	} else {
	    iw[k] = 0;
	    if (i__ >= 1 && i__ <= *n) {
		goto L100;
	    }
	}
	++info[2];
	info[1] = 1;
	iw[k] = 0;
	if (info[2] <= 1 && icntl[2] > 0) {
	    io___144.ciunit = icntl[2];
	    s_wsfe(&io___144);
	    do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer));
	    e_wsfe();
	}
	if (info[2] <= 10 && icntl[2] > 0) {
	    io___145.ciunit = icntl[2];
	    s_wsfe(&io___145);
	    do_fio(&c__1, (char *)&k, (ftnlen)sizeof(integer));
	    do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(integer));
	    do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
	    e_wsfe();
	}
	goto L100;
L80:
	if (perm[j] > perm[i__]) {
	    goto L90;
	}
	++iq[j];
	goto L100;
L90:
	++iq[i__];
L100:
	;
    }

/* ACCUMULATE ROW COUNTS TO GET POINTERS TO ROW ENDS */
/*     IN IPE. */
L110:
    *iwfr = 1;
    lbig = 0;
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	l = iq[i__];
	lbig = max(l,lbig);
	*iwfr += l;
	ipe[i__] = *iwfr - 1;
/* L120: */
    }

/* PERFORM IN-PLACE SORT */
    if (*nz == 0) {
	goto L250;
    }
    i__1 = *nz;
    for (k = 1; k <= i__1; ++k) {
	i__ = -iw[k];
	if (i__ <= 0) {
	    goto L160;
	}
	l = k;
	iw[k] = 0;
	i__2 = *nz;
	for (id = 1; id <= i__2; ++id) {
	    j = icn[l];
	    if (perm[i__] < perm[j]) {
		goto L130;
	    }
	    l = ipe[j];
	    ipe[j] = l - 1;
	    in = iw[l];
	    iw[l] = i__;
	    goto L140;
L130:
	    l = ipe[i__];
	    ipe[i__] = l - 1;
	    in = iw[l];
	    iw[l] = j;
L140:
	    i__ = -in;
	    if (i__ <= 0) {
		goto L160;
	    }
/* L150: */
	}
L160:
	;
    }

/* MAKE ROOM IN IW FOR ROW LENGTHS AND INITIALIZE FLAG. */
    k = *iwfr - 1;
    l = k + *n;
    *iwfr = l + 1;
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	flag__[i__] = 0;
	j = *n + 1 - i__;
	len = iq[j];
	if (len <= 0) {
	    goto L180;
	}
	i__2 = len;
	for (jdummy = 1; jdummy <= i__2; ++jdummy) {
	    iw[l] = iw[k];
	    --k;
	    --l;
/* L170: */
	}
L180:
	ipe[j] = l;
	--l;
/* L190: */
    }
    if (lbig >= icntl[4]) {
	goto L210;
    }

/* PLACE ROW LENGTHS IN IW */
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	k = ipe[i__];
	iw[k] = iq[i__];
	if (iq[i__] == 0) {
	    ipe[i__] = 0;
	}
/* L200: */
    }
    goto L250;


/* REMOVE DUPLICATE ENTRIES */
L210:
    *iwfr = 1;
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	k1 = ipe[i__] + 1;
	k2 = ipe[i__] + iq[i__];
	if (k1 <= k2) {
	    goto L220;
	}
	ipe[i__] = 0;
	goto L240;
L220:
	ipe[i__] = *iwfr;
	++(*iwfr);
	i__2 = k2;
	for (k = k1; k <= i__2; ++k) {
	    j = iw[k];
	    if (flag__[j] == i__) {
		goto L230;
	    }
	    iw[*iwfr] = j;
	    ++(*iwfr);
	    flag__[j] = i__;
L230:
	    ;
	}
	k = ipe[i__];
	iw[k] = *iwfr - k - 1;
L240:
	;
    }
L250:
    return 0;
} /* ma27jd_ */

/* Subroutine */ int ma27kd_(integer *n, integer *ipe, integer *iw, integer *
	lw, integer *iwfr, integer *ips, integer *ipv, integer *nv, integer *
	flag__, integer *ncmpa)
{
    /* System generated locals */
    integer i__1, i__2, i__3, i__4, i__5;

    /* Local variables */
    static integer i__, j, ie, je, me, ip, jp, ln, ml, js, ms, jp1, jp2, lwfr;
    extern /* Subroutine */ int ma27ud_(integer *, integer *, integer *, 
	    integer *, integer *, integer *);
    static integer minjs, kdummy;


/* USING A GIVEN PIVOTAL SEQUENCE AND A REPRESENTATION OF THE MATRIX THAT */
/*     INCLUDES ONLY NON-ZEROS OF THE STRICTLY UPPER-TRIANGULAR PART */
/*     OF THE PERMUTED MATRIX, CONSTRUCT TREE POINTERS. */

/* N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED. */
/* IPE(I) MUST BE SET TO POINT TO THE POSITION IN IW OF THE */
/*     START OF ROW I OR HAVE THE VALUE ZERO IF ROW I HAS NO OFF- */
/*     DIAGONAL NON-ZEROS. DURING EXECUTION IT IS USED AS FOLLOWS. */
/*     IF VARIABLE I IS ELIMINATED THEN IPE(I) POINTS TO THE LIST */
/*     OF VARIABLES FOR CREATED ELEMENT I. IF ELEMENT I IS */
/*     ABSORBED INTO NEWLY CREATED ELEMENT J THEN IPE(I)=-J. */
/* IW MUST BE SET ON ENTRY TO HOLD LISTS OF VARIABLES BY */
/*     ROWS, EACH LIST BEING HEADED BY ITS LENGTH. WHEN A VARIABLE */
/*     IS ELIMINATED ITS LIST IS REPLACED BY A LIST OF VARIABLES */
/*     IN THE NEW ELEMENT. */
/* LW MUST BE SET TO THE LENGTH OF IW. IT IS NOT ALTERED. */
/* IWFR MUST BE SET TO THE POSITION IN IW OF THE FIRST FREE VARIABLE. */
/*     IT IS REVISED DURING EXECUTION, CONTINUING TO HAVE THIS MEANING. */
/* IPS(I) MUST BE SET TO THE POSITION OF VARIABLE I IN THE REQUIRED */
/*     ORDERING. IT IS NOT ALTERED. */
/* IPV NEED NOT BE SET. IPV(K) IS SET TO HOLD THE K TH VARIABLE */
/*     IN PIVOT ORDER. */
/* NV NEED NOT BE SET. IF VARIABLE J HAS NOT BEEN ELIMINATED THEN */
/*     THE LAST ELEMENT WHOSE LEADING VARIABLE (VARIABLE EARLIEST */
/*     IN THE PIVOT SEQUENCE) IS J IS ELEMENT NV(J). IF ELEMENT J */
/*     EXISTS THEN THE LAST ELEMENT HAVING THE SAME LEADING */
/*     VARIABLE IS NV(J). IN BOTH CASES NV(J)=0 IF THERE IS NO SUCH */
/*     ELEMENT. IF ELEMENT J HAS BEEN MERGED INTO A LATER ELEMENT */
/*     THEN NV(J) IS THE DEGREE AT THE TIME OF ELIMINATION. */
/* FLAG IS USED AS WORKSPACE FOR VARIABLE FLAGS. */
/*     FLAG(JS)=ME IF JS HAS BEEN INCLUDED IN THE LIST FOR ME. */
/* NCMPA see INFO(11) in MA27A/AD. */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. External Subroutines .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */

/* INITIALIZATIONS */
    /* Parameter adjustments */
    --flag__;
    --nv;
    --ipv;
    --ips;
    --ipe;
    --iw;

    /* Function Body */
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	flag__[i__] = 0;
	nv[i__] = 0;
	j = ips[i__];
	ipv[j] = i__;
/* L10: */
    }
    *ncmpa = 0;

/* START OF MAIN LOOP */

    i__1 = *n;
    for (ml = 1; ml <= i__1; ++ml) {
/* ME=MS IS THE NAME OF THE VARIABLE ELIMINATED AND */
/*     OF THE ELEMENT CREATED IN THE MAIN LOOP. */
	ms = ipv[ml];
	me = ms;
	flag__[ms] = me;

/* MERGE ROW MS WITH ALL THE ELEMENTS HAVING MS AS LEADING VARIABLE. */
/* IP POINTS TO THE START OF THE NEW LIST. */
	ip = *iwfr;
/* MINJS IS SET TO THE POSITION IN THE ORDER OF THE LEADING VARIABLE */
/*     IN THE NEW LIST. */
	minjs = *n;
	ie = me;
	i__2 = *n;
	for (kdummy = 1; kdummy <= i__2; ++kdummy) {
/* SEARCH VARIABLE LIST OF ELEMENT IE. */
/* JP POINTS TO THE CURRENT POSITION IN THE LIST BEING SEARCHED. */
	    jp = ipe[ie];
/* LN IS THE LENGTH OF THE LIST BEING SEARCHED. */
	    ln = 0;
	    if (jp <= 0) {
		goto L60;
	    }
	    ln = iw[jp];

/* SEARCH FOR DIFFERENT VARIABLES AND ADD THEM TO LIST, */
/*     COMPRESSING WHEN NECESSARY */
	    i__3 = ln;
	    for (jp1 = 1; jp1 <= i__3; ++jp1) {
		++jp;
/* PLACE NEXT VARIABLE IN JS. */
		js = iw[jp];
/* JUMP IF VARIABLE HAS ALREADY BEEN INCLUDED. */
		if (flag__[js] == me) {
		    goto L50;
		}
		flag__[js] = me;
		if (*iwfr < *lw) {
		    goto L40;
		}
/* PREPARE FOR COMPRESSING IW BY ADJUSTING POINTER TO AND LENGTH OF */
/*     THE LIST FOR IE TO REFER TO THE REMAINING ENTRIES. */
		ipe[ie] = jp;
		iw[jp] = ln - jp1;
/* COMPRESS IW. */
		i__4 = ip - 1;
		ma27ud_(n, &ipe[1], &iw[1], &i__4, &lwfr, ncmpa);
/* COPY NEW LIST FORWARD */
		jp2 = *iwfr - 1;
		*iwfr = lwfr;
		if (ip > jp2) {
		    goto L30;
		}
		i__4 = jp2;
		for (jp = ip; jp <= i__4; ++jp) {
		    iw[*iwfr] = iw[jp];
		    ++(*iwfr);
/* L20: */
		}
L30:
		ip = lwfr;
		jp = ipe[ie];
/* ADD VARIABLE JS TO NEW LIST. */
L40:
		iw[*iwfr] = js;
/* Computing MIN */
		i__4 = minjs, i__5 = ips[js];
		minjs = min(i__4,i__5);
		++(*iwfr);
L50:
		;
	    }
/* RECORD ABSORPTION OF ELEMENT IE INTO NEW ELEMENT. */
L60:
	    ipe[ie] = -me;
/* PICK UP NEXT ELEMENT WITH LEADING VARIABLE MS. */
	    je = nv[ie];
/* STORE DEGREE OF IE. */
	    nv[ie] = ln + 1;
	    ie = je;
/* LEAVE LOOP IF THERE ARE NO MORE ELEMENTS. */
	    if (ie == 0) {
		goto L80;
	    }
/* L70: */
	}
L80:
	if (*iwfr > ip) {
	    goto L90;
	}
/* DEAL WITH NULL NEW ELEMENT. */
	ipe[me] = 0;
	nv[me] = 1;
	goto L100;
/* LINK NEW ELEMENT WITH OTHERS HAVING SAME LEADING VARIABLE. */
L90:
	minjs = ipv[minjs];
	nv[me] = nv[minjs];
	nv[minjs] = me;
/* MOVE FIRST ENTRY IN NEW LIST TO END TO ALLOW ROOM FOR LENGTH AT */
/*     FRONT. SET POINTER TO FRONT. */
	iw[*iwfr] = iw[ip];
	iw[ip] = *iwfr - ip;
	ipe[me] = ip;
	++(*iwfr);
L100:
	;
    }
    return 0;
} /* ma27kd_ */

/* Subroutine */ int ma27ld_(integer *n, integer *ipe, integer *nv, integer *
	ips, integer *ne, integer *na, integer *nd, integer *nsteps, integer *
	nemin)
{
    /* System generated locals */
    integer i__1, i__2;

    /* Local variables */
    static integer i__, k, l, ib, if__, il, is, nr, ison;


/* TREE SEARCH */

/* GIVEN SON TO FATHER TREE POINTERS, PERFORM DEPTH-FIRST */
/*     SEARCH TO FIND PIVOT ORDER AND NUMBER OF ELIMINATIONS */
/*     AND ASSEMBLIES AT EACH STAGE. */
/* N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED. */
/* IPE(I) MUST BE SET EQUAL TO -(FATHER OF NODE I) OR ZERO IF */
/*      NODE IS A ROOT. IT IS ALTERED TO POINT TO ITS NEXT */
/*      YOUNGER BROTHER IF IT HAS ONE, BUT OTHERWISE IS NOT */
/*      CHANGED. */
/* NV(I) MUST BE SET TO ZERO IF NO VARIABLES ARE ELIMINATED AT NODE */
/*      I AND TO THE DEGREE OTHERWISE. ONLY LEAF NODES CAN HAVE */
/*      ZERO VALUES OF NV(I). NV IS NOT ALTERED. */
/* IPS(I) NEED NOT BE SET. IT IS USED TEMPORARILY TO HOLD */
/*      -(ELDEST SON OF NODE I) IF IT HAS ONE AND 0 OTHERWISE. IT IS */
/*      EVENTUALLY SET TO HOLD THE POSITION OF NODE I IN THE ORDER. */
/* NE(IS) NEED NOT BE SET. IT IS SET TO THE NUMBER OF VARIABLES */
/*      ELIMINATED AT STAGE IS OF THE ELIMINATION. */
/* NA(IS) NEED NOT BE SET. IT IS SET TO THE NUMBER OF ELEMENTS */
/*      ASSEMBLED AT STAGE IS OF THE ELIMINATION. */
/* ND(IS) NEED NOT BE SET. IT IS SET TO THE DEGREE AT STAGE IS OF */
/*     THE ELIMINATION. */
/* NSTEPS NEED NOT BE SET. IT IS SET TO  THE NUMBER OF ELIMINATION */
/*      STEPS. */
/* NEMIN see ICNTL(5) in MA27A/AD. */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. Executable Statements .. */
/* INITIALIZE IPS AND NE. */
    /* Parameter adjustments */
    --nd;
    --na;
    --ne;
    --ips;
    --nv;
    --ipe;

    /* Function Body */
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	ips[i__] = 0;
	ne[i__] = 0;
/* L10: */
    }

/* SET IPS(I) TO -(ELDEST SON OF NODE I) AND IPE(I) TO NEXT YOUNGER */
/*     BROTHER OF NODE I IF IT HAS ONE. */
/* FIRST PASS IS FOR NODES WITHOUT ELIMINATIONS. */
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	if (nv[i__] > 0) {
	    goto L20;
	}
	if__ = -ipe[i__];
	is = -ips[if__];
	if (is > 0) {
	    ipe[i__] = is;
	}
	ips[if__] = -i__;
L20:
	;
    }
/* NR IS DECREMENTED FOR EACH ROOT NODE. THESE ARE STORED IN */
/*     NE(I),I=NR,N. */
    nr = *n + 1;
/* SECOND PASS TO ADD NODES WITH ELIMINATIONS. */
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	if (nv[i__] <= 0) {
	    goto L50;
	}
/* NODE IF IS THE FATHER OF NODE I. */
	if__ = -ipe[i__];
	if (if__ == 0) {
	    goto L40;
	}
	is = -ips[if__];
/* JUMP IF NODE IF HAS NO SONS YET. */
	if (is <= 0) {
	    goto L30;
	}
/* SET POINTER TO NEXT BROTHER */
	ipe[i__] = is;
/* NODE I IS ELDEST SON OF NODE IF. */
L30:
	ips[if__] = -i__;
	goto L50;
/* WE HAVE A ROOT */
L40:
	--nr;
	ne[nr] = i__;
L50:
	;
    }

/* DEPTH-FIRST SEARCH. */
/* IL HOLDS THE CURRENT TREE LEVEL. ROOTS ARE AT LEVEL N, THEIR SONS */
/*     ARE AT LEVEL N-1, ETC. */
/* IS HOLDS THE CURRENT ELIMINATION STAGE. WE ACCUMULATE THE NUMBER */
/*     OF ELIMINATIONS AT STAGE IS DIRECTLY IN NE(IS). THE NUMBER OF */
/*     ASSEMBLIES IS ACCUMULATED TEMPORARILY IN NA(IL), FOR TREE */
/*     LEVEL IL, AND IS TRANSFERED TO NA(IS) WHEN WE REACH THE */
/*     APPROPRIATE STAGE IS. */
    is = 1;
/* I IS THE CURRENT NODE. */
    i__ = 0;
    i__1 = *n;
    for (k = 1; k <= i__1; ++k) {
	if (i__ > 0) {
	    goto L60;
	}
/* PICK UP NEXT ROOT. */
	i__ = ne[nr];
	ne[nr] = 0;
	++nr;
	il = *n;
	na[*n] = 0;
/* GO TO SON FOR AS LONG AS POSSIBLE, CLEARING FATHER-SON POINTERS */
/*     IN IPS AS EACH IS USED AND SETTING NA(IL)=0 FOR ALL LEVELS */
/*     REACHED. */
L60:
	i__2 = *n;
	for (l = 1; l <= i__2; ++l) {
	    if (ips[i__] >= 0) {
		goto L80;
	    }
	    ison = -ips[i__];
	    ips[i__] = 0;
	    i__ = ison;
	    --il;
	    na[il] = 0;
/* L70: */
	}
/* RECORD POSITION OF NODE I IN THE ORDER. */
L80:
	ips[i__] = k;
	++ne[is];
/* JUMP IF NODE HAS NO ELIMINATIONS. */
	if (nv[i__] <= 0) {
	    goto L120;
	}
	if (il < *n) {
	    ++na[il + 1];
	}
	na[is] = na[il];
	nd[is] = nv[i__];
/* CHECK FOR STATIC CONDENSATION */
	if (na[is] != 1) {
	    goto L90;
	}
	if (nd[is - 1] - ne[is - 1] == nd[is]) {
	    goto L100;
	}
/* CHECK FOR SMALL NUMBERS OF ELIMINATIONS IN BOTH LAST TWO STEPS. */
L90:
	if (ne[is] >= *nemin) {
	    goto L110;
	}
	if (na[is] == 0) {
	    goto L110;
	}
	if (ne[is - 1] >= *nemin) {
	    goto L110;
	}
/* COMBINE THE LAST TWO STEPS */
L100:
	na[is - 1] = na[is - 1] + na[is] - 1;
	nd[is - 1] = nd[is] + ne[is - 1];
	ne[is - 1] = ne[is] + ne[is - 1];
	ne[is] = 0;
	goto L120;
L110:
	++is;
L120:
	ib = ipe[i__];
	if (ib >= 0) {
/* NODE I HAS A BROTHER OR IS A ROOT */
	    if (ib > 0) {
		na[il] = 0;
	    }
	    i__ = ib;
	} else {
/* GO TO FATHER OF NODE I */
	    i__ = -ib;
	    ++il;
	}
/* L160: */
    }
    *nsteps = is - 1;
    return 0;
} /* ma27ld_ */

/* Subroutine */ int ma27md_(integer *n, integer *nz, integer *irn, integer *
	icn, integer *perm, integer *na, integer *ne, integer *nd, integer *
	nsteps, integer *lstki, integer *lstkr, integer *iw, integer *info, 
	doublereal *ops)
{
    /* System generated locals */
    integer i__1, i__2, i__3;

    /* Local variables */
    static integer i__, k, nz1, nz2, nfr, iold, jold, iorg, jorg, inew, itop, 
	    lstk, nstk, irow;
    static doublereal delim;
    static integer nelim, itree, istki, nassr, istkr, nirnec, nrlnec, niradu, 
	    nrladu, numorg, nirtot, nrltot;


/* STORAGE AND OPERATION COUNT EVALUATION. */

/* EVALUATE NUMBER OF OPERATIONS AND SPACE REQUIRED BY FACTORIZATION */
/*     USING MA27B/BD.  THE VALUES GIVEN ARE EXACT ONLY IF NO NUMERICAL */
/*     PIVOTING IS PERFORMED AND THEN ONLY IF IRN(1) WAS NOT */
/*     EQUIVALENCED TO IW(1) BY THE USER BEFORE CALLING MA27A/AD.  IF */
/*     THE EQUIVALENCE HAS BEEN MADE ONLY AN UPPER BOUND FOR NIRNEC */
/*     AND NRLNEC CAN BE CALCULATED ALTHOUGH THE OTHER COUNTS WILL */
/*     STILL BE EXACT. */

/* N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED. */
/* NZ MUST BE SET TO THE NUMBER OF NON-ZEROS INPUT.  IT IS NOT ALTERED. */
/* IRN,ICN.  UNLESS IRN(1) HAS BEEN EQUIVALENCED TO IW(1) */
/*     IRN,ICN MUST BE SET TO THE ROW AND COLUMN INDICES OF THE */
/*     NON-ZEROS ON INPUT.  THEY ARE NOT ALTERED BY MA27M/MD. */
/* PERM MUST BE SET TO THE POSITION IN THE PIVOT ORDER OF EACH ROW. */
/*     IT IS NOT ALTERED. */
/* NA,NE,ND MUST BE SET TO HOLD, FOR EACH TREE NODE, THE NUMBER OF STACK */
/*     ELEMENTS ASSEMBLED, THE NUMBER OF ELIMINATIONS AND THE SIZE OF */
/*     THE ASSEMBLED FRONT MATRIX RESPECTIVELY.  THEY ARE NOT ALTERED. */
/* NSTEPS MUST BE SET TO HOLD THE NUMBER OF TREE NODES. IT IS NOT */
/*     ALTERED. */
/* LSTKI IS USED AS A WORK ARRAY BY MA27M/MD. */
/* LSTKR.  IF IRN(1) IS EQUIVALENCED TO IW(1)  THEN LSTKR(I) */
/*     MUST HOLD THE TOTAL NUMBER OF OFF-DIAGONAL ENTRIES (INCLUDING */
/*     DUPLICATES) IN ROW I (I=1,..,N) OF THE ORIGINAL MATRIX.  IT */
/*     IS USED AS WORKSPACE BY MA27M/MD. */
/* IW IS A WORKSPACE ARRAY USED BY OTHER SUBROUTINES AND PASSED TO THIS */
/*     SUBROUTINE ONLY SO THAT A TEST FOR EQUIVALENCE WITH IRN CAN BE */
/*     MADE. */

/* COUNTS FOR OPERATIONS AND STORAGE ARE ACCUMULATED IN VARIABLES */
/*     OPS,NRLTOT,NIRTOT,NRLNEC,NIRNEC,NRLADU,NRLNEC,NIRNEC. */
/* OPS NUMBER OF MULTIPLICATIONS AND ADDITIONS DURING FACTORIZATION. */
/* NRLADU,NIRADU REAL AND INTEGER STORAGE RESPECTIVELY FOR THE */
/*     MATRIX FACTORS. */
/* NRLTOT,NIRTOT REAL AND INTEGER STRORAGE RESPECTIVELY REQUIRED */
/*     FOR THE FACTORIZATION IF NO COMPRESSES ARE ALLOWED. */
/* NRLNEC,NIRNEC REAL AND INTEGER STORAGE RESPECTIVELY REQUIRED FOR */
/*     THE FACTORIZATION IF COMPRESSES ARE ALLOWED. */
/* INFO is an INTEGER array of length 20, see MA27A/AD. */
/* OPS ACCUMULATES THE NO. OF MULTIPLY/ADD PAIRS NEEDED TO CREATE THE */
/*     TRIANGULAR FACTORIZATION, IN THE DEFINITE CASE. */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */

    /* Parameter adjustments */
    --lstkr;
    --lstki;
    --perm;
    --irn;
    --icn;
    --nd;
    --ne;
    --na;
    --iw;
    --info;

    /* Function Body */
    if (*nz == 0) {
	goto L20;
    }
/* JUMP IF IW AND IRN HAVE NOT BEEN EQUIVALENCED. */
    if (irn[1] != iw[1]) {
	goto L20;
    }
/* RESET IRN(1). */
    irn[1] = iw[1] - 1;
/* THE TOTAL NUMBER OF OFF-DIAGONAL ENTRIES IS ACCUMULATED IN NZ2. */
/* LSTKI IS SET TO HOLD THE TOTAL NUMBER OF ENTRIES (INCUDING */
/*     THE DIAGONAL) IN EACH ROW IN PERMUTED ORDER. */
    nz2 = 0;
    i__1 = *n;
    for (iold = 1; iold <= i__1; ++iold) {
	inew = perm[iold];
	lstki[inew] = lstkr[iold] + 1;
	nz2 += lstkr[iold];
/* L10: */
    }
/* NZ1 IS THE NUMBER OF ENTRIES IN ONE TRIANGLE INCLUDING THE DIAGONAL. */
/* NZ2 IS THE TOTAL NUMBER OF ENTRIES INCLUDING THE DIAGONAL. */
    nz1 = nz2 / 2 + *n;
    nz2 += *n;
    goto L60;
/* COUNT (IN LSTKI) NON-ZEROS IN ORIGINAL MATRIX BY PERMUTED ROW (UPPER */
/*     TRIANGLE ONLY). INITIALIZE COUNTS. */
L20:
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	lstki[i__] = 1;
/* L30: */
    }
/* ACCUMULATE NUMBER OF NON-ZEROS WITH INDICES IN RANGE IN NZ1 */
/*     DUPLICATES ON THE DIAGONAL ARE IGNORED BUT NZ1 INCLUDES ANY */
/*     DIAGONALS NOT PRESENT ON INPUT. */
/* ACCUMULATE ROW COUNTS IN LSTKI. */
    nz1 = *n;
    if (*nz == 0) {
	goto L50;
    }
    i__1 = *nz;
    for (i__ = 1; i__ <= i__1; ++i__) {
	iold = irn[i__];
	jold = icn[i__];
/* JUMP IF INDEX IS OUT OF RANGE. */
	if (iold < 1 || iold > *n) {
	    goto L40;
	}
	if (jold < 1 || jold > *n) {
	    goto L40;
	}
	if (iold == jold) {
	    goto L40;
	}
	++nz1;
/* Computing MIN */
	i__2 = perm[iold], i__3 = perm[jold];
	irow = min(i__2,i__3);
	++lstki[irow];
L40:
	;
    }
L50:
    nz2 = nz1;
/* ISTKR,ISTKI CURRENT NUMBER OF STACK ENTRIES IN */
/*     REAL AND INTEGER STORAGE RESPECTIVELY. */
/* OPS,NRLADU,NIRADU,NIRTOT,NRLTOT,NIRNEC,NRLNEC,NZ2 ARE DEFINED ABOVE. */
/* NZ2 CURRENT NUMBER OF ORIGINAL MATRIX ENTRIES NOT YET PROCESSED. */
/* NUMORG CURRENT TOTAL NUMBER OF ROWS ELIMINATED. */
/* ITOP CURRENT NUMBER OF ELEMENTS ON THE STACK. */
L60:
    istki = 0;
    istkr = 0;
    *ops = 0.;
    nrladu = 0;
/*     ONE LOCATION IS NEEDED TO RECORD THE NUMBER OF BLOCKS */
/*     ACTUALLY USED. */
    niradu = 1;
    nirtot = nz1;
    nrltot = nz1;
    nirnec = nz2;
    nrlnec = nz2;
    numorg = 0;
    itop = 0;

/* EACH PASS THROUGH THIS LOOP PROCESSES A NODE OF THE TREE. */
    i__1 = *nsteps;
    for (itree = 1; itree <= i__1; ++itree) {
	nelim = ne[itree];
	delim = (doublereal) nelim;
	nfr = nd[itree];
	nstk = na[itree];
/* ADJUST STORAGE COUNTS ON ASSEMBLY OF CURRENT FRONTAL MATRIX. */
	nassr = nfr * (nfr + 1) / 2;
	if (nstk != 0) {
	    nassr = nassr - lstkr[itop] + 1;
	}
/* Computing MAX */
	i__2 = nrltot, i__3 = nrladu + nassr + istkr + nz1;
	nrltot = max(i__2,i__3);
/* Computing MAX */
	i__2 = nirtot, i__3 = niradu + nfr + 2 + istki + nz1;
	nirtot = max(i__2,i__3);
/* Computing MAX */
	i__2 = nrlnec, i__3 = nrladu + nassr + istkr + nz2;
	nrlnec = max(i__2,i__3);
/* Computing MAX */
	i__2 = nirnec, i__3 = niradu + nfr + 2 + istki + nz2;
	nirnec = max(i__2,i__3);
/* DECREASE NZ2 BY THE NUMBER OF ENTRIES IN ROWS BEING ELIMINATED AT */
/*     THIS STAGE. */
	i__2 = nelim;
	for (iorg = 1; iorg <= i__2; ++iorg) {
	    jorg = numorg + iorg;
	    nz2 -= lstki[jorg];
/* L70: */
	}
	numorg += nelim;
/* JUMP IF THERE ARE NO STACK ASSEMBLIES AT THIS NODE. */
	if (nstk <= 0) {
	    goto L90;
	}
/* REMOVE ELEMENTS FROM THE STACK.  THERE ARE ITOP ELEMENTS ON THE */
/*     STACK WITH THE APPROPRIATE ENTRIES IN LSTKR,LSTKI GIVING */
/*     THE REAL AND INTEGER STORAGE RESPECTIVELY FOR EACH STACK */
/*     ELEMENT. */
	i__2 = nstk;
	for (k = 1; k <= i__2; ++k) {
	    lstk = lstkr[itop];
	    istkr -= lstk;
	    lstk = lstki[itop];
	    istki -= lstk;
	    --itop;
/* L80: */
	}
/* ACCUMULATE NON-ZEROS IN FACTORS AND NUMBER OF OPERATIONS. */
L90:
	nrladu += nelim * ((nfr << 1) - nelim + 1) / 2;
	niradu = niradu + 2 + nfr;
	if (nelim == 1) {
	    --niradu;
	}
	*ops += (nfr * delim * (nfr + 1) - ((nfr << 1) + 1) * delim * (delim 
		+ 1) / 2 + delim * (delim + 1) * (delim * 2 + 1) / 6) / 2;
	if (itree == *nsteps) {
	    goto L100;
	}
/* JUMP IF ALL OF FRONTAL MATRIX HAS BEEN ELIMINATED. */
	if (nfr == nelim) {
	    goto L100;
	}
/* STACK REMAINDER OF ELEMENT. */
	++itop;
	lstkr[itop] = (nfr - nelim) * (nfr - nelim + 1) / 2;
	lstki[itop] = nfr - nelim + 1;
	istki += lstki[itop];
	istkr += lstkr[itop];
/* WE DO NOT NEED TO ADJUST THE COUNTS FOR THE REAL STORAGE BECAUSE */
/*     THE REMAINDER OF THE FRONTAL MATRIX IS SIMPLY MOVED IN THE */
/*     STORAGE FROM FACTORS TO STACK AND NO EXTRA STORAGE IS REQUIRED. */
/* Computing MAX */
	i__2 = nirtot, i__3 = niradu + istki + nz1;
	nirtot = max(i__2,i__3);
/* Computing MAX */
	i__2 = nirnec, i__3 = niradu + istki + nz2;
	nirnec = max(i__2,i__3);
L100:
	;
    }

/* ADJUST THE STORAGE COUNTS TO ALLOW FOR THE USE OF THE REAL AND */
/*     INTEGER STORAGE FOR PURPOSES OTHER THAN PURELY THE */
/*     FACTORIZATION ITSELF. */
/* THE SECOND TWO TERMS ARE THE MINUMUM AMOUNT REQUIRED BY MA27N/ND. */
/* Computing MAX */
    i__1 = nrlnec, i__2 = *n + max(*nz,nz1);
    nrlnec = max(i__1,i__2);
/* Computing MAX */
    i__1 = nrltot, i__2 = *n + max(*nz,nz1);
    nrltot = max(i__1,i__2);
    nrlnec = min(nrlnec,nrltot);
    nirnec = max(*nz,nirnec);
    nirtot = max(*nz,nirtot);
    nirnec = min(nirnec,nirtot);
    info[3] = nrltot;
    info[4] = nirtot;
    info[5] = nrlnec;
    info[6] = nirnec;
    info[7] = nrladu;
    info[8] = niradu;
    return 0;
} /* ma27md_ */

/* Subroutine */ int ma27nd_(integer *n, integer *nz, integer *nz1, 
	doublereal *a, integer *la, integer *irn, integer *icn, integer *iw, 
	integer *liw, integer *perm, integer *iw2, integer *icntl, integer *
	info)
{
    /* Format strings */
    static char fmt_40[] = "(\002 *** WARNING MESSAGE FROM SUBROUTINE MA27B"
	    "D\002,\002  *** INFO(1) =\002,i2)";
    static char fmt_50[] = "(i6,\002TH NON-ZERO (IN ROW\002,i6,\002 AND COLU"
	    "MN\002,i6,\002) IGNORED\002)";

    /* System generated locals */
    integer i__1, i__2;

    /* Builtin functions */
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);

    /* Local variables */
    static integer i__, k, j1, j2, ia, ii, jj, ich, iiw, iold, jold, inew;
    static doublereal anow;
    static integer jnew, ipos, jpos;
    static doublereal anext;

    /* Fortran I/O blocks */
    static cilist io___212 = { 0, 0, 0, fmt_40, 0 };
    static cilist io___213 = { 0, 0, 0, fmt_50, 0 };



/* SORT PRIOR TO FACTORIZATION USING MA27O/OD. */

/* THIS SUBROUTINE REORDERS THE USER'S INPUT SO THAT THE UPPER TRIANGLE */
/*     OF THE PERMUTED MATRIX, INCLUDING THE DIAGONAL, IS */
/*     HELD ORDERED BY ROWS AT THE END OF THE STORAGE FOR A AND IW. */
/*     IT IGNORES ENTRIES WITH ONE OR BOTH INDICES OUT OF RANGE AND */
/*     ACCUMULATES DIAGONAL ENTRIES. */
/*     IT ADDS EXPLICIT ZEROS ON THE DIAGONAL WHERE NECESSARY. */
/* N      - MUST BE SET TO THE ORDER OF THE MATRIX. */
/*          IT IS NOT ALTERED BY MA27N/ND. */
/* NZ     - ON ENTRY NZ MUST BE SET TO THE NUMBER */
/*          OF NON-ZEROS INPUT.  NOT ALTERED BY MA27N/ND. */
/* NZ1    - ON EXIT NZ1 WILL BE EQUAL TO THE NUMBER OF ENTRIES IN THE */
/*          SORTED MATRIX. */
/* A      - ON ENTRY A(I) MUST */
/*          HOLD THE VALUE OF THE ORIGINAL MATRIX ELEMENT IN POSITION */
/*          (IRN(I),ICN(I)),I=1,NZ.  ON EXIT A(LA-NZ1+I),I=1,NZ1 HOLDS */
/*          THE UPPER TRIANGLE OF THE PERMUTED MATRIX BY ROWS WITH */
/*          THE DIAGONAL ENTRY FIRST ALTHOUGH THERE IS NO FURTHER */
/*          ORDERING WITHIN THE ROWS THEMSELVES. */
/* LA     - LENGTH OF ARRAY A. MUST BE AT LEAST N+MAX(NZ,NZ1). */
/*          IT IS NOT ALTERED BY MA27N/ND. */
/* IRN    - IRN(I) MUST BE SET TO */
/*          HOLD THE ROW INDEX OF ENTRY A(I),I=1,NZ.  IRN WILL BE */
/*          UNALTERED BY MA27N/ND, UNLESS IT IS EQUIVALENCED WITH IW. */
/* ICN    - ICN(I) MUST BE SET TO */
/*          HOLD THE COLUMN INDEX OF ENTRY A(I),I=1,NZ.  ICN WILL BE */
/*          UNALTERED BY MA27N/ND, UNLESS IT IS EQUIVALENCED WITH IW. */
/* IW     - USED AS WORKSPACE AND ON */
/*          EXIT, ENTRIES IW(LIW-NZ1+I),I=1,NZ1 HOLD THE COLUMN INDICES */
/*          (THE ORIGINAL UNPERMUTED INDICES) OF THE CORRESPONDING ENTRY */
/*          OF A WITH THE FIRST ENTRY FOR EACH ROW FLAGGED NEGATIVE. */
/*          IRN(1) MAY BE EQUIVALENCED WITH IW(1) AND ICN(1) MAY BE */
/*          EQUIVALENCED WITH IW(K) WHERE K.GT.NZ. */
/* LIW    - LENGTH OF ARRAY IW. MUST BE AT LEAST AS */
/*          GREAT AS THE MAXIMUM OF NZ AND NZ1. */
/*          NOT ALTERED BY MA27N/ND. */
/* PERM   - PERM(I) HOLDS THE */
/*          POSITION IN THE TENTATIVE PIVOT ORDER OF ROW I IN THE */
/*          ORIGINAL MATRIX,I=1,N. NOT ALTERED BY MA27N/ND. */
/* IW2    - USED AS WORKSPACE. */
/*          SEE COMMENTS IN CODE IMMEDIATELY PRIOR TO */
/*          EACH USE. */
/* ICNTL is an INTEGER array of length 30, see MA27A/AD. */
/* INFO is an INTEGER array of length 20, see MA27A/AD. */
/*   INFO(1)  - ON EXIT FROM MA27N/ND, A ZERO VALUE OF */
/*          INFO(1) INDICATES THAT NO ERROR HAS BEEN DETECTED. */
/*          POSSIBLE NON-ZERO VALUES ARE .. */
/*          +1  WARNING.  INDICES OUT OF RANGE.  THESE ARE IGNORED, */
/*              THEIR NUMBER IS RECORDED IN INFO(2) OF MA27E/ED AND */
/*              MESSAGES IDENTIFYING THE FIRST TEN ARE OUTPUT ON UNIT */
/*              ICNTL(2). */
/*          -3  INTEGER ARRAY IW IS TOO SMALL. */
/*          -4  DOUBLE PRECISION ARRAY A IS TOO SMALL. */

/*     .. Parameters .. */
/*     .. */
/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */
    /* Parameter adjustments */
    --iw2;
    --perm;
    --a;
    --irn;
    --icn;
    --iw;
    --icntl;
    --info;

    /* Function Body */
    info[1] = 0;
/* INITIALIZE WORK ARRAY (IW2) IN PREPARATION FOR */
/*     COUNTING NUMBERS OF NON-ZEROS IN THE ROWS AND INITIALIZE */
/*     LAST N ENTRIES IN A WHICH WILL HOLD THE DIAGONAL ENTRIES */
    ia = *la;
    i__1 = *n;
    for (iold = 1; iold <= i__1; ++iold) {
	iw2[iold] = 1;
	a[ia] = 0.;
	--ia;
/* L10: */
    }
/* SCAN INPUT COPYING ROW INDICES FROM IRN TO THE FIRST NZ POSITIONS */
/*     IN IW.  THE NEGATIVE OF THE INDEX IS HELD TO FLAG ENTRIES FOR */
/*     THE IN-PLACE SORT.  ENTRIES IN IW CORRESPONDING TO DIAGONALS AND */
/*     ENTRIES WITH OUT-OF-RANGE INDICES ARE SET TO ZERO. */
/*     FOR DIAGONAL ENTRIES, REALS ARE ACCUMULATED IN THE LAST N */
/*     LOCATIONS OF A. */
/*     THE NUMBER OF ENTRIES IN EACH ROW OF THE PERMUTED MATRIX IS */
/*     ACCUMULATED IN IW2. */
/* INDICES OUT OF RANGE ARE IGNORED  AFTER BEING COUNTED AND */
/*     AFTER APPROPRIATE MESSAGES HAVE BEEN PRINTED. */
    info[2] = 0;
/* NZ1 IS THE NUMBER OF NON-ZEROS HELD AFTER INDICES OUT OF RANGE HAVE */
/*     BEEN IGNORED AND DIAGONAL ENTRIES ACCUMULATED. */
    *nz1 = *n;
    if (*nz == 0) {
	goto L80;
    }
    i__1 = *nz;
    for (k = 1; k <= i__1; ++k) {
	iold = irn[k];
	if (iold > *n || iold <= 0) {
	    goto L30;
	}
	jold = icn[k];
	if (jold > *n || jold <= 0) {
	    goto L30;
	}
	inew = perm[iold];
	jnew = perm[jold];
	if (inew != jnew) {
	    goto L20;
	}
	ia = *la - *n + iold;
	a[ia] += a[k];
	goto L60;
L20:
	inew = min(inew,jnew);
/* INCREMENT NUMBER OF ENTRIES IN ROW INEW. */
	++iw2[inew];
	iw[k] = -iold;
	++(*nz1);
	goto L70;
/* ENTRY OUT OF RANGE.  IT WILL BE IGNORED AND A FLAG SET. */
L30:
	info[1] = 1;
	++info[2];
	if (info[2] <= 1 && icntl[2] > 0) {
	    io___212.ciunit = icntl[2];
	    s_wsfe(&io___212);
	    do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer));
	    e_wsfe();
	}
	if (info[2] <= 10 && icntl[2] > 0) {
	    io___213.ciunit = icntl[2];
	    s_wsfe(&io___213);
	    do_fio(&c__1, (char *)&k, (ftnlen)sizeof(integer));
	    do_fio(&c__1, (char *)&irn[k], (ftnlen)sizeof(integer));
	    do_fio(&c__1, (char *)&icn[k], (ftnlen)sizeof(integer));
	    e_wsfe();
	}
L60:
	iw[k] = 0;
L70:
	;
    }
/* CALCULATE POINTERS (IN IW2) TO THE POSITION IMMEDIATELY AFTER THE END */
/*     OF EACH ROW. */
L80:
    if (*nz < *nz1 && *nz1 != *n) {
	goto L100;
    }
/* ROOM IS INCLUDED FOR THE DIAGONALS. */
    k = 1;
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	k += iw2[i__];
	iw2[i__] = k;
/* L90: */
    }
    goto L120;
/* ROOM IS NOT INCLUDED FOR THE DIAGONALS. */
L100:
    k = 1;
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	k = k + iw2[i__] - 1;
	iw2[i__] = k;
/* L110: */
    }
/* FAIL IF INSUFFICIENT SPACE IN ARRAYS A OR IW. */
L120:
    if (*nz1 > *liw) {
	goto L210;
    }
    if (*nz1 + *n > *la) {
	goto L220;
    }
/* NOW RUN THROUGH NON-ZEROS IN ORDER PLACING THEM IN THEIR NEW */
/* POSITION AND DECREMENTING APPROPRIATE IW2 ENTRY.  IF WE ARE */
/* ABOUT TO OVERWRITE AN ENTRY NOT YET MOVED, WE MUST DEAL WITH */
/* THIS AT THIS TIME. */
    if (*nz1 == *n) {
	goto L180;
    }
    i__1 = *nz;
    for (k = 1; k <= i__1; ++k) {
	iold = -iw[k];
	if (iold <= 0) {
	    goto L140;
	}
	jold = icn[k];
	anow = a[k];
	iw[k] = 0;
	i__2 = *nz;
	for (ich = 1; ich <= i__2; ++ich) {
	    inew = perm[iold];
	    jnew = perm[jold];
	    inew = min(inew,jnew);
	    if (inew == perm[jold]) {
		jold = iold;
	    }
	    jpos = iw2[inew] - 1;
	    iold = -iw[jpos];
	    anext = a[jpos];
	    a[jpos] = anow;
	    iw[jpos] = jold;
	    iw2[inew] = jpos;
	    if (iold == 0) {
		goto L140;
	    }
	    anow = anext;
	    jold = icn[jpos];
/* L130: */
	}
L140:
	;
    }
    if (*nz >= *nz1) {
	goto L180;
    }
/* MOVE UP ENTRIES TO ALLOW FOR DIAGONALS. */
    ipos = *nz1;
    jpos = *nz1 - *n;
    i__1 = *n;
    for (ii = 1; ii <= i__1; ++ii) {
	i__ = *n - ii + 1;
	j1 = iw2[i__];
	j2 = jpos;
	if (j1 > jpos) {
	    goto L160;
	}
	i__2 = j2;
	for (jj = j1; jj <= i__2; ++jj) {
	    iw[ipos] = iw[jpos];
	    a[ipos] = a[jpos];
	    --ipos;
	    --jpos;
/* L150: */
	}
L160:
	iw2[i__] = ipos + 1;
	--ipos;
/* L170: */
    }
/* RUN THROUGH ROWS INSERTING DIAGONAL ENTRIES AND FLAGGING BEGINNING */
/*     OF EACH ROW BY NEGATING FIRST COLUMN INDEX. */
L180:
    i__1 = *n;
    for (iold = 1; iold <= i__1; ++iold) {
	inew = perm[iold];
	jpos = iw2[inew] - 1;
	ia = *la - *n + iold;
	a[jpos] = a[ia];
	iw[jpos] = -iold;
/* L190: */
    }
/* MOVE SORTED MATRIX TO THE END OF THE ARRAYS. */
    ipos = *nz1;
    ia = *la;
    iiw = *liw;
    i__1 = *nz1;
    for (i__ = 1; i__ <= i__1; ++i__) {
	a[ia] = a[ipos];
	iw[iiw] = iw[ipos];
	--ipos;
	--ia;
	--iiw;
/* L200: */
    }
    goto L230;
/* **** ERROR RETURN **** */
L210:
    info[1] = -3;
    info[2] = *nz1;
    goto L230;
L220:
    info[1] = -4;
    info[2] = *nz1 + *n;

L230:
    return 0;
} /* ma27nd_ */

/* Subroutine */ int ma27od_(integer *n, integer *nz, doublereal *a, integer *
	la, integer *iw, integer *liw, integer *perm, integer *nstk, integer *
	nsteps, integer *maxfrt, integer *nelim, integer *iw2, integer *icntl,
	 doublereal *cntl, integer *info)
{
    /* Format strings */
    static char fmt_310[] = "(\002 *** WARNING MESSAGE FROM SUBROUTINE MA2"
	    "7BD\002,\002  *** INFO(1) =\002,i2,/,\002 PIVOT\002,i6,\002 HAS "
	    "DIFFERENT SIGN FROM THE PREVIOUS ONE\002)";

    /* System generated locals */
    integer i__1, i__2, i__3, i__4, i__5;
    doublereal d__1, d__2, d__3, d__4;

    /* Builtin functions */
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);

    /* Local variables */
    static integer i__, j, k, j1, j2, k1, k2, jj, kk, lt;
    static doublereal uu;
    static integer jj1, jjj, ifr, ibeg, iend, neig, iell;
    static doublereal amax;
    static integer jcol, nblk, iass, iorg, apos, astk, jmax, jnew;
    static doublereal rmax;
    static integer ipiv;
    static doublereal tmax;
    static integer ipos, istk, jpiv, jpos, kmax, irow, krow, nass, npiv, ntwo;
    static doublereal swop;
    static integer apos1, apos2, apos3, astk2, istk2, laell;
    extern /* Subroutine */ int ma27pd_(doublereal *, integer *, integer *, 
	    integer *, integer *, integer *, integer *, integer *);
    static integer iexch, liell, newel, jlast, azero, lnass;
    static doublereal amult;
    static integer ipmnp, jnext, lnpiv, iswop, iwpos, lapos2;
    static doublereal amult1, amult2;
    static integer npivp1, pospv1, pospv2, ncmpbi, posfac, ncmpbr, nirbdu, 
	    nrlbdu, ioldps;
    static doublereal detpiv, thresh;
    static integer ainput, jfirst, idummy, jdummy, jmxmip, kdummy, iinput, 
	    isnpiv, nfront, numass, numorg, numstk, pivsiz, ltopst, ntotpv;

    /* Fortran I/O blocks */
    static cilist io___280 = { 0, 0, 0, fmt_310, 0 };



/* FACTORIZATION SUBROUTINE */

/* THIS SUBROUTINE OPERATES ON THE INPUT MATRIX ORDERED BY MA27N/ND AND */
/*     PRODUCES THE FACTORS OF U AND D ('A'=UTRANSPOSE*D*U) FOR USE IN */
/*     SUBSEQUENT SOLUTIONS.  GAUSSIAN ELIMINATION IS USED WITH PIVOTS */
/*     CHOSEN FROM THE DIAGONAL.  TO ENSURE STABILITY, BLOCK PIVOTS OF */
/*     ORDER 2 WILL BE USED IF THE DIAGONAL ENTRY IS NOT LARGE ENOUGH. */

/* N      - MUST BE SET TO THE ORDER OF THE MATRIX. IT IS NOT ALTERED. */
/* NZ     - MUST BE SET TO THE NUMBER OF NON-ZEROS IN UPPER TRIANGLE OF */
/*          PERMUTED MATRIX.  NOT ALTERED BY MA27O/OD. */
/* A      - MUST BE SET ON INPUT TO MATRIX HELD BY ROWS REORDERED BY */
/*          PERMUTATION FROM MA27A/AD IN A(LA-NZ+I),I=1,NZ.   ON */
/*          EXIT FROM MA27O/OD, THE FACTORS OF U AND D ARE HELD IN */
/*          POSITIONS 1 TO POSFAC-1. */
/* LA     - LENGTH OF ARRAY A.  A VALUE FOR LA */
/*          SUFFICIENT FOR DEFINITE SYSTEMS */
/*          WILL HAVE BEEN PROVIDED BY MA27A/AD. NOT ALTERED BY MA27O/OD. */
/* IW     - MUST BE SET SO THAT,ON INPUT, IW(LIW-NZ+I),I=1,NZ */
/*          HOLDS THE COLUMN INDEX OF THE ENTRY IN A(LA-NZ+I).  ON EXIT, */
/*          IW HOLDS INTEGER INDEXING INFORMATION ON THE FACTORS. */
/*          THE ABSOLUTE VALUE OF THE FIRST ENTRY IN IW WILL BE SET TO */
/*          THE NUMBER OF BLOCK PIVOTS ACTUALLY USED.  THIS MAY BE */
/*          DIFFERENT FROM NSTEPS SINCE NUMERICAL CONSIDERATIONS */
/*          MAY PREVENT US CHOOSING A PIVOT AT EACH STAGE.  IF THIS ENTRY */
/*          IN IW IS NEGATIVE, THEN AT LEAST ONE TWO BY TWO */
/*          PIVOT HAS BEEN USED DURING THE DECOMPOSITION. */
/*          INTEGER INFORMATION ON EACH BLOCK PIVOT ROW FOLLOWS.  FOR */
/*          EACH BLOCK PIVOT ROW THE COLUMN INDICES ARE PRECEDED BY A */
/*          COUNT OF THE NUMBER OF ROWS AND COLUMNS IN THE BLOCK PIVOT */
/*          WHERE, IF ONLY ONE ROW IS PRESENT, ONLY THE NUMBER OF */
/*          COLUMNS TOGETHER WITH A NEGATIVE FLAG IS HELD.  THE FIRST */
/*          COLUMN INDEX FOR A TWO BY TWO PIVOT IS FLAGGED NEGATIVE. */
/* LIW    - LENGTH OF ARRAY IW.  A VALUE FOR LIW SUFFICIENT FOR */
/*          DEFINITE SYSTEMS */
/*          WILL HAVE BEEN PROVIDED BY MA27A/AD.  NOT ALTERED BY MA27O/OD */
/* PERM   - PERM(I) MUST BE SET TO POSITION OF ROW/COLUMN I IN THE */
/*          TENTATIVE PIVOT ORDER GENERATED BY MA27A/AD. */
/*          IT IS NOT ALTERED BY MA27O/OD. */
/* NSTK   - MUST BE LEFT UNCHANGED SINCE OUTPUT FROM MA27A/AD. NSTK(I) */
/*          GIVES THE NUMBER OF GENERATED STACK ELEMENTS ASSEMBLED AT */
/*          STAGE I.  IT IS NOT ALTERED BY MA27O/OD. */
/* NSTEPS - LENGTH OF ARRAYS NSTK AND NELIM. VALUE IS GIVEN ON OUTPUT */
/*          FROM MA27A/AD (WILL NEVER EXCEED N). IT IS NOT ALTERED BY */
/*          MA27O/OD. */
/* MAXFRT - NEED NOT BE SET ON INPUT.  ON OUTPUT */
/*          MAXFRT WILL BE SET TO THE MAXIMUM FRONT SIZE ENCOUNTERED */
/*          DURING THE DECOMPOSITION. */
/* NELIM  - MUST BE UNCHANGED SINCE OUTPUT FROM MA27A/AD. NELIM(I) */
/*          GIVES THE NUMBER OF ORIGINAL ROWS ASSEMBLED AT STAGE I. */
/*          IT IS NOT ALTERED BY MA27O/OD. */
/* IW2    - INTEGER ARRAY OF LENGTH N. USED AS WORKSPACE BY MA27O/OD. */
/*          ALTHOUGH WE COULD HAVE USED A SHORT WORD INTEGER IN THE IBM */
/*          VERSION, WE HAVE NOT DONE SO BECAUSE THERE IS A SPARE */
/*          FULL INTEGER ARRAY (USED IN THE SORT BEFORE MA27O/OD) */
/*          AVAILABLE WHEN MA27O/OD IS CALLED FROM MA27B/BD. */
/* ICNTL is an INTEGER array of length 30, see MA27A/AD. */
/* CNTL is a DOUBLE PRECISION array of length 5, see MA27A/AD. */
/* INFO is an INTEGER array of length 20, see MA27A/AD. */
/*   INFO(1)  - INTEGER VARIABLE.  DIAGNOSTIC FLAG.  A ZERO VALUE ON EXIT */
/*          INDICATES SUCCESS.  POSSIBLE NEGATIVE VALUES ARE ... */
/*          -3  INSUFFICIENT STORAGE FOR IW. */
/*          -4  INSUFFICIENT STORAGE FOR A. */
/*          -5  ZERO PIVOT FOUND IN FACTORIZATION OF DEFINITE MATRIX. */

/*     .. Parameters .. */
/*     .. */
/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. External Subroutines .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Statement Functions .. */
/*     .. */
/*     .. Statement Function definitions .. */
/* THE FOLLOWING ARITHMETIC FUNCTION GIVES THE DISPLACEMENT FROM */
/*     THE START OF THE ASSEMBLED MATRIX(OF ORDER IX) OF THE DIAGONAL */
/*     ENTRY IN ITS ROW IY. */
/*     .. */
/*     .. Executable Statements .. */
/* INITIALIZATION. */
/* NBLK IS THE NUMBER OF BLOCK PIVOTS USED. */
    /* Parameter adjustments */
    --iw2;
    --perm;
    --a;
    --iw;
    --nelim;
    --nstk;
    --icntl;
    --cntl;
    --info;

    /* Function Body */
    nblk = 0;
    ntwo = 0;
    neig = 0;
    ncmpbi = 0;
    ncmpbr = 0;
    *maxfrt = 0;
    nrlbdu = 0;
    nirbdu = 0;
/* A PRIVATE VARIABLE UU IS SET TO CNTL(1), SO THAT CNTL(1) WILL REMAIN */
/* UNALTERED. */
    uu = min(cntl[1],.5);
    uu = max(uu,-.5);
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	iw2[i__] = 0;
/* L10: */
    }
/* IWPOS IS POINTER TO FIRST FREE POSITION FOR FACTORS IN IW. */
/* POSFAC IS POINTER FOR FACTORS IN A. AT EACH PASS THROUGH THE */
/*     MAJOR LOOP POSFAC INITIALLY POINTS TO THE FIRST FREE LOCATION */
/*     IN A AND THEN IS SET TO THE POSITION OF THE CURRENT PIVOT IN A. */
/* ISTK IS POINTER TO TOP OF STACK IN IW. */
/* ISTK2 IS POINTER TO BOTTOM OF STACK IN IW (NEEDED BY COMPRESS). */
/* ASTK IS POINTER TO TOP OF STACK IN A. */
/* ASTK2 IS POINTER TO BOTTOM OF STACK IN A (NEEDED BY COMPRESS). */
/* IINPUT IS POINTER TO CURRENT POSITION IN ORIGINAL ROWS IN IW. */
/* AINPUT IS POINTER TO CURRENT POSITION IN ORIGINAL ROWS IN A. */
/* AZERO IS POINTER TO LAST POSITION ZEROED IN A. */
/* NTOTPV IS THE TOTAL NUMBER OF PIVOTS SELECTED. THIS IS USED */
/*     TO DETERMINE WHETHER THE MATRIX IS SINGULAR. */
    iwpos = 2;
    posfac = 1;
    istk = *liw - *nz + 1;
    istk2 = istk - 1;
    astk = *la - *nz + 1;
    astk2 = astk - 1;
    iinput = istk;
    ainput = astk;
    azero = 0;
    ntotpv = 0;
/* NUMASS IS THE ACCUMULATED NUMBER OF ROWS ASSEMBLED SO FAR. */
    numass = 0;

/* EACH PASS THROUGH THIS MAIN LOOP PERFORMS ALL THE OPERATIONS */
/*     ASSOCIATED WITH ONE SET OF ASSEMBLY/ELIMINATIONS. */
    i__1 = *nsteps;
    for (iass = 1; iass <= i__1; ++iass) {
/* NASS WILL BE SET TO THE NUMBER OF FULLY ASSEMBLED VARIABLES IN */
/*     CURRENT NEWLY CREATED ELEMENT. */
	nass = nelim[iass];
/* NEWEL IS A POINTER INTO IW TO CONTROL OUTPUT OF INTEGER INFORMATION */
/*     FOR NEWLY CREATED ELEMENT. */
	newel = iwpos + 1;
/* SYMBOLICALLY ASSEMBLE INCOMING ROWS AND GENERATED STACK ELEMENTS */
/*     ORDERING THE RESULTANT ELEMENT ACCORDING TO PERMUTATION PERM.  WE */
/*     ASSEMBLE THE STACK ELEMENTS FIRST BECAUSE THESE WILL ALREADY BE */
/*     ORDERED. */
/* SET HEADER POINTER FOR MERGE OF INDEX LISTS. */
	jfirst = *n + 1;
/* INITIALIZE NUMBER OF VARIABLES IN CURRENT FRONT. */
	nfront = 0;
	numstk = nstk[iass];
	ltopst = 1;
	lnass = 0;
/* JUMP IF NO STACK ELEMENTS ARE BEING ASSEMBLED AT THIS STAGE. */
	if (numstk == 0) {
	    goto L80;
	}
	j2 = istk - 1;
	lnass = nass;
	ltopst = (iw[istk] + 1) * iw[istk] / 2;
	i__2 = numstk;
	for (iell = 1; iell <= i__2; ++iell) {
/* ASSEMBLE ELEMENT IELL PLACING */
/*     THE INDICES INTO A LINKED LIST IN IW2 ORDERED */
/*     ACCORDING TO PERM. */
	    jnext = jfirst;
	    jlast = *n + 1;
	    j1 = j2 + 2;
	    j2 = j1 - 1 + iw[j1 - 1];
/* RUN THROUGH INDEX LIST OF STACK ELEMENT IELL. */
	    i__3 = j2;
	    for (jj = j1; jj <= i__3; ++jj) {
		j = iw[jj];
/* JUMP IF ALREADY ASSEMBLED */
		if (iw2[j] > 0) {
		    goto L60;
		}
		jnew = perm[j];
/* IF VARIABLE WAS PREVIOUSLY FULLY SUMMED BUT WAS NOT PIVOTED ON */
/*     EARLIER BECAUSE OF NUMERICAL TEST, INCREMENT NUMBER OF FULLY */
/*     SUMMED ROWS/COLUMNS IN FRONT. */
		if (jnew <= numass) {
		    ++nass;
		}
/* FIND POSITION IN LINKED LIST FOR NEW VARIABLE.  NOTE THAT WE START */
/*     FROM WHERE WE LEFT OFF AFTER ASSEMBLY OF PREVIOUS VARIABLE. */
		i__4 = *n;
		for (idummy = 1; idummy <= i__4; ++idummy) {
		    if (jnext == *n + 1) {
			goto L30;
		    }
		    if (perm[jnext] > jnew) {
			goto L30;
		    }
		    jlast = jnext;
		    jnext = iw2[jlast];
/* L20: */
		}
L30:
		if (jlast != *n + 1) {
		    goto L40;
		}
		jfirst = j;
		goto L50;
L40:
		iw2[jlast] = j;
L50:
		iw2[j] = jnext;
		jlast = j;
/* INCREMENT NUMBER OF VARIABLES IN THE FRONT. */
		++nfront;
L60:
		;
	    }
/* L70: */
	}
	lnass = nass - lnass;
/* NOW INCORPORATE ORIGINAL ROWS.  NOTE THAT THE COLUMNS IN THESE ROWS */
/*     NEED NOT BE IN ORDER. WE ALSO PERFORM */
/*     A SWOP SO THAT THE DIAGONAL ENTRY IS THE FIRST IN ITS */
/*     ROW.  THIS ALLOWS US TO AVOID STORING THE INVERSE OF ARRAY PERM. */
L80:
	numorg = nelim[iass];
	j1 = iinput;
	i__2 = numorg;
	for (iorg = 1; iorg <= i__2; ++iorg) {
	    j = -iw[j1];
	    i__3 = *liw;
	    for (idummy = 1; idummy <= i__3; ++idummy) {
		jnew = perm[j];
/* JUMP IF VARIABLE ALREADY INCLUDED. */
		if (iw2[j] > 0) {
		    goto L130;
		}
/* HERE WE MUST ALWAYS START OUR SEARCH AT THE BEGINNING. */
		jlast = *n + 1;
		jnext = jfirst;
		i__4 = *n;
		for (jdummy = 1; jdummy <= i__4; ++jdummy) {
		    if (jnext == *n + 1) {
			goto L100;
		    }
		    if (perm[jnext] > jnew) {
			goto L100;
		    }
		    jlast = jnext;
		    jnext = iw2[jlast];
/* L90: */
		}
L100:
		if (jlast != *n + 1) {
		    goto L110;
		}
		jfirst = j;
		goto L120;
L110:
		iw2[jlast] = j;
L120:
		iw2[j] = jnext;
/* INCREMENT NUMBER OF VARIABLES IN FRONT. */
		++nfront;
L130:
		++j1;
		if (j1 > *liw) {
		    goto L150;
		}
		j = iw[j1];
		if (j < 0) {
		    goto L150;
		}
/* L140: */
	    }
L150:
	    ;
	}
/* NOW RUN THROUGH LINKED LIST IW2 PUTTING INDICES OF VARIABLES IN NEW */
/*     ELEMENT INTO IW AND SETTING IW2 ENTRY TO POINT TO THE RELATIVE */
/*     POSITION OF THE VARIABLE IN THE NEW ELEMENT. */
	if (newel + nfront < istk) {
	    goto L160;
	}
/* COMPRESS IW. */
	ma27pd_(&a[1], &iw[1], &istk, &istk2, &iinput, &c__2, &ncmpbr, &
		ncmpbi);
	if (newel + nfront < istk) {
	    goto L160;
	}
	info[2] = *liw + 1 + newel + nfront - istk;
	goto L770;
L160:
	j = jfirst;
	i__2 = nfront;
	for (ifr = 1; ifr <= i__2; ++ifr) {
	    ++newel;
	    iw[newel] = j;
	    jnext = iw2[j];
	    iw2[j] = newel - (iwpos + 1);
	    j = jnext;
/* L170: */
	}

/* ASSEMBLE REALS INTO FRONTAL MATRIX. */
	*maxfrt = max(*maxfrt,nfront);
	iw[iwpos] = nfront;
/* FIRST ZERO OUT FRONTAL MATRIX AS APPROPRIATE FIRST CHECKING TO SEE */
/*     IF THERE IS SUFFICIENT SPACE. */
	laell = (nfront + 1) * nfront / 2;
	apos2 = posfac + laell - 1;
	if (numstk != 0) {
	    lnass = lnass * ((nfront << 1) - lnass + 1) / 2;
	}
	if (posfac + lnass - 1 >= astk) {
	    goto L180;
	}
	if (apos2 < astk + ltopst - 1) {
	    goto L190;
	}
/* COMPRESS A. */
L180:
	ma27pd_(&a[1], &iw[1], &astk, &astk2, &ainput, &c__1, &ncmpbr, &
		ncmpbi);
	if (posfac + lnass - 1 >= astk) {
	    goto L780;
	}
	if (apos2 >= astk + ltopst - 1) {
	    goto L780;
	}
L190:
	if (apos2 <= azero) {
	    goto L220;
	}
	apos = azero + 1;
/* Computing MIN */
	i__2 = apos2, i__3 = astk - 1;
	lapos2 = min(i__2,i__3);
	if (lapos2 < apos) {
	    goto L210;
	}
	i__2 = lapos2;
	for (k = apos; k <= i__2; ++k) {
	    a[k] = 0.;
/* L200: */
	}
L210:
	azero = apos2;
/* JUMP IF THERE ARE NO STACK ELEMENTS TO ASSEMBLE. */
L220:
	if (numstk == 0) {
	    goto L260;
	}
/* PLACE REALS CORRESPONDING TO STACK ELEMENTS IN CORRECT POSITIONS IN A. */
	i__2 = numstk;
	for (iell = 1; iell <= i__2; ++iell) {
	    j1 = istk + 1;
	    j2 = istk + iw[istk];
	    i__3 = j2;
	    for (jj = j1; jj <= i__3; ++jj) {
		irow = iw[jj];
		irow = iw2[irow];
		apos = posfac + (irow - 1) * (2 * nfront - irow + 2) / 2;
		i__4 = j2;
		for (jjj = jj; jjj <= i__4; ++jjj) {
		    j = iw[jjj];
		    apos2 = apos + iw2[j] - irow;
		    a[apos2] += a[astk];
		    a[astk] = 0.;
		    ++astk;
/* L230: */
		}
/* L240: */
	    }
/* INCREMENT STACK POINTER. */
	    istk = j2 + 1;
/* L250: */
	}
/* INCORPORATE REALS FROM ORIGINAL ROWS. */
L260:
	i__2 = numorg;
	for (iorg = 1; iorg <= i__2; ++iorg) {
	    j = -iw[iinput];
/* WE CAN DO THIS BECAUSE THE DIAGONAL IS NOW THE FIRST ENTRY. */
	    irow = iw2[j];
	    apos = posfac + (irow - 1) * (2 * nfront - irow + 2) / 2;
/* THE FOLLOWING LOOP GOES FROM 1 TO NZ BECAUSE THERE MAY BE DUPLICATES. */
	    i__3 = *nz;
	    for (idummy = 1; idummy <= i__3; ++idummy) {
		apos2 = apos + iw2[j] - irow;
		a[apos2] += a[ainput];
		++ainput;
		++iinput;
		if (iinput > *liw) {
		    goto L280;
		}
		j = iw[iinput];
		if (j < 0) {
		    goto L280;
		}
/* L270: */
	    }
L280:
	    ;
	}
/* RESET IW2 AND NUMASS. */
	numass += numorg;
	j1 = iwpos + 2;
	j2 = iwpos + nfront + 1;
	i__2 = j2;
	for (k = j1; k <= i__2; ++k) {
	    j = iw[k];
	    iw2[j] = 0;
/* L290: */
	}
/* PERFORM PIVOTING ON ASSEMBLED ELEMENT. */
/* NPIV IS THE NUMBER OF PIVOTS SO FAR SELECTED. */
/* LNPIV IS THE NUMBER OF PIVOTS SELECTED AFTER THE LAST PASS THROUGH */
/*     THE THE FOLLOWING LOOP. */
	lnpiv = -1;
	npiv = 0;
	i__2 = nass;
	for (kdummy = 1; kdummy <= i__2; ++kdummy) {
	    if (npiv == nass) {
		goto L660;
	    }
	    if (npiv == lnpiv) {
		goto L660;
	    }
	    lnpiv = npiv;
	    npivp1 = npiv + 1;
/* JPIV IS USED AS A FLAG TO INDICATE WHEN 2 BY 2 PIVOTING HAS OCCURRED */
/*     SO THAT IPIV IS INCREMENTED CORRECTLY. */
	    jpiv = 1;
/* NASS IS MAXIMUM POSSIBLE NUMBER OF PIVOTS. */
/* WE EITHER TAKE THE DIAGONAL ENTRY OR THE 2 BY 2 PIVOT WITH THE */
/*     LARGEST OFF-DIAGONAL AT EACH STAGE. */
/* EACH PASS THROUGH THIS LOOP TRIES TO CHOOSE ONE PIVOT. */
	    i__3 = nass;
	    for (ipiv = npivp1; ipiv <= i__3; ++ipiv) {
		--jpiv;
/* JUMP IF WE HAVE JUST PROCESSED A 2 BY 2 PIVOT. */
		if (jpiv == 1) {
		    goto L640;
		}
		i__4 = nfront - npiv;
		i__5 = ipiv - npiv;
		apos = posfac + (i__5 - 1) * (2 * i__4 - i__5 + 2) / 2;
/* IF THE USER HAS INDICATED THAT THE MATRIX IS DEFINITE, WE */
/*     DO NOT NEED TO TEST FOR STABILITY BUT WE DO CHECK TO SEE IF THE */
/*     PIVOT IS NON-ZERO OR HAS CHANGED SIGN. */
/*     IF IT IS ZERO, WE EXIT WITH AN ERROR. IF IT HAS CHANGED SIGN */
/*     AND U WAS SET NEGATIVE, THEN WE AGAIN EXIT IMMEDIATELY.  IF THE */
/*     PIVOT CHANGES SIGN AND U WAS ZERO, WE CONTINUE WITH THE */
/*     FACTORIZATION BUT PRINT A WARNING MESSAGE ON UNIT ICNTL(2). */
/* ISNPIV HOLDS A FLAG FOR THE SIGN OF THE PIVOTS TO DATE SO THAT */
/*     A SIGN CHANGE WHEN DECOMPOSING AN ALLEGEDLY DEFINITE MATRIX CAN */
/*     BE DETECTED. */
		if (uu > 0.) {
		    goto L320;
		}
		if ((d__1 = a[apos], abs(d__1)) <= cntl[3]) {
		    goto L790;
		}
/* JUMP IF THIS IS NOT THE FIRST PIVOT TO BE SELECTED. */
		if (ntotpv > 0) {
		    goto L300;
		}
/* SET ISNPIV. */
		if (a[apos] > 0.) {
		    isnpiv = 1;
		}
		if (a[apos] < 0.) {
		    isnpiv = -1;
		}
L300:
		if (a[apos] > 0. && isnpiv == 1) {
		    goto L560;
		}
		if (a[apos] < 0. && isnpiv == -1) {
		    goto L560;
		}
		if (info[1] != 2) {
		    info[2] = 0;
		}
		++info[2];
		info[1] = 2;
		i__ = ntotpv + 1;
		if (icntl[2] > 0 && info[2] <= 10) {
		    io___280.ciunit = icntl[2];
		    s_wsfe(&io___280);
		    do_fio(&c__1, (char *)&info[1], (ftnlen)sizeof(integer));
		    do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(integer));
		    e_wsfe();
		}
		isnpiv = -isnpiv;
		if (uu == 0.) {
		    goto L560;
		}
		goto L800;
L320:
		amax = 0.;
		tmax = amax;
/* FIND LARGEST ENTRY TO RIGHT OF DIAGONAL IN ROW OF PROSPECTIVE PIVOT */
/*     IN THE FULLY-SUMMED PART.  ALSO RECORD COLUMN OF THIS LARGEST */
/*     ENTRY. */
		j1 = apos + 1;
		j2 = apos + nass - ipiv;
		if (j2 < j1) {
		    goto L340;
		}
		i__4 = j2;
		for (jj = j1; jj <= i__4; ++jj) {
		    if ((d__1 = a[jj], abs(d__1)) <= amax) {
			goto L330;
		    }
		    jmax = ipiv + jj - j1 + 1;
		    amax = (d__1 = a[jj], abs(d__1));
L330:
		    ;
		}
/* DO SAME AS ABOVE FOR NON-FULLY-SUMMED PART ONLY HERE WE DO NOT NEED */
/*     TO RECORD COLUMN SO LOOP IS SIMPLER. */
L340:
		j1 = j2 + 1;
		j2 = apos + nfront - ipiv;
		if (j2 < j1) {
		    goto L360;
		}
		i__4 = j2;
		for (jj = j1; jj <= i__4; ++jj) {
/* Computing MAX */
		    d__2 = (d__1 = a[jj], abs(d__1));
		    tmax = max(d__2,tmax);
/* L350: */
		}
/* NOW CALCULATE LARGEST ENTRY IN OTHER PART OF ROW. */
L360:
		rmax = max(tmax,amax);
		apos1 = apos;
		kk = nfront - ipiv;
		lt = ipiv - (npiv + 1);
		if (lt == 0) {
		    goto L380;
		}
		i__4 = lt;
		for (k = 1; k <= i__4; ++k) {
		    ++kk;
		    apos1 -= kk;
/* Computing MAX */
		    d__2 = rmax, d__3 = (d__1 = a[apos1], abs(d__1));
		    rmax = max(d__2,d__3);
/* L370: */
		}
/* JUMP IF STABILITY TEST SATISFIED. */
L380:
/* Computing MAX */
		d__2 = cntl[3], d__3 = uu * rmax;
		if ((d__1 = a[apos], abs(d__1)) > max(d__2,d__3)) {
		    goto L450;
		}
/* CHECK BLOCK PIVOT OF ORDER 2 FOR STABILITY. */
		if (abs(amax) <= cntl[3]) {
		    goto L640;
		}
		i__4 = nfront - npiv;
		i__5 = jmax - npiv;
		apos2 = posfac + (i__5 - 1) * (2 * i__4 - i__5 + 2) / 2;
		detpiv = a[apos] * a[apos2] - amax * amax;
		thresh = abs(detpiv);
/* SET THRESH TO U TIMES THE RECIPROCAL OF THE MAX-NORM OF THE INVERSE */
/*     OF THE PROSPECTIVE BLOCK. */
/* Computing MAX */
		d__3 = (d__1 = a[apos], abs(d__1)) + amax, d__4 = (d__2 = a[
			apos2], abs(d__2)) + amax;
		thresh /= uu * max(d__3,d__4);
/* CHECK 2 BY 2 PIVOT FOR STABILITY. */
/* FIRST CHECK AGAINST ROW IPIV. */
		if (thresh <= rmax) {
		    goto L640;
		}
/* FIND LARGEST ENTRY IN ROW JMAX. */
/* FIND MAXIMUM TO THE RIGHT OF THE DIAGONAL. */
		rmax = 0.;
		j1 = apos2 + 1;
		j2 = apos2 + nfront - jmax;
		if (j2 < j1) {
		    goto L400;
		}
		i__4 = j2;
		for (jj = j1; jj <= i__4; ++jj) {
/* Computing MAX */
		    d__2 = rmax, d__3 = (d__1 = a[jj], abs(d__1));
		    rmax = max(d__2,d__3);
/* L390: */
		}
/* NOW CHECK TO THE LEFT OF THE DIAGONAL. */
/* WE USE TWO LOOPS TO AVOID TESTING FOR ROW IPIV INSIDE THE LOOP. */
L400:
		kk = nfront - jmax + 1;
		apos3 = apos2;
		jmxmip = jmax - ipiv - 1;
		if (jmxmip == 0) {
		    goto L420;
		}
		i__4 = jmxmip;
		for (k = 1; k <= i__4; ++k) {
		    apos2 -= kk;
		    ++kk;
/* Computing MAX */
		    d__2 = rmax, d__3 = (d__1 = a[apos2], abs(d__1));
		    rmax = max(d__2,d__3);
/* L410: */
		}
L420:
		ipmnp = ipiv - npiv - 1;
		if (ipmnp == 0) {
		    goto L440;
		}
		apos2 -= kk;
		++kk;
		i__4 = ipmnp;
		for (k = 1; k <= i__4; ++k) {
		    apos2 -= kk;
		    ++kk;
/* Computing MAX */
		    d__2 = rmax, d__3 = (d__1 = a[apos2], abs(d__1));
		    rmax = max(d__2,d__3);
/* L430: */
		}
L440:
		if (thresh <= rmax) {
		    goto L640;
		}
		pivsiz = 2;
		goto L460;
L450:
		pivsiz = 1;
L460:
		irow = ipiv - npiv;

/* PIVOT HAS BEEN CHOSEN.  IF BLOCK PIVOT OF ORDER 2, PIVSIZ IS EQUAL TO */
/*     TWO OTHERWISE PIVSIZ EQUALS ONE.. */
/* THE FOLLOWING LOOP MOVES THE PIVOT BLOCK TO THE TOP LEFT HAND CORNER */
/*     OF THE FRONTAL MATRIX. */
		i__4 = pivsiz;
		for (krow = 1; krow <= i__4; ++krow) {
/* WE JUMP IF SWOP IS NOT NECESSARY. */
		    if (irow == 1) {
			goto L530;
		    }
		    j1 = posfac + irow;
		    j2 = posfac + nfront - (npiv + 1);
		    if (j2 < j1) {
			goto L480;
		    }
		    apos2 = apos + 1;
/* SWOP PORTION OF ROWS WHOSE COLUMN INDICES ARE GREATER THAN LATER ROW. */
		    i__5 = j2;
		    for (jj = j1; jj <= i__5; ++jj) {
			swop = a[apos2];
			a[apos2] = a[jj];
			a[jj] = swop;
			++apos2;
/* L470: */
		    }
L480:
		    j1 = posfac + 1;
		    j2 = posfac + irow - 2;
		    apos2 = apos;
		    kk = nfront - (irow + npiv);
		    if (j2 < j1) {
			goto L500;
		    }
/* SWOP PORTION OF ROWS/COLUMNS WHOSE INDICES LIE BETWEEN THE TWO ROWS. */
		    i__5 = j2;
		    for (jjj = j1; jjj <= i__5; ++jjj) {
			jj = j2 - jjj + j1;
			++kk;
			apos2 -= kk;
			swop = a[apos2];
			a[apos2] = a[jj];
			a[jj] = swop;
/* L490: */
		    }
L500:
		    if (npiv == 0) {
			goto L520;
		    }
		    apos1 = posfac;
		    ++kk;
		    apos2 -= kk;
/* SWOP PORTION OF COLUMNS WHOSE INDICES ARE LESS THAN EARLIER ROW. */
		    i__5 = npiv;
		    for (jj = 1; jj <= i__5; ++jj) {
			++kk;
			apos1 -= kk;
			apos2 -= kk;
			swop = a[apos2];
			a[apos2] = a[apos1];
			a[apos1] = swop;
/* L510: */
		    }
/* SWOP DIAGONALS AND INTEGER INDEXING INFORMATION */
L520:
		    swop = a[apos];
		    a[apos] = a[posfac];
		    a[posfac] = swop;
		    ipos = iwpos + npiv + 2;
		    iexch = iwpos + irow + npiv + 1;
		    iswop = iw[ipos];
		    iw[ipos] = iw[iexch];
		    iw[iexch] = iswop;
L530:
		    if (pivsiz == 1) {
			goto L550;
		    }
/* SET VARIABLES FOR THE SWOP OF SECOND ROW OF BLOCK PIVOT. */
		    if (krow == 2) {
			goto L540;
		    }
		    irow = jmax - (npiv + 1);
		    jpos = posfac;
		    posfac = posfac + nfront - npiv;
		    ++npiv;
		    apos = apos3;
		    goto L550;
/* RESET VARIABLES PREVIOUSLY SET FOR SECOND PASS. */
L540:
		    --npiv;
		    posfac = jpos;
L550:
		    ;
		}

		if (pivsiz == 2) {
		    goto L600;
		}
/* PERFORM THE ELIMINATION USING ENTRY (IPIV,IPIV) AS PIVOT. */
/* WE STORE U AND DINVERSE. */
L560:
		a[posfac] = 1. / a[posfac];
		if (a[posfac] < 0.) {
		    ++neig;
		}
		j1 = posfac + 1;
		j2 = posfac + nfront - (npiv + 1);
		if (j2 < j1) {
		    goto L590;
		}
		ibeg = j2 + 1;
		i__4 = j2;
		for (jj = j1; jj <= i__4; ++jj) {
		    amult = -a[jj] * a[posfac];
		    iend = ibeg + nfront - (npiv + jj - j1 + 2);
/* THE FOLLOWING SPECIAL COMMENT FORCES VECTORIZATION ON THE CRAY-1. */
/* DIR$ IVDEP */
		    i__5 = iend;
		    for (irow = ibeg; irow <= i__5; ++irow) {
			jcol = jj + irow - ibeg;
			a[irow] += amult * a[jcol];
/* L570: */
		    }
		    ibeg = iend + 1;
		    a[jj] = amult;
/* L580: */
		}
L590:
		++npiv;
		++ntotpv;
		jpiv = 1;
		posfac = posfac + nfront - npiv + 1;
		goto L640;
/* PERFORM ELIMINATION USING BLOCK PIVOT OF ORDER TWO. */
/* REPLACE BLOCK PIVOT BY ITS INVERSE. */
/* SET FLAG TO INDICATE USE OF 2 BY 2 PIVOT IN IW. */
L600:
		ipos = iwpos + npiv + 2;
		++ntwo;
		iw[ipos] = -iw[ipos];
		pospv1 = posfac;
		pospv2 = posfac + nfront - npiv;
		swop = a[pospv2];
		if (detpiv < 0.) {
		    ++neig;
		}
		if (detpiv > 0. && swop < 0.) {
		    neig += 2;
		}
		a[pospv2] = a[pospv1] / detpiv;
		a[pospv1] = swop / detpiv;
		a[pospv1 + 1] = -a[pospv1 + 1] / detpiv;
		j1 = pospv1 + 2;
		j2 = pospv1 + nfront - (npiv + 1);
		if (j2 < j1) {
		    goto L630;
		}
		jj1 = pospv2;
		ibeg = pospv2 + nfront - (npiv + 1);
		i__4 = j2;
		for (jj = j1; jj <= i__4; ++jj) {
		    ++jj1;
		    amult1 = -(a[pospv1] * a[jj] + a[pospv1 + 1] * a[jj1]);
		    amult2 = -(a[pospv1 + 1] * a[jj] + a[pospv2] * a[jj1]);
		    iend = ibeg + nfront - (npiv + jj - j1 + 3);
/* THE FOLLOWING SPECIAL COMMENT FORCES VECTORIZATION ON THE CRAY-1. */
/* DIR$ IVDEP */
		    i__5 = iend;
		    for (irow = ibeg; irow <= i__5; ++irow) {
			k1 = jj + irow - ibeg;
			k2 = jj1 + irow - ibeg;
			a[irow] = a[irow] + amult1 * a[k1] + amult2 * a[k2];
/* L610: */
		    }
		    ibeg = iend + 1;
		    a[jj] = amult1;
		    a[jj1] = amult2;
/* L620: */
		}
L630:
		npiv += 2;
		ntotpv += 2;
		jpiv = 2;
		posfac = pospv2 + nfront - npiv + 1;
L640:
		;
	    }
/* L650: */
	}
/* END OF MAIN ELIMINATION LOOP. */

L660:
	if (npiv != 0) {
	    ++nblk;
	}
	ioldps = iwpos;
	iwpos = iwpos + nfront + 2;
	if (npiv == 0) {
	    goto L690;
	}
	if (npiv > 1) {
	    goto L680;
	}
	iw[ioldps] = -iw[ioldps];
	i__2 = nfront;
	for (k = 1; k <= i__2; ++k) {
	    j1 = ioldps + k;
	    iw[j1] = iw[j1 + 1];
/* L670: */
	}
	--iwpos;
	goto L690;
L680:
	iw[ioldps + 1] = npiv;
/* COPY REMAINDER OF ELEMENT TO TOP OF STACK */
L690:
	liell = nfront - npiv;
	if (liell == 0 || iass == *nsteps) {
	    goto L750;
	}
	if (iwpos + liell < istk) {
	    goto L700;
	}
	ma27pd_(&a[1], &iw[1], &istk, &istk2, &iinput, &c__2, &ncmpbr, &
		ncmpbi);
L700:
	istk = istk - liell - 1;
	iw[istk] = liell;
	j1 = istk;
	kk = iwpos - liell - 1;
/* THE FOLLOWING SPECIAL COMMENT FORCES VECTORIZATION ON THE CRAY-1. */
/* DIR$ IVDEP */
	i__2 = liell;
	for (k = 1; k <= i__2; ++k) {
	    ++j1;
	    ++kk;
	    iw[j1] = iw[kk];
/* L710: */
	}
/* WE COPY IN REVERSE DIRECTION TO AVOID OVERWRITE PROBLEMS. */
	laell = (liell + 1) * liell / 2;
	kk = posfac + laell;
	if (kk != astk) {
	    goto L720;
	}
	astk -= laell;
	goto L740;
/* THE MOVE AND ZEROING OF ARRAY A IS PERFORMED WITH TWO LOOPS SO */
/* THAT THE CRAY-1 WILL VECTORIZE THEM SAFELY. */
L720:
	kmax = kk - 1;
/* THE FOLLOWING SPECIAL COMMENT FORCES VECTORIZATION ON THE CRAY-1. */
/* DIR$ IVDEP */
	i__2 = laell;
	for (k = 1; k <= i__2; ++k) {
	    --kk;
	    --astk;
	    a[astk] = a[kk];
/* L730: */
	}
/* Computing MIN */
	i__2 = kmax, i__3 = astk - 1;
	kmax = min(i__2,i__3);
	i__2 = kmax;
	for (k = kk; k <= i__2; ++k) {
	    a[k] = 0.;
/* L735: */
	}
L740:
/* Computing MIN */
	i__2 = azero, i__3 = astk - 1;
	azero = min(i__2,i__3);
L750:
	if (npiv == 0) {
	    iwpos = ioldps;
	}
/* L760: */
    }

/* END OF LOOP ON TREE NODES. */

    iw[1] = nblk;
    if (ntwo > 0) {
	iw[1] = -nblk;
    }
    nrlbdu = posfac - 1;
    nirbdu = iwpos - 1;
    if (ntotpv == *n) {
	goto L810;
    }
    info[1] = 3;
    info[2] = ntotpv;
    goto L810;
/* **** ERROR RETURNS **** */
L770:
    info[1] = -3;
    goto L810;
L780:
    info[1] = -4;
/* Computing MAX */
    i__1 = posfac + lnass, i__2 = apos2 - ltopst + 2;
    info[2] = *la + max(i__1,i__2) - astk;
    goto L810;
L790:
    info[1] = -5;
    info[2] = ntotpv + 1;
    goto L810;
L800:
    info[1] = -6;
    info[2] = ntotpv + 1;
L810:
    info[9] = nrlbdu;
    info[10] = nirbdu;
    info[12] = ncmpbr;
    info[13] = ncmpbi;
    info[14] = ntwo;
    info[15] = neig;
    return 0;
} /* ma27od_ */

/* Subroutine */ int ma27pd_(doublereal *a, integer *iw, integer *j1, integer 
	*j2, integer *itop, integer *ireal, integer *ncmpbr, integer *ncmpbi)
{
    /* System generated locals */
    integer i__1;

    /* Local variables */
    static integer jj, jjj, ipos;

/* THIS SUBROUTINE PERFORMS A VERY SIMPLE COMPRESS (BLOCK MOVE). */
/*     ENTRIES J1 TO J2 (INCL.) IN A OR IW AS APPROPRIATE ARE MOVED TO */
/*     OCCUPY THE POSITIONS IMMEDIATELY PRIOR TO POSITION ITOP. */
/* A/IW HOLD THE ARRAY BEING COMPRESSED. */
/* J1/J2 DEFINE THE ENTRIES BEING MOVED. */
/* ITOP DEFINES THE POSITION IMMEDIATELY AFTER THE POSITIONS TO WHICH */
/*     J1 TO J2 ARE MOVED. */
/* IREAL MUST BE SET BY THE USER TO 2 IF THE MOVE IS ON ARRAY IW, */
/*     ANY OTHER VALUE WILL PERFORM THE MOVE ON A. */
/* NCMPBR and NCMPBI, see INFO(12) and INFO(13) in MA27A/AD (ACCUMULATE */
/*     THE NUMBER OF COMPRESSES OF THE REALS AND INTEGERS PERFORMED BY */
/*     MA27B/BD. */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. Executable Statements .. */
    /* Parameter adjustments */
    --iw;
    --a;

    /* Function Body */
    ipos = *itop - 1;
    if (*j2 == ipos) {
	goto L50;
    }
    if (*ireal == 2) {
	goto L20;
    }
    ++(*ncmpbr);
    if (*j1 > *j2) {
	goto L40;
    }
    i__1 = *j2;
    for (jjj = *j1; jjj <= i__1; ++jjj) {
	jj = *j2 - jjj + *j1;
	a[ipos] = a[jj];
	--ipos;
/* L10: */
    }
    goto L40;
L20:
    ++(*ncmpbi);
    if (*j1 > *j2) {
	goto L40;
    }
    i__1 = *j2;
    for (jjj = *j1; jjj <= i__1; ++jjj) {
	jj = *j2 - jjj + *j1;
	iw[ipos] = iw[jj];
	--ipos;
/* L30: */
    }
L40:
    *j2 = *itop - 1;
    *j1 = ipos + 1;
L50:
    return 0;
} /* ma27pd_ */

/* Subroutine */ int ma27qd_(integer *n, doublereal *a, integer *la, integer *
	iw, integer *liw, doublereal *w, integer *maxfnt, doublereal *rhs, 
	integer *iw2, integer *nblk, integer *latop, integer *icntl)
{
    /* System generated locals */
    integer i__1, i__2, i__3;

    /* Local variables */
    static integer j, k, j1, j2, j3, k1, k2, k3;
    static doublereal w1, w2;
    static integer jj, ifr, ist, iblk, apos, irhs, ilvl, ipiv, jpiv, ipos, 
	    npiv, irow, liell;

/* THIS SUBROUTINE PERFORMS FORWARD ELIMINATION */
/*     USING THE FACTOR U TRANSPOSE STORED IN A/IA AFTER MA27B/BD. */

/* N      - MUST BE SET TO THE ORDER OF THE MATRIX. NOT ALTERED */
/*          BY MA27Q/QD. */
/* A      - MUST BE SET TO HOLD THE REAL VALUES */
/*          CORRESPONDING TO THE FACTORS OF DINVERSE AND U.  THIS MUST BE */
/*          UNCHANGED SINCE THE PRECEDING CALL TO MA27B/BD.  NOT ALTERED */
/*          BY MA27Q/QD. */
/* LA     - LENGTH OF ARRAY A.  NOT ALTERED BY MA27Q/QD. */
/* IW     - HOLDS THE INTEGER INDEXING */
/*          INFORMATION FOR THE MATRIX FACTORS IN A.  THIS MUST BE */
/*          UNCHANGED SINCE THE PRECEDING CALL TO MA27B/BD.  NOT ALTERED */
/*          BY MA27Q/QD. */
/* LIW    - LENGTH OF ARRAY IW.  NOT ALTERED BY MA27Q/QD. */
/* W      - USED */
/*          AS WORKSPACE BY MA27Q/QD TO HOLD THE COMPONENTS OF THE RIGHT */
/*          HAND SIDES CORRESPONDING TO CURRENT BLOCK PIVOTAL ROWS. */
/* MAXFNT - MUST BE SET TO THE LARGEST NUMBER OF */
/*          VARIABLES IN ANY BLOCK PIVOT ROW.  THIS VALUE WILL HAVE */
/*          BEEN OUTPUT BY MA27B/BD.  NOT ALTERED BY MA27Q/QD. */
/* RHS    - ON INPUT, */
/*          MUST BE SET TO HOLD THE RIGHT HAND SIDES FOR THE EQUATIONS */
/*          WHICH THE USER DESIRES TO SOLVE.  ON OUTPUT, RHS WILL HOLD */
/*          THE MODIFIED VECTORS CORRESPONDING TO PERFORMING FORWARD */
/*          ELIMINATION ON THE RIGHT HAND SIDES. */
/* IW2    - NEED NOT BE SET ON ENTRY. ON EXIT IW2(I) (I = 1,NBLK) */
/*          WILL HOLD POINTERS TO THE */
/*          BEGINNING OF EACH BLOCK PIVOT IN ARRAY IW. */
/* NBLK   - NUMBER OF BLOCK PIVOT ROWS. NOT ALTERED BY MA27Q/QD. */
/* LATOP  - NEED NOT BE SET ON ENTRY. ON EXIT, IT IS THE POSITION IN */
/*          A OF THE LAST ENTRY IN THE  FACTORS. IT MUST BE PASSED */
/*          UNCHANGED TO MA27R/RD. */
/* ICNTL is an INTEGER array of length 30, see MA27A/AD. */
/*   ICNTL(IFRLVL+I) I=1,20 IS USED TO CONTROL WHETHER DIRECT OR INDIRECT */
/*     ACCESS IS USED BY MA27C/CD.  INDIRECT ACCESS IS EMPLOYED */
/*     IN FORWARD AND BACK SUBSTITUTION RESPECTIVELY IF THE SIZE OF */
/*     A BLOCK IS LESS THAN ICNTL(IFRLVL+MIN(10,NPIV)) AND */
/*     ICNTL(IFRLVL+10+MIN(10,NPIV)) RESPECTIVELY, WHERE NPIV IS THE */
/*     NUMBER OF PIVOTS IN THE BLOCK. */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */
/* APOS. RUNNING POINTER TO CURRENT PIVOT POSITION IN ARRAY A. */
/* IPOS. RUNNING POINTER TO BEGINNING OF BLOCK PIVOT ROW IN IW. */
    /* Parameter adjustments */
    --rhs;
    --a;
    --iw;
    --w;
    --iw2;
    --icntl;

    /* Function Body */
    apos = 1;
    ipos = 1;
    j2 = 0;
    iblk = 0;
    npiv = 0;
    i__1 = *n;
    for (irow = 1; irow <= i__1; ++irow) {
	if (npiv > 0) {
	    goto L90;
	}
	++iblk;
	if (iblk > *nblk) {
	    goto L150;
	}
	ipos = j2 + 1;
/* SET UP POINTER IN PREPARATION FOR BACK SUBSTITUTION. */
	iw2[iblk] = ipos;
/* ABS(LIELL) IS NUMBER OF VARIABLES (COLUMNS) IN BLOCK PIVOT ROW. */
	liell = -iw[ipos];
/* NPIV IS NUMBER OF PIVOTS (ROWS) IN BLOCK PIVOT. */
	npiv = 1;
	if (liell > 0) {
	    goto L10;
	}
	liell = -liell;
	++ipos;
	npiv = iw[ipos];
L10:
	j1 = ipos + 1;
	j2 = ipos + liell;
	ilvl = min(npiv,10);
	if (liell < icntl[ilvl + 5]) {
	    goto L90;
	}

/* PERFORM OPERATIONS USING DIRECT ADDRESSING. */

/* LOAD APPROPRIATE COMPONENTS OF RIGHT HAND SIDES INTO ARRAY W. */
	ifr = 0;
	i__2 = j2;
	for (jj = j1; jj <= i__2; ++jj) {
	    j = (i__3 = iw[jj], abs(i__3));
	    ++ifr;
	    w[ifr] = rhs[j];
/* L20: */
	}
/* JPIV IS USED AS A FLAG SO THAT IPIV IS INCREMENTED CORRECTLY AFTER */
/* THE USE OF A 2 BY 2 PIVOT. */
	jpiv = 1;
	j3 = j1;
/* PERFORM OPERATIONS. */
	i__2 = npiv;
	for (ipiv = 1; ipiv <= i__2; ++ipiv) {
	    --jpiv;
	    if (jpiv == 1) {
		goto L70;
	    }
/* JUMP IF WE HAVE A 2 BY 2 PIVOT. */
	    if (iw[j3] < 0) {
		goto L40;
	    }
/* PERFORM FORWARD SUBSTITUTION USING 1 BY 1 PIVOT. */
	    jpiv = 1;
	    ++j3;
	    ++apos;
	    ist = ipiv + 1;
	    if (liell < ist) {
		goto L70;
	    }
	    w1 = w[ipiv];
	    k = apos;
	    i__3 = liell;
	    for (j = ist; j <= i__3; ++j) {
		w[j] += a[k] * w1;
		++k;
/* L30: */
	    }
	    apos = apos + liell - ist + 1;
	    goto L70;
/* PERFORM OPERATIONS WITH 2 BY 2 PIVOT. */
L40:
	    jpiv = 2;
	    j3 += 2;
	    apos += 2;
	    ist = ipiv + 2;
	    if (liell < ist) {
		goto L60;
	    }
	    w1 = w[ipiv];
	    w2 = w[ipiv + 1];
	    k1 = apos;
	    k2 = apos + liell - ipiv;
	    i__3 = liell;
	    for (j = ist; j <= i__3; ++j) {
		w[j] = w[j] + w1 * a[k1] + w2 * a[k2];
		++k1;
		++k2;
/* L50: */
	    }
L60:
	    apos = apos + (liell - ist + 1 << 1) + 1;
L70:
	    ;
	}
/* RELOAD W BACK INTO RHS. */
	ifr = 0;
	i__2 = j2;
	for (jj = j1; jj <= i__2; ++jj) {
	    j = (i__3 = iw[jj], abs(i__3));
	    ++ifr;
	    rhs[j] = w[ifr];
/* L80: */
	}
	npiv = 0;
	goto L140;

/* PERFORM OPERATIONS USING INDIRECT ADDRESSING. */

/* JUMP IF WE HAVE A 2 BY 2 PIVOT. */
L90:
	if (iw[j1] < 0) {
	    goto L110;
	}
/* PERFORM FORWARD SUBSTITUTION USING 1 BY 1 PIVOT. */
	--npiv;
	++apos;
	++j1;
	if (j1 > j2) {
	    goto L140;
	}
	irhs = iw[j1 - 1];
	w1 = rhs[irhs];
	k = apos;
	i__2 = j2;
	for (j = j1; j <= i__2; ++j) {
	    irhs = (i__3 = iw[j], abs(i__3));
	    rhs[irhs] += a[k] * w1;
	    ++k;
/* L100: */
	}
	apos = apos + j2 - j1 + 1;
	goto L140;
/* PERFORM OPERATIONS WITH 2 BY 2 PIVOT */
L110:
	npiv += -2;
	j1 += 2;
	apos += 2;
	if (j1 > j2) {
	    goto L130;
	}
	irhs = -iw[j1 - 2];
	w1 = rhs[irhs];
	irhs = iw[j1 - 1];
	w2 = rhs[irhs];
	k1 = apos;
	k3 = apos + j2 - j1 + 2;
	i__2 = j2;
	for (j = j1; j <= i__2; ++j) {
	    irhs = (i__3 = iw[j], abs(i__3));
	    rhs[irhs] = rhs[irhs] + w1 * a[k1] + w2 * a[k3];
	    ++k1;
	    ++k3;
/* L120: */
	}
L130:
	apos = apos + (j2 - j1 + 1 << 1) + 1;
L140:
	;
    }
L150:
    *latop = apos - 1;
    return 0;
} /* ma27qd_ */

/* Subroutine */ int ma27rd_(integer *n, doublereal *a, integer *la, integer *
	iw, integer *liw, doublereal *w, integer *maxfnt, doublereal *rhs, 
	integer *iw2, integer *nblk, integer *latop, integer *icntl)
{
    /* System generated locals */
    integer i__1, i__2, i__3;

    /* Local variables */
    static integer j, k, j1, j2;
    static doublereal w1, w2;
    static integer jj, jj1, jj2, ifr, ist, iblk, apos, irhs, ilvl, ipiv, jpiv,
	     loop, ipos, jpos, npiv, apos2, i1rhs, i2rhs, liell, iirhs, iipiv;

/* THIS SUBROUTINE PERFORMS BACKWARD ELIMINATION OPERATIONS */
/*     USING THE FACTORS DINVERSE AND U */
/*     STORED IN A/IW AFTER MA27B/BD. */

/* N      - MUST BE SET TO THE ORDER OF THE MATRIX. NOT ALTERED */
/*          BY MA27R/RD. */
/* A      - MUST BE SET TO HOLD THE REAL VALUES CORRESPONDING */
/*          TO THE FACTORS OF DINVERSE AND U.  THIS MUST BE */
/*          UNCHANGED SINCE THE PRECEDING CALL TO MA27B/BD.  NOT ALTERED */
/*          BY MA27R/RD. */
/* LA     - LENGTH OF ARRAY A. NOT ALTERED BY MA27R/RD. */
/* IW     - HOLDS THE INTEGER INDEXING */
/*          INFORMATION FOR THE MATRIX FACTORS IN A.  THIS MUST BE */
/*          UNCHANGED SINCE THE PRECEDING CALL TO MA27B/BD.  NOT ALTERED */
/*          BY MA27R/RD. */
/* LIW    - LENGTH OF ARRAY IW.  NOT ALTERED BY MA27R/RD. */
/* W      - USED */
/*          AS WORKSPACE BY MA27R/RD TO HOLD THE COMPONENTS OF THE RIGHT */
/*          HAND SIDES CORRESPONDING TO CURRENT BLOCK PIVOTAL ROWS. */
/* MAXFNT - INTEGER VARIABLE.  MUST BE SET TO THE LARGEST NUMBER OF */
/*          VARIABLES IN ANY BLOCK PIVOT ROW.  THIS VALUE WAS GIVEN */
/*          ON OUTPUT FROM MA27B/BD.  NOT ALTERED BY MA27R/RD. */
/* RHS    - ON INPUT, */
/*          MUST BE SET TO HOLD THE RIGHT HAND SIDE MODIFIED BY THE */
/*          FORWARD SUBSTITUTION OPERATIONS.  ON OUTPUT, RHS WILL HOLD */
/*          THE SOLUTION VECTOR. */
/* IW2    - ON ENTRY IW2(I) (I = 1,NBLK) */
/*          MUST HOLD POINTERS TO THE */
/*          BEGINNING OF EACH BLOCK PIVOT IN ARRAY IW, AS SET BY */
/*          MA27Q/QD. */
/* NBLK   - NUMBER OF BLOCK PIVOT ROWS. NOT ALTERED BY MA27R/RD. */
/* LATOP  - IT IS THE POSITION IN */
/*          A OF THE LAST ENTRY IN THE  FACTORS. IT MUST BE UNCHANGED */
/*          SINCE THE CALL TO MA27Q/QD.  IT IS NOT ALTERED BY MA27R/RD. */
/* ICNTL is an INTEGER array of length 30, see MA27A/AD. */
/*   ICNTL(IFRLVL+I) I=1,20 IS USED TO CONTROL WHETHER DIRECT OR INDIRECT */
/*     ACCESS IS USED BY MA27C/CD.  INDIRECT ACCESS IS EMPLOYED */
/*     IN FORWARD AND BACK SUBSTITUTION RESPECTIVELY IF THE SIZE OF */
/*     A BLOCK IS LESS THAN ICNTL(IFRLVL+MIN(10,NPIV)) AND */
/*     ICNTL(IFRLVL+10+MIN(10,NPIV)) RESPECTIVELY, WHERE NPIV IS THE */
/*     NUMBER OF PIVOTS IN THE BLOCK. */


/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */
/* APOS. RUNNING POINTER TO CURRENT PIVOT POSITION IN ARRAY A. */
/* IPOS. RUNNING POINTER TO BEGINNING OF CURRENT BLOCK PIVOT ROW. */
    /* Parameter adjustments */
    --rhs;
    --a;
    --iw;
    --w;
    --iw2;
    --icntl;

    /* Function Body */
    apos = *latop + 1;
    npiv = 0;
    iblk = *nblk + 1;
/* RUN THROUGH BLOCK PIVOT ROWS IN THE REVERSE ORDER. */
    i__1 = *n;
    for (loop = 1; loop <= i__1; ++loop) {
	if (npiv > 0) {
	    goto L110;
	}
	--iblk;
	if (iblk < 1) {
	    goto L190;
	}
	ipos = iw2[iblk];
/* ABS(LIELL) IS NUMBER OF VARIABLES (COLUMNS) IN BLOCK PIVOT ROW. */
	liell = -iw[ipos];
/* NPIV IS NUMBER OF PIVOTS (ROWS) IN BLOCK PIVOT. */
	npiv = 1;
	if (liell > 0) {
	    goto L10;
	}
	liell = -liell;
	++ipos;
	npiv = iw[ipos];
L10:
	jpos = ipos + npiv;
	j2 = ipos + liell;
	ilvl = min(10,npiv) + 10;
	if (liell < icntl[ilvl + 5]) {
	    goto L110;
	}

/* PERFORM OPERATIONS USING DIRECT ADDRESSING. */

	j1 = ipos + 1;
/* LOAD APPROPRIATE COMPONENTS OF RHS INTO W. */
	ifr = 0;
	i__2 = j2;
	for (jj = j1; jj <= i__2; ++jj) {
	    j = (i__3 = iw[jj], abs(i__3));
	    ++ifr;
	    w[ifr] = rhs[j];
/* L20: */
	}
/* JPIV IS USED AS A FLAG SO THAT IPIV IS INCREMENTED CORRECTLY AFTER */
/*     THE USE OF A 2 BY 2 PIVOT. */
	jpiv = 1;
/* PERFORM ELIMINATIONS. */
	i__2 = npiv;
	for (iipiv = 1; iipiv <= i__2; ++iipiv) {
	    --jpiv;
	    if (jpiv == 1) {
		goto L90;
	    }
	    ipiv = npiv - iipiv + 1;
	    if (ipiv == 1) {
		goto L30;
	    }
/* JUMP IF WE HAVE A 2 BY 2 PIVOT. */
	    if (iw[jpos - 1] < 0) {
		goto L60;
	    }
/* PERFORM BACK-SUBSTITUTION USING 1 BY 1 PIVOT. */
L30:
	    jpiv = 1;
	    apos -= liell + 1 - ipiv;
	    ist = ipiv + 1;
	    w1 = w[ipiv] * a[apos];
	    if (liell < ist) {
		goto L50;
	    }
	    jj1 = apos + 1;
	    i__3 = liell;
	    for (j = ist; j <= i__3; ++j) {
		w1 += a[jj1] * w[j];
		++jj1;
/* L40: */
	    }
L50:
	    w[ipiv] = w1;
	    --jpos;
	    goto L90;
/* PERFORM BACK-SUBSTITUTION OPERATIONS WITH 2 BY 2 PIVOT */
L60:
	    jpiv = 2;
	    apos2 = apos - (liell + 1 - ipiv);
	    apos = apos2 - (liell + 2 - ipiv);
	    ist = ipiv + 1;
	    w1 = w[ipiv - 1] * a[apos] + w[ipiv] * a[apos + 1];
	    w2 = w[ipiv - 1] * a[apos + 1] + w[ipiv] * a[apos2];
	    if (liell < ist) {
		goto L80;
	    }
	    jj1 = apos + 2;
	    jj2 = apos2 + 1;
	    i__3 = liell;
	    for (j = ist; j <= i__3; ++j) {
		w1 += w[j] * a[jj1];
		w2 += w[j] * a[jj2];
		++jj1;
		++jj2;
/* L70: */
	    }
L80:
	    w[ipiv - 1] = w1;
	    w[ipiv] = w2;
	    jpos += -2;
L90:
	    ;
	}
/* RELOAD WORKING VECTOR INTO SOLUTION VECTOR. */
	ifr = 0;
	i__2 = j2;
	for (jj = j1; jj <= i__2; ++jj) {
	    j = (i__3 = iw[jj], abs(i__3));
	    ++ifr;
	    rhs[j] = w[ifr];
/* L100: */
	}
	npiv = 0;
	goto L180;

/* PERFORM OPERATIONS USING INDIRECT ADDRESSING. */

L110:
	if (npiv == 1) {
	    goto L120;
	}
/* JUMP IF WE HAVE A 2 BY 2 PIVOT. */
	if (iw[jpos - 1] < 0) {
	    goto L150;
	}
/* PERFORM BACK-SUBSTITUTION USING 1 BY 1 PIVOT. */
L120:
	--npiv;
	apos -= j2 - jpos + 1;
	iirhs = iw[jpos];
	w1 = rhs[iirhs] * a[apos];
	j1 = jpos + 1;
	if (j1 > j2) {
	    goto L140;
	}
	k = apos + 1;
	i__2 = j2;
	for (j = j1; j <= i__2; ++j) {
	    irhs = (i__3 = iw[j], abs(i__3));
	    w1 += a[k] * rhs[irhs];
	    ++k;
/* L130: */
	}
L140:
	rhs[iirhs] = w1;
	--jpos;
	goto L180;
/* PERFORM OPERATIONS WITH 2 BY 2 PIVOT */
L150:
	npiv += -2;
	apos2 = apos - (j2 - jpos + 1);
	apos = apos2 - (j2 - jpos + 2);
	i1rhs = -iw[jpos - 1];
	i2rhs = iw[jpos];
	w1 = rhs[i1rhs] * a[apos] + rhs[i2rhs] * a[apos + 1];
	w2 = rhs[i1rhs] * a[apos + 1] + rhs[i2rhs] * a[apos2];
	j1 = jpos + 1;
	if (j1 > j2) {
	    goto L170;
	}
	jj1 = apos + 2;
	jj2 = apos2 + 1;
	i__2 = j2;
	for (j = j1; j <= i__2; ++j) {
	    irhs = (i__3 = iw[j], abs(i__3));
	    w1 += rhs[irhs] * a[jj1];
	    w2 += rhs[irhs] * a[jj2];
	    ++jj1;
	    ++jj2;
/* L160: */
	}
L170:
	rhs[i1rhs] = w1;
	rhs[i2rhs] = w2;
	jpos += -2;
L180:
	;
    }
L190:
    return 0;
} /* ma27rd_ */