/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:46 */
/*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 "sherql.h"
#include <stdlib.h>
void /*FUNCTION*/ sherql(
float *ar,
float *ai,
long lda,
long n,
float eval[],
float *vr,
float *vi,
float work[],
long *ierr)
{
#define AR(I_,J_) (*(ar+(I_)*(lda)+(J_)))
#define AI(I_,J_) (*(ai+(I_)*(lda)+(J_)))
#define VR(I_,J_) (*(vr+(I_)*(lda)+(J_)))
#define VI(I_,J_) (*(vi+(I_)*(lda)+(J_)))
long int i, j, k, l, n1, n2;
float f, fi, g, gi, h, hh, s, scale, si;
/* 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 SHERQL Krogh Changes to use M77CON
*>> 1994-04-20 SHERQL CLL Edited comment to make DP & SP files similar.
*>> 1992-04-23 CLL Declaring all variables.
*>> 1992-04-08 SHERQL Removed unused label 220.
*>> 1991-10-23 SHERQL Krogh Initial version, converted from EISPACK.
*
* This subroutine reduces a complex hermitian matrix to a real
* symmetric tridiagonal matrix using unitary similarity
* transformations, computes the eigenvalues and eigenvectors of the
* tridiagonal matrix, and then computes the eigenvectors of the
* original matrix.
*
* This subroutine is a concatenation of a slight modification of the
* EISPACK subroutine HTRIDI, a call to SIMQL which is a slight
* modification of the EISPACK subroutine IMTQL2, and a slight
* modification of the EISPACK subroutine HTRIBK. These subroutines
* in turn are a translation of a complex analogue of the ALGOL
* procedure TRED1, Num. Math. 11, 181-195(1968) by Martin, Reinsch,
* and Wilkinson. Handbook for Auto. Comp., Vol.ii-Linear Algebra,
* 212-226(1971), and the complex analogue of the ALGOL procedure
* TRBAK1, Num. Math. 11, 181-195(1968) by Martin, Reinsch, and
* Wilkinson. Handbook for Auto. Comp., vol.ii-Linear Algebra,
* 212-226(1971).
* EISPACK routines are dated August 1983.
*
* On input
*
* AR and AI contain the real and imaginary parts,
* respectively, of the complex hermitian input matrix.
* only the lower triangle of the matrix need be supplied.
*
* LDA must be set to the row dimension of two-dimensional
* array parameters as declared in the calling program
* dimension statement.
*
* N is the order of the matrix.
*
* On output
*
* AR and AI contain information about the unitary transformations
* used in the reduction in their full lower triangles. Their
* strict upper triangles and the diagonal of AR are unaltered.
*
* 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, and no eigenvectors are
* computed.
*
* VR and VI contain the the real and imaginary parts of the
* eigenvectors with column k corresponding to the eigenvalue in
* EVAL(k). Note that the last component of each returned vector
* is real and that vector euclidean norms are preserved.
*
* 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 "?": ?HERQL, ?IMQL
* Both versions use IERM1
* ------------------------------------------------------------------ */
/* ------------------------------------------------------------------ */
n1 = n;
n2 = n1 + n1;
Work[n2] = 1.0e0;
Work[n2 + n1] = 0.0e0;
for (i = 1; i <= n1; i++)
{
Eval[i] = AR(i - 1,i - 1);
}
if (n1 <= 1)
{
if (n1 == 1)
{
AR(0,0) = 1.0e0;
AI(0,0) = 0.0e0;
}
else
{
ierm1( "SHERQL", -1, 0, "Require N > 0.", "N", n1, '.' );
*ierr = -1;
}
return;
}
for (i = n1; i >= 1; i--)
{
l = i - 1;
h = 0.0e0;
scale = 0.0e0;
if (l < 1)
goto L_130;
/* .......... scale row (algol tol then not needed) .......... */
for (k = 1; k <= l; k++)
{
scale += fabsf( AR(k - 1,i - 1) ) + fabsf( AI(k - 1,i - 1) );
}
if (scale != 0.0e0)
goto L_140;
Work[n1 + l] = 1.0e0;
Work[n2 + l] = 0.0e0;
L_130:
Work[i] = 0.0e0;
goto L_290;
L_140:
for (k = 1; k <= l; k++)
{
AR(k - 1,i - 1) /= scale;
AI(k - 1,i - 1) /= scale;
h += AR(k - 1,i - 1)*AR(k - 1,i - 1) + AI(k - 1,i - 1)*
AI(k - 1,i - 1);
}
g = sqrtf( h );
Work[i] = scale*g;
f = fabsf( AR(l - 1,i - 1) );
si = fabsf( AI(l - 1,i - 1) );
if (f > si)
{
f *= sqrtf( 1.0e0 + powif(si/f,2) );
}
else if (si != 0.0e0)
{
f = si*sqrtf( 1.0e0 + powif(f/si,2) );
}
else
{
Work[n1 + l] = -Work[n1 + i];
si = Work[n2 + i];
AR(l - 1,i - 1) = g;
goto L_170;
}
/* .......... Form next diagonal element of matrix t .......... */
Work[n1 + l] = (AI(l - 1,i - 1)*Work[n2 + i] - AR(l - 1,i - 1)*
Work[n1 + i])/f;
si = (AR(l - 1,i - 1)*Work[n2 + i] + AI(l - 1,i - 1)*Work[n1 + i])/
f;
h += f*g;
g = 1.0e0 + g/f;
AR(l - 1,i - 1) *= g;
AI(l - 1,i - 1) *= g;
if (l == 1)
goto L_270;
f = 0.0e0;
L_170:
for (j = 1; j <= l; j++)
{
g = 0.0e0;
gi = 0.0e0;
/* .......... form element of a*u .......... */
for (k = 1; k <= j; k++)
{
g += AR(k - 1,j - 1)*AR(k - 1,i - 1) + AI(k - 1,j - 1)*
AI(k - 1,i - 1);
gi += -AR(k - 1,j - 1)*AI(k - 1,i - 1) + AI(k - 1,j - 1)*
AR(k - 1,i - 1);
}
for (k = j + 1; k <= l; k++)
{
g += AR(j - 1,k - 1)*AR(k - 1,i - 1) - AI(j - 1,k - 1)*
AI(k - 1,i - 1);
gi += -AR(j - 1,k - 1)*AI(k - 1,i - 1) - AI(j - 1,k - 1)*
AR(k - 1,i - 1);
}
/* .......... form element of p .......... */
Work[j] = g/h;
Work[n2 + j] = gi/h;
f += Work[j]*AR(j - 1,i - 1) - Work[n2 + j]*AI(j - 1,i - 1);
}
hh = f/(h + h);
/* .......... form reduced a .......... */
for (j = 1; j <= l; j++)
{
f = AR(j - 1,i - 1);
g = Work[j] - hh*f;
Work[j] = g;
fi = -AI(j - 1,i - 1);
gi = Work[n2 + j] - hh*fi;
Work[n2 + j] = -gi;
for (k = 1; k <= j; k++)
{
AR(k - 1,j - 1) += -f*Work[k] - g*AR(k - 1,i - 1) +
fi*Work[n2 + k] + gi*AI(k - 1,i - 1);
AI(k - 1,j - 1) += -f*Work[n2 + k] - g*AI(k - 1,i - 1) -
fi*Work[k] - gi*AR(k - 1,i - 1);
}
}
L_270:
for (k = 1; k <= l; k++)
{
AR(k - 1,i - 1) *= scale;
AI(k - 1,i - 1) *= scale;
}
Work[n2 + l] = -si;
L_290:
hh = Eval[i];
Eval[i] = AR(i - 1,i - 1);
AR(i - 1,i - 1) = hh;
AI(i - 1,i - 1) = scale*sqrtf( h );
}
/* ------------------------------------------------------------------
* */
for (i = 1; i <= n1; i++)
{
for (j = 1; j <= n1; j++)
{
VR(j - 1,i - 1) = 0.0e0;
}
VR(i - 1,i - 1) = 1.0e0;
}
simql( vr, lda, n1, eval, work, ierr );
if (*ierr != 0)
return;
/* ------------------------------------------------------------------
*
* .......... transform the eigenvectors of the real symmetric
* tridiagonal matrix to those of the hermitian
* tridiagonal matrix. .......... */
for (k = 1; k <= n1; k++)
{
for (j = 1; j <= n1; j++)
{
VI(j - 1,k - 1) = -VR(j - 1,k - 1)*Work[n2 + k];
VR(j - 1,k - 1) *= Work[n1 + k];
}
}
/* .......... Recover and apply the Householder matrices .......... */
for (i = 2; i <= n1; i++)
{
l = i - 1;
h = AI(i - 1,i - 1);
if (h != 0.0e0)
{
for (j = 1; j <= n1; j++)
{
s = 0.0e0;
si = 0.0e0;
for (k = 1; k <= l; k++)
{
s += AR(k - 1,i - 1)*VR(j - 1,k - 1) - AI(k - 1,i - 1)*
VI(j - 1,k - 1);
si += AR(k - 1,i - 1)*VI(j - 1,k - 1) + AI(k - 1,i - 1)*
VR(j - 1,k - 1);
}
/* .......... double divisions avoid possible underflow .......... */
s = (s/h)/h;
si = (si/h)/h;
for (k = 1; k <= l; k++)
{
VR(j - 1,k - 1) += -s*AR(k - 1,i - 1) - si*AI(k - 1,i - 1);
VI(j - 1,k - 1) += -si*AR(k - 1,i - 1) + s*AI(k - 1,i - 1);
}
}
}
}
return;
#undef VI
#undef VR
#undef AI
#undef AR
} /* end of function */