symbfct.c

斯坦福大学Grant和Boyd教授等开发的凸优化matlab工具箱

/* symfct.f -- translated by f2c (version 19951025).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

typedef int integer;                  /* removed "long" */
#if !defined(max)
#define  max(A, B)   ((A) > (B) ? (A) : (B))
#endif

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

/*   Version:        0.3 */
/*   Last modified:  January 12, 1995 */
/*   Authors:        Esmond G. Ng and Barry W. Peyton */

/*   Mathematical Sciences Section, Oak Ridge National Laboratory */

/* *********************************************************************** */
/* *********************************************************************** */
/* **************    SFINIT  ..... SET UP FOR SYMB. FACT.     ************ */
/* *********************************************************************** */
/* *********************************************************************** */

/*   PURPOSE: */
/*       THIS SUBROUTINE COMPUTES THE STORAGE REQUIREMENTS AND SETS UP */
/*       PRELIMINARY DATA STRUCTURES FOR THE SYMBOLIC FACTORIZATION. */

/*   NOTE: */
/*       THIS VERSION PRODUCES THE MAXIMAL SUPERNODE PARTITION (I.E., */
/*       THE ONE WITH THE FEWEST POSSIBLE SUPERNODES). */

/*   INPUT PARAMETERS: */
/*       NEQNS       -   NUMBER OF EQUATIONS. */
/*       NNZA        -   LENGTH OF ADJACENCY STRUCTURE. */
/*       XADJ(*)     -   ARRAY OF LENGTH NEQNS+1, CONTAINING POINTERS */
/*                       TO THE ADJACENCY STRUCTURE. */
/*       ADJNCY(*)   -   ARRAY OF LENGTH XADJ(NEQNS+1)-1, CONTAINING */
/*                       THE ADJACENCY STRUCTURE. */
/*       IWSIZ       -   SIZE OF INTEGER WORKING STORAGE. */

/*   UPDATED PARAMETERS: */
/* (* JFS Sept 1, 1998: moved "PERM" and "INVP" to updated section *) */
/*       (PERM,INVP)     -   ON INPUT, THE GIVEN PERM AND INVERSE PERM */
/*                           VECTORS.  ON OUTPUT, THE NEW PERM AND */
/*                           INVERSE PERM VECTORS OF THE EQUIVALENT */
/*                           ORDERING. */

/*   OUTPUT PARAMETERS: */
/*       COLCNT(*)   -   ARRAY OF LENGTH NEQNS, CONTAINING THE NUMBER */
/*                       OF NONZEROS IN EACH COLUMN OF THE FACTOR, */
/*                       INCLUDING THE DIAGONAL ENTRY. */
/*       NNZL        -   NUMBER OF NONZEROS IN THE FACTOR, INCLUDING */
/*                       THE DIAGONAL ENTRIES. */
/*       NSUB        -   NUMBER OF SUBSCRIPTS. */
/*       NSUPER      -   NUMBER OF SUPERNODES (<= NEQNS). */
/*       SNODE(*)    -   ARRAY OF LENGTH NEQNS FOR RECORDING */
/*                       SUPERNODE MEMBERSHIP. */
/*       XSUPER(*)   -   ARRAY OF LENGTH NEQNS+1, CONTAINING THE */
/*                       SUPERNODE PARTITIONING. */
/*       IFLAG(*)    -   ERROR FLAG. */
/*                          0: SUCCESSFUL SF INITIALIZATION. */
/*                         -1: INSUFFICENT WORKING STORAGE */
/*                             [IWORK(*)]. */

/*   WORK PARAMETERS: */
/*       IWORK(*)    -   INTEGER WORK ARRAY OF LENGTH 7*NEQNS+3. */

/*   FIRST CREATED ON    NOVEMBER 14, 1994. */
/*   LAST UPDATED ON     January 12, 1995. */

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

/* Subroutine */ int sfinit_(neqns, nnza, xadj, adjncy, perm, invp, colcnt, 
	nnzl, nsub, nsuper, snode, xsuper, iwsiz, iwork, iflag)
integer *neqns, *nnza, *xadj, *adjncy, *perm, *invp, *colcnt, *nnzl, *nsub, *
	nsuper, *snode, *xsuper, *iwsiz, *iwork, *iflag;
{
    /* System generated locals */
    integer i__1;

    /* Local variables */
    extern /* Subroutine */ int fsup1_(), fsup2_();
    static integer i__;
    extern /* Subroutine */ int fcnthn_(), chordr_(), etordr_();


/*       ----------- */
/*       PARAMETERS. */
/*       ----------- */

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

/*       -------------------------------------------------------- */
/*       RETURN IF THERE IS INSUFFICIENT INTEGER WORKING STORAGE. */
/*       -------------------------------------------------------- */
    /* Parameter adjustments */
    --iwork;
    --xsuper;
    --snode;
    --colcnt;
    --invp;
    --perm;
    --xadj;
    --adjncy;

    /* Function Body */
    *iflag = 0;
    if (*iwsiz < *neqns * 7 + 3) {
	*iflag = -1;
	return 0;
    }
/* (* JFS: */
/* ------------------------------------------------------------ */
/* Handle the case of diagonal matrices separately, to avoid */
/* an unsolved BUG in FCNTHN. */
/* This patch is due to Jos F. Sturm, Sept 1, 1998. */
/* ------------------------------------------------------------ */
    if (xadj[*neqns + 1] - 1 == 0) {
	i__1 = *neqns;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    colcnt[i__] = 1;
	    snode[i__] = i__;
	    xsuper[i__] = i__;
/* L10: */
	}
	xsuper[*neqns + 1] = *neqns + 1;
	*nnzl = *neqns;
	*nsub = *neqns;
	*nsuper = *neqns;
	*iflag = 0;
	return 0;
    }
/* ------------------------------------------------------------ */
/* end of patch JFS *) */

/*       ------------------------------------------ */
/*       COMPUTE ELIMINATION TREE AND POSTORDERING. */
/*       ------------------------------------------ */
    etordr_(neqns, &xadj[1], &adjncy[1], &perm[1], &invp[1], &iwork[1], &
	    iwork[*neqns + 1], &iwork[(*neqns << 1) + 1], &iwork[*neqns * 3 + 
	    1]);

/*       --------------------------------------------- */
/*       COMPUTE ROW AND COLUMN FACTOR NONZERO COUNTS. */
/*       --------------------------------------------- */
    fcnthn_(neqns, nnza, &xadj[1], &adjncy[1], &perm[1], &invp[1], &iwork[1], 
	    &snode[1], &colcnt[1], nnzl, &iwork[*neqns + 1], &iwork[(*neqns <<
	     1) + 1], &xsuper[1], &iwork[*neqns * 3 + 1], &iwork[(*neqns << 2)
	     + 2], &iwork[*neqns * 5 + 3], &iwork[*neqns * 6 + 4]);

/*       --------------------------------------------------------- */
/*       REARRANGE CHILDREN SO THAT THE LAST CHILD HAS THE MAXIMUM */
/*       NUMBER OF NONZEROS IN ITS COLUMN OF L. */
/*       --------------------------------------------------------- */
    chordr_(neqns, &xadj[1], &adjncy[1], &perm[1], &invp[1], &colcnt[1], &
	    iwork[1], &iwork[*neqns + 1], &iwork[(*neqns << 1) + 1], &iwork[*
	    neqns * 3 + 1]);

/*       ---------------- */
/*       FIND SUPERNODES. */
/*       ---------------- */
    fsup1_(neqns, &iwork[1], &colcnt[1], nsub, nsuper, &snode[1]);
    fsup2_(neqns, nsuper, &iwork[1], &snode[1], &xsuper[1]);

    return 0;
} /* sfinit_ */

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

/*   Version:        0.3 */
/*   Last modified:  December 27, 1994 */
/*   Authors:        Joseph W.H. Liu */

/*   Mathematical Sciences Section, Oak Ridge National Laboratory */

/* *********************************************************************** */
/* *********************************************************************** */
/* **********     ETORDR ..... ELIMINATION TREE REORDERING     *********** */
/* *********************************************************************** */
/* *********************************************************************** */

/*   WRITTEN BY JOSEPH LIU (JUL 17, 1985) */

/*   PURPOSE: */
/*       TO DETERMINE AN EQUIVALENT REORDERING BASED ON THE STRUCTURE OF */
/*       THE ELIMINATION TREE.  A POSTORDERING OF THE GIVEN ELIMINATION */
/*       TREE IS RETURNED. */

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

/*   UPDATED PARAMETERS: */
/*       (PERM,INVP)     -   ON INPUT, THE GIVEN PERM AND INVERSE PERM */
/*                           VECTORS.  ON OUTPUT, THE NEW PERM AND */
/*                           INVERSE PERM VECTORS OF THE EQUIVALENT */
/*                           ORDERING. */

/*   OUTPUT PARAMETERS: */
/*       PARENT          -   THE PARENT VECTOR OF THE ELIMINATION TREE */
/*                           ASSOCIATED WITH THE NEW ORDERING. */

/*   WORKING PARAMETERS: */
/*       FSON            -   THE FIRST SON VECTOR. */
/*       BROTHR          -   THE BROTHER VECTOR. */
/*       INVPOS          -   THE INVERSE PERM VECTOR FOR THE */
/*                           POSTORDERING. */

/*   PROGRAM SUBROUTINES: */
/*       BETREE, ETPOST, ETREE , INVINV. */

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

/* Subroutine */ int etordr_(neqns, xadj, adjncy, perm, invp, parent, fson, 
	brothr, invpos)
integer *neqns, *xadj, *adjncy, *perm, *invp, *parent, *fson, *brothr, *
	invpos;
{
    extern /* Subroutine */ int etree_(), betree_(), invinv_(), etpost_();


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



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

/*       ----------------------------- */
/*       COMPUTE THE ELIMINATION TREE. */
/*       ----------------------------- */
    /* Parameter adjustments */
    --invpos;
    --brothr;
    --fson;
    --parent;
    --invp;
    --perm;
    --adjncy;
    --xadj;

    /* Function Body */
    etree_(neqns, &xadj[1], &adjncy[1], &perm[1], &invp[1], &parent[1], &
	    invpos[1]);

/*       -------------------------------------------------------- */
/*       COMPUTE A BINARY REPRESENTATION OF THE ELIMINATION TREE. */
/*       -------------------------------------------------------- */
    betree_(neqns, &parent[1], &fson[1], &brothr[1]);

/*       ------------------------------- */
/*       POSTORDER THE ELIMINATION TREE. */
/*       ------------------------------- */
    etpost_(neqns, &fson[1], &brothr[1], &invpos[1], &parent[1], &perm[1]);

/*       -------------------------------------------------------- */
/*       COMPOSE THE ORIGINAL ORDERING WITH THE NEW POSTORDERING. */
/*       -------------------------------------------------------- */
    invinv_(neqns, &invp[1], &invpos[1], &perm[1]);

    return 0;
} /* etordr_ */

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

/*   Version:        0.3 */
/*   Last modified:  December 27, 1994 */
/*   Authors:        Joseph W.H. Liu */

/*   Mathematical Sciences Section, Oak Ridge National Laboratory */

/* *********************************************************************** */
/* *********************************************************************** */
/* ****************     ETREE ..... ELIMINATION TREE     ***************** */
/* *********************************************************************** */
/* *********************************************************************** */

/*   WRITTEN BY JOSEPH LIU (JUL 17, 1985) */

/*   PURPOSE: */
/*       TO DETERMINE THE ELIMINATION TREE FROM A GIVEN ORDERING AND */
/*       THE ADJACENCY STRUCTURE.  THE PARENT VECTOR IS RETURNED. */

/*   INPUT PARAMETERS: */
/*       NEQNS           -   NUMBER OF EQUATIONS. */
/*       (XADJ,ADJNCY)   -   THE ADJACENCY STRUCTURE. */
/*       (PERM,INVP)     -   PERMUTATION AND INVERSE PERMUTATION VECTORS */

/*   OUTPUT PARAMETERS: */
/*       PARENT          -   THE PARENT VECTOR OF THE ELIMINATION TREE. */

/*   WORKING PARAMETERS: */
/*       ANCSTR          -   THE ANCESTOR VECTOR. */

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

/* Subroutine */ int etree_(neqns, xadj, adjncy, perm, invp, parent, ancstr)
integer *neqns, *xadj, *adjncy, *perm, *invp, *parent, *ancstr;
{
    /* System generated locals */
    integer i__1, i__2;

    /* Local variables */
    static integer node, next, i__, j, jstop, jstrt, nbr;


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



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


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

    /* Parameter adjustments */
    --ancstr;
    --parent;
    --invp;
    --perm;
    --adjncy;
    --xadj;

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

    i__1 = *neqns;
    for (i__ = 1; i__ <= i__1; ++i__) {
	parent[i__] = 0;
	ancstr[i__] = 0;
	node = perm[i__];

	jstrt = xadj[node];
	jstop = xadj[node + 1] - 1;
	if (jstrt <= jstop) {
	    i__2 = jstop;
	    for (j = jstrt; j <= i__2; ++j) {
		nbr = adjncy[j];
		nbr = invp[nbr];
		if (nbr < i__) {
/*                       --------------------------------
----------- */
/*                       FOR EACH NBR, FIND THE ROOT OF IT
S CURRENT */
/*                       ELIMINATION TREE.  PERFORM PATH C
OMPRESSION */
/*                       AS THE SUBTREE IS TRAVERSED. */
/*                       --------------------------------
----------- */
L100:
		    if (ancstr[nbr] == i__) {
			goto L300;
		    }
		    if (ancstr[nbr] > 0) {
			next = ancstr[nbr];
			ancstr[nbr] = i__;
			nbr = next;
			goto L100;
		    }
/*                       --------------------------------
------------ */
/*                       NOW, NBR IS THE ROOT OF THE SUBTR
EE.  MAKE I */
/*                       THE PARENT NODE OF THIS ROOT. */
/*                       --------------------------------
------------ */
		    parent[nbr] = i__;
		    ancstr[nbr] = i__;
		}
L300:
		;
	    }
	}
/* L400: */
    }