/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:48 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "ssymql.h"
#include <stdlib.h>
void /*FUNCTION*/ ssymql(
float *a,
long lda,
long n,
float eval[],
float work[],
long *ierr)
{
#define A(I_,J_)	(*(a+(I_)*(lda)+(J_)))
	long int i, j, k, l;
	float f, g, h, hh, scale;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const Eval = &eval[0] - 1;
	float *const Work = &work[0] - 1;
		/* end of OFFSET VECTORS */
 
	/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
	 * ALL RIGHTS RESERVED.
	 * Based on Government Sponsored Research NAS7-03001.
	 *>> 1994-10-20 SSYMQL  Krogh  Changes to use M77CON
	 *>> 1992-04-24 SSYMQL  CLL  Minor edits.
	 *>> 1992-04-08 SSYMQL  Krogh Unused label 130 removed.
	 *>> 1991-10-23 SSYMQL  Krogh Initial version, converted from EISPACK.
	 *
	 *     This subroutine is a slight modification of the EISPACK subroutine
	 *     TRED2 (dated August 1983) which reduces a real matrix to a
	 *     symmetric tridiagonal matrix using and accumulating orthogonal
	 *     similarity transformations.  This subroutine then calls SIMQL to
	 *     compute the eigenvalues and eigenvectors of the matrix.
	 *     TRED2 was in turn a translation of the ALGOL procedure, TRED2,
	 *     Num. Math. 11, 181-195(1968) by Martin, Reinsch, and Wilkinson.
	 *
	 *     On input
	 *
	 *     LDA  must be set to the row dimension of the two-dimensional array
	 *          A as declared in the calling program dimension statement.
	 *
	 *     N    is the order of the matrix.
	 *
	 *     A    contains the real symmetric input matrix.  only the
	 *          lower triangle of the matrix need be supplied.
	 *
	 *     On output
	 *
	 *     A    contains the eigenvectors.
	 *
	 *     EVAL contains the eigenvalues in ascending order.  If an
	 *          error exit is made, the eigenvalues are correct but
	 *          unordered for indices 1,2,...,IERR-1.
	 *
	 *     WORK used for working storage.
	 *
	 *     IERR is set to:
	 *          zero  for normal return,
	 *          J     if the J-th eigenvalue has not been determined after 30
	 *                iterations.
	 *     ------------------------------------------------------------------
	 *--S replaces "?": ?SYMQL, ?IMQL
	 *     ------------------------------------------------------------------ */
 
	/*     ------------------------------------------------------------------
	 * */
	for (i = 1; i <= n; i++)
	{
		Eval[i] = A(i - 1,n - 1);
	}
 
	if (n <= 1)
	{
		if (n == 1)
		{
			A(0,0) = 1.0e0;
		}
		else
		{
			ierm1( "SSYMQL", -1, 0, "Require N > 0.", "N", n, '.' );
			*ierr = -1;
		}
		return;
	}
	for (i = n; i >= 2; i--)
	{
		l = i - 1;
		h = 0.0e0;
		scale = 0.0e0;
		if (l >= 2)
		{
			/*     .......... scale row (algol tol then not needed) .......... */
			for (k = 1; k <= l; k++)
			{
				scale += fabsf( Eval[k] );
			}
			if (scale != 0.0e0)
				goto L_140;
		}
 
		Work[i] = Eval[l];
 
		for (j = 1; j <= l; j++)
		{
			Eval[j] = A(j - 1,l - 1);
			A(j - 1,i - 1) = 0.0e0;
			A(i - 1,j - 1) = 0.0e0;
		}
 
		goto L_290;
 
L_140:
		for (k = 1; k <= l; k++)
		{
			Eval[k] /= scale;
			h += Eval[k]*Eval[k];
		}
 
		f = Eval[l];
		g = -signf( sqrtf( h ), f );
		Work[i] = scale*g;
		h += -f*g;
		Eval[l] = f - g;
		/*     .......... form A*u .......... */
		for (j = 1; j <= l; j++)
		{
			Work[j] = 0.0e0;
		}
 
		for (j = 1; j <= l; j++)
		{
			f = Eval[j];
			A(i - 1,j - 1) = f;
			g = Work[j] + A(j - 1,j - 1)*f;
 
			for (k = j + 1; k <= l; k++)
			{
				g += A(j - 1,k - 1)*Eval[k];
				Work[k] += A(j - 1,k - 1)*f;
			}
 
			Work[j] = g;
		}
		/*     .......... form P .......... */
		f = 0.0e0;
 
		for (j = 1; j <= l; j++)
		{
			Work[j] /= h;
			f += Work[j]*Eval[j];
		}
 
		hh = f/(h + h);
		/*     .......... form q .......... */
		for (j = 1; j <= l; j++)
		{
			Work[j] += -hh*Eval[j];
		}
		/*     .......... form reduced A .......... */
		for (j = 1; j <= l; j++)
		{
			f = Eval[j];
			g = Work[j];
 
			for (k = j; k <= l; k++)
			{
				A(j - 1,k - 1) += -f*Work[k] - g*Eval[k];
			}
 
			Eval[j] = A(j - 1,l - 1);
			A(j - 1,i - 1) = 0.0e0;
		}
 
L_290:
		Eval[i] = h;
	}
	/*     .......... accumulation of transformation matrices .......... */
	for (i = 2; i <= n; i++)
	{
		l = i - 1;
		A(l - 1,n - 1) = A(l - 1,l - 1);
		A(l - 1,l - 1) = 1.0e0;
		h = Eval[i];
		if (h != 0.0e0)
		{
 
			for (k = 1; k <= l; k++)
			{
				Eval[k] = A(i - 1,k - 1)/h;
			}
 
			for (j = 1; j <= l; j++)
			{
				g = 0.0e0;
 
				for (k = 1; k <= l; k++)
				{
					g += A(i - 1,k - 1)*A(j - 1,k - 1);
				}
 
				for (k = 1; k <= l; k++)
				{
					A(j - 1,k - 1) += -g*Eval[k];
				}
			}
		}
 
		for (k = 1; k <= l; k++)
		{
			A(i - 1,k - 1) = 0.0e0;
		}
 
	}
 
	for (i = 1; i <= n; i++)
	{
		Eval[i] = A(i - 1,n - 1);
		A(i - 1,n - 1) = 0.0e0;
	}
 
	A(n - 1,n - 1) = 1.0e0;
	simql( a, lda, n, eval, work, ierr );
	return;
#undef	A
} /* end of function */