/*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_devvun s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#include "p_devvun.h"
/*     program DRDEVVUN
 *>> 1996-05-28 DRDEVVUN  Krogh Added external statement.
 *>> 1994-10-19 DRDEVVUN  Krogh  Changes to use M77CON
 *>> 1994-09-22 DRDEVVUN  CLL
 *>> 1992-04-23 CLL
 *>> 1992-03-04 DRDEVVUN  Krogh Initial version.
 *     Demonstrate unsymmetric eigenvalue/eigenvector subroutine DEVVUN.
 *     ------------------------------------------------------------------
 *--D replaces "?": DR?EVVUN, ?EVVUN, ?VECP, ?MATP, ?DOT
 *     ------------------------------------------------------------------ */
		/* PARAMETER translations */
#define	LDA	3
		/* end of PARAMETER translations */
 
 
int main( )
{
	long int i, iflag[LDA], j, _i, _r;
	double asav[LDA][LDA], d[LDA][LDA], vec[LDA][LDA], vi[LDA], vr[LDA],
	 work[LDA];
	static double a[LDA][LDA];
	static double anorm = 30.0e0;
	static long n = LDA;
	static int _aini = 1;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	long *const Iflag = &iflag[0] - 1;
	double *const Vi = &vi[0] - 1;
	double *const Vr = &vr[0] - 1;
	double *const Work = &work[0] - 1;
		/* end of OFFSET VECTORS */
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp0[] = {8.0e0,-1.0e0,-5.0e0};
		for (i = 1, _r = 0; i <= 3; i++)
		{
			a[i - 1][0] = _itmp0[_r++];
			}
		}
		{ static double _itmp1[] = {-4.0e0,4.0e0,-2.0e0};
		for (i = 1, _r = 0; i <= 3; i++)
		{
			a[i - 1][1] = _itmp1[_r++];
			}
		}
		{ static double _itmp2[] = {18.0e0,-5.0e0,-7.0e0};
		for (i = 1, _r = 0; i <= 3; i++)
		{
			a[i - 1][2] = _itmp2[_r++];
			}
		}
		_aini = 0;
	}
 
	/*     ------------------------------------------------------------------ */
	printf("DRDEVVUN..  Demo driver for DEVVUN.\n");
 
	/*     First copy A() to ASAV() for later residual check.
	 * */
	for (j = 1; j <= n; j++)
	{
		for (i = 1; i <= n; i++)
		{
			asav[j - 1][i - 1] = a[j - 1][i - 1];
		}
	}
 
	devvun( &a[0][0], LDA, n, vr, vi, (double*)vec, iflag, work );
 
	printf(" IFLAG(1) =%2ld\n", Iflag[1]);
	if (Iflag[1] <= 2)
	{
		dvecp( vr, n, " Real part of the eigenvalues" );
		dvecp( vi, n, " Imaginary part of the eigenvalues" );
		dmatp( (double*)vec, LDA, n, n, "0 Eigenvectors as columns or pairs of columns"
		  );
 
		/*           As a check compute D = (ASAV*VEC - VEC*EVAL) / ANORM.
		 *           Expect D to be close to the machine precision.
		 * */
		for (j = 1; j <= n; j++)
		{
			if (Vi[j] == 0.0e0)
			{
				/*                 Compute residual for a real eigenvalue and eigenvector */
				for (i = 1; i <= n; i++)
				{
					d[j - 1][i - 1] = (ddot( n, &asav[0][i - 1], LDA,
					 &vec[j - 1][0], 1 ) - vec[j - 1][i - 1]*Vr[j])/
					 anorm;
				}
			}
			else if (Vi[j] > 0.0e0)
			{
				for (i = 1; i <= n; i++)
				{
					d[j - 1][i - 1] = (ddot( n, &asav[0][i - 1], LDA,
					 &vec[j - 1][0], 1 ) - vec[j - 1][i - 1]*Vr[j] +
					 vec[j][i - 1]*Vi[j])/anorm;
					d[j][i - 1] = (ddot( n, &asav[0][i - 1], LDA,
					 &vec[j][0], 1 ) - vec[j - 1][i - 1]*Vi[j] - vec[j][i - 1]*
					 Vr[j])/anorm;
				}
			}
		}
		dmatp( (double*)d, LDA, n, n, "0 Packed residual matrix D = (A*EVEC - EVEC*EVAL) / ANORM"
		  );
	}
	else
	{
		printf("\n Failure in DEVVUN. \n");
	}
	exit(0);
} /* end of function */