/*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 "sevvun.h" #include void /*FUNCTION*/ sevvun( float *a, long lda, long n, float evalr[], float evali[], float *vec, long iflag[], float work[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) #define VEC(I_,J_) (*(vec+(I_)*(lda)+(J_))) LOGICAL32 notlas; long int en, enm2, i, igh, ii, itn, its, j, k, l, low, ltype, m, mp, mp2, na; float norm, p, q, r, ra, s, sa, t, tst1, tst2, vi, vr, 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; 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. * This file contains SEVVUN and CDIV. *>> 1996-03-30 DEVUNN Krogh MIN0 => MIN, added external statement. *>> 1994-11-02 SEVVUN Krogh Changes to use M77CON *>> 1992-04-23 SEVVUN CLL Declaring all variables. *>> 1992-04-22 SEVVUN Krogh Removed unused labels 330 and 1220. *>> 1992-03-05 SEVVUN Krogh Initial version, converted from EISPACK. * * This subroutine uses slight modifcations of EISPACK routines * BALANC, ELMHES, ELTRAN, HQR2 and BALBAK to get the eigenvalues and * eigenvectors of a general real matrix by the QR method. The first * two are are encapsulated in SEVBH. The following three are given * inline here. * * ELTRAN is a translation of the algol procedure ELMTRANS, * Num. Math. 16, 181-204(1970) by Peters and Wilkinson. * Handbook for Auto. Comp., Vol.ii-Linear Algebra, 372-395(1971). * * HQR2 is a translation of the ALGOL procedure HQR2, * Num. Math. 16, 181-204(1970) by Peters and Wilkinson. * Handbook for Auto. Comp., Vol.ii-Linear Algebra, 372-395(1971). * * BALBAK is a translation of the ALGOL procedure BALBAK, * Num. Math. 13, 293-304(1969) by Parlett and Reinsch. * Handbook for Auto. Comp., Vol.ii-Linear Algebra, 315-326(1971). * * This subroutine finds the eigenvalues and eigenvectors * of a real general matrix by the QR method. * * On input * * 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. * A contains the input matrix whose eigenvalues and eigenvectors * are desired. * * 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. * VEC contains the real and imaginary parts of the eigenvectors. * if the I-th eigenvalue is real, the I-th column of VEC * contains its eigenvector. If the I-th eigenvalue is complex * with positive imaginary part, the I-th and (I+1)-th * columns of VEC contain the real and imaginary parts of its * eigenvector. The eigenvectors are unnormalized. If an * error exit is made, none of the eigenvectors has been found. * 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 "?": ?EVVUN, ?EVBH, ?NRM2, ?SCAL, ?CDIV * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ * */ sevbh( a, lda, n, &low, &igh, iflag, work ); ltype = 1; /* -------------------------------- ELTRAN --------------------------- * * Accumulate the stabilized elementary similarity transformations * used in the preceding reduction of the matrix to upper Hessenberg * form. * * .......... initialize VEC to identity matrix .......... */ for (j = 1; j <= n; j++) { for (i = 1; i <= n; i++) { VEC(j - 1,i - 1) = 0.0e0; } VEC(j - 1,j - 1) = 1.0e0; } if (igh > low + 1) { for (mp = igh - 1; mp >= (low + 1); mp--) { for (i = mp + 1; i <= igh; i++) { VEC(mp - 1,i - 1) = A(mp - 2,i - 1); } i = Iflag[mp]; if (i != mp) { for (j = mp; j <= igh; j++) { VEC(j - 1,mp - 1) = VEC(j - 1,i - 1); VEC(j - 1,i - 1) = 0.0e0; } VEC(mp - 1,i - 1) = 1.0e0; } } } /* -------------------------------- HQR2 ----------------------------- * * Find the eigenvalues of a real upper Hessenberg matrix by the QR * method. * */ Iflag[1] = 0; norm = 0.0e0; k = 1; /* .......... store roots isolated by balanc * 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_80: if (en < low) goto L_340; its = 0; na = en - 1; enm2 = na - 1; /* .......... look for single small sub-diagonal element */ L_90: 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_110; } /* .......... form shift .......... */ L_110: 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( "SEVVUN", 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; } /* .......... accumulate transformations .......... */ for (i = low; i <= igh; i++) { p = x*VEC(k - 1,i - 1) + y*VEC(k,i - 1) + zz*VEC(k + 1,i - 1); VEC(k - 1,i - 1) -= p; VEC(k,i - 1) += -p*q; VEC(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; } /* .......... accumulate transformations .......... */ for (i = low; i <= igh; i++) { p = x*VEC(k - 1,i - 1) + y*VEC(k,i - 1); VEC(k - 1,i - 1) -= p; VEC(k,i - 1) += -p*q; } } L_260: ; } goto L_90; /* .......... one root found .......... */ L_270: A(en - 1,en - 1) = x + t; Evalr[en] = A(en - 1,en - 1); Evali[en] = 0.0e0; en = na; goto L_80; /* .......... two roots found .......... */ L_280: p = (y - x)/2.0e0; q = p*p + w; zz = sqrtf( fabsf( q ) ); A(en - 1,en - 1) = x + t; x = A(en - 1,en - 1); A(na - 1,na - 1) = y + 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; x = A(na - 1,en - 1); s = fabsf( x ) + fabsf( zz ); p = x/s; q = zz/s; r = sqrtf( p*p + q*q ); p /= r; q /= r; /* .......... row modification .......... */ for (j = na; j <= n; j++) { zz = A(j - 1,na - 1); A(j - 1,na - 1) = q*zz + p*A(j - 1,en - 1); A(j - 1,en - 1) = q*A(j - 1,en - 1) - p*zz; } /* .......... column modification .......... */ for (i = 1; i <= en; i++) { zz = A(na - 1,i - 1); A(na - 1,i - 1) = q*zz + p*A(en - 1,i - 1); A(en - 1,i - 1) = q*A(en - 1,i - 1) - p*zz; } /* .......... accumulate transformations .......... */ for (i = low; i <= igh; i++) { zz = VEC(na - 1,i - 1); VEC(na - 1,i - 1) = q*zz + p*VEC(en - 1,i - 1); VEC(en - 1,i - 1) = q*VEC(en - 1,i - 1) - p*zz; } } else { /* .......... complex pair .......... */ ltype = 2; Evalr[na] = x + p; Evalr[en] = x + p; Evali[na] = zz; Evali[en] = -zz; } en = enm2; goto L_80; /* .......... all roots found. Backsubstitute to find * vectors of upper triangular form .......... */ L_340: if (norm == 0.0e0) goto L_1000; for (en = n; en >= 1; en--) { p = Evalr[en]; q = Evali[en]; na = en - 1; switch (SARITHIF(q)) { case -1: goto L_710; case 0: goto L_600; case 1: goto L_800; } /* .......... real vector .......... */ L_600: m = en; A(en - 1,en - 1) = 1.0e0; for (i = na; i >= 1; i--) { w = A(i - 1,i - 1) - p; r = 0.0e0; for (j = m; j <= en; j++) { r += A(j - 1,i - 1)*A(en - 1,j - 1); } if (Evali[i] < 0.0e0) { zz = w; s = r; } else { m = i; if (Evali[i] == 0.0e0) { t = w; if (t == 0.0e0) { tst1 = norm; t = tst1; L_640: t *= 0.01e0; tst2 = norm + t; if (tst2 > tst1) goto L_640; } A(en - 1,i - 1) = -r/t; } else { /* .......... solve real equations .......... */ x = A(i,i - 1); y = A(i - 1,i); q = (Evalr[i] - p)*(Evalr[i] - p) + Evali[i]*Evali[i]; t = (x*s - zz*r)/q; A(en - 1,i - 1) = t; if (fabsf( x ) > fabsf( zz )) { A(en - 1,i) = (-r - w*t)/x; } else { A(en - 1,i) = (-s - y*t)/zz; } } /* .......... overflow control .......... */ t = fabsf( A(en - 1,i - 1) ); if (t != 0.0e0) { tst1 = t; tst2 = tst1 + 1.0e0/tst1; if (tst2 <= tst1) { for (j = i; j <= en; j++) { A(en - 1,j - 1) /= t; } } } } } /* .......... end real vector .......... */ goto L_800; /* .......... complex vector .......... */ L_710: m = na; /* .......... last vector component chosen imaginary so that * eigenvector matrix is triangular .......... */ if (fabsf( A(na - 1,en - 1) ) > fabsf( A(en - 1,na - 1) )) { A(na - 1,na - 1) = q/A(na - 1,en - 1); A(en - 1,na - 1) = -(A(en - 1,en - 1) - p)/A(na - 1,en - 1); } else { scdiv( 0.0e0, -A(en - 1,na - 1), A(na - 1,na - 1) - p, q, &A(na - 1,na - 1), &A(en - 1,na - 1) ); } A(na - 1,en - 1) = 0.0e0; A(en - 1,en - 1) = 1.0e0; enm2 = na - 1; for (i = enm2; i >= 1; i--) { w = A(i - 1,i - 1) - p; ra = 0.0e0; sa = 0.0e0; for (j = m; j <= en; j++) { ra += A(j - 1,i - 1)*A(na - 1,j - 1); sa += A(j - 1,i - 1)*A(en - 1,j - 1); } if (Evali[i] < 0.0e0) { zz = w; r = ra; s = sa; } else { m = i; if (Evali[i] == 0.0e0) { scdiv( -ra, -sa, w, q, &A(na - 1,i - 1), &A(en - 1,i - 1) ); } else { /* .......... solve complex equations .......... */ x = A(i,i - 1); y = A(i - 1,i); vr = (Evalr[i] - p)*(Evalr[i] - p) + Evali[i]* Evali[i] - q*q; vi = (Evalr[i] - p)*2.0e0*q; if (vr == 0.0e0 && vi == 0.0e0) { tst1 = norm*(fabsf( w ) + fabsf( q ) + fabsf( x ) + fabsf( y ) + fabsf( zz )); vr = tst1; L_740: vr *= 0.01e0; tst2 = tst1 + vr; if (tst2 > tst1) goto L_740; } scdiv( x*r - zz*ra + q*sa, x*s - zz*sa - q*ra, vr, vi, &A(na - 1,i - 1), &A(en - 1,i - 1) ); if (fabsf( x ) > fabsf( zz ) + fabsf( q )) { A(na - 1,i) = (-ra - w*A(na - 1,i - 1) + q* A(en - 1,i - 1))/x; A(en - 1,i) = (-sa - w*A(en - 1,i - 1) - q* A(na - 1,i - 1))/x; } else { scdiv( -r - y*A(na - 1,i - 1), -s - y*A(en - 1,i - 1), zz, q, &A(na - 1,i), &A(en - 1,i) ); } } /* .......... overflow control .......... */ t = fmaxf( fabsf( A(na - 1,i - 1) ), fabsf( A(en - 1,i - 1) ) ); if (t != 0.0e0) { tst1 = t; tst2 = tst1 + 1.0e0/tst1; if (tst2 <= tst1) { for (j = i; j <= en; j++) { A(na - 1,j - 1) /= t; A(en - 1,j - 1) /= t; } } } } } /* .......... end complex vector .......... */ L_800: ; } /* .......... end back substitution. * vectors of isolated roots .......... */ for (i = 1; i <= n; i++) { if (i >= low && i <= igh) goto L_840; for (j = i; j <= n; j++) { VEC(j - 1,i - 1) = A(j - 1,i - 1); } L_840: ; } /* .......... multiply by transformation matrix to give * vectors of original full matrix. */ for (j = n; j >= low; j--) { m = min( j, igh ); for (i = low; i <= igh; i++) { zz = 0.0e0; for (k = low; k <= m; k++) { zz += VEC(k - 1,i - 1)*A(j - 1,k - 1); } VEC(j - 1,i - 1) = zz; } } L_1000: ; /* ------------------------------ BALBAK ----------------------------- * * Form eigenvectors of a real general matrix by back transforming * those of the corresponding balanced matrix determined by BALANC. * */ if (igh != low) { for (i = low; i <= igh; i++) { s = Work[i]; /* .......... left hand eigenvectors are back transformed * if the foregoing statement is replaced by * S=1.0E0/WORK(I). .......... */ for (j = 1; j <= n; j++) { VEC(j - 1,i - 1) *= s; } } } /* ......... for I=LOW-1 step -1 until 1, * IGH+1 step 1 until N do -- .......... */ for (ii = 1; ii <= n; ii++) { i = ii; if ((i < low) || (i > igh)) { if (i < low) i = low - ii; k = Work[i]; if (k != i) { for (j = 1; j <= n; j++) { s = VEC(j - 1,i - 1); VEC(j - 1,i - 1) = VEC(j - 1,k - 1); VEC(j - 1,k - 1) = s; } } } } /* Normalize the eigenvectors */ for (i = 1; i <= n; i++) { p = snrm2( n, &VEC(i - 1,0), 1 ); if (Evali[i] == 0.e0) { sscal( n, signf( 1.e0, VEC(i - 1,0) )/p, &VEC(i - 1,0), 1 ); } else if (Evali[i] > 0.e0) { q = snrm2( n, &VEC(i,0), 1 ); if (p > q) { p *= sqrtf( 1.e0 + powif(q/p,2) ); } else if (q != 0.e0) { p = q*sqrtf( 1.e0 + powif(p/q,2) ); } else { goto L_1220; } if (VEC(i,0) == 0.e0) { p = signf( 1.e0, VEC(i,0) )/p; q = 0.e0; } else if (fabsf( VEC(i - 1,0) ) > fabsf( VEC(i,0) )) { p = 1.e0/(p*sqrtf( 1.e0 + powif(VEC(i,0)/VEC(i - 1,0),2) )); q = p*(VEC(i,0)/fabsf( VEC(i - 1,0) )); p = signf( p, VEC(i - 1,0) ); } else { q = 1.e0/(p*sqrtf( 1.e0 + powif(VEC(i - 1,0)/VEC(i,0),2) )); p = q*(VEC(i - 1,0)/fabsf( VEC(i,0) )); q = signf( q, VEC(i,0) ); } for (j = 1; j <= n; j++) { r = p*VEC(i - 1,j - 1) + q*VEC(i,j - 1); VEC(i,j - 1) = p*VEC(i,j - 1) - q*VEC(i - 1,j - 1); VEC(i - 1,j - 1) = r; } } L_1220: ; } /*-- 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]; for (j = 1; j <= n; j++) { Work[j] = VEC(i - 1,j - 1); } L_2120: k = Iflag[m]; Iflag[m] = m; if (k != l) { Evalr[m] = Evalr[k]; Evali[m] = Evali[k]; for (j = 1; j <= n; j++) { VEC(m - 1,j - 1) = VEC(k - 1,j - 1); } m = k; goto L_2120; } else { Evalr[m] = p; Evali[m] = q; for (j = 1; j <= n; j++) { VEC(m - 1,j - 1) = Work[j]; } goto L_2110; } } } Iflag[1] = ltype; return; #undef VEC #undef A } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ scdiv( float ar, float ai, float br, float bi, float *cr, float *ci) { float ais, ars, bis, brs, s; /* complex division, (CR,CI) = (AR,AI)/(BR,BI) * ------------------------------------------------------------------ */ s = fabsf( br ) + fabsf( bi ); ars = ar/s; ais = ai/s; brs = br/s; bis = bi/s; s = SQ(brs) + SQ(bis); *cr = (ars*brs + ais*bis)/s; *ci = (ais*brs - ars*bis)/s; return; } /* end of function */