/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:41 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dherql.h" #include void /*FUNCTION*/ dherql( double *ar, double *ai, long lda, long n, double eval[], double *vr, double *vi, double 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; double f, fi, g, gi, h, hh, s, scale, si; /* 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 DHERQL Krogh Changes to use M77CON *>> 1994-04-20 DHERQL CLL Edited comment to make DP & SP files similar. *>> 1992-04-23 CLL Declaring all variables. *>> 1992-04-08 DHERQL Removed unused label 220. *>> 1991-10-23 DHERQL 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 DIMQL 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. * * ------------------------------------------------------------------ *--D 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( "DHERQL", -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 += fabs( AR(k - 1,i - 1) ) + fabs( 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 = sqrt( h ); Work[i] = scale*g; f = fabs( AR(l - 1,i - 1) ); si = fabs( AI(l - 1,i - 1) ); if (f > si) { f *= sqrt( 1.0e0 + powi(si/f,2) ); } else if (si != 0.0e0) { f = si*sqrt( 1.0e0 + powi(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*sqrt( 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; } dimql( 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 */