mmd.c

LU分解求解矩阵方程组的解

#ifdef _CRAY
#define int short
typedef short shortint;
#elif defined (_LONGINT)
#define int long
typedef long shortint;
#else
typedef int shortint;
#endif

/* *************************************************************** */
/* *************************************************************** */
/* ****     GENMMD ..... MULTIPLE MINIMUM EXTERNAL DEGREE     **** */
/* *************************************************************** */
/* *************************************************************** */

/*     AUTHOR - JOSEPH W.H. LIU */
/*              DEPT OF COMPUTER SCIENCE, YORK UNIVERSITY. */

/*     PURPOSE - THIS ROUTINE IMPLEMENTS THE MINIMUM DEGREE */
/*        ALGORITHM.  IT MAKES USE OF THE IMPLICIT REPRESENTATION */
/*        OF ELIMINATION GRAPHS BY QUOTIENT GRAPHS, AND THE */
/*        NOTION OF INDISTINGUISHABLE NODES.  IT ALSO IMPLEMENTS */
/*        THE MODIFICATIONS BY MULTIPLE ELIMINATION AND MINIMUM */
/*        EXTERNAL DEGREE. */
/*        --------------------------------------------- */
/*        CAUTION - THE ADJACENCY VECTOR ADJNCY WILL BE */
/*        DESTROYED. */
/*        --------------------------------------------- */

/*     INPUT PARAMETERS - */
/*        NEQNS  - NUMBER OF EQUATIONS. */
/*        (XADJ,ADJNCY) - THE ADJACENCY STRUCTURE. */
/*        DELTA  - TOLERANCE VALUE FOR MULTIPLE ELIMINATION. */
/*        MAXINT - MAXIMUM MACHINE REPRESENTABLE (SHORT) INTEGER */
/*                 (ANY SMALLER ESTIMATE WILL DO) FOR MARKING */
/*                 NODES. */

/*     OUTPUT PARAMETERS - */
/*        PERM   - THE MINIMUM DEGREE ORDERING. */
/*        INVP   - THE INVERSE OF PERM. */
/*        NOFSUB - AN UPPER BOUND ON THE NUMBER OF NONZERO */
/*                 SUBSCRIPTS FOR THE COMPRESSED STORAGE SCHEME. */

/*     WORKING PARAMETERS - */
/*        DHEAD  - VECTOR FOR HEAD OF DEGREE LISTS. */
/*        INVP   - USED TEMPORARILY FOR DEGREE FORWARD LINK. */
/*        PERM   - USED TEMPORARILY FOR DEGREE BACKWARD LINK. */
/*        QSIZE  - VECTOR FOR SIZE OF SUPERNODES. */
/*        LLIST  - VECTOR FOR TEMPORARY LINKED LISTS. */
/*        MARKER - A TEMPORARY MARKER VECTOR. */

/*     PROGRAM SUBROUTINES - */
/*        MMDELM, MMDINT, MMDNUM, MMDUPD. */

/* *************************************************************** */

/* Subroutine */ int genmmd_dist_(int *neqns, int *xadj, shortint *adjncy, 
	shortint *invp, shortint *perm, int *delta, shortint *dhead, 
	shortint *qsize, shortint *llist, shortint *marker, int *maxint, 
	int *nofsub)
{
    /* System generated locals */
    int i__1;

    /* Local variables */
    static int mdeg, ehead, i, mdlmt, mdnode;
    extern /* Subroutine */ int mmdelm_dist(int *, int *, shortint *, 
	    shortint *, shortint *, shortint *, shortint *, shortint *, 
	    shortint *, int *, int *), mmdupd_dist(int *, int *, 
	    int *, shortint *, int *, int *, shortint *, shortint 
	    *, shortint *, shortint *, shortint *, shortint *, int *, 
	    int *), mmdint_dist(int *, int *, shortint *, shortint *, 
	    shortint *, shortint *, shortint *, shortint *, shortint *), 
	    mmdnum_dist(int *, shortint *, shortint *, shortint *);
    static int nextmd, tag, num;


/* *************************************************************** */


/* *************************************************************** */

    /* Parameter adjustments */
    --marker;
    --llist;
    --qsize;
    --dhead;
    --perm;
    --invp;
    --adjncy;
    --xadj;

    /* Function Body */
    if (*neqns <= 0) {
	return 0;
    }

/*        ------------------------------------------------ */
/*        INITIALIZATION FOR THE MINIMUM DEGREE ALGORITHM. */
/*        ------------------------------------------------ */
    *nofsub = 0;
    mmdint_dist(neqns, &xadj[1], &adjncy[1], &dhead[1], &invp[1], &perm[1],
		&qsize[1], &llist[1], &marker[1]);

/*        ---------------------------------------------- */
/*        NUM COUNTS THE NUMBER OF ORDERED NODES PLUS 1. */
/*        ---------------------------------------------- */
    num = 1;

/*        ----------------------------- */
/*        ELIMINATE ALL ISOLATED NODES. */
/*        ----------------------------- */
    nextmd = dhead[1];
L100:
    if (nextmd <= 0) {
	goto L200;
    }
    mdnode = nextmd;
    nextmd = invp[mdnode];
    marker[mdnode] = *maxint;
    invp[mdnode] = -num;
    ++num;
    goto L100;

L200:
/*        ---------------------------------------- */
/*        SEARCH FOR NODE OF THE MINIMUM DEGREE. */
/*        MDEG IS THE CURRENT MINIMUM DEGREE; */
/*        TAG IS USED TO FACILITATE MARKING NODES. */
/*        ---------------------------------------- */
    if (num > *neqns) {
	goto L1000;
    }
    tag = 1;
    dhead[1] = 0;
    mdeg = 2;
L300:
    if (dhead[mdeg] > 0) {
	goto L400;
    }
    ++mdeg;
    goto L300;
L400:
/*            ------------------------------------------------- */
/*            USE VALUE OF DELTA TO SET UP MDLMT, WHICH GOVERNS */
/*            WHEN A DEGREE UPDATE IS TO BE PERFORMED. */
/*            ------------------------------------------------- */
    mdlmt = mdeg + *delta;
    ehead = 0;

L500:
    mdnode = dhead[mdeg];
    if (mdnode > 0) {
	goto L600;
    }
    ++mdeg;
    if (mdeg > mdlmt) {
	goto L900;
    }
    goto L500;
L600:
/*                ---------------------------------------- */
/*                REMOVE MDNODE FROM THE DEGREE STRUCTURE. */
/*                ---------------------------------------- */
    nextmd = invp[mdnode];
    dhead[mdeg] = nextmd;
    if (nextmd > 0) {
	perm[nextmd] = -mdeg;
    }
    invp[mdnode] = -num;
    *nofsub = *nofsub + mdeg + qsize[mdnode] - 2;
    if (num + qsize[mdnode] > *neqns) {
	goto L1000;
    }
/*                ---------------------------------------------- */
/*                ELIMINATE MDNODE AND PERFORM QUOTIENT GRAPH */
/*                TRANSFORMATION.  RESET TAG VALUE IF NECESSARY. */
/*                ---------------------------------------------- */
    ++tag;
    if (tag < *maxint) {
	goto L800;
    }
    tag = 1;
    i__1 = *neqns;
    for (i = 1; i <= i__1; ++i) {
	if (marker[i] < *maxint) {
	    marker[i] = 0;
	}
/* L700: */
    }
L800:
    mmdelm_dist(&mdnode, &xadj[1], &adjncy[1], &dhead[1], &invp[1], &perm[1],
		&qsize[1], &llist[1], &marker[1], maxint, &tag);
    num += qsize[mdnode];
    llist[mdnode] = ehead;
    ehead = mdnode;
    if (*delta >= 0) {
	goto L500;
    }
L900:
/*            ------------------------------------------- */
/*            UPDATE DEGREES OF THE NODES INVOLVED IN THE */
/*            MINIMUM DEGREE NODES ELIMINATION. */
/*            ------------------------------------------- */
    if (num > *neqns) {
	goto L1000;
    }
    mmdupd_dist(&ehead, neqns, &xadj[1], &adjncy[1], delta, &mdeg, &dhead[1],
		&invp[1], &perm[1], &qsize[1], &llist[1], &marker[1], maxint,
		&tag);
    goto L300;

L1000:
    mmdnum_dist(neqns, &perm[1], &invp[1], &qsize[1]);
    return 0;

} /* genmmd_dist_ */

/* *************************************************************** */
/* *************************************************************** */
/* ***     MMDINT ..... MULT MINIMUM DEGREE INITIALIZATION     *** */
/* *************************************************************** */
/* *************************************************************** */

/*     AUTHOR - JOSEPH W.H. LIU */
/*              DEPT OF COMPUTER SCIENCE, YORK UNIVERSITY. */

/*     PURPOSE - THIS ROUTINE PERFORMS INITIALIZATION FOR THE */
/*        MULTIPLE ELIMINATION VERSION OF THE MINIMUM DEGREE */
/*        ALGORITHM. */

/*     INPUT PARAMETERS - */
/*        NEQNS  - NUMBER OF EQUATIONS. */
/*        (XADJ,ADJNCY) - ADJACENCY STRUCTURE. */

/*     OUTPUT PARAMETERS - */
/*        (DHEAD,DFORW,DBAKW) - DEGREE DOUBLY LINKED STRUCTURE. */
/*        QSIZE  - SIZE OF SUPERNODE (INITIALIZED TO ONE). */
/*        LLIST  - LINKED LIST. */
/*        MARKER - MARKER VECTOR. */

/* *************************************************************** */

/* Subroutine */ int mmdint_dist(int *neqns, int *xadj, shortint *adjncy, 
	shortint *dhead, shortint *dforw, shortint *dbakw, shortint *qsize, 
	shortint *llist, shortint *marker)
{
    /* System generated locals */
    int i__1;

    /* Local variables */
    static int ndeg, node, fnode;


/* *************************************************************** */


/* *************************************************************** */

    /* Parameter adjustments */
    --marker;
    --llist;
    --qsize;
    --dbakw;
    --dforw;
    --dhead;
    --adjncy;
    --xadj;

    /* Function Body */
    i__1 = *neqns;
    for (node = 1; node <= i__1; ++node) {
	dhead[node] = 0;
	qsize[node] = 1;
	marker[node] = 0;
	llist[node] = 0;
/* L100: */
    }
/*        ------------------------------------------ */
/*        INITIALIZE THE DEGREE DOUBLY LINKED LISTS. */
/*        ------------------------------------------ */
    i__1 = *neqns;
    for (node = 1; node <= i__1; ++node) {
	ndeg = xadj[node + 1] - xadj[node] + 1;
	fnode = dhead[ndeg];
	dforw[node] = fnode;
	dhead[ndeg] = node;
	if (fnode > 0) {
	    dbakw[fnode] = node;
	}
	dbakw[node] = -ndeg;
/* L200: */
    }
    return 0;

} /* mmdint_dist */

/* *************************************************************** */
/* *************************************************************** */
/* **     MMDELM ..... MULTIPLE MINIMUM DEGREE ELIMINATION     *** */
/* *************************************************************** */
/* *************************************************************** */

/*     AUTHOR - JOSEPH W.H. LIU */
/*              DEPT OF COMPUTER SCIENCE, YORK UNIVERSITY. */

/*     PURPOSE - THIS ROUTINE ELIMINATES THE NODE MDNODE OF */
/*        MINIMUM DEGREE FROM THE ADJACENCY STRUCTURE, WHICH */
/*        IS STORED IN THE QUOTIENT GRAPH FORMAT.  IT ALSO */
/*        TRANSFORMS THE QUOTIENT GRAPH REPRESENTATION OF THE */
/*        ELIMINATION GRAPH. */

/*     INPUT PARAMETERS - */
/*        MDNODE - NODE OF MINIMUM DEGREE. */
/*        MAXINT - ESTIMATE OF MAXIMUM REPRESENTABLE (SHORT) */
/*                 INT. */
/*        TAG    - TAG VALUE. */

/*     UPDATED PARAMETERS - */
/*        (XADJ,ADJNCY) - UPDATED ADJACENCY STRUCTURE. */
/*        (DHEAD,DFORW,DBAKW) - DEGREE DOUBLY LINKED STRUCTURE. */
/*        QSIZE  - SIZE OF SUPERNODE. */
/*        MARKER - MARKER VECTOR. */
/*        LLIST  - TEMPORARY LINKED LIST OF ELIMINATED NABORS. */

/* *************************************************************** */

/* Subroutine */ int mmdelm_dist(int *mdnode, int *xadj, shortint *adjncy,
	 shortint *dhead, shortint *dforw, shortint *dbakw, shortint *qsize, 
 	 shortint *llist, shortint *marker, int *maxint, int *tag)
{
    /* System generated locals */
    int i__1, i__2;

    /* Local variables */
    static int node, link, rloc, rlmt, i, j, nabor, rnode, elmnt, xqnbr, 
	    istop, jstop, istrt, jstrt, nxnode, pvnode, nqnbrs, npv;


/* *************************************************************** */


/* *************************************************************** */

/*        ----------------------------------------------- */
/*        FIND REACHABLE SET AND PLACE IN DATA STRUCTURE. */
/*        ----------------------------------------------- */
    /* Parameter adjustments */
    --marker;
    --llist;
    --qsize;
    --dbakw;
    --dforw;
    --dhead;
    --adjncy;
    --xadj;

    /* Function Body */
    marker[*mdnode] = *tag;
    istrt = xadj[*mdnode];
    istop = xadj[*mdnode + 1] - 1;
/*        ------------------------------------------------------- */
/*        ELMNT POINTS TO THE BEGINNING OF THE LIST OF ELIMINATED */
/*        NABORS OF MDNODE, AND RLOC GIVES THE STORAGE LOCATION */
/*        FOR THE NEXT REACHABLE NODE. */
/*        ------------------------------------------------------- */
    elmnt = 0;
    rloc = istrt;
    rlmt = istop;
    i__1 = istop;
    for (i = istrt; i <= i__1; ++i) {
	nabor = adjncy[i];
	if (nabor == 0) {
	    goto L300;
	}
	if (marker[nabor] >= *tag) {
	    goto L200;
	}
	marker[nabor] = *tag;
	if (dforw[nabor] < 0) {
	    goto L100;
	}
	adjncy[rloc] = nabor;
	++rloc;
	goto L200;
L100:
	llist[nabor] = elmnt;
	elmnt = nabor;
L200:
	;
    }
L300:
/*            ----------------------------------------------------- */
/*            MERGE WITH REACHABLE NODES FROM GENERALIZED ELEMENTS. */
/*            ----------------------------------------------------- */
    if (elmnt <= 0) {
	goto L1000;
    }
    adjncy[rlmt] = -elmnt;
    link = elmnt;
L400:
    jstrt = xadj[link];
    jstop = xadj[link + 1] - 1;
    i__1 = jstop;
    for (j = jstrt; j <= i__1; ++j) {
	node = adjncy[j];
	link = -node;
	if (node < 0) {
	    goto L400;
	} else if (node == 0) {
	    goto L900;
	} else {
	    goto L500;
	}
L500:
	if (marker[node] >= *tag || dforw[node] < 0) {
	    goto L800;
	}
	marker[node] = *tag;
/*                            --------------------------------- */
/*                            USE STORAGE FROM ELIMINATED NODES */
/*                            IF NECESSARY. */
/*                            --------------------------------- */
L600:
	if (rloc < rlmt) {
	    goto L700;
	}
	link = -adjncy[rlmt];
	rloc = xadj[link];
	rlmt = xadj[link + 1] - 1;
	goto L600;
L700:
	adjncy[rloc] = node;
	++rloc;
L800:
	;
    }
L900:
    elmnt = llist[elmnt];
    goto L300;
L1000:
    if (rloc <= rlmt) {
	adjncy[rloc] = 0;
    }
/*        -------------------------------------------------------- */
/*        FOR EACH NODE IN THE REACHABLE SET, DO THE FOLLOWING ... */
/*        -------------------------------------------------------- */
    link = *mdnode;
L1100:
    istrt = xadj[link];
    istop = xadj[link + 1] - 1;
    i__1 = istop;
    for (i = istrt; i <= i__1; ++i) {
	rnode = adjncy[i];
	link = -rnode;
	if (rnode < 0) {
	    goto L1100;
	} else if (rnode == 0) {
	    goto L1800;
	} else {
	    goto L1200;
	}
L1200:
/*                -------------------------------------------- */
/*                IF RNODE IS IN THE DEGREE LIST STRUCTURE ... */
/*                -------------------------------------------- */
	pvnode = dbakw[rnode];
	if (pvnode == 0 || pvnode == -(*maxint)) {
	    goto L1300;
	}
/*                    ------------------------------------- */
/*                    THEN REMOVE RNODE FROM THE STRUCTURE. */
/*                    ------------------------------------- */
	nxnode = dforw[rnode];
	if (nxnode > 0) {
	    dbakw[nxnode] = pvnode;
	}
	if (pvnode > 0) {
	    dforw[pvnode] = nxnode;
	}
	npv = -pvnode;
	if (pvnode < 0) {
	    dhead[npv] = nxnode;
	}
L1300:
/*                ---------------------------------------- */
/*                PURGE INACTIVE QUOTIENT NABORS OF RNODE. */
/*                ---------------------------------------- */
	jstrt = xadj[rnode];
	jstop = xadj[rnode + 1] - 1;
	xqnbr = jstrt;
	i__2 = jstop;
	for (j = jstrt; j <= i__2; ++j) {
	    nabor = adjncy[j];
	    if (nabor == 0) {
		goto L1500;
	    }
	    if (marker[nabor] >= *tag) {
		goto L1400;
	    }
	    adjncy[xqnbr] = nabor;
	    ++xqnbr;
L1400:
	    ;
	}
L1500:
/*                ---------------------------------------- */
/*                IF NO ACTIVE NABOR AFTER THE PURGING ... */
/*                ---------------------------------------- */