/*-*- mode: C; kept-old-versions: 12;  kept-new-versions: 20; -*-
 *
 * srqfnc.f -- translated by f2c (version 20031025) and by
 * $Id: f2c-clean,v 1.10 2002/03/28 16:37:27 maechler Exp $
 * plus extended manual code cleaning by Martin Maechler, ETH Zurich
 */

#include <Rmath.h>

/* Table of constant values */

static int c__1 = 1;
static double c_9995 = .9995;


/* BLAS : */
#include <R_ext/BLAS.h>
/* double ddot_(int *, double *, int *, double *, int *); */
/* int   daxpy_(int *, double*, double *, int *, double *, int *); */

#define DDOT(_n_, _X_, _Y_)  F77_CALL(ddot)(_n_, _X_, &c__1, _Y_, &c__1)

/* SparseM -- ./sparskit2.f : */
#include "sparseM.h"

/* Cholesky related : */
#include "cholesky.h"

/* advance declaration : */
static void
slpfnc(int *n1, int *m,
       int *nnza1, double *a1, int *ja1, int *ia1,
       double *ao1, int *jao1, int *iao1,
       int *n2,
       int *nnza2, double *a2, int *ja2, int *ia2,
       double *ao2, int *jao2, int *iao2,
       int *nnzdmax, double *d, int *jd, int *id,
       double *dsub, int *jdsub, int *nsubmax, int *lindx,
       int *xlindx, int *nnzlmax, double *lnz, int *xlnz,
       int *invp, int *perm, int *iwmax, int *iwork, int *colcnt,
       int *snode, int *xsuper, int *split,
       int *tmpmax, double *tmpvec, double *rhs, double *newrhs,
       int *cachsz, int *level, double *x1, double *x2,
       double *s, double *u, double *c1, double *c2,
       double *y, double *b, double *r2, double *z1,
       double *z2, double *w, double *q1, double *q2,
       int *nnzemax, double *e, int *je, int *ie,
       int *nnzgmax, double *g, int *jg, int *ig,
       int *nnzhmax, double *h, int *jh, int *ih,
       double *dy, double *dx1, double *dx2, double *ds,
       double *dz1, double * dz2, double *dw,
       double *dxdz1, double *dxdz2, double *dsdw,
       double *xi1, double *xi2, int *maxn1n2,
       double *ww1, double *ww2, double *ww3,
       double *small, int *ierr, int *maxit, double *timewd);

static void
boundc(double *x1, double *dx1, double *x2,
       double *dx2, double *s, double *ds, double *z1,
       double *dz1, double *z2, double *dz2, double *w,
       double *dw, int *n1, int *n2, double *beta,
       double *deltap, double *deltad);

/* This is called from R : */
int
F77_SUB(srqfnc)(int *n1, int *m, int *nnza1,
	double *a1, int *ja1, int *ia1,
	double *ao1, int *jao1, int *iao1,
	int *n2, int *nnza2,
	double *a2, int *ja2, int *ia2,
	double *ao2, int *jao2, int *iao2,
	int *nnzdmax, double *d, int *jd, int *id,
	double *dsub, int *jdsub, int *nnzemax, double *e,
	int *je, int *ie, int *nnzgmax,
	double *g, int *jg, int *ig, int *nnzhmax,
	double *h, int *jh, int *ih,
	int *nsubmax, int *lindx, int *xlindx,
	int *nnzlmax, double *lnz, int *xlnz, int *iw,
	int *iwmax, int *iwork, int *xsuper, int *tmpmax,
	double *tmpvec, int *maxn1n2,
	double *ww1, double *wwm, double *wwn1, double *wwn2,
	int *cachsz, int *level, double *x1, double *x2,
	double *s, double *u,
	double *c1, double *c2, double *y, double *small,
	int *ierr, int *maxit, double *timewd)
{
    /* System generated locals */
    int iw_dim1, wwm_dim1, wwn1_dim1, wwn2_dim1;

    /* Parameter adjustments */
    wwm_dim1  = *m;	wwm -= wwm_dim1;
    wwn1_dim1 = *n1;	wwn1 -= wwn1_dim1;
    wwn2_dim1 = *n2;	wwn2 -= wwn2_dim1;
    iw_dim1   = *m;	iw -= iw_dim1;

    /* Function Body */
    slpfnc(n1, m, nnza1, a1, ja1, ia1, ao1, jao1, iao1,
	    n2, nnza2, a2, ja2, ia2, ao2, jao2, iao2, nnzdmax,
	    d, jd, id, dsub, jdsub,
	    nsubmax, lindx, xlindx, nnzlmax,
	    lnz, xlnz, &iw[iw_dim1], &iw[(iw_dim1 << 1)],
	    iwmax, iwork, &iw[iw_dim1 * 3], &iw[(iw_dim1 << 2)],
	    xsuper, &iw[iw_dim1 * 5], tmpmax, tmpvec,
	    &wwm[(wwm_dim1 << 1)], &wwm[wwm_dim1 * 3], cachsz, level,
	    x1, x2, s, u,
	    c1, c2, y, &wwm[wwm_dim1], &wwn2[wwn2_dim1],
	    &wwn1[wwn1_dim1], &wwn2[(wwn2_dim1 << 1)],
	    &wwn1[(wwn1_dim1 << 1)], &wwn1[wwn1_dim1 * 3],
	    &wwn2[wwn2_dim1 * 3], nnzemax, e, je, ie,
	    nnzgmax, g, jg, ig, nnzhmax, h, jh, ih,
	    &wwm[(wwm_dim1 << 2)],
	    &wwn1[(wwn1_dim1 << 2)], &wwn2[(wwn2_dim1 << 2)],
	    &wwn1[wwn1_dim1 * 5], &wwn1[wwn1_dim1 * 6],
	    &wwn2[wwn2_dim1 * 5],
	    &wwn1[wwn1_dim1 * 7], &wwn1[(wwn1_dim1 << 3)],
	    &wwn2[wwn2_dim1 * 6], &wwn1[wwn1_dim1 * 9],
	    &wwn1[wwn1_dim1 * 10], &wwn2[wwn2_dim1 * 7], maxn1n2, ww1,
	    &wwm[wwm_dim1 * 5], &wwm[wwm_dim1 * 6], small, ierr, maxit, timewd);
    return 0;
} /* srqfnc_ */

static void
slpfnc(int *n1, int *m,
       int *nnza1, double *a1, int *ja1, int *ia1,
       double *ao1, int *jao1, int *iao1,
       int *n2,
       int *nnza2, double *a2, int *ja2, int *ia2,
       double *ao2, int *jao2, int *iao2,
       int *nnzdmax, double *d, int *jd, int *id,
       double *dsub, int *jdsub, int *nsubmax, int *lindx,
       int *xlindx, int *nnzlmax, double *lnz, int *xlnz,
       int *invp, int *perm, int *iwmax, int *iwork, int *colcnt,
       int *snode, int *xsuper, int *split,
       int *tmpmax, double *tmpvec, double *rhs, double *newrhs,
       int *cachsz, int *level, double *x1, double *x2,
       double *s, double *u, double *c1, double *c2,
       double *y, double *b, double *r2, double *z1,
       double *z2, double *w, double *q1, double *q2,
       int *nnzemax, double *e, int *je, int *ie,
       int *nnzgmax, double *g, int *jg, int *ig,
       int *nnzhmax, double *h, int *jh, int *ih,
       double *dy, double *dx1, double *dx2, double *ds,
       double *dz1, double * dz2, double *dw,
       double *dxdz1, double *dxdz2, double *dsdw,
       double *xi1, double *xi2, int *maxn1n2,
       double *ww1, double *ww2, double *ww3,
       double *small, int *ierr, int *maxit, double *timewd)
{

/* Sparse implentation of LMS's interior point method via
    Ng-Peyton's sparse Cholesky factorization for sparse
    symmetric positive definite

 INPUT:
     n1 -- the number of row in the coefficient matrix A1'
     m -- the number of column in the coefficient matrix A1'
     nnza1 -- the number of non-zero elements in A'
     a1 -- an nnza1-vector of non-zero values of the design
          matrix (A1') stored in csr format
     ja1 -- an nnza1-vector of indices of the non-zero elements of
           the coefficient matrix
     ia1 -- an (n1+1)-vector of pointers to the begining of each
           row in a1 and ja1
     ao1 -- an nnza1-vector of work space for the transpose of
           the design matrix stored in csr format or the
           design matrix stored in csc format
     jao1 -- an nnza1-vector of work space for the indices of the
            transpose of the design matrix
     iao1 -- an (n1+1)-vector of pointers to the begining of each
            column in ao1 and jao1
     n2 -- the number of row in the constraint matrix A2'
     nnza2 -- the number of non-zero elements in A2'
     a2 -- an nnza2-vector of non-zero values of the contraint
          matrix (A2') stored in csr format
     ja2 -- an nnza2-vector of indices of the non-zero elements of
           the constraint matrix
     ia2 -- an (n2+1)-vector of pointers to the begining of each
           row in a2 and ja2
     ao2 -- an nnza2-vector of work space for the transpose of
           the constraint matrix stored in csr format or the
           constraint matrix stored in csc format
     jao2 -- an nnza2-vector of work space for the indices of the
            transpose of the constraint matrix
     iao2 -- an (n2+1)-vector of pointers to the begining of each
            column in ao2 and jao2
     nnzdmax -- upper bound of the non-zero elements in A1A1'
     d -- an nnzdmax-vector of non-zero values used to store
          the transpose of the design matrix multiplied by the design
          matrix (A1A1') stored in csr format;
          also used to store A1Q1^(-1) and A2Q2^(-1) later
     jd -- an nnzdmax-vector of indices in d
     id -- an (m+1)-vector of pointers to the begining of each
           row in d and jd
     dsub -- the values of d excluding the diagonal elements
     jdsub -- the indices to dsub
     nsubmax -- upper bound of the dimension of lindx
     lindx -- an nsub-vector of integer which contains, in
           column major order, the row subscripts of the nonzero
           entries in L in a compressed storage format
     xlindx -- an (m+1)-vector of int of pointers for lindx
     nnzlmax -- the upper bound of the non-zero entries in
                L stored in lnz, including the diagonal entries
     lnz -- First contains the non-zero entries of d; later
            contains the entries of the Cholesky factor
     xlnz -- column pointer for L stored in lnz
     invp -- an n1-vector of int of inverse permutation
             vector
     perm -- an n1-vector of int of permutation vector
     iw -- int work array of length m
     iwmax -- upper bound of the general purpose int
              working storage iwork; set at 7*m+3
     iwork -- an iwsiz-vector of int as work space
     colcnt -- array of length m, containing the number of
               non-zeros in each column of the factor, including
               the diagonal entries
     snode -- array of length m for recording supernode
              membership
     xsuper -- array of length m+1 containing the supernode
               partitioning
     split -- an m-vector with splitting of supernodes so that
              they fit into cache
     tmpmax -- upper bound of the dimension of tmpvec
     tmpvec -- a tmpmax-vector of temporary vector
     rhs -- m-vector to store the rhs
     newrhs -- extra work vector for right-hand side and
               solution
     cachsz -- size of the cache (in kilobytes) on the target
               machine
     level -- level of loop unrolling while performing numerical
              factorization
     x1 -- an n1-vector, the initial feasible solution for the primal
           solution that corresponds to the design matrix A1'
     x2 -- an n2-vector, the initial feasible solution for the primal
           solution that corresponds to the constraint matrix A2'
     s -- an n1-vector
     u -- an n1-vector of the upper bound for x1
     c1 -- an n1-vector in the primal; negative response in the
           regression quantile setting
     c2 -- an n2-vector, the negative rhs of the inequality constraint
     y -- an m-vector, the initial dual solution
     b -- an n1-vector, usualy the rhs of the equality constraint
          X'a = (1-tau)X'e in the rq setting
     r2 -- an n2-vector of residuals
     z1 -- an n1-vector of the dual slack variable
     z2 -- an n2-vector
     w -- an n-vector
     q1 -- an n1-vector of work array containing the diagonal
          elements of the Q1^(-1) matrix
     q2 -- an n2-vector of work array containing the diagonal
          elements of the Q2^(-1) matrix
     e -- an nnzdmax-vector containing the non-zero entries of
          A1Q1^(-1)A1' stored in csr format
     je -- an nnzdmax-vector of indices for e
     ie -- an (m+1)-vector of pointers to the begining of each
           row in e and je
     nnzgmax -- upper bound of the non-zero elements in g,jg
     g -- an nnzgmax-vector containing the non-zero entries of
          A2Q2^(-1)A2' stored in csr format
     jg -- an nnzgmax-vector of indices for g
     ig -- an (m+1)-vector of pointers to the begining of each
           row in g and jg
     nnzhmax -- upper bound of the non-zero elements in h,jh
     h -- an nnzhmax-vector containing the non-zero entries of
          AQ^(-1)A' stored in csr format
     jh -- an nnzhmax-vector of indices for h
     ih -- an (m+1)-vector of pointers to the begining of each
           row in h and jh
     dy -- an m-vector of work array
     dx1 -- an n1-vector of work array
     dx2 -- an n2-vector of work array
     ds -- an n1-vector of work array
     dz1 -- an n1-vector of work array
     dz2 -- an n2-vector of work array
     dw -- an n1-vector of work array
     dxdz1 -- an n1-vector of work array
     dxdz2 -- an n2-vector of work array
     dsdw -- an n1-vector of work arry
     xi1 -- an n1-vector of work array
     xi2 -- an n2-vector of work array
     xinv1 -- an n1-vector of work array
     xinv2 -- an n2-vector of work array
     sinv -- work array
     maxn1n2 -- max(n1,n2)
     ww1 -- an maxn1n2-vector of work array
     ww2 -- an m-vector of work array
     ww3 -- an m-vector of work array
     small -- convergence tolerance for interior algorithm
     ierr -- error flag :
       1 -- insufficient storage (work space) when calling extract;
       3 -- insufficient storage in iwork when calling ordmmd;
       4 -- insufficient storage in iwork when calling sfinit;
       5 -- nnzl > nnzlmax when calling sfinit
       6 -- nsub > nsubmax when calling sfinit
       7 -- insufficient work space in iwork when calling symfct
       8 -- inconsistancy in input when calling symfct
       9 -- tmpsiz > tmpmax when calling symfct; increase tmpmax
       10 -- nonpositive diagonal encountered    when calling blkfct
       11 -- insufficient work storage in tmpvec when calling blkfct
       12 -- insufficient work storage in iwork  when calling blkfct
       13 -- nnzd > nnzdmax in e,je when calling amub
       14 -- nnzd > nnzdmax in g,jg when calling amub
       15 -- nnzd > nnzdmax in h,jh when calling aplb
       17 -- tiny diagonals replaced with Inf when calling blkfct()
       		{new error code, was confounded with '10'}

     maxit -- maximal number of iterations; on return: the number of iterations

     timewd -- vector of length 7: [7]: amount of time for this subroutine
               [1:6] time info for the phases of cholfct() only.

 OUTPUT:
     y -- an m-vector of primal solution
*/

    /* Local variables */
    int i, it, nnzd, nsuper, nnzdsub;
    double mu, gap, deltad, deltap;
    double timbeg, timend;

/* Parameter adjustments */
    --dsdw;
    --dxdz1;
    --dxdz2;
    --dz1; --dz2;
    --dx1; --dx2;
    --ds;
    --dw;
    --dy;
    --c1;  --c2;
    --q1;  --q2;
    --x1;  --x2;
    --z1;  --z2;
    --r2;
    --xi1; --xi2;
    --b;
    --s;
    --u;
    --w;
    --y;

    --ww1;
    --ww2;
    --ww3;

    --newrhs;
    --rhs;

    --perm;
    --invp;

    --dsub;
    --timewd;

    /* Function Body */
    for (i = 1; i <= 7; ++i)
	timewd[i] = 0.;

/*  Compute the initial gap */

    gap =
	DDOT(n1, &z1[1], &x1[1]) +
	DDOT(n2, &z2[1], &x2[1]) +
	DDOT(n1, &w[1],  &s[1]);

/*  Start iteration */
    it = 0;
    while(gap >= *small && it <= *maxit) {

	++it;

/*  Create the diagonal matrix Q1^(-1) stored in q1 as an n1-vector,
    the diagonal matrix Q2^(-1) stored in q2 as an n2-vector,
    and store the residuals in r1 in ds, and r3 in dy temporarily,
    and r2 in r2 permanently
*/

	/* amux: obtain A1x1 and store the value in ww2 */
	F77_CALL(amux)(m, &x1[1], &ww2[1], ao1, jao1, iao1);

	/* amux: obtain A2x2 and store the value in ww3 */
	F77_CALL(amux)(m, &x2[1], &ww3[1], ao2, jao2, iao2);

	/*    obtain  A2'y and store it temporarily in r2 */
	F77_CALL(amux)(n2, &y[1], &r2[1], a2, ja2, ia2);

	for (i = 1; i <= *n1; ++i) {
	    q1[i] = 1. / (z1[i] / x1[i] + w[i] / s[i]);
	    ds[i] = z1[i] - w[i];
	}
	for (i = 1; i <= *n2; ++i) {
	    q2[i] = x2[i] / z2[i];
	    r2[i] = c2[i] - r2[i];
	}
	for (i = 1; i <= *m; ++i) {
	    dy[i] = b[i] - ww2[i] - ww3[i];
	}

/*  Obtain AQA = A1Q1^(-1)A1' + A2Q2^(-1)A2' in 5 steps :

 *  Step1: Obtain A1Q1^(-1) and store the values in d,jd,id in csr format
 *         Also compute A1Q1^(-1)r1 and store the values in ww2 to be used
 *         to generate r3;
 *  Step2: Compute A1Q1^(-1)A1' and store the values in e,je,ie
 *  Step3: Obtain A2Q2^(-1) and store the values in d,jd,id in csr format
 *         Also compute A2Q2^(-1)r2 and store the values in in ww3 to
 *         be used to generate r3;
 *  Step4: Compute A2Q2^(-1)A2' and store the value in g,jg,ig
 *  Step5: Compute AQA and store the values in h,jh,ih
 */

/*   Step 1 */

	F77_CALL(amudia)(m, &c__1, ao1, jao1, iao1, &q1[1], d, jd, id);
	F77_CALL(amux)(m, &ds[1], &ww2[1], d, jd, id);

/*   Step 2 */

	F77_CALL(amub)(m, m, &c__1, d, jd, id, a1, ja1, ia1,
		       e, je, ie, nnzemax, iwork, ierr);
	if (*ierr) { *ierr = 13; goto L100; }

/*   Step 3 */

	F77_CALL(amudia)(m, &c__1, ao2, jao2, iao2, &q2[1], d, jd, id);
	F77_CALL(amux)(m, &r2[1], &ww3[1], d, jd, id);

/*   Step 4 */

	F77_CALL(amub)(m, m, &c__1, d, jd, id, a2, ja2, ia2,
		       g, jg, ig, nnzgmax, iwork, ierr);
	if (*ierr) { *ierr = 14; goto L100; }

/*   Step 5 */

	F77_CALL(aplb)(m, m, &c__1, e, je, ie, g, jg, ig,
		       h, jh, ih, nnzhmax, iwork, ierr);
	if (*ierr) { *ierr = 15; goto L100; }

	/* Generate rhs = r3 + A1Q1^(-1) r1 + A2Q2^(-1) r2 : */

	for (i = 1; i <= *m; ++i) {
	    rhs[i] = dy[i] + ww2[i] + ww3[i];
	}

/*  Extract the non-diagonal structure of h,jh,ih and store in dsub,jdsub */

	nnzd = ih[*m] - 1;
	nnzdsub = nnzd - *m;
	i = *nnzhmax + 1;
	F77_CALL(extract)(h, jh, ih, &dsub[1], jdsub, m, nnzhmax, &i, ierr);

	if (*ierr) { *ierr = 1; goto L100; }

/*  Compute dy = (AQ^(-1)A')^(-1)rhs; result returned via dy

Call chlfct to perform Cholesky's decomposition of h,jh,ih */

	F77_CALL(chlfct)(m, xlindx, lindx, &invp[1], &perm[1], iwork, &nnzdsub,
			 jdsub, colcnt, &nsuper, snode, xsuper, nnzlmax,
			 nsubmax, xlnz, lnz, ih, jh, h, cachsz,
			 tmpmax, level, tmpvec, split, ierr, &it, &timewd[1]);

	if (*ierr) { goto L100; }

/* Call blkslv: Numerical solution for the new rhs stored in rhs */

	for (i = 1; i <= *m; ++i) {
	    newrhs[i] = rhs[perm[i]];
	}
	timbeg = F77_CALL(gtimer)();
	F77_CALL(blkslv)(&nsuper, xsuper, xlindx, lindx, xlnz, lnz, &newrhs[1]);
	timend = F77_CALL(gtimer)();
	timewd[7] = timewd[7] + timend - timbeg;
	for (i = 1; i <= *m; ++i) {
	    dy[i] = newrhs[invp[i]];
	}

/*  Compute dx1 = Q1^(-1)(A1'dy - r1), ds = -dx1, dz1, dz2  and dw */

	F77_CALL(amux)(n1, &dy[1], &dx1[1], a1, ja1, ia1);
	F77_CALL(amux)(n2, &dy[1], &dx2[1], a2, ja2, ia2);
	for (i = 1; i <= *n1; ++i) {
	    dx1[i] = q1[i] * (dx1[i] - ds[i]);
	    ds[i] = -dx1[i];
	    dz1[i] = -z1[i] * (dx1[i] / x1[i] + 1.);
	    dw[i] = -w[i] * (ds[i] / s[i] + 1.);
	}
	for (i = 1; i <= *n2; ++i) {
	    dx2[i] = q2[i] * (dx2[i] - r2[i]);
	    dz2[i] = -z2[i] * (dx2[i] / x2[i] + 1.);
	}

/* Compute the maximum allowable step lengths */

	boundc(&x1[1], &dx1[1], &x2[1], &dx2[1],
	       &s[1],  &ds[1],  &z1[1], &dz1[1],
	       &z2[1], &dz2[1], &w[1], &dw[1],
	       n1, n2, &c_9995, &deltap, &deltad);

	if (deltap * deltad < 1.) {

	    /* Update mu */
	    double g1;

	    mu = DDOT(n1, &z1[1], &x1[1]) +
		DDOT(n2, &z2[1], &x2[1]) +
		DDOT(n1, &w[1], &s[1]);

	    g1 = mu +
		deltap *          DDOT(n1, &z1[1], &dx1[1]) +
		deltad *          DDOT(n1, &dz1[1], &x1[1]) +
		deltad * deltap * DDOT(n1, &dz1[1], &dx1[1]) +
		deltap *          DDOT(n2, &z2[1], &dx2[1]) +
		deltad *          DDOT(n2, &dz2[1], &x2[1]) +
		deltad * deltap * DDOT(n2, &dz2[1], &dx2[1]) +
		deltap *          DDOT(n1, &w[1], &ds[1]) +
		deltad *          DDOT(n1, &dw[1], &s[1]) +
		deltad * deltap * DDOT(n1, &dw[1], &ds[1]);

	    g1 /= mu;
	    mu *=  g1 * g1 * g1 / (2. * *n1 + *n2);

/* Compute dx1dz1, dx2dz2 and dsdw */

	    for (i = 1; i <= *n1; ++i) {
		dxdz1[i] = dx1[i] * dz1[i];
		dsdw[i] = ds[i] * dw[i];
		xi1[i] = dxdz1[i] / x1[i] - dsdw[i] / s[i] -
		     mu * (    1. / x1[i] -      1. / s[i]);
		ww1[i] = q1[i] * xi1[i];
	    }

/* Compute A1Q1^(-1)(X1^(-1)*dx1dz1 - S^(-1)*dsdw - mu(X1^(-1) - S^(-1))) and
   store it in ww2 temporarily */

	    F77_CALL(amux)(m, &ww1[1], &ww2[1], ao1, jao1, iao1);
	    for (i = 1; i <= *n2; ++i) {
		dxdz2[i] = dx2[i] * dz2[i];
		xi2[i] = (dxdz2[i] - mu) / x2[i];
		ww1[i] = q2[i] * xi2[i];
	    }

/* Compute A2Q2^(-1)(X2^(-1)*dx2dz2 - mu X2^(-1)) and store it in ww3 temporarily */

	    F77_CALL(amux)(m, &ww1[1], &ww3[1], ao2, jao2, iao2);
	    for (i = 1; i <= *m; ++i) {
		rhs[i] +=  ww2[i] + ww3[i];
	    }


/* Compute (AQ^(-1)A')^(-1)rhs and return the result in dy

Call blkslv: Numerical solution for the new rhs stored in rhs */

	    for (i = 1; i <= *m; ++i) {
		newrhs[i] = rhs[perm[i]];
	    }
	    timbeg = F77_CALL(gtimer)();
	    F77_CALL(blkslv)(&nsuper, xsuper, xlindx, lindx, xlnz, lnz, &newrhs[1]);
	    timend = F77_CALL(gtimer)();
	    timewd[7] = timewd[7] + timend - timbeg;
	    for (i = 1; i <= *m; ++i) {
		dy[i] = newrhs[invp[i]];
	    }

/*  Compute dx1=Q1^(-1)(A1'dy-X1^(-1)*dx1dz1-S^(-1)*dsdw
    -mu*(X1^(-1)-S^(-1))-r1), ds = -dx1, dz1, dz2  and dw */

	    F77_CALL(amux)(n1, &dy[1], &dx1[1], a1, ja1, ia1);
	    F77_CALL(amux)(n2, &dy[1], &dx2[1], a2, ja2, ia2);
	    for (i = 1; i <= *n1; ++i) {
		dx1[i] = q1[i] * (dx1[i] - xi1[i] - z1[i] + w[i]);
		ds[i] = -dx1[i];
		dz1[i] = -z1[i] + (mu - z1[i] * dx1[i] - dxdz1[i]) / x1[i];
		dw[i] = -w[i] + (mu - w[i] * ds[i] - dsdw[i]) / s[i];
	    }
	    for (i = 1; i <= *n2; ++i) {
		dx2[i] = q2[i] * (dx2[i] - xi2[i] - r2[i]);
		dz2[i] = -z2[i] + (mu - z2[i] * dx2[i] - dxdz2[i]) / x2[i];
	    }

/* Compute the maximum allowable step lengths */

	    boundc(&x1[1], &dx1[1], &x2[1], &dx2[1], &s[1], &ds[1],
		   &z1[1], &dz1[1], &z2[1], &dz2[1], &w[1], &dw[1],
		   n1, n2, &c_9995, &deltap, &deltad);
	}

/* Take the step */

	F77_CALL(daxpy)(n1, &deltap, &dx1[1], &c__1, &x1[1], &c__1);
	F77_CALL(daxpy)(n2, &deltap, &dx2[1], &c__1, &x2[1], &c__1);
	F77_CALL(daxpy)(n1, &deltap, &ds[1],  &c__1, &s[1],  &c__1);
	F77_CALL(daxpy)(n1, &deltad, &dw[1],  &c__1, &w[1],  &c__1);
	F77_CALL(daxpy)(n1, &deltad, &dz1[1], &c__1, &z1[1], &c__1);
	F77_CALL(daxpy)(n2, &deltad, &dz2[1], &c__1, &z2[1], &c__1);
	F77_CALL(daxpy)(m,  &deltad, &dy[1],  &c__1, &y[1],  &c__1);
	gap =
	    DDOT(n1, &z1[1], &x1[1]) +
	    DDOT(n2, &z2[1], &x2[1]) +
	    DDOT(n1, &w[1],  &s[1] );

    } /* end {while} iterations */

L100:
    *maxit = it;
    return;
} /* slpfnc */

static void
boundc(double *x1, double *dx1, double *x2,
       double *dx2, double *s, double *ds, double *z1,
       double *dz1, double *z2, double *dz2, double *w,
       double *dw, int *n1, int *n2, double *beta,
       double *deltap, double *deltad)
{
    int i;

    *deltap = R_PosInf;
    *deltad = R_PosInf;

    for (i = 0; i < *n1; ++i) {
	if(dx1[i] < 0) *deltap = fmin2(*deltap, -x1[i] / dx1[i]);
	if(ds[i]  < 0) *deltap = fmin2(*deltap, -s[i] / ds[i]);
	if(dz1[i] < 0) *deltad = fmin2(*deltad, -z1[i] / dz1[i]);
	if(dw[i]  < 0) *deltad = fmin2(*deltad, -w[i] / dw[i]);
    }
    for (i = 0; i < *n2; ++i) {
	if(dx2[i] < 0) *deltap = fmin2(*deltap, -x2[i] / dx2[i]);
	if(dz2[i] < 0) *deltad = fmin2(*deltad, -z2[i] / dz2[i]);
    }
    *deltap = fmin2(1., *beta * *deltap);
    *deltad = fmin2(1., *beta * *deltad);
    return;
} /* boundc */