/* linmin.f -- translated by f2c (version 19950314).
   You must link the resulting object file with the libraries:
	-lF77 -lI77 -lm   (in that order)
*/

#include "f2c.h"

/* Common Block Declarations */

struct {
    char argz[512];
} argz_;

#define argz_1 argz_

struct {
    doublereal cosine;
} gravec_;

#define gravec_1 gravec_

struct {
    integer numcal;
} numcal_;

#define numcal_1 numcal_

struct {
    char keywrd[80];
} keywrd_;

#define keywrd_1 keywrd_

struct {
    integer numat, nat[86], nfirst[86], nmidle[86], nlast[86], norbs, nelecs, 
	    nalpha, nbeta, nclose, nopen, ndumy;
    doublereal fract;
} molkst_;

#define molkst_1 molkst_

/* Table of constant values */

static integer c__10 = 10;
static integer c__1 = 1;
static logical c_true = TRUE_;
static logical c_false = FALSE_;

/* Subroutine */ int linmin_(doublereal *xparam, doublereal *step, doublereal 
	*pvect, integer *nvar, doublereal *funct, logical *okf, logical *okc)
{
    /* Initialized data */

    static integer icalcn = 0;

    /* Format strings */
    static char fmt_230[] = "(\002 ---QLINMN \002,/5x,\002LEFT   ...\002,f17"
	    ".8,f17.11/5x,\002CENTER ...\002,f17.8,f17.11,/5x,\002RIGHT  .."
	    ".\002,f17.8,f17.11,/)";
    static char fmt_240[] = "(5x,\002LEFT    ...\002,f17.8,f17.11,/5x,\002CE"
	    "NTER  ...\002,f17.8,f17.11,/5x,\002RIGHT   ...\002,f17.8,f17.11,"
	    "/5x,\002NEW     ...\002,f17.8,f17.11,/)";

    /* System generated locals */
    integer i__1, i__2;
    doublereal d__1, d__2;

    /* Builtin functions */
    integer i_indx(char *, char *, ftnlen, ftnlen), pow_ii(integer *, integer 
	    *), s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(
	    void);
    double d_sign(doublereal *, doublereal *);

    /* Local variables */
    static doublereal aabs, beta, grad, pabs, alfs, fmin;
    static integer left;
    static doublereal fmax;
    static integer ictr;
    static doublereal drop, tiny;
    static integer i;
    extern doublereal reada_(char *, integer *, ftnlen);
    static logical halfe;
    static doublereal gamma, s, alpha, angle;
    static integer right;
    static doublereal xminm, xcrit, xmaxm;
    static integer iquit;
    static logical print;
    static doublereal estor, xstor[258], vt[3];
    extern /* Subroutine */ int compfg_(doublereal *, logical *, doublereal *,
	     logical *, doublereal *, logical *), exchng_(doublereal *, 
	    doublereal *, doublereal *, doublereal *, doublereal *, 
	    doublereal *, doublereal *, doublereal *, integer *);
    static integer center;
    static logical fulscf, askful;
    static integer maxlin;
    static doublereal energy, funold, stlast, ymaxst, ssqlst, sqstor, tee, 
	    fin, phi[3], eps, xxm;

    /* Fortran I/O blocks */
    static cilist io___15 = { 0, 6, 0, "('  FULL SCF CALCULATIONS:',3X,L1)", 
	    0 };
    static cilist io___35 = { 0, 6, 0, fmt_230, 0 };
    static cilist io___44 = { 0, 6, 0, fmt_240, 0 };
    static cilist io___46 = { 0, 6, 0, "('  TINY',F14.9)", 0 };


/* COMDECK SIZES */
/************************************************************************
****/
/*  THIS FILE CONTAINS ALL THE ARRAY SIZES FOR USE IN MOPAC.              
**/
/*                                                                        
**/
/*    THERE ARE ONLY  PARAMETERS THAT THE PROGRAMMER NEED SET:            
**/
/*    MAXHEV = MAXIMUM NUMBER OF HEAVY ATOMS (HEAVY: NON-HYDROGEN ATOMS)  
**/
/*    MAXLIT = MAXIMUM NUMBER OF HYDROGEN ATOMS.                          
**/
/*    MAXTIM = DEFAULT TIME FOR A JOB. (SECONDS)                          
**/
/*    MAXDMP = DEFAULT TIME FOR AUTOMATIC RESTART FILE GENERATION (SECS)  
**/
/*                                                                        
**/
/*                                                                        
**/
/************************************************************************
****/
/*                                                                        
**/
/*  THE FOLLOWING CODE DOES NOT NEED TO BE ALTERED BY THE PROGRAMMER      
**/
/*                                                                        
**/
/************************************************************************
****/
/*                                                                        
**/
/*   ALL OTHER PARAMETERS ARE DERIVED FUNCTIONS OF THESE TWO PARAMETERS   
**/
/*                                                                        
**/
/*     NAME                   DEFINITION                                  
**/
/*    NUMATM         MAXIMUM NUMBER OF ATOMS ALLOWED.                     
**/
/*    MAXORB         MAXIMUM NUMBER OF ORBITALS ALLOWED.                  
**/
/*    MAXPAR         MAXIMUM NUMBER OF PARAMETERS FOR OPTIMISATION.       
**/
/*    N2ELEC         MAXIMUM NUMBER OF TWO ELECTRON INTEGRALS ALLOWED.    
**/
/*    MPACK          AREA OF LOWER HALF TRIANGLE OF DENSITY MATRIX.       
**/
/*    MORB2          SQUARE OF THE MAXIMUM NUMBER OF ORBITALS ALLOWED.    
**/
/*    MAXHES         AREA OF HESSIAN MATRIX                               
**/
/************************************************************************
****/
/************************************************************************
****/
/*  FOR SHORT VERSION USE LINE WITH NMECI=1, FOR LONG VERSION USE LINE    
**/
/*  WITH NMECI=10                                                         
**/
/************************************************************************
****/
/*     PARAMETER (NMECI=1,   NPULAY=1) */
/************************************************************************
****/
/* DECK MOPAC */
/* next line added for Unix implementation for command line arguments */

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

/*  LINMIN DOES A LINE MINIMISATION. */

/*  ON INPUT:  XPARAM = STARTING COORDINATE OF SEARCH. */
/*             STEP   = STEP SIZE FOR INITIATING SEARCH. */
/*             PVECT  = DIRECTION OF SEARCH. */
/*             NVAR   = NUMBER OF VARIABLES IN XPARAM. */
/*             FUNCT  = INITIAL VALUE OF THE FUNCTION TO BE MINIMIZED. */
/*             ISOK   = NOT IMPORTANT. */
/*             COSINE = COSINE OF ANGLE OF CURRENT AND PREVIOUS GRADIENT. 
*/

/*  ON OUTPUT: XPARAM = COORDINATE OF MINIMUM OF FUNCTI0N. */
/*             STEP   = NEW STEP SIZE, USED IN NEXT CALL OF LINMIN. */
/*             PVECT  = UNCHANGED, OR NEGATED, DEPENDING ON STEP. */
/*             FUNCT  = FINAL, MINIMUM VALUE OF THE FUNCTION. */
/*             OKF    = TRUE IF LINMIN IMPROVED FUNCT, FALSE OTHERWISE. */
/*             OKC    = TRUE IF LINMIN FOUND THE MINIMUM, FALSE OTHERWISE 
*/

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

/*  THE FOLLOWING COMMON IS USED TO FIND OUT IF A NON-VARIATIONALLY */
/*  OPTIMIZED WAVE-FUNCTION IS BEING USED. */

    /* Parameter adjustments */
    --pvect;
    --xparam;

    /* Function Body */
    if (icalcn != numcal_1.numcal) {
	halfe = i_indx(keywrd_1.keywrd, "C.I.", 80L, 4L) != 0 || 
		molkst_1.nclose != molkst_1.nopen;
	askful = i_indx(keywrd_1.keywrd, "FULSCF", 80L, 6L) != 0;
	drop = 2e-5;
	if (i_indx(keywrd_1.keywrd, "PREC", 80L, 4L) != 0) {
	    drop *= .01;
	}
	if (i_indx(keywrd_1.keywrd, "GNORM", 80L, 5L) != 0) {
/* Computing MIN */
	    i__1 = i_indx(keywrd_1.keywrd, "GNORM", 80L, 5L);
	    d__1 = reada_(keywrd_1.keywrd, &i__1, 80L);
	    drop *= min(d__1,1.);
	}
	xmaxm = .4;
	i = 2;
	*step = 1.;
	maxlin = 15;
	xcrit = 1e-4;
	if (i_indx(keywrd_1.keywrd, "FORCE", 80L, 5L) != 0) {
	    i = 3;
	    xcrit = 1e-5;
	}
	angle = .8;
	if (halfe) {
	    angle = -2.;
	}
	gravec_1.cosine = 99.99;

/*  ANGLE IS USED TO DECIDE IF P IS TO BE UPDATED AS CALCULATION */
/*        PROCEEDS. */

	ymaxst = .4;
	i__1 = -i;
	eps = (doublereal) pow_ii(&c__10, &i__1);
	tee = eps;
	print = i_indx(keywrd_1.keywrd, "LINMIN", 80L, 6L) != 0;
	icalcn = numcal_1.numcal;
    }
    fulscf = askful || gravec_1.cosine > angle;
    if (print) {
	s_wsfe(&io___15);
	do_fio(&c__1, (char *)&fulscf, (ftnlen)sizeof(logical));
	e_wsfe();
    }
    xmaxm = 0.;
    i__1 = *nvar;
    for (i = 1; i <= i__1; ++i) {
	pabs = (d__1 = pvect[i], abs(d__1));
/* L10: */
	xmaxm = max(xmaxm,pabs);
    }
    xminm = xmaxm;
    xmaxm = ymaxst / xmaxm;
    fin = *funct;
    ssqlst = *funct;
    iquit = 0;
    phi[0] = *funct;
    vt[0] = 0.;
    vt[1] = *step / 4.;
    if (vt[1] > xmaxm) {
	vt[1] = xmaxm;
    }
    fmax = *funct;
    fmin = *funct;
    *step = vt[1];
    i__1 = *nvar;
    for (i = 1; i <= i__1; ++i) {
/* L20: */
	xparam[i] += *step * pvect[i];
    }
    compfg_(&xparam[1], &c_true, &phi[1], &fulscf, &grad, &c_false);
    if (phi[1] > fmax) {
	fmax = phi[1];
    }
    if (phi[1] < fmin) {
	fmin = phi[1];
    }
    exchng_(&phi[1], &sqstor, &energy, &estor, &xparam[1], xstor, step, &alfs,
	     nvar);
    if (phi[0] <= phi[1]) {
	goto L30;
    }
    goto L40;
L30:
    vt[2] = -vt[1];
    left = 3;
    center = 1;
    right = 2;
    goto L50;
L40:
    vt[2] = vt[1] * 2.;
    left = 1;
    center = 2;
    right = 3;
L50:
    stlast = vt[2];
    *step = stlast - *step;
    i__1 = *nvar;
    for (i = 1; i <= i__1; ++i) {
/* L60: */
	xparam[i] += *step * pvect[i];
    }
    compfg_(&xparam[1], &c_true, funct, &fulscf, &grad, &c_false);
    if (*funct > fmax) {
	fmax = *funct;
    }
    if (*funct < fmin) {
	fmin = *funct;
    }
    if (*funct < sqstor) {
	exchng_(funct, &sqstor, &energy, &estor, &xparam[1], xstor, step, &
		alfs, nvar);
    }
    if (*funct < fin) {
	iquit = 1;
    }
    phi[2] = *funct;
    if (print) {
	s_wsfe(&io___35);
	do_fio(&c__1, (char *)&vt[0], (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&phi[0], (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&vt[1], (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&phi[1], (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&vt[2], (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&phi[2], (ftnlen)sizeof(doublereal));
	e_wsfe();
    }
    *okc = TRUE_;
    i__1 = maxlin;
    for (ictr = 3; ictr <= i__1; ++ictr) {
	alpha = vt[1] - vt[2];
	beta = vt[2] - vt[0];
	gamma = vt[0] - vt[1];
	if ((d__1 = alpha * beta * gamma, abs(d__1)) > 1e-20) {
	    alpha = -(phi[0] * alpha + phi[1] * beta + phi[2] * gamma) / (
		    alpha * beta * gamma);
	} else {
	    goto L190;
	}
	beta = (phi[0] - phi[1]) / gamma - alpha * (vt[0] + vt[1]);
	if (alpha <= 0.) {
	    goto L70;
	} else {
	    goto L100;
	}
L70:
	if (phi[right - 1] > phi[left - 1]) {
	    goto L80;
	}
	*step = vt[right - 1] * 3. - vt[center - 1] * 2.;
	goto L90;
L80:
	*step = vt[left - 1] * 3. - vt[center - 1] * 2.;
L90:
	s = *step - stlast;
	if (abs(s) > xmaxm) {
	    s = d_sign(&xmaxm, &s) * (xmaxm / s * .01f + 1);
	}
	*step = s + stlast;
	goto L110;
L100:
	*step = -beta / (alpha * 2.);
	s = *step - stlast;
	xxm = xmaxm * 2.;
	if (abs(s) > xxm) {
	    s = d_sign(&xxm, &s) * (xxm / s * .01f + 1);
	}
	*step = s + stlast;
L110:
	if (ictr <= 3) {
	    goto L120;
	}
	aabs = (d__1 = s * xminm, abs(d__1));
	if (aabs < xcrit) {
	    goto L190;
	}
L120:
	i__2 = *nvar;
	for (i = 1; i <= i__2; ++i) {
/* L130: */
	    xparam[i] += s * pvect[i];
	}
	funold = *funct;
	compfg_(&xparam[1], &c_true, funct, &fulscf, &grad, &c_false);
	if (*funct > fmax) {
	    fmax = *funct;
	}
	if (*funct < fmin) {
	    fmin = *funct;
	}
	if (*funct < sqstor) {
	    exchng_(funct, &sqstor, &energy, &estor, &xparam[1], xstor, step, 
		    &alfs, nvar);
	}
	if (*funct < fin) {
	    iquit = 1;
	}
	if (print) {
	    s_wsfe(&io___44);
	    do_fio(&c__1, (char *)&vt[left - 1], (ftnlen)sizeof(doublereal));
	    do_fio(&c__1, (char *)&phi[left - 1], (ftnlen)sizeof(doublereal));
	    do_fio(&c__1, (char *)&vt[center - 1], (ftnlen)sizeof(doublereal))
		    ;
	    do_fio(&c__1, (char *)&phi[center - 1], (ftnlen)sizeof(doublereal)
		    );
	    do_fio(&c__1, (char *)&vt[right - 1], (ftnlen)sizeof(doublereal));
	    do_fio(&c__1, (char *)&phi[right - 1], (ftnlen)sizeof(doublereal))
		    ;
	    do_fio(&c__1, (char *)&(*step), (ftnlen)sizeof(doublereal));
	    do_fio(&c__1, (char *)&(*funct), (ftnlen)sizeof(doublereal));
	    e_wsfe();
	}

/* TEST TO EXIT FROM LINMIN IF NOT DROPPING IN VALUE OF FUNCTION FAST.
 */

/* Computing MAX */
	d__1 = (ssqlst - fmin) * .2;
	tiny = max(d__1,drop);
	tiny = min(tiny,.5);
	if (print) {
	    s_wsfe(&io___46);
	    do_fio(&c__1, (char *)&tiny, (ftnlen)sizeof(doublereal));
	    e_wsfe();
	}
	if ((d__1 = funold - *funct, abs(d__1)) < tiny && iquit == 1) {
	    goto L190;
	}
	if ((d__1 = *step - stlast, abs(d__1)) <= eps * (d__2 = *step + 
		stlast, abs(d__2)) + tee && iquit == 1) {
	    goto L190;
	}
	stlast = *step;
	if (*step > vt[right - 1] || *step > vt[center - 1] && *funct < phi[
		center - 1] || *step > vt[left - 1] && *step < vt[center - 1] 
		&& *funct > phi[center - 1]) {
	    goto L140;
	}
	vt[right - 1] = *step;
	phi[right - 1] = *funct;
	goto L150;
L140:
	vt[left - 1] = *step;
	phi[left - 1] = *funct;
L150:
	if (vt[center - 1] < vt[right - 1]) {
	    goto L160;
	}
	i = center;
	center = right;
	right = i;
L160:
	if (vt[left - 1] < vt[center - 1]) {
	    goto L170;
	}
	i = left;
	left = center;
	center = i;
L170:
	if (vt[center - 1] < vt[right - 1]) {
	    goto L180;
	}
	i = center;
	center = right;
	right = i;
L180:
	;
    }
    *okc = FALSE_;
L190:
    exchng_(&sqstor, funct, &estor, &energy, xstor, &xparam[1], &alfs, step, 
	    nvar);
    *okf = *funct < ssqlst;
    if (*funct >= ssqlst) {
	return 0;
    }
    if (*step >= 0.) {
	goto L220;
    } else {
	goto L200;
    }
L200:
    *step = -(*step);
    i__1 = *nvar;
    for (i = 1; i <= i__1; ++i) {
/* L210: */
	pvect[i] = -pvect[i];
    }
L220:
    return 0;


} /* linmin_ */