/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:10 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_dherql s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#include "p_dherql.h"
/*     program DRDHERQL
 *>> 1996-05-28 DRDHERQL  Krogh Added external statement.
 *>> 1994-10-19 DRDHERQL  Krogh  Changes to use M77CON
 *>> 1994-09-23 DRDHERQL  CLL
 *>> 1992-04-23 CLL
 *>> 1992-03-04 DRDHERQL  Krogh Initial version.
 *     Demonstrate Hermitian eigenvalue/eigenvector subroutine DHERQL.
 *     ------------------------------------------------------------------
 *--D replaces "?": DR?HERQL, ?HERQL, ?VECP, ?MATP, ?DOT
 *     ------------------------------------------------------------------ */
		/* PARAMETER translations */
#define	LDA	4
#define	LDA3	(3*LDA)
		/* end of PARAMETER translations */
 
 
int main( )
{
	long int i, ierr, j, _i, _r;
	double aisav[LDA][LDA], arsav[LDA][LDA], di[LDA][LDA], dr[LDA][LDA],
	 eval[LDA], vi[LDA][LDA], vr[LDA][LDA], work[LDA3];
	static double ai[LDA][LDA], ar[LDA][LDA];
	static double anorm = 46.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! */
		ar[0][0] = 25.0e0;
		{ static double _itmp0[] = {-3.0e0,25.0e0};
		for (i = 1, _r = 0; i <= 2; i++)
		{
			ar[i - 1][1] = _itmp0[_r++];
			}
		}
		{ static double _itmp1[] = {-8.0e0,0.0e0,25.0e0};
		for (i = 1, _r = 0; i <= 3; i++)
		{
			ar[i - 1][2] = _itmp1[_r++];
			}
		}
		{ static double _itmp2[] = {0.0e0,-8.0e0,3.0e0,25.0e0};
		for (i = 1, _r = 0; i <= 4; i++)
		{
			ar[i - 1][3] = _itmp2[_r++];
			}
		}
		ai[0][0] = 0.0e0;
		{ static double _itmp3[] = {4.0e0,0.0e0};
		for (i = 1, _r = 0; i <= 2; i++)
		{
			ai[i - 1][1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {-6.0e0,0.0e0,0.0e0};
		for (i = 1, _r = 0; i <= 3; i++)
		{
			ai[i - 1][2] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {0.0e0,6.0e0,4.0e0,0.0e0};
		for (i = 1, _r = 0; i <= 4; i++)
		{
			ai[i - 1][3] = _itmp5[_r++];
			}
		}
		_aini = 0;
	}
 
	/*     ------------------------------------------------------------------ */
	printf("DRDHERQL..  Demo driver for DHERQL.\n");
 
	/*     First copy AR() and AI() to ARSAV() and AISAV() for later
	 *     residual check.
	 * */
	for (i = 1; i <= n; i++)
	{
		for (j = 1; j <= i; j++)
		{
			arsav[j - 1][i - 1] = ar[j - 1][i - 1];
			arsav[i - 1][j - 1] = arsav[j - 1][i - 1];
			aisav[j - 1][i - 1] = ai[j - 1][i - 1];
			aisav[i - 1][j - 1] = -aisav[j - 1][i - 1];
		}
	}
 
	dherql( (double*)ar, (double*)ai, LDA, n, eval, (double*)vr, (double*)vi,
	 work, &ierr );
	if (ierr == 0)
	{
		dvecp( eval, n, "0 Eigenvalues" );
		dmatp( (double*)vr, LDA, n, n, "0 Real parts of eigenvectors as column vectors"
		  );
		dmatp( (double*)vi, LDA, n, n, "0 Imaginary parts of eigenvectors as column vectors"
		  );
 
		/*        As a check compute D = (ASAV*EVEC - EVEC*EVAL) / ANORM.
		 *        Expect D to be close to the machine precision.
		 * */
		for (j = 1; j <= LDA; j++)
		{
			for (i = 1; i <= LDA; i++)
			{
				dr[j - 1][i - 1] = (ddot( LDA, &arsav[0][i - 1], LDA,
				 &vr[j - 1][0], 1 ) - ddot( LDA, &aisav[0][i - 1],
				 LDA, &vi[j - 1][0], 1 ) - vr[j - 1][i - 1]*Eval[j])/
				 anorm;
				di[j - 1][i - 1] = (ddot( LDA, &arsav[0][i - 1], LDA,
				 &vi[j - 1][0], 1 ) + ddot( LDA, &aisav[0][i - 1],
				 LDA, &vr[j - 1][0], 1 ) - vi[j - 1][i - 1]*Eval[j])/
				 anorm;
			}
		}
		dmatp( (double*)dr, LDA, n, n, "0Real part of residual matrix D = (A*EVEC - EVEC*EVAL) / ANORM"
		  );
		dmatp( (double*)di, LDA, n, n, "0Imag part of residual matrix D = (A*EVEC - EVEC*EVAL) / ANORM"
		  );
	}
	else
	{
		printf("\n Convergence failure in DHERQL, IERR =%5ld\n", ierr);
	}
	exit(0);
} /* end of function */