/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:13 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_dsymql s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#include "p_dsymql.h"
/*     program DRDSYMQL
 *>> 1996-05-28 DRDSYMQL  Krogh Added external statement.
 *>> 1994-10-19 DRDSYMQL  Krogh  Changes to use M77CON
 *>> 1994-09-22 DRDSYMQL  CLL
 *>> 1992-04-23 CLL
 *>> 1992-03-04 DRDSYMQL  Krogh Initial version.
 *     Demonstrate symmetric eigenvalue/eigenvector subroutine DSYMQL.
 *     ------------------------------------------------------------------
 *--D replaces "?": DR?SYMQL, ?SYMQL, ?VECP, ?MATP, ?DOT
 *     ------------------------------------------------------------------ */
		/* PARAMETER translations */
#define	LDA	4
		/* end of PARAMETER translations */
 
 
int main( )
{
	long int i, ierr, j, _i, _r;
	double asav[LDA][LDA], d[LDA][LDA], eval[LDA], work[LDA];
	static double a[LDA][LDA];
	static double anorm = 11.0e0;
	static long n = LDA;
	static int _aini = 1;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	double *const Eval = &eval[0] - 1;
	double *const Work = &work[0] - 1;
		/* end of OFFSET VECTORS */
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		a[0][0] = 5.0e0;
		{ static double _itmp0[] = {4.0e0,5.0e0};
		for (j = 1, _r = 0; j <= 2; j++)
		{
			a[j - 1][1] = _itmp0[_r++];
			}
		}
		{ static double _itmp1[] = {1.0e0,1.0e0,4.0e0};
		for (j = 1, _r = 0; j <= 3; j++)
		{
			a[j - 1][2] = _itmp1[_r++];
			}
		}
		{ static double _itmp2[] = {1.0e0,1.0e0,2.0e0,4.0e0};
		for (j = 1, _r = 0; j <= 4; j++)
		{
			a[j - 1][3] = _itmp2[_r++];
			}
		}
		_aini = 0;
	}
 
	/*     ------------------------------------------------------------------ */
	printf("DRDSYMQL..  Demo driver for DSYMQL.\n");
 
	/*     First copy A() to ASAV() for later residual check.
	 * */
	for (i = 1; i <= n; i++)
	{
		for (j = 1; j <= i; j++)
		{
			asav[j - 1][i - 1] = a[j - 1][i - 1];
			asav[i - 1][j - 1] = asav[j - 1][i - 1];
		}
	}
	dsymql( (double*)a, LDA, n, eval, work, &ierr );
	if (ierr == 0)
	{
		dvecp( eval, n, "0 Eigenvalues" );
		dmatp( (double*)a, LDA, n, n, "0 Eigenvectors as column vectors"
		  );
 
		/*        As a check compute D = (ASAV*EVEC - EVEC*EVAL) / ANORM.
		 *        The EVEC's are in the array A().
		 *        Expect D to be close to machine precision.
		 * */
		for (j = 1; j <= n; j++)
		{
			for (i = 1; i <= n; i++)
			{
				d[j - 1][i - 1] = (ddot( n, &asav[0][i - 1], LDA,
				 &a[j - 1][0], 1 ) - a[j - 1][i - 1]*Eval[j])/anorm;
			}
		}
		dmatp( (double*)d, LDA, n, n, "0 Residual matrix D = (A*EVEC - EVEC*EVAL) / ANORM"
		  );
 
	}
	else
	{
		printf("\n Convergence failure in DSYMQL, IERR =%5ld\n", ierr);
	}
	exit(0);
} /* end of function */