/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:19 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "sevun.h" #include void /*FUNCTION*/ sevun( float *a, long lda, long n, float evalr[], float evali[], long iflag[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) LOGICAL32 notlas; long int en, enm2, i, igh, itn, its, j, k, l, low, ltype, m, mp2, na; float norm, p, q, r, s, t, tst1, tst2, w, x, y, zz; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Evali = &evali[0] - 1; float *const Evalr = &evalr[0] - 1; long *const Iflag = &iflag[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. *>> 1996-03-30 SEVUN Krogh MIN0 => MIN *>> 1994-10-20 SEVUN Krogh Changes to use M77CON *>> 1992-04-24 SEVUN CLL Minor edits. *>> 1992-04-23 SEVUN Krogh Made DP version compatible with SP version. *>> 1991-10-25 SEVUN Krogh Initial version, converted from EISPACK. * * This subroutine uses slight modifcations of EISPACK routines * BALANC, ELMHES, and HQR to get the eigenvalues of a general real * matrix by the QR method. The first two are are encapsulated in * SEVBH. HQR, which is inline here, is a translation of the ALGOL * procedure HQR, Num. Math. 14, 219-231(1970) by Martin, Peters, and * Wilkinson. * Handbook for Auto. Comp., Vol.ii-Linear Algebra, 359-371(1971). * * On input * A contains the input matrix whose eigenvalues are desired. * 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 * A has been destroyed. * EVALR and EVALI contain real and imaginary parts, respectively, of * the eigenvalues. The eigenvalues are given in order of * increasing real parts. When real parts are equal they are * given in order of increasing absolute complex part. Complex * conjugate pairs of values appear consecutively with * the eigenvalue having the positive imaginary part first. If * an error exit is made, the eigenvalues should be correct * (but unordered) for indices IFLAG(1)+1,...,N. * IFLAG(1) is set to * 1 If all eigenvalues are real. * 2 If some eigenvalues are complex. * 3 If N < 1 on the initial entry. * 4 If the limit of 30*N iterations is exhausted. * ------------------------------------------------------------------ *--S replaces "?": ?EVUN, ?EVBH * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ * */ sevbh( a, lda, n, &low, &igh, iflag, evalr ); ltype = 1; norm = 0.0e0; k = 1; /* .......... store roots isolated by SEVBH * and compute matrix norm .......... */ for (i = 1; i <= n; i++) { for (j = k; j <= n; j++) { norm += fabsf( A(j - 1,i - 1) ); } k = i; if ((i < low) || (i > igh)) { Evalr[i] = A(i - 1,i - 1); Evali[i] = 0.0e0; } } en = igh; t = 0.0e0; itn = 30*n; /* .......... search for next eigenvalues .......... */ L_60: if (en < low) goto L_300; its = 0; na = en - 1; enm2 = na - 1; /* .......... look for single small sub-diagonal element */ L_70: for (l = en; l >= (low + 1); l--) { s = fabsf( A(l - 2,l - 2) ) + fabsf( A(l - 1,l - 1) ); if (s == 0.0e0) s = norm; tst1 = s; tst2 = tst1 + fabsf( A(l - 2,l - 1) ); if (tst2 == tst1) goto L_100; } /* .......... form shift .......... */ L_100: x = A(en - 1,en - 1); if (l == en) goto L_270; y = A(na - 1,na - 1); w = A(na - 1,en - 1)*A(en - 1,na - 1); if (l == na) goto L_280; if (itn <= 0) { /* .......... set error -- all eigenvalues have not * converged after 30*N iterations .......... */ ermsg( "SEVUN", en, 0, "ERROR NO. is index of eigenvalue causing convergence failure." , '.' ); Iflag[1] = 4; if (en <= 0) Iflag[1] = 3; return; } if (its != 10 && its != 20) goto L_130; /* .......... form exceptional shift .......... */ t += x; for (i = low; i <= en; i++) { A(i - 1,i - 1) -= x; } s = fabsf( A(na - 1,en - 1) ) + fabsf( A(enm2 - 1,na - 1) ); x = 0.75e0*s; y = x; w = -0.4375e0*s*s; L_130: its += 1; itn -= 1; /* .......... look for two consecutive small sub-diagonal elements. */ for (m = enm2; m >= l; m--) { zz = A(m - 1,m - 1); r = x - zz; s = y - zz; p = (r*s - w)/A(m - 1,m) + A(m,m - 1); q = A(m,m) - zz - r - s; r = A(m,m + 1); s = fabsf( p ) + fabsf( q ) + fabsf( r ); p /= s; q /= s; r /= s; if (m == l) goto L_150; tst1 = fabsf( p )*(fabsf( A(m - 2,m - 2) ) + fabsf( zz ) + fabsf( A(m,m) )); tst2 = tst1 + fabsf( A(m - 2,m - 1) )*(fabsf( q ) + fabsf( r )); if (tst2 == tst1) goto L_150; } L_150: mp2 = m + 2; for (i = mp2; i <= en; i++) { A(i - 3,i - 1) = 0.0e0; if (i != mp2) A(i - 4,i - 1) = 0.0e0; } /* .......... double QR step involving rows L to EN and * columns M to EN .......... */ for (k = m; k <= na; k++) { notlas = k != na; if (k != m) { p = A(k - 2,k - 1); q = A(k - 2,k); r = 0.0e0; if (notlas) r = A(k - 2,k + 1); x = fabsf( p ) + fabsf( q ) + fabsf( r ); if (x == 0.0e0) goto L_260; p /= x; q /= x; r /= x; } s = signf( sqrtf( p*p + q*q + r*r ), p ); if (k != m) { A(k - 2,k - 1) = -s*x; } else { if (l != m) A(k - 2,k - 1) = -A(k - 2,k - 1); } p += s; x = p/s; y = q/s; zz = r/s; q /= p; r /= p; if (notlas) { /* .......... row modification .......... */ for (j = k; j <= n; j++) { p = A(j - 1,k - 1) + q*A(j - 1,k) + r*A(j - 1,k + 1); A(j - 1,k - 1) += -p*x; A(j - 1,k) += -p*y; A(j - 1,k + 1) += -p*zz; } j = min( en, k + 3 ); /* .......... column modification .......... */ for (i = 1; i <= j; i++) { p = x*A(k - 1,i - 1) + y*A(k,i - 1) + zz*A(k + 1,i - 1); A(k - 1,i - 1) -= p; A(k,i - 1) += -p*q; A(k + 1,i - 1) += -p*r; } } else { /* .......... row modification .......... */ for (j = k; j <= n; j++) { p = A(j - 1,k - 1) + q*A(j - 1,k); A(j - 1,k - 1) += -p*x; A(j - 1,k) += -p*y; } j = min( en, k + 3 ); /* .......... column modification .......... */ for (i = 1; i <= j; i++) { p = x*A(k - 1,i - 1) + y*A(k,i - 1); A(k - 1,i - 1) -= p; A(k,i - 1) += -p*q; } } L_260: ; } goto L_70; /* .......... one root found .......... */ L_270: Evalr[en] = x + t; Evali[en] = 0.0e0; en = na; goto L_60; /* .......... two roots found .......... */ L_280: p = (y - x)/2.0e0; q = p*p + w; zz = sqrtf( fabsf( q ) ); x += t; if (q >= 0.0e0) { /* .......... real pair .......... */ zz = p + signf( zz, p ); Evalr[na] = x + zz; Evalr[en] = Evalr[na]; if (zz != 0.0e0) Evalr[en] = x - w/zz; Evali[na] = 0.0e0; Evali[en] = 0.0e0; } else { /* .......... complex pair .......... */ ltype = 2; Evalr[na] = x + p; Evalr[en] = x + p; Evali[na] = zz; Evali[en] = -zz; } en = enm2; goto L_60; L_300: ; /*-- Begin mask code changes * Set up for Shell sort * Sort so real parts are algebraically increasing * For = real parts, so abs(imag. parts) are increasing * For both =, sort on index -- preserves complex pair order *-- End mask code changes */ for (i = 1; i <= n; i++) { Iflag[i] = i; } l = 1; for (k = 1; k <= n; k++) { l = 3*l + 1; if (l >= n) goto L_2020; } L_2020: l = max( 1, (l - 4)/9 ); L_2030: for (j = l + 1; j <= n; j++) { k = Iflag[j]; p = Evalr[k]; i = j - l; L_2040: switch (SARITHIF(p - Evalr[Iflag[i]])) { case -1: goto L_2070; case 0: goto L_2050; case 1: goto L_2080; } L_2050: switch (SARITHIF(fabsf( Evali[k] ) - fabsf( Evali[Iflag[i]] ))) { case -1: goto L_2070; case 0: goto L_2060; case 1: goto L_2080; } L_2060: if (k > Iflag[i]) goto L_2080; L_2070: Iflag[i + l] = Iflag[i]; i -= l; if (i > 0) goto L_2040; L_2080: Iflag[i + l] = k; } l = (l - 1)/3; if (l != 0) goto L_2030; /* Indices in IFLAG now give the desired order -- * Move entries to get this order. */ L_2110: for (i = l + 1; i <= n; i++) { if (Iflag[i] != i) { l = i; m = i; p = Evalr[i]; q = Evali[i]; L_2120: k = Iflag[m]; Iflag[m] = m; if (k != l) { Evalr[m] = Evalr[k]; Evali[m] = Evali[k]; m = k; goto L_2120; } else { Evalr[m] = p; Evali[m] = q; goto L_2110; } } } Iflag[1] = ltype; return; #undef A } /* end of function */