/* diag.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 fmo[46225], xdumy[20339];
} scrach_;

#define scrach_1 scrach_

/* Subroutine */ int diag_(doublereal *fao, doublereal *vector, integer *nocc,
	 doublereal *eig, integer *mdim, integer *n)
{
    /* Initialized data */

    static logical first = TRUE_;

    /* System generated locals */
    integer vector_dim1, vector_offset, i__1, i__2, i__3;
    doublereal d__1;

    /* Builtin functions */
    double sqrt(doublereal), d_sign(doublereal *, doublereal *);

    /* Local variables */
    static doublereal beta;
    static integer lumo;
    static doublereal tiny, a, b, c, d, e;
    static integer i, j, k, m;
    static doublereal alpha;
    static integer j1, k2, ij, kk;
    static doublereal bigeps, ws[215];
    extern /* Subroutine */ int epseta_(doublereal *, doublereal *);
    static doublereal eta, eps, sum;

/* 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 */

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

/*   "FAST" DIAGONALISATION PROCEDURE. */

/*    ON INPUT FAO CONTAINS THE LOWER HALF TRIANGLE OF THE MATRIX TO BE */
/*                         DIAGONALISED, PACKED. */
/*             VECTOR  CONTAINS THE OLD EIGENVECTORS ON INPUT, THE NEW */
/*             VECTORS ON EXITING. */
/*             NOCC = NUMBER OF OCCUPIED MOLECULAR ORBITALS. */
/*             EIG  = EIGENVALUES FROM AN EXACT DIAGONALISATION */
/*             MDIM = DECLARED SIZE OF MATRIX "C". */
/*             N = NUMBER OF ATOMIC ORBITALS IN BASIS SET */

/*  DIAG IS A PSEUDO-DIAGONALISATION PROCEDURE, IN THAT THE VECTORS THAT 
*/
/*       ARE GENERATED BY IT ARE MORE NEARLY ABLE TO BLOCK-DIAGONALISE */
/*       THE FOCK MATRIX OVER MOLECULAR ORBITALS THAN THE STARTING */
/*       VECTORS. IT MUST BE CONSIDERED PSEUDO FOR SEVERAL REASONS: */
/*       (A) IT DOES NOT GENERATE EIGENVECTORS - THE SECULAR DETERMINANT 
*/
/*           IS NOT DIAGONALISED, ONLY THE OCCUPIED-VIRTUAL INTERSECTION. 
*/
/*       (B) MANY SMALL ELEMENTS IN THE SEC.DET. ARE IGNORED AS BEING TOO 
*/
/*           SMALL COMPARED WITH THE LARGEST ELEMENT. */
/*       (C) WHEN ELEMENTS ARE ELIMINATED BY ROTATION, THE REST OF THE */
/*           SEC. DET. IS ASSUMED NOT TO CHANGE, I.E. ELEMENTS CREATED */
/*           ARE IGNORED. */
/*       (D) THE ROTATION REQUIRED TO ELIMINATE THOSE ELEMENTS CONSIDERED 
*/
/*           SIGNIFICANT IS APPROXIMATED TO USING THE EIGENVALUES OF THE 
*/
/*           EXACT DIAGONALISATION THROUGHOUT THE REST OF THE ITERATIVE */
/*           PROCEDURE. */

/*  (NOTE:- IN AN ITERATIVE PROCEDURE ALL THE APPROXIMATIONS PRESENT IN */
/*          DIAG BECOME VALID AT SELF-CONSISTENCY, SELF-CONSISTENCY IS */
/*          NOT SLOWED DOWN BY USE OF THESE APPROXIMATIONS) */

/*    REFERENCE: */
/*             "FAST SEMIEMPIRICAL CALCULATIONS", */
/*             STEWART. J.J.P., CSASZAR, P., PULAY, P., J. COMP. CHEM., */
/*             3, 227, (1982) */

/* ***********************************************************************
 */
/*             FMO  IS A WORK-SPACE OF SIZE (N-NOCC)*NOCC, IT WILL HOLD */
/*                  THE FOCK MOLECULAR ORBITAL INTERACTION MATRIX. */

/*  FIRST, CONSTRUCT THAT PART OF A SECULAR DETERMINANT OVER MOLECULAR */
/*  ORBITALS WHICH CONNECTS THE OCCUPIED AND VIRTUAL SETS. */

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

    /* Parameter adjustments */
    --fao;
    --eig;
    vector_dim1 = *mdim;
    vector_offset = vector_dim1 + 1;
    vector -= vector_offset;

    /* Function Body */
    if (first) {
	first = FALSE_;

/*   EPS IS THE SMALLEST NUMBER WHICH, WHEN ADDED TO 1.D0, IS NOT */
/*   EQUAL TO 1.D0 */
	epseta_(&eps, &eta);

/*   INCREASE EPS TO ALLOW FOR A LOT OF ROUND-OFF */

	bigeps = sqrt(eps) * 10.;
    }
    tiny = 0.;
    lumo = *nocc + 1;
    ij = 0;
    i__1 = *n;
    for (i = lumo; i <= i__1; ++i) {
	kk = 0;
	i__2 = *n;
	for (j = 1; j <= i__2; ++j) {
	    sum = 0.;
	    i__3 = j;
	    for (k = 1; k <= i__3; ++k) {
		++kk;
/* L10: */
		sum += fao[kk] * vector[k + i * vector_dim1];
	    }
	    if (j == *n) {
		goto L30;
	    }
	    j1 = j + 1;
	    k2 = kk;
	    i__3 = *n;
	    for (k = j1; k <= i__3; ++k) {
		k2 = k2 + k - 1;
/* L20: */
		sum += fao[k2] * vector[k + i * vector_dim1];
	    }
L30:
	    ws[j - 1] = sum;
	}
	i__2 = *nocc;
	for (j = 1; j <= i__2; ++j) {
	    ++ij;
	    sum = 0.;
	    i__3 = *n;
	    for (k = 1; k <= i__3; ++k) {
/* L40: */
		sum += ws[k - 1] * vector[k + j * vector_dim1];
	    }
	    if (tiny < abs(sum)) {
		tiny = abs(sum);
	    }
/* L50: */
	    scrach_1.fmo[ij - 1] = sum;
	}
/* L60: */
    }
    tiny *= .04;
/* ***********************************************************************
 */

/*   NOW DO A CRUDE 2 BY 2 ROTATION TO "ELIMINATE" SIGNIFICANT ELEMENTS */

/* ***********************************************************************
 */
    ij = 0;
    i__1 = *n;
    for (i = lumo; i <= i__1; ++i) {
	i__2 = *nocc;
	for (j = 1; j <= i__2; ++j) {
	    ++ij;
	    if ((d__1 = scrach_1.fmo[ij - 1], abs(d__1)) < tiny) {
		goto L80;
	    }

/*      BEGIN 2 X 2 ROTATIONS */

	    a = eig[j];
	    b = eig[i];
	    c = scrach_1.fmo[ij - 1];
	    d = a - b;

/*    USE BIGEPS TO DETERMINE WHETHER TO DO A 2 BY 2 ROTATION */

	    if ((d__1 = c / d, abs(d__1)) < bigeps) {
		goto L80;
	    }

/*  AT THIS POINT WE KNOW THAT */
	    d__1 = sqrt(c * 4. * c + d * d);
	    e = d_sign(&d__1, &d);
	    alpha = sqrt((d / e + 1.) * .5);
	    d__1 = sqrt(1. - alpha * alpha);
	    beta = -d_sign(&d__1, &c);

/*      ROTATION OF PSEUDO-EIGENVECTORS */

	    i__3 = *n;
	    for (m = 1; m <= i__3; ++m) {
		a = vector[m + j * vector_dim1];
		b = vector[m + i * vector_dim1];
		vector[m + j * vector_dim1] = alpha * a + beta * b;
		vector[m + i * vector_dim1] = alpha * b - beta * a;
/* L70: */
	    }
L80:
	    ;
	}
/* L90: */
    }
    return 0;
} /* diag_ */