/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:15 */
/*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 "dimql.h"
#include <stdlib.h>
void /*FUNCTION*/ dimql(
double *a,
long lda,
long n,
double eval[],
double work[],
long *ierr)
{
#define A(I_,J_)	(*(a+(I_)*(lda)+(J_)))
	long int i, j, k, l, m;
	double b, c, f, g, p, r, s, tst1, tst2;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	double *const Eval = &eval[0] - 1;
	double *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 DIMQL  Krogh  Changes to use M77CON
	 *>> 1992-04-24 DIMQL  CLL   Minor edits
	 *>> 1991-10-23 DIMQL  Krogh Initial version, converted from EISPACK.
	 *
	 *     This subroutine is a slight modification of the EISPACK subroutine
	 *     IMTQL2 (dated August 1983) which computes the eigenvalues and
	 *     eigenvectors of a real symmetric tridiagonal matrix.
	 *     IMTQL2 was in turn a translation of the ALGOL procedure, IMTQL2,
	 *     Num. Math. 12, 377-383(1968) by Martin and Wilkinson as modified
	 *     in Num. Math. 15, 450(1970) by Dubrelle.  Handbook for Auto.
	 *     Comp., vol.ii-linear algebra, 241-248(1971).
	 *
	 *     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.
	 *     EVAL contains the diagonal elements of the input matrix.
	 *     WORK contains the subdiagonal elements of the input matrix
	 *          in its last n-1 positions.  e(1) is arbitrary.
	 *     A    contains the transformation matrix produced in the
	 *          reduction to tridiagonal form.  if the eigenvectors
	 *          of the tridiagonal matrix are desired, A must contain
	 *          the identity matrix.
	 *
	 *     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.
	 *     ------------------------------------------------------------------
	 *--D replaces "?": ?IMQL
	 *     ------------------------------------------------------------------ */
 
	/*     ------------------------------------------------------------------
	 * */
	*ierr = 0;
 
	for (i = 2; i <= n; i++)
	{
		Work[i - 1] = Work[i];
	}
 
	Work[n] = 0.0e0;
 
	for (l = 1; l <= n; l++)
	{
		j = 0;
		/*     .......... look for small sub-diagonal element .......... */
L_605:
		for (m = l; m <= n; m++)
		{
			if (m == n)
				goto L_620;
			tst1 = fabs( Eval[m] ) + fabs( Eval[m + 1] );
			tst2 = tst1 + fabs( Work[m] );
			if (tst2 == tst1)
				goto L_620;
		}
 
L_620:
		p = Eval[l];
		if (m != l)
		{
			if (j == 30)
			{
				/*                   .......... set error -- no convergence to an
				 *                              eigenvalue after 30 iterations .......... */
				ermsg( "DIMQL", l, 0, "ERROR NO. is index of eigenvalue causing convergence failure."
				 , '.' );
				*ierr = l;
				return;
			}
			j += 1;
			/*        .......... form shift .......... */
			g = (Eval[l + 1] - p)/(2.0e0*Work[l]);
			if (fabs( g ) > 1.0e0)
			{
				r = g*sqrt( 1.0e0 + 1.0e0/SQ(g) );
			}
			else
			{
				r = sign( sqrt( 1.0e0 + SQ(g) ), g );
			}
			g = Eval[m] - p + Work[l]/(g + r);
			s = 1.0e0;
			c = 1.0e0;
			p = 0.0e0;
			for (i = m - 1; i >= l; i--)
			{
				f = s*Work[i];
				b = c*Work[i];
				if (fabs( f ) > fabs( g ))
				{
					r = fabs( f )*sqrt( 1.e0 + powi(g/f,2) );
				}
				else if (g != 0.0e0)
				{
					r = fabs( g )*sqrt( 1.e0 + powi(f/g,2) );
				}
				else
				{
					/*        .......... recover from underflow .......... */
					Work[i + 1] = 0.0e0;
					Eval[i + 1] -= p;
					Work[m] = 0.0e0;
					goto L_605;
				}
				Work[i + 1] = r;
				s = f/r;
				c = g/r;
				g = Eval[i + 1] - p;
				r = (Eval[i] - g)*s + 2.0e0*c*b;
				p = s*r;
				Eval[i + 1] = g + p;
				g = c*r - b;
				/*        .......... form vector .......... */
				for (k = 1; k <= n; k++)
				{
					f = A(i,k - 1);
					A(i,k - 1) = s*A(i - 1,k - 1) + c*f;
					A(i - 1,k - 1) = c*A(i - 1,k - 1) - s*f;
				}
 
			}
 
			Eval[l] -= p;
			Work[l] = g;
			Work[m] = 0.0e0;
			goto L_605;
		}
	}
	/*     .......... order eigenvalues and eigenvectors .......... */
	for (i = 1; i <= (n - 1); i++)
	{
		k = i;
		p = Eval[i];
 
		for (j = i + 1; j <= n; j++)
		{
			if (Eval[j] < p)
			{
				k = j;
				p = Eval[j];
			}
		}
 
		if (k != i)
		{
			Eval[k] = Eval[i];
			Eval[i] = p;
 
			for (j = 1; j <= n; j++)
			{
				p = A(i - 1,j - 1);
				A(i - 1,j - 1) = A(k - 1,j - 1);
				A(k - 1,j - 1) = p;
			}
		}
 
	}
 
	return;
#undef	A
} /* end of function */