/*****************************************************************************
*              RESOLUTION OF A SET OF NON-LINEAR EQUATIONS                   *
* -------------------------------------------------------------------------- *
* COPYRIGHT 1991, BY R.S. BAIN                                               *
* Originally wrapped by bain@luther.che.wisc.edu on Tue Mar 30 16:27:48 1993 *
*                                                                            *
* As of this writing, we do not know how to contact Rod Bain, the            * 
* author of NNES.                                                            *
*                                                                            *
* -- David M. Gay (dmg@bell-labs.com)                                        *
*    Bell Labs, Murray Hill                                                  *
*    8 February 1999                                                         *
*                                                                            *
*                                  Visual C++ Release By J-P Moreau, Paris.  *
*                                               (PART 2/3)                   *
*                                           (www.jpmoreau.fr)                *
*****************************************************************************/
#include <basis.h>
#include <vmblock.h>

    //Headers of functions used below
    void ataov (bool *, int, int, int, int, REAL **, REAL **, REAL *);
	void atvov (bool *, int, int, int, int, REAL **, REAL *, REAL *);
	void avmul (int, int, int, int, REAL **, REAL *, REAL *);
    void bakdif(bool *, int, int, REAL *, REAL, REAL *, REAL *, REAL **, REAL *);
    void bnddif(bool *, int, int, REAL, REAL *, REAL *, REAL *, REAL *, REAL **, REAL *);
    void deufls(bool *, bool, bool *, bool, bool *, bool *, bool *, bool *, bool, bool *, 
		        int *, int *, int *, int *, int *, int, int, int, int *, int, int, int, 
			    int *, REAL, REAL, REAL, REAL, REAL, REAL *, REAL, REAL *, REAL, REAL *, 
				REAL, REAL **, REAL **, REAL *, REAL *, REAL *, REAL *, REAL *, REAL **, 
				REAL *, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *);
 	REAL Dot_Product(int, REAL *, REAL *);
	void fcn1(bool *, int, REAL *, REAL *);
	void fcnevl(bool *, int, int, int, int, REAL, REAL *, REAL *, REAL *);
    void fordif(bool *, int, int, REAL, REAL *, REAL *, REAL **, REAL *);
   	void innerp (bool, bool *, int, int, int, int, int, int, REAL *, REAL *, REAL *);
	void jacob(bool *, int, REAL **, REAL *);
	void Line0(FILE *);
	void Line1(FILE *);
	void matcop (int, int, int, int, int, int, REAL **, REAL **);
    REAL MAX(REAL, REAL);
    REAL MIN(REAL, REAL);
	void mmulv (int, REAL **, REAL *, REAL *);
    void mmulvt(int, REAL **, REAL *, REAL *);
	void Msg(FILE *, char *);
	void qrupda (bool *, int, int, REAL, REAL **, REAL **, REAL *, REAL *);
    REAL Sign(REAL, REAL);
	void twonrm (bool *, int, int, REAL, REAL *, REAL *);
	void uvmul (int, int, int, int, REAL **, REAL *, REAL *);


void dogleg (bool *frstdg, bool *newtkn, bool overch, bool *overfl, int *maxexp, int n, 
			 int *notrst, int nunit, int output, REAL *beta, REAL caulen, REAL *delta, 
			 REAL etafac, REAL newlen, REAL sqrtz, REAL *delf, REAL *s, REAL *scalex, 
			 REAL *sn, REAL *ssdhat, REAL *vhat)  {
/*-------------------------------------------------------------------------
!    FEB. 23, 1992
!
!    THIS SUBROUTINE FINDS A TRUST REGION STEP USING THE
!    (DOUBLE) DOGLEG METHOD.
!------------------------------------------------------------------------*/
  REAL delfts, eta, factor, gamma, lambda, temp, tempv, zeta;
  int  i;
  char ss[80];

  extern FILE *fp;
  extern bool wrnsup;

  *overfl=FALSE;
  eta = ONE;

  if (output > 3) {
    Line0(fp);
    sprintf(ss,"    TRUST REGION STEP: %6d  TRUST REGION LENGTH, DELTA:%11.3f", *notrst, delta);
	Msg(fp,ss);
    Line0(fp);
    sprintf(ss,"    LENGTH OF NEWTON STEP, NEWLEN: %11.3f", newlen);
	Msg(fp,ss);
  }

// CHECK FOR NEWTON STEP WITHIN TRUST REGION - IF SO USE NEWTON STEP.

  if (newlen <= *delta) {
    for (i=1; i<=n; i++) s[i]=sn[i];
    *newtkn=TRUE;
    temp=*delta;
    *delta=newlen;
    if (output > 3) {
      Line0(fp);
      Msg(fp,"       NEWTON STEP WITHIN ACCEPTABLE RANGE ( <= THAN DELTA)");
      if (temp == *delta) {
        sprintf(ss,"   *       DELTA STAYS AT LENGTH OF NEWTON STEP: %11.3f", delta);
		Msg(fp,ss);
      }
      else
        Msg(fp,"       DELTA SET TO LENGTH OF NEWTON STEP.");
      Line0(fp);
      Msg(fp,"       FULL NEWTON STEP ATTEMPTED.");
    }
    return;
  }
  else {

/*  NEWTON STEP NOT WITHIN TRUST REGION - APPLY [DOUBLE)
    DOGLEG PROCEDURE.  (IF ETAFAC EQUALS 1.0 THEN THE SINGLE
    DOGLEG PROCEDURE IS BEING APPLIED).  */

    *newtkn=FALSE;
    if (*frstdg) {

/*  SPECIAL SECTION FOR FIRST DOGLEG STEP - CALCULATES
    CAUCHY POINT (MINIMIZER OF MODEL FUNCTION IN STEEPEST
    DESCENT DIRECTION OF THE OBJECTIVE FUNCTION).  */

      *frstdg=FALSE;
      if (output > 4) {
        Line0(fp);
        if (etafac == ONE)
          Msg(fp,"       FIRST SINGLE DOGLEG STEP");
        else
          Msg(fp,"       FIRST DOUBLE DOGLEG STEP");
        Line0(fp);
        Msg(fp,"          SCALED CAUCHY STEP.");
        Line0(fp);
      }

//  NOTE: BETA AND SQRTZ WERE CALCULATED IN SUBROUTINE DELCAU

      zeta=sqrtz*sqrtz;

//  FIND STEP TO CAUCHY POINT

      factor=-(zeta/(*beta));
      for (i=1; i<=n; i++) {
        ssdhat[i]=factor*(delf[i]/scalex[i]);
        if (output > 4) {
          sprintf(ss,"           SSDHAT(%3d) = %12.3f", i, ssdhat[i]);
          Msg(fp,ss);
        }
      }
      innerp (overch, overfl, *maxexp, n, n, n, nunit, output, &delfts, delf, sn);
      *overfl=FALSE;

//  PROTECT AGAINST (RARE) CASE WHEN CALCULATED DIRECTIONAL DERIVATIVE EQUALS ZERO.

      if (delfts != ZERO) {

//  STANDARD EXECUTION
        gamma=(zeta/ABS(delfts))*(zeta/(*beta));
        eta=etafac + (ONE-etafac)*gamma;
      }
      else {
		if (output > 1 && (!wrnsup)) {
          Line0(fp);
          Msg(fp,"    WARNING: DELFTS=0,  ETA SET TO 1.0 TO AVOID DIVISION BY ZERO");
        }
        eta=ONE;
      }
      if (output > 4) {
        Line0(fp);
        sprintf(ss,"          ETA = %11.3f", eta); Msg(fp,ss);
        Line0(fp);
        Msg(fp,"        VHAT VECTOR     VHAT(I) = ETA*SN(I)*SCALEX(I) - SSDHAT(I)");
        Line0(fp);
      }
      for (i=1; i<=n; i++) {
        vhat[i]=eta*scalex[i]*sn[i] - ssdhat[i];
        if (output > 4) {
          sprintf(ss,"            VHAT(%3d) = %12.3f", i, vhat[i]);
          Msg(fp,ss);
        }
      }
	} // if frstdg...
  } // if newlen...

// ETA*NEWLEN <= DELTA MEANS TAKE STEP IN NEWTON DIRECTION TO TRUST REGION BOUNDARY

  if (eta*newlen <= *delta) {
	if (output > 4) {
      Line0(fp);
      Msg(fp,"          ETA*NEWLEN <= DELTA     S(I) = (DELTA/NEWLEN)*SN(I)");
    }
    for (i=1; i<=n; i++)  s[i]=(*delta/newlen)*sn[i];
  }
  else {

//  DISTANCE TO CAUCHY POINT >= DELTA MEANS TAKE STEP IN
//  STEEPEST DESCENT DIRECTION TO TRUST REGION BOUNDARY.
	if (caulen >= *delta) {
	  if (output > 4) {
        Line0(fp);
        Msg(fp,"          CAULEN >= DELTA,  S(I)=(DELTA/CAULEN)*(SSDHAT(I)/SCALEX(I))");
      }
      for (i=1; i<=n; i++)  s[i] = (*delta/caulen)*(ssdhat[i]/scalex[i]);
    }
    else {

//    TAKE (DOUBLE) DOGLEG STEP

      innerp (overch, overfl, *maxexp, n, n, n, nunit, output, &temp, ssdhat, vhat);
      innerp (overch, overfl, *maxexp, n, n, n, nunit, output, &tempv, vhat, vhat);
      *overfl=FALSE;
      lambda=(-temp + SQRT(temp*temp - tempv*(caulen*caulen-(*delta)*(*delta))))/tempv;
      if (output > 4) {
        Line0(fp);
        Msg(fp,"          S(I) = (SSDHAT(I)+LAMBDA*VHAT(I))/SCALEX(I)");
        Line0(fp);
        sprintf(ss,"          WHERE LAMBDA = %12.3f", lambda); Msg(fp,ss);
      }
      for (i=1; i<=n; i++)  s[i]=(ssdhat[i] + lambda*vhat[i])/scalex[i];
    }
  }
  if (output > 3) {
    Line0(fp); Line0(fp);
    Msg(fp,"          REVISED STEP FROM SUBROUTINE DOGLEG:");
    Line0(fp);
    for (i=1; i<=n; i++) {
      sprintf(ss,"           S(%3d) = %12.3f", i, s[i]); Msg(fp,ss);
	}
  }
} // dogleg()


void gradf (bool overch, bool *overfl, bool restrt, bool sclfch, bool sclxch,
            int jupdm, int maxexp, int n, int nunit, int output, int qnupdm,
			REAL *delf, REAL *fvecc, REAL **jac, REAL *scalef, REAL *scalex) {
/*---------------------------------------------------------------------------
!    FEB. 23, 1992
!
!    THIS SUBROUTINE COMPUTES THE GRADIENT OF THE FUNCTION
!
!           F=1/2{SCALEF*FVECC)^(SCALEF*FVECC)
!
!    WHICH IS USED AS THE OBJECTIVE FUNCTION FOR MINIMIZATION.
!
!    NOTE: WHEN THE FACTORED FORM OF THE JACOBIAN IS UPDATED IN QUASI-NEWTON
!          METHODS THE GRADIENT IS UPDATED AS WELL IN THE SAME SUBROUTINE -
!          IT IS PRINTED HERE THOUGH.  IN THESE CASES QNUPDM > 0.
!--------------------------------------------------------------------------*/
  int i; REAL wv1[10]; char s[80];

  extern FILE *fp;
 
  if ((restrt) || jupdm == 0 || qnupdm == 0) {

//  GRADIENT NOT ALREADY UPDATED:    FIND DELF = J^F

    for (i=1; i<=n; i++)
      if (sclfch)
        wv1[i]=fvecc[i]*scalef[i]*scalef[i];
      else
        wv1[i]=fvecc[i];

	mmulvt(n,jac,wv1,delf);  // modified 04/05/07

//  PRINT GRADIENT VECTOR, DELF

    if (output > 3) {
	  if (! sclxch) {
        Line0(fp);
        if (! sclfch)
          Msg(fp,"    GRADIENT OF OBJECTIVE FUNCTION:");
        else                                                                            
          Msg(fp,"    GRADIENT OF SCALED OBJECTIVE FUNCTION:");
        Line0(fp);
        for (i=1; i<=n; i++) {
          sprintf(s,"    DELF(%3d) = %e", i, delf[i]); Msg(fp,s);  
		}
      }
	  else {
        Line0(fp);
        Msg(fp,"    GRADIENT OF OBJECTIVE FUNCTION IN SCALED X UNITS:");
        Line0(fp);
        for (i=1; i<=n; i++) {
          sprintf(s,"    DELF(%3d) = %e          DELF(%3d) = %12.3f",
                  i, delf[i], i, scalef[i]*scalef[i]*delf[i]/scalex[i]);
          Msg(fp,s);
        }
	  }
	}
  }

} // gradf()

void initch (bool *instop, bool *linesr, bool *newton, bool *overfl, bool *sclfch, 
             bool *sclxch, int acptcr, int *contyp, int jactyp, int jupdm, int maxexp, 
			 int n, int nunit, int output, int qnupdm, int stopcr, int trupdm, 
			 REAL *epsmch, REAL *fcnold, REAL ftol, REAL *boundl, REAL *boundu,
			 REAL *fvecc, REAL *scalef, REAL *scalex, REAL *xc)  {
/*----------------------------------------------------------------------------
!    AUG. 27, 1991
!
!    THIS SUBROUTINE FIRST CHECKS TO SEE IF N IS WITHIN THE ACCEPTABLE RANGE.
!
!    THE SECOND CHECK IS TO SEE IF THE INITIAL ESTIMATE IS
!    ALREADY A SOLUTION BY THE FUNCTION VALUE CRITERION, FTOL.
!
!    THE THIRD CHECK IS MADE TO SEE IF THE NEWTON OPTION IS BEING
!    USED WITH THE LINE SEARCH.  IF NOT A WARNING IS GIVEN AND
!    THE LINE SEARCH OPTION IS INVOKED.
!
!    THE FOURTH CHECK IS TO ENSURE APPLICABILITY OF SELECTED
!    VALUES FOR INTEGER CONSTANTS.
!
!    THE FIFTH CHECK IS TO WARN THE USER IF INITIAL ESTIMATES ARE NOT WITHIN
!    THE RANGES SET BY THE BOUNDL AND BOUNDU VECTORS.
!    CONTYP IS CHANGED FROM 0 TO 1 IF ANY BOUND HAS BEEN SET BY THE USER
!
!    THE SIXTH CHECK ENSURES BOUNDL[I) < BOUNDU[I) FOR ALL I.
!--------------------------------------------------------------------------*/
// Label: e130
  bool frster;
  REAL temp1, temp2;
  int i; char s[80];

  extern FILE *fp;

  *instop=FALSE;
  temp1=-pow(TEN,maxexp);
  temp2=pow(TEN,maxexp);

// CHECK FOR N IN RANGE

  if (n <= 0) {
    *instop=TRUE;
    Line1(fp);(fp); Line0(fp);
    Msg(fp,"  N IS OUT OF RANGE - RESET TO POSITIVE INTEGER.");
    Line0(fp); Line1(fp);(fp);
  }

// CHECK FOR SCALING FACTORS POSITIVE

  frster=TRUE;
  *sclfch=FALSE;
  *sclxch=FALSE;
  for (i=1; i<=n; i++) {
	if (scalef[i] <= ZERO) {
	  if (frster) {
        *instop=TRUE;
        frster=FALSE;
        Line1(fp);(fp);
      }
      Line0(fp);
      sprintf(s,"       SCALEF(%3d) = %12.3f    SHOULD BE POSITIVE", i, scalef[i]);
	  Msg(fp,s);
    }
    if (scalef[i] != ONE)  *sclfch=TRUE;
    if (scalex[i] <= ZERO) {
	  if (frster) {
        *instop=TRUE;
        frster=FALSE;
        Line1(fp);(fp);
      }
      Line0(fp);
      sprintf(s,"       SCALEX(%3d) = %12.3f    SHOULD BE POSITIVE", i, scalex[i]);
	  Msg(fp,s);
    }
    if (scalex[i] != ONE)  *sclxch=TRUE;
  }
  if (!frster) {
    Line0(fp); Line1(fp);(fp);
  }

// EVALUATE INITIAL RESIDUAL VECTOR AND OBJECTIVE FUNCTION AND
// CHECK TO SEE IF THE INITIAL GUESS IS ALREADY A SOLUTION.

  fcn1(overfl, n, fvecc, xc);

// NOTE: NUMBER OF LINE SEARCH FUNCTION EVALUATIONS, NFUNC,
// INITIALIZED AT 1 WHICH REPRESENTS THIS EVALUATION.

  if (*overfl) {
    Line1(fp);(fp); Line0(fp);
    Msg(fp,"       OVERFLOW IN INITIAL FUNCTION VECTOR EVALUATION.");
    Line0(fp); Line1(fp);(fp);
    *instop=TRUE;
    return;
  }

  fcnevl(overfl, maxexp, n, nunit, output, *epsmch, fcnold, fvecc, scalef);

  if (*overfl) {
    Line1(fp);(fp); Line0(fp);
    Msg(fp,"       OVERFLOW IN INITIAL OBJECTIVE FUNCTION EVALUATION.");
    Line0(fp); Line1(fp);(fp);
    *instop=TRUE;
    return;
  }

// CHECK FOR SOLUTION USING SECOND STOPPING CRITERION

  for (i=1; i<=n; i++)
    if (ABS(fvecc[i]) > ftol) goto e130;

  *instop=TRUE;
  Line1(fp);(fp); Line0(fp);
  Msg(fp,"  WARNING: THIS IS ALREADY A SOLUTION BY THE CRITERIA OF THE SOLVER.");
  Line0(fp);

// IF THE PROBLEM IS BADLY SCALED THE OBJECTIVE FUNCTION MAY MEET THE
// TOLERANCE ALTHOUGH THE INITIAL ESTIMATE IS NOT THE SOLUTION.

  Msg(fp,"  THIS MAY POSSIBLY BE ALLEVIATED BY RESCALING THE PROBLEM IF THE");
  Msg(fp,"  INITIAL ESTIMATE IS KNOWN NOT TO BE A SOLUTION.");
  Line0(fp); Line1(fp);(fp);

// CHECK FOR NEWTON'S METHOD REQUESTED BUT LINE SEARCH NOT BEING USED.

e130:;
  if ((*newton) && (!(*linesr))) {
    *linesr=TRUE;
    Line1(fp);(fp); Line0(fp);
    Msg(fp,"  WARNING: INCOMPATIBLE OPTIONS: NEWTON=TRUE AND LINESR=FALSE.");
    Msg(fp,"  LINESR SET TO TRUE; EXECUTION OF NEWTON METHOD CONTINUING.");
    Line0(fp); Line1(fp);(fp);
  }

// CHECK INTEGER CONSTANTS

  if (acptcr != 1 && acptcr != 12) {
    Line1(fp); Line0(fp);
    sprintf(s,"  ACPTCR NOT AN ACCEPTABLE VALUE: %5d", acptcr); Msg(fp,s);
    *instop=TRUE;
    Line0(fp); Line1(fp);(fp);
  }
  if (jactyp < 0 || jactyp > 3) {
    Line1(fp); Line0(fp);
    sprintf(s,"  JACTYP: %5d - NOT IN PROPER RANGE.", jactyp); Msg(fp,s);
    *instop=TRUE;
    Line0(fp); Line1(fp);(fp);
  }
  if (stopcr != 1 && stopcr != 12 && stopcr != 2 && stopcr == 3) {
    Line1(fp); Line0(fp);
    sprintf(s,"  STOPCR NOT AN ACCEPTABLE VALUE: %5d", stopcr); Msg(fp,s);
    *instop=TRUE;
    Line0(fp); Line1(fp);
  }
  if (qnupdm < 0 || qnupdm > 1) {
    Line1(fp); Line0(fp);
    sprintf(s,"  QNUPDM:%5d - NOT IN PROPER RANGE.", qnupdm); Msg(fp,s);
    *instop=TRUE;
    Line0(fp); Line1(fp);
  }
  if (trupdm < 0 || trupdm > 1) {
    Line1(fp); Line0(fp);
    sprintf(s,"  TRUPDM: %5d - NOT IN PROPER RANGE.", trupdm); Msg(fp,s);
    *instop=TRUE;
    Line0(fp); Line1(fp);
  }
  if (jupdm < 0 || jupdm > 2) {
    Line1(fp); Line0(fp);
    sprintf(s,"  JUPDM: %5d - NOT IN PROPER RANGE.", jupdm); Msg(fp,s);
    *instop=TRUE;
    Line0(fp); Line1(fp);
  }

// CHECK FOR INITIAL ESTIMATES NOT WITHIN SPECIFIED BOUNDS AND
// SET CONTYP TO 1 => AT LEAST ONE BOUND IS IN EFFECT.

  *contyp=0;
  for (i=1; i<=n; i++)
    if (boundl[i] != temp1 || boundu[i] != temp2) {
      *contyp=1;
      return;
    }

  frster=TRUE;
  if (*contyp != 0) {
    for (i=1; i<=n; i++) {

// CHECK FOR INITIAL ESTIMATES OUT OF RANGE AND LOWER
// BOUND GREATER THAN OR EQUAL TO THE UPPER BOUND.

	  if (xc[i] < boundl[i] || xc[i] > boundu[i]) {
	    if (frster) {
          *instop=TRUE;
          frster=FALSE;
          Line1(fp); Line0(fp);
          Msg(fp,"       COMPONENTS MUST BE WITHIN BOUNDS:");
          Line0(fp);
          Msg(fp,"      NO         XC               BOUNDL          BOUNDU");
          Line0(fp);
		}
        sprintf(s,"     %3d   %12.3f      %12.3f    %12.3f", i, xc[i], boundl[i], boundu[i]);
		Msg(fp,s);
      }
    }

    if (!frster) {
      Line0(fp); Line1(fp);
    }

    frster=TRUE;
    for (i=1; i<=n; i++)
	  if (boundl[i] >= boundu[i]) {
	    if (frster) {
          frster=FALSE;
          Line1(fp); Line0(fp);
          Msg(fp,"       LOWER BOUND MUST BE LESS THAN UPPER BOUND - VIOLATIONS LISTED.");
          Line0(fp);
        }
        sprintf(s,"       BOUNDL(%3d) = %12.3f    BOUNDU(%3d) = %12.3f", i, boundl[i], i, boundu[i]);
		Msg(fp,s);
      }

	if (!frster) {
      Line0(fp); Line1(fp);
    }
  }

} // initch()


void jaccd (int n, int nunit, int output, REAL epsmch, REAL *fvecj1, 
            REAL *fvecj2, REAL **jacfdm, REAL *scalex, REAL *xc)  {
/*---------------------------------------------------------------------
!    FEB. 11, 1991
!
!    THIS SUBROUTINE EVALUATES THE JACOBIAN USING CENTRAL DIFFERENCES.
!
!    FVECJ1 AND FVECJ2 ARE TEMPORARY VECTORS TO HOLD THE
!    RESIDUAL VECTORS FOR THE CENTRAL DIFFERENCE CALCULATION.
!--------------------------------------------------------------------*/
  bool overfl;
  int j,k;
  REAL curtep, deltaj, tempj;

  extern bool wrnsup;
  extern FILE *fp;

  overfl=FALSE;
  curtep=pow(epsmch, 0.33);

  for (j=1; j<=n; j++) {
    deltaj=curtep*Sign((MAX(ABS(xc[j]), ONE/scalex[j])),xc[j]);
    tempj=xc[j];
    xc[j]=xc[j] + deltaj;

//  NOTE: THIS STEP IS FOR FLOATING POINT ACCURACY ONLY.

    deltaj=xc[j] - tempj;

    fcn1(&overfl, n, fvecj1, xc);
    xc[j]=tempj-deltaj;
    fcn1(&overfl, n, fvecj2, xc);
    if ((overfl) && output > 2 && (!wrnsup)) {
      Line0(fp);
      Msg(fp,"    WARNING: OVERFLOW IN FUNCTION VECTOR IN <JACCD>.");
    }
    for (k=1; k<=n; k++)
      jacfdm[k][j]=(fvecj1[k] - fvecj2[k])/(TWO*deltaj);
    xc[j]=tempj;
  }
} // jaccd()

void matprt (int, int, REAL **);

void jacfd (int jactyp, int n, int nunit, int output, REAL epsmch,
            REAL *boundl, REAL *boundu, REAL *fvecc, REAL *fvecj1, 
            REAL **jacfdm, REAL *scalex, REAL *xc)  {
/*----------------------------------------------------------------------------
!   FEB. 15, 1991
!
!   THIS SUBROUTINE EVALUATES THE JACOBIAN USING ONE-SIDED FINITE DIFFERENCES.
!
!   JACTYP "1" SIGNIFIES FORWARD DIFFERENCES
!   JACTYP "2" SIGNIFIES BACKWARD DIFFERENCES
!
!   FVECJ1 IS A TEMPORARY VECTOR WHICH STORES THE RESIDUAL
!   VECTOR FOR THE FINITE DIFFERENCE CALCULATION.
!---------------------------------------------------------------------------*/
  REAL deltaj, sqrtep, tempj;
  bool overfl;
  int  j, k;
  char s[80];

  extern bool wrnsup;
  extern FILE *fp;

  sqrtep=SQRT(epsmch);

// FINITE-DIFFERENCE CALCULATION BY COLUMNS

  for (j=1; j<=n; j++) {

//  DELTAJ IS THE STEP SIZE - IT IS ALWAYS POSITIVE

    deltaj=sqrtep*MAX(ABS(xc[j]), ONE/scalex[j]);
    tempj=xc[j];
//  temporary storage of xc(j)
    if (jactyp == 1) {
	  if (xc[j]+deltaj <= boundu[j]) {
//    STEP WITHIN BOUNDS - COMPLETE FORWARD DIFFERENCE
        xc[j] += deltaj;
        deltaj=xc[j] - tempj;
        fordif (&overfl, j, n, deltaj, fvecc, fvecj1, jacfdm, xc);
      }
      else {

//     STEP WOULD VIOLATE BOUNDU - TRY BACKWARD DIFFERENCE

		if (xc[j]-deltaj >= boundl[j]) {
          xc[j] -= deltaj;
          bakdif (&overfl, j, n, &deltaj, tempj, fvecc, fvecj1, jacfdm, xc);
        }
        else {
/*      STEP WOULD ALSO VIOLATE BOUNDL - IF THE DIFFERENCE IN THE
        BOUNDS, (BOUNDU-BOUNDL), IS GREATER THAN DELTAJ CALCULATE
        THE FUNCTION VECTOR AT EACH BOUND AND USE THIS DIFFERENCE -
        THIS REQUIRES ONE EXTRA FUNCTION EVALUATION.
        THE CURRENT FVECC IS STORED IN WV3, THEN REPLACED.  */
		  if (boundu[j]-boundl[j] >= deltaj) {
            bnddif (&overfl, j, n, epsmch, boundl, boundu, fvecc, fvecj1, jacfdm, xc);
            if (output > 2 && (!wrnsup) && (!overfl)) {
              Line0(fp);
              Msg(fp,"    WARNING: BOUNDS TOO CLOSE FOR 1-SIDED FINITE-DIFFERENCES.");
              sprintf(s,"             LOWER AND UPPER BOUNDS USED FOR JACOBIAN COLUMN: %3d", j);
              Msg(fp,s);
              Msg(fp,"           THIS REQUIRED ONE EXTRA FUNCTION EVALUATION.");
            }
          }
          else {

//        BOUNDS ARE EXTREMELY CLOSE (BUT NOT EQUAL OR
//        THE PROGRAM WOULD HAVE STOPPED IN INITCH).

			if (output > 2 && (!wrnsup) && (!overfl)) {
              Line0(fp);
              Msg(fp,"    WARNING: BOUNDS TOO CLOSE FOR 1-SIDED FINITE-DIFFERENCES.");
              sprintf(s,"           BOUNDS ARE EXTREMELY CLOSE FOR COMPONENT: %3d", j);
			  Msg(fp,s);
              Msg(fp,"           FINITE DIFFERENCE JACOBIAN IS UNRELIABLE.");
            }
            bnddif (&overfl, j, n, epsmch, boundl, boundu, fvecc, fvecj1, jacfdm, xc);
          }
        }
      }
    }
    else  { //jactyp<>1
	  if (xc[j]-deltaj >= boundl[j]) {
        xc[j]=xc[j] - deltaj;
        bakdif (&overfl, j, n, &deltaj, tempj, fvecc, fvecj1, jacfdm, xc);
      }
      else {
		if (xc[j]+deltaj <= boundu[j]) {
          xc[j]=xc[j] + deltaj;
          fordif (&overfl, j, n, deltaj, fvecc, fvecj1, jacfdm, xc);
        }
        else {
		  if (boundu[j]-boundl[j] >= deltaj) {
            bnddif (&overfl, j, n, epsmch, boundl, boundu, fvecc, fvecj1, jacfdm, xc);
            if (output > 2 && (!wrnsup) && (!overfl)) {
              Line0(fp);
              Msg(fp,"    WARNING: BOUNDS TOO CLOSE FOR 1-SIDED FINITE-DIFFERENCES.");
              sprintf(s,"           LOWER AND UPPER BOUNDS USED FOR JACOBIAN COLUMN: %3d", j);
			  Msg(fp,s);
              Msg(fp,"           THIS REQUIRED ONE EXTRA FUNCTION EVALUATION.");
            }
          }
          else {
            bnddif (&overfl, j, n, epsmch, boundl, boundu, fvecc, fvecj1, jacfdm, xc);
            for (k=1; k<=n; k++)  jacfdm[k][j]=ZERO;
            if (output > 2 && (!wrnsup) && (!overfl)) {
              Line0(fp);
              Msg(fp,"    WARNING: BOUNDS TOO CLOSE FOR 1-SIDED FINITE-DIFFERENCES.");
              sprintf(s,"           BOUNDS ARE EXTREMELY CLOSE FOR COMPONENT: %3d", j);
			  Msg(fp,s);
              Msg(fp,"           FINITE DIFFERENCE JACOBIAN IS UNRELIABLE.");
            }
		  }
        }
      }
    }
    if ((overfl) && (output > 2) && (!wrnsup)) {
      Line0(fp);
      Msg(fp,"    WARNING: OVERFLOW IN FUNCTION VECTOR IN SUBROUTINE JACFD.");
      sprintf(s,"           BOUNDS ARE EXTREMELY CLOSE FOR COMPONENT: %3d", j);
	  Msg(fp,s);
      for (k=1; k<=n; k++) jacfdm[k][j]=ZERO;
    }
    xc[j]=tempj;
  }
} // jacfd()


void jacobi (bool checkj, bool *jacerr, bool *overfl, int jactyp, int n, int nunit, 
             int output, REAL epsmch, REAL *fdtolj, REAL *boundl, REAL *boundu,
	  	     REAL *fvecc, REAL *fvecj1, REAL *fvecj2, REAL **jac, REAL **jacfdm, 
			 REAL *scalex, REAL *xc)  {
/*------------------------------------------------------------------------
!    APR. 13, 1991
!
!    THIS SUBROUTINE EVALUATES THE JACOBIAN.  IF CHECKJ IS TRUE
!    THEN THE ANALYTICAL JACOBIAN IS CHECKED NUMERICALLY.
!
!    JACOB IS A USER-SUPPLIED ANALYTICAL JACOBIAN USED ONLY IF JACTYP=0.
!    THE JACOBIAN NAME MAY BE CHANGED BY USING THE EXTERNAL STATEMENT
!    IN THE MAIN DRIVER (NOT USED HERE).
!
!    JACFD ESTIMATES THE JACOBIAN USING FINITE DIFFERENCES:
!    FORWARD IF JACTYP=1 OR BACKWARD IF JACTYP=2.
!
!    JACCD ESTIMATES THE JACOBIAN USING CENTRAL DIFFERENCES.
!
!    IF THE ANALYTICAL JACOBIAN IS CHECKED THE FINITE DIFFERENCE
!    JACOBIAN IS STORED IN "JACFDM" AND THEN COMPARED.
!
!    FRSTER  INDICATES FIRST ERROR - USED ONLY TO SET BORDERS FOR OUTPUT
!    JACERR  FLAG TO INDICATE TO THE CALLING PROGRAM AN ERROR
!            IN THE ANALYTICAL JACOBIAN
!-----------------------------------------------------------------------*/
  bool frster;
  int i,j;
  char s[80];

  extern FILE *fp;

  frster=TRUE;
  *jacerr=FALSE;
  *overfl=FALSE;

  if (jactyp == 0)
    jacob (overfl, n, jac, xc);
  else if (jactyp == 1 || jactyp == 2)
    jacfd (jactyp, n, nunit, output, epsmch, boundl, boundu, fvecc, fvecj1, jac, scalex, xc);
  else
    jaccd (n, nunit, output, epsmch, fvecj1, fvecj2, jac, scalex, xc);

  if (jactyp == 0 && (checkj)) {

// NOTE: JACTYP=0 SENT TO JACFD PRODUCES A FORWARD DIFFERENCE ESTIMATE OF THE JACOBIAN.

    jacfd (jactyp, n, nunit, output, epsmch, boundl, boundu, fvecc, fvecj1, jacfdm, scalex, xc);

    for (j=1; j<=n; j++)
      for (i=1; i<=n; i++) {
	    if (ABS((jacfdm[i][j]-jac[i][j])/MAX(MAX(ABS(jac[i][j]),ABS(jacfdm[i][j])),ONE)) > *fdtolj) {
          *jacerr=TRUE;
          if (output >= 0) {
		    if (frster) {
              frster=FALSE;
              Line1(fp);
            }
            Line0(fp);
            sprintf(s,"    CHECK JACOBIAN TERM (%3d,%3d):", i, j); Msg(fp,s);
            Line0(fp);
            sprintf(s,"    ANALYTICAL DERIVATIVE IS %12.3f", jac[i][j]); Msg(fp,s);
            sprintf(s,"     NUMERICAL DERIVATIVE IS %12.3f", jacfdm[i][j]); Msg(fp,s);
            Line0(fp); Line1(fp);
          }
        }
      }
  }

} // jacobi()

void line (bool *abort, bool *absnew, bool deuflh, bool *geoms, bool newton, bool overch, 
		   bool *overfl, bool *qnfail, bool *qrsing, bool *restrt, bool *sclfch, bool *sclxch,
		   int *acpcod, int acptcr, int *contyp, int isejac, int *itnum, int jupdm, int maxexp, 
		   int *maxlin, int mgll, int mnew, int n, int narmij, int *nfunc, int nunit, int output, 
		   int qnupdm, int *stopcr, int *trmcod, REAL *alpha, REAL *confac, REAL epsmch, REAL *fcnmax, 
		   REAL *fcnnew, REAL fcnold, REAL lam0, REAL maxstp, REAL newlen, REAL *sbrnrm, REAL *sigma,
           REAL **a, REAL *boundl, REAL *boundu, REAL *delf, REAL *ftrack, REAL *fvec, REAL **h, 
		   REAL *hhpi, REAL **jac, REAL *rdiag, REAL *rhs, REAL *s, REAL *sbar, REAL *scalef, 
		   REAL *scalex, REAL *sn, REAL *strack, REAL *xc, REAL *xplus)  {
/*-------------------------------------------------------------
!    SEPT. 9, 1991
!------------------------------------------------------------*/
// Label: e90

  REAL delfts, lambda, mu, newmax, norm;
  bool convio;
  int  i,j,k;
  REAL wv2[10];
  char ss[80];

  extern bool wrnsup;
  extern FILE *fp;
  extern REAL nrmpre;

  convio=FALSE;
  *overfl=FALSE;
  *abort = FALSE;
  *qnfail = FALSE;

  if ((newton) || (*absnew)) {
    if (newton)

//  FIND NEXT ITERATE FOR PURE NEWTON'S METHOD
      for (k=1; k<=n; k++)  xplus[k]=xc[k] + sn[k];
    else

/*  FIND NEXT ITERATE FOR "ABSOLUTE" NEWTON'S METHOD.
    IF COMPONENT I WOULD BE OUTSIDE ITS BOUND THEN TAKE ABSOLUTE VALUE OF
    THE VIOLATION AND GO THIS DISTANCE INTO THE FEASIBLE REGION.
    ENSURE THAT THIS REFLECTION OFF ONE BOUND DOES NOT VIOLATE THE OTHER */

      for (i=1; i<=n; i++) {
        wv2[i]=ZERO;
        if (sn[i] >= ZERO)
		  if (xc[i]+sn[i] > boundu[i]) {
            convio=TRUE;
            wv2[i]=TWO;
            xplus[i]=MAX(TWO*boundu[i]-xc[i]-sn[i], boundl[i]);
          }
          else
            xplus[i]=xc[i] + sn[i];
        else
		  if (xc[i]+sn[i] < boundl[i]) {
            convio=TRUE;
            wv2[i]=-TWO;
            xplus[i]=MIN(TWO*boundl[i]-xc[i]-sn[i], boundu[i]);
          }
          else
            xplus[i]=xc[i] + sn[i];
      }

	if ((convio) && output > 4) {
      Line0(fp);
      Msg(fp,"       CONSTRAINT VIOLATORS IN ABSOLUTE NEWTON'S METHOD");
      Line0(fp);
      Msg(fp,"       COMPONENT  PROPOSED POINT  VIOLATED BOUND  FEASIBLE VALUE");
      Line0(fp);
      for (i=1; i<=n; i++) {
	    if (wv2[i] > ZERO) {
          sprintf(ss,"     %6d     %12.3f     %12.3f    %12.3f", i, (xc[i]+sn[i]), boundu[i],
                  xplus[i]); Msg(fp,ss);
        }
		else if (wv2[i] < ZERO) {
           sprintf(ss,"     %6d     %12.3f     %12.3f    %12.3f", i, (xc[i]+sn[i]), boundl[i],
                  xplus[i]); Msg(fp,ss);
        }
      }
    }
    fcn1(overfl, n, fvec, xplus);
    nfunc++;

    if (*overfl) {
      *overfl=FALSE;
      *fcnnew=pow(TEN, maxexp);
      if (output > 2) {
        Line0(fp);
        Msg(fp,"       POTENTIAL OVERFLOW DETECTED IN FUNCTION EVALUATION.");
      }
      return;
	}

    fcnevl (overfl, maxexp, n, nunit, output, epsmch, fcnnew, fvec, scalef);

//  RETURN FROM PURE NEWTON'S METHOD - OTHERWISE CONDUCT LINE SEARCH

  }

  if (output > 3) {
    Line0(fp); Line0(fp);
    Msg(fp,"    SUMMARY OF LINE SEARCH");
    Line0(fp);
  }

// SHORTEN NEWTON STEP IF LENGTH IS GREATER THAN MAXSTP

  if (newlen > maxstp) {
	if (output > 3) {
      Line0(fp);
      Msg(fp,"       LENGTH OF NEWTON STEP SHORTENED TO MAXSTP.");
    }
    for (k=1; k<=n; k++) sn[k] *= maxstp/newlen;
  }

// CHECK DIRECTIONAL DERIVATIVE (MAGNITUDE AND SIGN).

  innerp (overch, overfl, maxexp, n, n, n, nunit, output, &delfts, delf, sn);

  if (*overfl) {
    *overfl=FALSE;
    if (output > 2 && (!wrnsup)) {
      Line0(fp);
      sprintf(ss,"    WARNING: DIRECTIONAL DERIVATIVE, DELFTS, SET TO %12.3f", delfts);
	  Msg(fp,ss);
    }
  }

// REVERSE SEARCH DIRECTION IF DIRECTIONAL DERIVATIVE IS POSITIVE

  if (delfts > ZERO) {
    for (k=1; k<=n; k++)  sn[k]=-sn[k];
    delfts=-delfts;
    if (output > 2 && (!wrnsup)) {
      Line0(fp);
      Msg(fp,"    WARNING: DIRECTIONAL DERIVATIVE IS POSITIVE: DIRECTION REVERSED.");
    }
  }

// OUTPUT INFORMATION

  if (output > 3) {
    if (ABS(delfts)<100000)
      sprintf(ss,"       INNER PRODUCT OF DELF AND SN, DELFTS: .......%12.3f", delfts);
	else
      sprintf(ss,"       INNER PRODUCT OF DELF AND SN, DELFTS: .......%e", delfts);
	Msg(fp,ss);
    if (!(*sclxch)) {
      sprintf(ss,"       LENGTH OF NEWTON STEP, NEWLEN: ..............%12.3f", newlen);
	  Msg(fp,ss);
    }
    else {
      sprintf(ss,"       LENGTH OF SCALED NEWTON STEP, NEWLEN: .......%12.3f", newlen);
      Msg(fp,ss);
    }
    sprintf(ss,"       MAXIMUM STEP LENGTH ALLOWED, MAXSTP: ........%12.3f", maxstp);
  }

// ESTABLISH INITIAL RELAXATION FACTOR

  if (deuflh)

/*   AT FIRST STEP IN DAMPED NEWTON OR AFTER EXPLICIT JACOBIAN EVALUATION
     IN QUASI-NEWTON OR IF THE STEP SIZE IS WITHIN STOPPING TOLERANCE BUT
     STOPCR=3, THE LINE SEARCH IS STARTED AT LAMBDA=1.  */

    if ((isejac == 1 || *trmcod == 1) && *stopcr == 3)
      lambda=lam0;
    else {

      for (k=1; k<=n; k++)  wv2[k]=(sbar[k] - sn[k])*scalex[k];

      twonrm (overfl, maxexp, n, epsmch, &norm, wv2);

//    PREVENT DIVIDE BY ZERO IF NORM IS ZERO (UNDERFLOWS)

      if (norm < epsmch)
//    START LINE SEARCH AT LAMBDA=LAM0, USE DUMMY MU
        mu=TEN;
      else
        mu=nrmpre*lambda/norm;

      if (output > 4) {
        Line0(fp);
        sprintf(ss,"       DEUFLHARD TEST RATIO, MU: %11.3f", mu);
      }

//    SET INITIAL LAMBDA DEPENDING ON MU.   THIS IS A MODIFICATION OF
//    DEUFLHARD"S METHOD WHERE THE CUTOFF VALUE WOULD BE 0.7 FOR LAM0=1.0

      if (mu > lam0/TEN)
        lambda=lam0;
      else
        lambda=lam0/TEN;
	}
  else
    lambda=lam0;

//  STORE LENGTH OF NEWTON STEP.  IF NEWTON STEP LENGTH WAS
//  GREATER THAN MAXSTP IT WAS SHORTENED TO MAXSTP.

  nrmpre=MIN(maxstp,newlen);

//  ESTABLISH FCNMAX AND NEWMAX FOR NONMONOTONIC LINE SEARCH

  newmax=newlen;
  *fcnmax=fcnold;
  if (isejac > narmij)
    if (isejac < narmij+mgll)
	  for (j=1; j<=mnew; j++) {
        *fcnmax=MAX(*fcnmax,ftrack[j-1]);
        newmax=MAX(newmax,strack[j-1]);
      }
    else
      for (j=1; j<=mnew; j++) {
        *fcnmax=MAX(*fcnmax, ftrack[j]);
        newmax=MAX(newmax, strack[j]);
      }

  if (output > 3) {
    Line0(fp); Line0(fp);
    if (!(*sclxch))
      Msg(fp,"       LINE SEARCH");
    else
      Msg(fp,"       LINE SEARCH (X'S GIVEN IN UNSCALED UNITS)");
  }

// *** CONDUCT LINE SEARCH ***

  deufls (abort, deuflh, geoms, overch, overfl, qnfail, qrsing, restrt,
          sclfch, sclxch, acpcod, &acptcr, contyp, itnum, &jupdm, maxexp,
          *maxlin, n, nfunc, nunit, output, qnupdm, stopcr, *alpha, *confac,
          delfts, epsmch, *fcnmax, fcnnew, fcnold, &lambda, newmax, sbrnrm,
          *sigma, a, h, boundl, boundu, delf, fvec, hhpi, jac, rdiag, rhs,
          s, sbar, scalef, scalex, sn, xc, xplus);

} // line()

void llfa (bool overch, bool *overfl, bool *sclfch, bool *sclxch, int isejac, int maxexp, 
		   int n, int nunit, int output, REAL epsmch, REAL *omega, REAL **a, REAL *delf, 
		   REAL *fvec, REAL *fvecc, REAL **jac, REAL **plee, REAL *rdiag, REAL *s, 
		   REAL *scalef, REAL *scalex, REAL *xc, REAL *xplus)  {
/*--------------------------------------------------------
!    FEB. 23, 1991
!
!    THE LEE AND LEE QUASI-NEWTON METHOD IS APPLIED TO
!    THE FACTORED FORM OF THE JACOBIAN.
!
!    NOTE: T AND W ARE TEMPORARY WORKING VECTORS ONLY.
!-------------------------------------------------------*/
//Label: e10
  REAL denom, sqrtep, sum;
  REAL *t, *w, *wv3;
  REAL **temp1, **temp2;
  bool skipup;
  int  i,j,k;
  void *vmblock = NULL;

  extern FILE *fp;

  vmblock = vminit();
  t       = (REAL *) vmalloc(vmblock, VEKTOR,  n+1, 0);
  w       = (REAL *) vmalloc(vmblock, VEKTOR,  n+1, 0);  
  wv3     = (REAL *) vmalloc(vmblock, VEKTOR,  n+1, 0);
  temp1   = (REAL **) vmalloc(vmblock, MATRIX,  n+1, n+1);
  temp2   = (REAL **) vmalloc(vmblock, MATRIX,  n+1, n+1);
  if (! vmcomplete(vmblock))
  {
    LogError ("No Memory", 0, __FILE__, __LINE__);
    return;
  }    
  
  *overfl=FALSE;
  sqrtep=SQRT(epsmch);
  skipup=TRUE;

  for (i=1; i<=n; i++) {
    a[i][i]=rdiag[i];
    s[i]=xplus[i] - xc[i];
  }

/*  R IS IN THE UPPER TRIANGLE OF A.

        T=RS   */

  uvmul(n, n, n, n, a, s, t);

// FORM PART OF NUMERATOR AND CHECK TO SEE IF A SIGNIFICANT
// CHANGE WOULD BE MADE TO THE JACOBIAN.

  for (i=1; i<=n; i++) {
    for (k=1; k<=n; k++)  w[k]=jac[1][k];
    innerp (overch, overfl, maxexp, n, n, n, nunit, output, &sum, w, t);
    w[i]=scalef[i]*(fvec[i] - fvecc[i]) - sum;

// TEST TO ENSURE VECTOR W IS NONZERO.  IF W[I]=0 FOR
// ALL I THEN THE UPDATE IS SKIPPED - SKIPUP IS TRUE.

    if (ABS(w[i]) > sqrtep*scalef[i]*(ABS(fvec[i]) + ABS(fvecc[i])))
      skipup=FALSE;
// update to be performed
    else
      w[i]=ZERO;
  }

  if (!skipup) {

//  T=Q^W   Q^ IS STORED IN JAC

    avmul (n, n, n, n, jac, w, t);

//  FIND DENOMINATOR; FORM W=S^P (P IS SYMMETRIC SO PS IS FOUND)

    avmul (n, n, n, n, plee, s, w);
    if ((*sclxch))
//  SCALE W TO FIND DENOMINATOR
      for (k=1; k<=n; k++)  wv3[k]=w[k]*scalex[k]*scalex[k];
	else {
	  for (k= 1; k<=n; k++) {
        temp1[k][1]=w[k];
        temp2[k][1]=wv3[k];
      }
      matcop (n, n, 1, 1, n, 1, temp1, temp2);
      for (k=1; k<=n; k++)  wv3[k]=temp2[k][1];
    }

    innerp (overch, overfl, maxexp, n, n, n, nunit, output, &denom, wv3, s);

//  IF OVERFLOW WOULD OCCUR MAKE NO CHANGE TO JACOBIAN
    if (*overfl) {
	  if (output > 3) {
        Line0(fp);
        Msg(fp,"    WARNING: JACOBIAN NOT UPDATED TO AVOID OVERFLOW IN DENOMINATOR OF");
        Msg(fp,"    LEE AND LEE UPDATE.");
      }
      goto e10;     //free memory and return
    }

//  IF DENOM IS ZERO THE SOLVER IS PROBABLY NEAR SOLUTION -
//  AVOID OVERFLOW AND CONTINUE WITH SAME JACOBIAN.

    if (denom == ZERO) goto e10;

//  THE SCALED VERSION OF S IS TAKEN TO THE UPDATE

    for (k=1; k<=n; k++)  w[k] *= scalex[k]*scalex[k]/denom;

//  UPDATE THE QR DECOMPOSITION USING A SERIES OF GIVENS (JACOBI) ROTATIONS

    qrupda (overfl, maxexp, n, epsmch, a, jac, t, w);

//  RESET RDIAG AS DIAGONAL OF CURRENT R WHICH IS IN THE UPPER TRIANGLE OF A
    for (i=1; i<=n; i++)  rdiag[i]=a[i][i];

//  UPDATE P MATRIX
    denom=pow(*omega,(isejac+2)) + denom;
    plee[1][1] -=  wv3[1]*wv3[1]/denom;
    for (j=2; j<=n; j++) {
	  for (i=1; i<j; i++) {
        plee[i][j] -= wv3[i]*wv3[j]/denom;
        plee[j][i]=plee[i][j];
      }
      plee[j][j] -= wv3[j]*wv3[j]/denom;
    }    
  }

//  UPDATE THE GRADIENT VECTOR, DELF.
//  DELF = (QR)^F = R^Q^F = R^JAC F

  if (*sclfch)
    for (k=1; k<=n; k++)  w[k]=fvec[k]*scalef[k];
  else {
    for (k= 1; k<=n; k++) {
      temp1[k][1]=fvec[k];
      temp2[k][1]=w[k];
    }
    matcop (n, n, 1, 1, n, 1, temp1, temp2);
    for (k=1; k<=n; k++)  w[k]=temp2[k][1];
  }
  avmul (n, n, n, n, jac, w, t);
  mmulv (n, a, t, delf);

e10: vmfree(vmblock);  //free memory
} // llfa()


void llun (bool overch, bool *overfl, int isejac, int maxexp, int n, int nunit,
           int output, REAL epsmch, REAL omega, REAL *fvec, REAL *fvecc, REAL **jac, 
           REAL **plee, REAL *s, REAL *scalex, REAL *xc, REAL *xplus)  {
/*-------------------------------------------------------
!
!    FEB. 13, 1991
!
!    UPDATE THE JACOBIAN USING THE LEE AND LEE METHOD.
!------------------------------------------------------*/
//Label: e10
  REAL denom, sqrtep, sum, tempi;
  REAL *wv1,*temp,*temp1;
  int  i,j;
  void *vmblock = NULL;

  extern FILE *fp;

  vmblock = vminit();
  temp    = (REAL *) vmalloc(vmblock, VEKTOR,  n+1, 0);
  temp1   = (REAL *) vmalloc(vmblock, VEKTOR,  n+1, 0);  
  wv1     = (REAL *) vmalloc(vmblock, VEKTOR,  n+1, 0);
  if (! vmcomplete(vmblock))
  {
    LogError ("No Memory", 0, __FILE__, __LINE__);
    return;
  }    

  sqrtep=SQRT(epsmch);
  for (i=1; i<=n; i++)  s[i]=(xplus[i] - xc[i])*scalex[i];
  for (i=1; i<=n; i++) {
    for (j=1; j<=n; j++)  temp[j]=plee[j][i]; 
    wv1[i] = Dot_Product(n, s, temp);
  }

  innerp (overch, overfl, maxexp, n, n, n, nunit, output, &denom, wv1, s);

//IF OVERFLOW WOULD OCCUR MAKE NO CHANGE TO JACOBIAN

  if (*overfl) {
    if (output > 3) {
      Line0(fp);
      Msg(fp,"    WARNING: JACOBIAN NOT UPDATED TO AVOID OVERFLOW IN DENOMINATOR OF");
      Msg(fp,"    LEE AND LEE UPDATE.");
    }
    goto e10;
  }

// IF DENOM IS ZERO THE SOLVER MUST BE VERY NEAR SOLUTION -
// AVOID OVERFLOW AND CONTINUE WITH SAME JACOBIAN

  if (denom == ZERO)  goto e10;
  for (i=1; i<=n; i++) {
    for (j=1; j<=n; j++) {
      temp[i]=jac[i][j];
      temp1[i]=xplus[j] - xc[j];
    }
    sum = Dot_Product(n, temp, temp1);
    tempi=fvec[i] - fvecc[i] - sum;
    if (ABS(tempi) >= sqrtep*(ABS(fvec[i]) + ABS(fvecc[i]))) {
      tempi /= denom;
      for (j=1; j<=n; j++) jac[i][j] += tempi*wv1[j]*scalex[j];
    }
  }

// UPDATE P MATRIX

  denom=pow(omega,isejac+2) + denom;
  plee[1][1] -= wv1[1]*wv1[1]/denom;
  for (j=2; j<=n; j++) {
    for (i=1; i<j; i++) {
      plee[i][j] -= wv1[i]*wv1[j]/(denom*scalex[i]*scalex[j]);
      plee[j][i]=plee[i][j];
    }
    plee[j][j] -= wv1[j]*wv1[j]/denom;
  }

e10: free(vmblock);
} // llun()


void matprt (int nrowpr, int ncolpr, REAL **a)  {
/*----------------------------------------------------------------------------
!    FEB. 6, 1991
!
!    THIS SUBROUTINE PRINTS RECTANGULAR BLOCKS STARTING WITH ELEMENT A[1,1]
!    OF SIZE NROWPR BY NCOLPR FOR MATRIX A [WHICH HAS DECLARED SIZE NROWA BY
!    NCOLA].  THE MATRIX IS PRINTED AS A BLOCK FOR SIZES UP TO 5X5 OR BY
!    COLUMNS IF IT IS LARGER.
!
!    NROWPR IS THE NUMBER OF ROWS TO BE PRINTED
!    NCOLPR IS THE NUMBER OF COLUMNS TO BE PRINTED
!
!    IF MATRIX PRINTING IS TO BE SUPPRESSED THEN LOGICAL
!    VARIABLE MATSUP MUST BE SET TO TRUE BEFORE THE CALL TO NNES.
!---------------------------------------------------------------------------*/
  int i, j, k, limit;
  char s[80]; char s1[13]; char s2[4];

  extern bool matsup;
  extern FILE *fp;

  if (matsup) {
    Line0(fp);
    Msg(fp,"       MATRIX PRINTING SUPPRESSED.");
    return;
  }

// FOR NCOLPR <= 5 WRITE MATRIX AS A WHOLE
  strcpy(s,""); strcpy(s1,""); strcpy(s2,"");
  Line0(fp);
  if (ncolpr <= 5) {
    strcpy(s,"           ");
    for (k= 1; k<=ncolpr; k++) {
      sprintf(s2,"%3d", k); strcat(s,s2); strcat(s,"        ");
    }
    Msg(fp,s); Line0(fp);
    for (i=1; i<=nrowpr; i++) {
      strcpy(s,"");
      sprintf(s2,"%3d", i); strcat(s,s2); strcat(s,"     ");
      for (k= 1; k<=ncolpr; k++) {
        sprintf(s1,"%10.3f", a[i][k]); strcat(s,s1); strcat(s," ");
      }
      Msg(fp,s);
    }
    Line0(fp);
  }
  else {

//  LIMIT IS THE NUMBER OF GROUPS OF 5 COLUMNS
    limit=ncolpr / 5;

//  WRITE COMPLETE BLOCKS FIRST (LEFTOVERS LATER)
    for (j=1; j<=limit; j++) {
      strcpy(s,"           ");
      for (k= 1+(j-1)*5; k<=5+(j-1)*5; k++) {
        sprintf(s2,"%3d", k); strcat(s,s2); strcat(s,"        ");
      }
      Msg(fp,s); Line0(fp);
      for (i=1; i<=nrowpr; i++) {
        strcpy(s,"");
        sprintf(s2,"%3d", i); strcat(s,s2); strcat(s,"     ");
		for (k= 1+(j-1)*5; k<= 5+(j-1)*5; k++) {
          sprintf(s1,"%10.3f", a[i][k]); strcat(s,s1); strcat(s," ");
        }
        Msg(fp,s);
      }
      Line0(fp);
    }

//  WRITE REMAINING ELEMENTS
    strcpy(s,"           ");
    for (k= 5*limit+1; k<=ncolpr; k++) {
      sprintf(s2,"%3d", k); strcat(s,s2); strcat(s,"        ");
    }
    Msg(fp,s); Line0(fp);
    for (i=1; i<=nrowpr; i++) {
      strcpy(s,"");
      sprintf(s1,"%3d", i); strcat(s,s1); strcat(s,"     ");
      for (k= 5*limit+1; k<=ncolpr; k++) {
        sprintf(s1,"%10.3f", a[i][k]); strcat(s,s1); strcat(s," ");
      }
      Msg(fp,s);
    }
    Line0(fp);
  }

} // matprt()


void maxst (bool *overfl, int maxexp, int n, int nunit, int output,
            REAL epsmch, REAL *maxstp, REAL mstpf, REAL *scalex, REAL *xc) {
/*-------------------------------------------------------------------------
!
!    FEB. 11, 1991
!
!    THIS SUBROUTINE ESTABLISHES A MAXIMUM STEP LENGTH BASED ON THE 2-NORMS
!    OF THE INITIAL ESTIMATES AND THE SCALING FACTORS MULTIPLIED BY A
!    FACTOR MSTPF.
!
!         MAXSTP=MSTPF*MAX(NORM1,NORM2)
!
!              WHERE   MSTPF  USER-CHOSEN FACTOR (DEFAULT: 1000)
!                      NORM1  2-NORM OF SCALED STARTING ESTIMATES
!                      NORM2  2-NORM OF COMPONENT SCALING FACTORS
!------------------------------------------------------------------------*/
  REAL norm1, norm2;
  REAL wv1[10];
  int k;
  char s[80];

  extern FILE *fp;
  extern bool wrnsup;

  *overfl=FALSE;
  for (k=1; k<=n; k++)  wv1[k]=scalex[k]*xc[k];

  twonrm (overfl, maxexp, n, epsmch, &norm1, wv1);

  if (*overfl) {
    *maxstp=pow(TEN, maxexp);
    if (output > 2 && (!wrnsup)) {
      Line0(fp);
      sprintf(s,"    WARNING: NORM OF SCALED INITIAL ESTIMATE SET TO %12.3f", norm1); Msg(fp,s);
      Line0(fp);
      sprintf(s,"       MAXIMUM STEP SIZE, MAXSTP, SET TO %12.3f", maxstp); Msg(fp,s);
    }
    return;
  }

  twonrm (overfl, maxexp, n, epsmch, &norm2, scalex);

  if (*overfl) {
    *maxstp=pow(TEN, maxexp);
    if (output > 2 && (!wrnsup)) {
      Line0(fp);
      sprintf(s,"    WARNING: NORM OF SCALING FACTORS SET TO %12.3f", norm2); Msg(fp,s);
      Line0(fp);
      sprintf(s,"       MAXIMUM STEP SIZE, MAXSTP, SET TO %12.3f", *maxstp); Msg(fp,s);
    }
    return;
  }

  *maxstp = mstpf*MAX(norm1,norm2);

} // maxst()

void nersl (bool newton, bool restrt, bool sclfch, bool sclxch, int acpcod, 
			int jupdm, int n, int nunit, int output, REAL fcnnew, REAL *fvec, REAL *xplus) {
/*------------------------------------------------------
!    SEPT. 2, 1991
!
!    THE RESULTS OF EACH ITERATION ARE PRINTED.
!-----------------------------------------------------*/
  int i; char s[80];

  extern FILE *fp;

  Line0(fp);
  Line0(fp);
  if (!sclxch)
    Msg(fp,"    SUMMARY OF ITERATION RESULTS");
  else
    Msg(fp,"    SUMMARY OF ITERATION RESULTS (X'S GIVEN IN UNSCALED UNITS)");
  Line0(fp);
  Msg(fp,"        UPDATED ESTIMATES               UPDATED FUNCTION VALUES");
  Line0(fp);
  for (i=1; i<=n; i++) {
    if(ABS(fvec[i])<100000)
      sprintf(s,"      X(%3d) = %12.3f             F(%3d) = %12.3f", i, xplus[i], i, fvec[i]);
	else
      sprintf(s,"      X(%3d) = %12.3f             F(%3d) = %e", i, xplus[i], i, fvec[i]);
    Msg(fp,s);
  }
  Line0(fp);
  if (!sclfch) {
    sprintf(s,"      OBJECTIVE FUNCTION VALUE: %e", fcnnew); Msg(fp,s);
  }
  else {
    sprintf(s,"      SCALED OBJECTIVE FUNCTION VALUE: %e", fcnnew); Msg(fp,s);
  }
  Line0(fp);
  if (output > 3 && (!newton)) {
    sprintf(s,"      STEP ACCEPTANCE CODE, ACPCOD: %9d", acpcod); Msg(fp,s);
  }
  if ((restrt) && output >= 3 &&  jupdm != 0) {
    if (output > 3)  Line0(fp);
    Msg(fp,"      NOTE: JACOBIAN EVALUATED EXPLICITLY AT THIS STEP.");
    Line0(fp);
  }
} // nersl()


void nestop (bool *absnew, bool linesr, bool *newton, bool *sclfch, bool *sclxch,
             int acptcr, int *itnum, int n, int *nac1, int *nac2, int *nac12, int *nfunc, 
             int *njetot, int nunit, int output, int stopcr, int *trmcod, REAL *fcnnew,
			 REAL ftol, REAL nsttol, REAL *stpmax, REAL stptol, REAL *fvec, REAL *scalef, 
			 REAL *scalex, REAL *xc, REAL *xplus)  {
/*---------------------------------------------------------------------------
!      FEB. 23, 1992
!
! *** THIS SUBROUTINE CHECKS TO SEE IF THE CONVERGENCE CRITERIA HAVE BEEN MET
!     AND PRINTS MAIN RESULTS TO OUTPUT FILE.  ***
!--------------------------------------------------------------------------*/
// Label: e100
  REAL MAX1, max2, ratio1, sum;
  int i; char s[80];

  extern FILE *fp;
  extern nfetot;

  if (output > 3) {
    Line0(fp);
    Msg(fp,"    CONVERGENCE TESTING");
    Line0(fp);
  }

// IF THE NEWTON STEP WAS WITHIN TOLERANCE THEN TRMCOD IS 1

  if (*trmcod == 1) {
	if (output > 3) {
      Line0(fp);
      if (!(*sclxch)) {
        sprintf(s,"    MAXIMUM NEWTON STEP LENGTH STPMAX: %12.3f", stpmax);
		Msg(fp,s);
      }
      else {
        sprintf(s,"    MAXIMUM SCALED NEWTON STEP LENGTH STPMAX: %12.3f", stpmax);
        Msg(fp,s);
      }
      sprintf(s,"      FIRST CONVERGENCE CRITERION MET;  NSTTOL IS: %12.3f", nsttol);
	  Msg(fp,s);
      Line0(fp);

//    SKIP CHECKING OTHER STEP SIZE CRITERION AS TRMCOD IS ALREADY 1
      goto e100;
    }
  }

/*  IF THE NEWTON STEP WAS NOT WITHIN TOLERANCE THEN, IF
    STOPCR IS NOT EQUAL TO 2, THE SECOND STEP SIZE STOPPING
    CRITERION MUST BE CHECKED.  */

  if (stopcr != 2 && *trmcod != 1) {
    MAX1=ZERO;
    for (i=1; i<=n; i++) {
      ratio1=(ABS(xplus[i]-xc[i]))/MAX(ABS(xplus[i]), ONE/scalex[i]);
      MAX1=MAX(MAX1, ratio1);
      if (output > 4) {
        sprintf(s,"    RELATIVE STEP SIZE (%3d) = %12.3f", i, ratio1);
        Msg(fp,s);
      }
    }
    if (output > 4)  Line0(fp);
    if (output > 3) {
	  if (!(*sclxch)) {
        sprintf(s,"    MAXIMUM STEP SIZE: %12.3f   STPTOL: %11.3f", MAX1, stptol);
        Msg(fp,s);
	  }
      else {
        sprintf(s,"    MAXIMUM RELATIVE STEP SIZE: %12.3f   STPTOL: %11.3f", MAX1, stptol);
        Msg(fp,s);
      }
      Line0(fp);
    }
    if (MAX1 < stptol)  *trmcod=1;
  }

/*  NOTE: CONTINUATION AT LABEL 100 MEANS THAT TRMCOD WAS 1 ON ENTRY
          SO THE STEP SIZE CRITERION ABOVE DID NOT NEED TO BE CHECKED.

    THE SECOND STOPPING CRITERION IS CHECKED IF NEEDED.  */

e100:if (stopcr == 2 || stopcr == 12 || (stopcr == 3 && *trmcod == 1)) {
    
    max2=ZERO;

    for (i=1; i<=n; i++) {
      max2=MAX(max2,scalef[i]*ABS(fvec[i]));
      if (output > 4)
	    if (!(*sclfch)) {
          sprintf(s,"   ABSOLUTE FUNCTION VECTOR (%3d) = %12.3f", i, ABS(fvec[i]));
          Msg(fp,s);
        }
        else {
          sprintf(s,"   ABSOLUTE SCALED FUNCTION VECTOR (%3d) = %12.3f", i, (scalef[i]*ABS(fvec[i])));
          Msg(fp,s);
        }
    }

    if (output > 4)  Line0(fp);
    if (output > 3) {
	  if (!(*sclfch)) {
        sprintf(s,"    MAXIMUM ABSOLUTE FUNCTION: %e     FTOL: %e", max2, ftol); Msg(fp,s);
      }
      else {
        sprintf(s,"    MAX ABSOLUTE SCALED FUNCTION: %e     FTOL: %e", max2, ftol); Msg(fp,s);
      }
      Line0(fp);
    }

    if (max2 < ftol)
      if (stopcr == 3 && *trmcod == 1)

//    BOTH NEEDED STOPPING CRITERIA HAVE BEEN MET
        *trmcod=3;

      else if (stopcr == 12 && *trmcod == 1)

//    BOTH STOPPING CRITERIA HAVE BEEN MET ALTHOUGH
//    EITHER ONE WOULD BE SATISFACTORY.
        *trmcod=12;

      else if (stopcr == 2 || stopcr == 12)
        *trmcod=2;
    else if (stopcr == 3)

//    ONLY THE FIRST STOPPING CRITERION WAS MET - TRMCOD
//    MUST BE RESET FROM 1 BACK TO 0.
      *trmcod=0;

 }

//   PRINT FINAL RESULTS IF CONVERGENCE REACHED

  if (*trmcod > 0) {
    if (output > 0) {
      if (output == 1)  Line1(fp);;
      Line0(fp);
      sprintf(s,"    CONVERGENCE REACHED; TERMINATION CODE: ...............%6d", *trmcod); Msg(fp,s);

//   ITERATION RESULTS NOT PRINTED IN NERSL

      Line0(fp); Line0(fp);
      if ((*sclfch) || (*sclxch)) {
        Msg(fp,"                           UNSCALED RESULTS");
        Line0(fp);
      }
      Msg(fp,"          FINAL ESTIMATES                 FINAL FUNCTION VALUES");
      Line0(fp);
      for (i=1; i<=n; i++) {
        sprintf(s,"    X(%3d) = %15.6f            F(%3d) = %15.6f", i, xplus[i], i, fvec[i]);
        Msg(fp,s); 
      }
      Line0(fp);

      if (sclfch) {
//      NEED UNSCALED OBJECTIVE FUNCTION
        sum = Dot_Product(n, fvec, fvec);
        sprintf(s,"    FINAL OBJECTIVE FUNCTION VALUE: %e", sum/TWO); Msg(fp,s);
      }
      else {
        sprintf(s,"    FINAL OBJECTIVE FUNCTION VALUE: %e", fcnnew); Msg(fp,s);
      }
      if ((*sclfch) || (*sclxch)) {
        Line0(fp); Line0(fp);
        Msg(fp,"                           SCALED RESULTS");
        Line0(fp);
        Msg(fp,"          FINAL ESTIMATES                  FINAL FUNCTION VALUES");
        Line0(fp);
        for (i=1; i<=n; i++) {
          sprintf(s,"    X(%3d) = %15.6f            F(%3d) = %15.6f", i, scalex[i]*xplus[i], i, scalef[i]*fvec[i]);
		  Msg(fp,s);
        }
        Line0(fp);
        sprintf(s,"    FINAL OBJECTIVE FUNCTION VALUE: %e", fcnnew); Msg(fp,s);
      }
	}
  }
  else {
    if (output > 3) {
      Msg(fp,"    CONVERGENCE NOT REACHED."); Line0(fp);
    }
    return;
  }

// TERMINATION HAS BEEN REACHED

  if (output > 0) {
    Line0(fp); Line0(fp);
    sprintf(s,"    TOTAL NUMBER OF ITERATIONS: ..........................%6d", *itnum); Msg(fp,s);
    if (!(*newton) && !(*absnew))
	  if (linesr) {
        sprintf(s,"    TOTAL NUMBER OF LINE SEARCH FUNCTION EVALUATIONS: ....%6d", *nfunc);
        Msg(fp,s);
      }
      else {
        sprintf(s,"    TOTAL NUMBER OF TRUST REGION FUNCTION EVALUATIONS: ...%6d", *nfunc);
        Msg(fp,s);
      }

    sprintf(s,"    TOTAL NUMBER OF EXPLICIT JACOBIAN EVALUATIONS: .......%6d", njetot); Msg(fp,s);
    sprintf(s,"    TOTAL NUMBER OF FUNCTION EVALUATIONS: ................%6d", nfetot); Msg(fp,s);
    Line0(fp);

    if ((!newton) && (!absnew) && acptcr != 1 && output > 2) {
      sprintf(s,"    NUMBER OF STEPS ACCEPTED BY FUNCTION VALUE ONLY: .....%6d", nac1); Msg(fp,s);
      sprintf(s,"    NUMBER OF STEPS ACCEPTED BY STEP SIZE VALUE ONLY: ....%6d", nac2); Msg(fp,s);
      sprintf(s,"    NUMBER OF STEPS ACCEPTED BY EITHER CRITERION: ........%6d", nac12); Msg(fp,s);
      Line0(fp);
    }
    Line1(fp);
  }

} //nestop()

// end of file unnes1.cpp