/* hqrii.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_

/* Table of constant values */

static integer c__1 = 1;

/* Subroutine */ int hqrii_(doublereal *a, integer *n, integer *m, doublereal 
	*e, doublereal *v)
{
    /* System generated locals */
    integer v_dim1, v_offset, i__1, i__2, i__3;
    doublereal d__1, d__2, d__3, d__4;

    /* Builtin functions */
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
    /* Subroutine */ int s_stop(char *, ftnlen);
    double sqrt(doublereal), d_sign(doublereal *, doublereal *);

    /* Local variables */
    static integer iord;
    static doublereal seps, zero, sinv, summ, c, h;
    static integer i, j, k, l;
    static doublereal r, s, t, u, w[1290]	/* was [5][258] */, z;
    static integer itere;
    static doublereal ee, ff;
    static integer ig, ii;
    static doublereal ra, fn;
    static integer kk, ll;
    static doublereal rn, vn, gersch, ww;
    static integer im1, ip1, nm1, nm2, kp1;
    static doublereal sorter, del, eps, sum, eps1, eps2, eps3;

    /* Fortran I/O blocks */
    static cilist io___1 = { 0, 6, 0, "(////10X,'IN HQRII, N =',I4,' M =',I4)"
	    , 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 */

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

/* HQRII IS A DIAGONALISATION ROUTINE, WRITTEN BY YOSHITAKA BEPPU OF */
/*       NAGOYA UNIVERSITY, JAPAN. */
/*       FOR DETAILS SEE 'COMPUTERS & CHEMISTRY' VOL.6 1982. PAGE 000. */

/* ON INPUT    A       = MATRIX TO BE DIAGONALIZED */
/*             N       = SIZE OF MATRIX TO BE DIAGONALIZED. */
/*             M       = NUMBER OF EIGENVECTORS NEEDED. */
/*             E       = ARRAY OF SIZE AT LEAST N */
/*             V       = ARRAY OF SIZE AT LEAST NMAX*M */

/* ON OUTPUT   E       = EIGENVALUES */
/*             V       = EIGENVECTORS IN ARRAY OF SIZE NMAX*M */

/* ***********************************************************************
 */
    /* Parameter adjustments */
    --a;
    --e;
    v_dim1 = *n;
    v_offset = v_dim1 + 1;
    v -= v_offset;

    /* Function Body */
    if (*n <= 1 || *m <= 1 || *m > *n) {
	if (*n == 1 && *m == 1) {
	    e[1] = a[1];
	    v[v_dim1 + 1] = 1.;
	    return 0;
	}
	s_wsfe(&io___1);
	do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
	do_fio(&c__1, (char *)&(*m), (ftnlen)sizeof(integer));
	e_wsfe();
	s_stop("", 0L);
    }

/* EPS3 AND EPS ARE MACHINE-PRECISION DEPENDENT */

    eps3 = 1e-30;
    zero = 0.;
    ll = *n * (*n + 1) / 2 + 1;
    eps = 1e-8;
    iord = -1;
    nm1 = *n - 1;
    if (*n == 2) {
	goto L80;
    }
    nm2 = *n - 2;
/*     HOUSEHOLDER TRANSFORMATION */
    i__1 = nm2;
    for (k = 1; k <= i__1; ++k) {
	kp1 = k + 1;
	w[k * 5 - 4] = a[k * (k + 1) / 2];
	sum = 0.f;
	i__2 = *n;
	for (j = kp1; j <= i__2; ++j) {
	    w[j * 5 - 4] = a[j * (j - 1) / 2 + k];
/* L10: */
/* Computing 2nd power */
	    d__1 = w[j * 5 - 4];
	    sum = d__1 * d__1 + sum;
	}
	d__1 = sqrt(sum);
	s = d_sign(&d__1, &w[kp1 * 5 - 4]);
	w[k * 5 - 5] = -s;
	w[kp1 * 5 - 4] += s;
	a[k + kp1 * (kp1 - 1) / 2] = w[kp1 * 5 - 4];
	h = w[kp1 * 5 - 4] * s;
	if (abs(h) < 1e-35) {
	    goto L70;
	}
/* #      IF(H.EQ.0.D0) GOTO 70 */
	summ = 0.;
	i__2 = *n;
	for (i = kp1; i <= i__2; ++i) {
	    sum = 0.;
	    i__3 = i;
	    for (j = kp1; j <= i__3; ++j) {
/* L20: */
		sum += a[j + i * (i - 1) / 2] * w[j * 5 - 4];
	    }
	    if (i >= *n) {
		goto L40;
	    }
	    ip1 = i + 1;
	    i__3 = *n;
	    for (j = ip1; j <= i__3; ++j) {
/* L30: */
		sum += a[i + j * (j - 1) / 2] * w[j * 5 - 4];
	    }
L40:
	    w[i * 5 - 5] = sum / h;
/* L50: */
	    summ = w[i * 5 - 5] * w[i * 5 - 4] + summ;
	}
	u = summ * .5 / h;
	i__2 = *n;
	for (j = kp1; j <= i__2; ++j) {
	    w[j * 5 - 5] = w[j * 5 - 4] * u - w[j * 5 - 5];
	    i__3 = j;
	    for (i = kp1; i <= i__3; ++i) {
/* L60: */
		a[i + j * (j - 1) / 2] = w[i * 5 - 5] * w[j * 5 - 4] + w[j * 
			5 - 5] * w[i * 5 - 4] + a[i + j * (j - 1) / 2];
	    }
	}
L70:
	a[k * (k + 1) / 2] = h;
    }
L80:
    w[nm1 * 5 - 4] = a[nm1 * (nm1 + 1) / 2];
    w[*n * 5 - 4] = a[*n * (*n + 1) / 2];
    w[nm1 * 5 - 5] = a[nm1 + *n * (*n - 1) / 2];
    w[*n * 5 - 5] = 0.;
    gersch = abs(w[1]) + abs(w[0]);
    i__1 = nm1;
    for (i = 1; i <= i__1; ++i) {
/* L90: */
/* Computing MAX */
	d__4 = (d__1 = w[(i + 1) * 5 - 4], abs(d__1)) + (d__2 = w[i * 5 - 5], 
		abs(d__2)) + (d__3 = w[(i + 1) * 5 - 5], abs(d__3));
	gersch = max(d__4,gersch);
    }
    del = eps * gersch;
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	w[i * 5 - 3] = w[i * 5 - 5];
	e[i] = w[i * 5 - 4];
/* L100: */
	v[i + *m * v_dim1] = e[i];
    }
    if (del == zero) {
	goto L210;
    }
/*     QR-METHOD WITH ORIGIN SHIFT */
    k = *n;
L110:
    l = k;
L120:
    if ((d__1 = w[(l - 1) * 5 - 3], abs(d__1)) < del) {
	goto L130;
    }
    --l;
    if (l > 1) {
	goto L120;
    }
L130:
    if (l == k) {
	goto L160;
    }
    ww = (e[k - 1] + e[k]) * .5;
    r = e[k] - ww;
/* Computing 2nd power */
    d__2 = w[(k - 1) * 5 - 3];
    d__1 = sqrt(d__2 * d__2 + r * r);
    z = d_sign(&d__1, &r) + ww;
    ee = e[l] - z;
    e[l] = ee;
    ff = w[l * 5 - 3];
    r = sqrt(ee * ee + ff * ff);
    j = l;
    goto L150;
L140:
/* Computing 2nd power */
    d__1 = e[j];
/* Computing 2nd power */
    d__2 = w[j * 5 - 3];
    r = sqrt(d__1 * d__1 + d__2 * d__2);
    w[(j - 1) * 5 - 3] = s * r;
    ee = e[j] * c;
    ff = w[j * 5 - 3] * c;
L150:
    r += 1e-15;
    c = e[j] / r;
    s = w[j * 5 - 3] / r;
    ww = e[j + 1] - z;
    e[j] = (ff * c + ww * s) * s + ee + z;
    e[j + 1] = c * ww - s * ff;
    ++j;
    if (j < k) {
	goto L140;
    }
    w[(k - 1) * 5 - 3] = e[k] * s;
    e[k] = e[k] * c + z;
    goto L110;
L160:
    --k;
    if (k > 1) {
	goto L110;
    }
/*    *    *    *    *    *    *    *    *    *    *    *    * */

/*   AT THIS POINT THE ARRAY 'E' CONTAINS THE UN-ORDERED EIGENVALUES */

/*    *    *    *    *    *    *    *    *    *    *    *    * */
/*     STRAIGHT SELECTION SORT OF EIGENVALUES */
    sorter = 1.;
    if (iord < 0) {
	sorter = -1.;
    }
    j = *n;
L170:
    l = 1;
    ii = 1;
    ll = 1;
    i__1 = j;
    for (i = 2; i <= i__1; ++i) {
	if ((e[i] - e[l]) * sorter > 0.) {
	    goto L180;
	}
	l = i;
	goto L190;
L180:
	ii = i;
	ll = l;
L190:
	;
    }
    if (ii == ll) {
	goto L200;
    }
    ww = e[ll];
    e[ll] = e[ii];
    e[ii] = ww;
L200:
    j = ii - 1;
    if (j >= 2) {
	goto L170;
    }
L210:
    if (*m == 0) {
	return 0;
    }
/* ************** */
/*  ORDERING OF EIGENVALUES COMPLETE. */
/* ************** */
/*      INVERSE-ITERATION FOR EIGENVECTORS */
    fn = (real) (*n);
    eps1 = 1e-5;
    seps = sqrt(eps);
    eps2 = .05;
    rn = 0.;
    ra = eps * .6180339887485;
/*    0.618... IS THE FIBONACCI NUMBER (-1+SQRT(5))/2. */
    ig = 1;
    i__1 = *m;
    for (i = 1; i <= i__1; ++i) {
	im1 = i - 1;
	i__3 = *n;
	for (j = 1; j <= i__3; ++j) {
	    w[j * 5 - 3] = 0.;
	    w[j * 5 - 2] = w[j * 5 - 5];
	    w[j * 5 - 1] = v[j + *m * v_dim1] - e[i];
	    rn += ra;
	    if (rn >= eps) {
		rn -= eps;
	    }
/* L220: */
	    v[j + i * v_dim1] = rn;
	}
	i__3 = nm1;
	for (j = 1; j <= i__3; ++j) {
	    if ((d__1 = w[j * 5 - 1], abs(d__1)) >= (d__2 = w[j * 5 - 5], abs(
		    d__2))) {
		goto L230;
	    }
	    w[j * 5 - 4] = -w[j * 5 - 1] / w[j * 5 - 5];
	    w[j * 5 - 1] = w[j * 5 - 5];
	    t = w[(j + 1) * 5 - 1];
	    w[(j + 1) * 5 - 1] = w[j * 5 - 2];
	    w[j * 5 - 2] = t;
	    w[j * 5 - 3] = w[(j + 1) * 5 - 2];
	    if (w[j * 5 - 3] == zero) {
		w[j * 5 - 3] = del;
	    }
	    w[(j + 1) * 5 - 2] = 0.;
	    goto L240;
L230:
	    if (w[j * 5 - 1] == zero) {
		w[j * 5 - 1] = del;
	    }
	    w[j * 5 - 4] = -w[j * 5 - 5] / w[j * 5 - 1];
L240:
	    w[(j + 1) * 5 - 2] = w[j * 5 - 3] * w[j * 5 - 4] + w[(j + 1) * 5 
		    - 2];
/* L250: */
	    w[(j + 1) * 5 - 1] = w[j * 5 - 2] * w[j * 5 - 4] + w[(j + 1) * 5 
		    - 1];
	}
	if ((d__1 = w[*n * 5 - 1], abs(d__1)) < eps3) {
	    w[*n * 5 - 1] = del;
	}
	for (itere = 1; itere <= 5; ++itere) {
	    if (itere == 1) {
		goto L270;
	    }
	    i__3 = nm1;
	    for (j = 1; j <= i__3; ++j) {
		if (w[j * 5 - 3] == zero) {
		    goto L260;
		}
		t = v[j + i * v_dim1];
		v[j + i * v_dim1] = v[j + 1 + i * v_dim1];
		v[j + 1 + i * v_dim1] = t;
L260:
		v[j + 1 + i * v_dim1] = v[j + i * v_dim1] * w[j * 5 - 4] + v[
			j + 1 + i * v_dim1];
	    }
L270:
	    v[*n + i * v_dim1] /= w[*n * 5 - 1];
	    v[nm1 + i * v_dim1] = (v[nm1 + i * v_dim1] - v[*n + i * v_dim1] * 
		    w[nm1 * 5 - 2]) / w[nm1 * 5 - 1];
/* Computing MAX */
	    d__3 = (d__1 = v[*n + i * v_dim1], abs(d__1)), d__4 = (d__2 = v[
		    nm1 + i * v_dim1], abs(d__2)), d__3 = max(d__3,d__4);
	    vn = max(d__3,1e-20);
	    if (*n == 2) {
		goto L290;
	    }
	    k = nm2;
L280:
	    v[k + i * v_dim1] = (v[k + i * v_dim1] - v[k + 1 + i * v_dim1] * 
		    w[k * 5 - 2] - v[k + 2 + i * v_dim1] * w[k * 5 - 3]) / w[
		    k * 5 - 1];
/* Computing MAX */
	    d__2 = (d__1 = v[k + i * v_dim1], abs(d__1)), d__2 = max(d__2,vn);
	    vn = max(d__2,1e-20);
	    --k;
	    if (k >= 1) {
		goto L280;
	    }
L290:
	    s = eps1 / vn;
	    i__3 = *n;
	    for (j = 1; j <= i__3; ++j) {
/* L300: */
		v[j + i * v_dim1] *= s;
	    }
	    if (itere > 1 && vn > 1.) {
		goto L320;
	    }
/* L310: */
	}
/*     TRANSFORMATION OF EIGENVECTORS */
L320:
	if (*n == 2) {
	    goto L360;
	}
	i__3 = nm2;
	for (j = 1; j <= i__3; ++j) {
	    k = *n - j - 1;
	    if (a[k * (k + 1) / 2] == zero) {
		goto L350;
	    }
	    kp1 = k + 1;
	    sum = 0.;
	    i__2 = *n;
	    for (kk = kp1; kk <= i__2; ++kk) {
/* L330: */
		sum += a[k + kk * (kk - 1) / 2] * v[kk + i * v_dim1];
	    }
	    s = -sum / a[k * (k + 1) / 2];
	    i__2 = *n;
	    for (kk = kp1; kk <= i__2; ++kk) {
/* L340: */
		v[kk + i * v_dim1] = a[k + kk * (kk - 1) / 2] * s + v[kk + i *
			 v_dim1];
	    }
L350:
	    ;
	}
L360:
	i__3 = i;
	for (j = ig; j <= i__3; ++j) {
	    if ((d__1 = e[j] - e[i], abs(d__1)) < eps2) {
		goto L380;
	    }
/* L370: */
	}
	j = i;
L380:
	ig = j;
	if (ig == i) {
	    goto L410;
	}
/*     RE-ORTHOGONALISATION */
	i__3 = im1;
	for (k = ig; k <= i__3; ++k) {
	    sum = 0.;
	    i__2 = *n;
	    for (j = 1; j <= i__2; ++j) {
/* L390: */
		sum = v[j + k * v_dim1] * v[j + i * v_dim1] + sum;
	    }
	    s = -sum;
	    i__2 = *n;
	    for (j = 1; j <= i__2; ++j) {
/* L400: */
		v[j + i * v_dim1] = v[j + k * v_dim1] * s + v[j + i * v_dim1];
	    }
	}
/*     NORMALISATION */
L410:
	sum = 1e-24;
	i__2 = *n;
	for (j = 1; j <= i__2; ++j) {
/* L420: */
/* Computing 2nd power */
	    d__1 = v[j + i * v_dim1];
	    sum += d__1 * d__1;
	}
	sinv = 1. / sqrt(sum);
	i__2 = *n;
	for (j = 1; j <= i__2; ++j) {
/* L430: */
	    v[j + i * v_dim1] *= sinv;
	}
    }
    return 0;
} /* hqrii_ */