/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:54 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "scomqr.h" #include #include void /*FUNCTION*/ scomqr( long nm, long n, long low, long igh, float *hr, float *hi, float z[], long *ierr) { #define HR(I_,J_) (*(hr+(I_)*(nm)+(J_))) #define HI(I_,J_) (*(hi+(I_)*(nm)+(J_))) long int _l0, en, enm1, i, its, j, l, ll, lp1; float machep, norm, si, sr, ti, tr, xi, xr, yi, yr, z1[2], z2[2], zdenom[2], zin[2], znum[2], zout[2], zquot[2], zzi, zzr; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Z = &z[0] - 1; float *const Z1 = &z1[0] - 1; float *const Z2 = &z2[0] - 1; float *const Zdenom = &zdenom[0] - 1; float *const Zin = &zin[0] - 1; float *const Znum = &znum[0] - 1; float *const Zout = &zout[0] - 1; float *const Zquot = &zquot[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. *>> 2001-01-24 SCOMQR Krogh ZSQRT -> CSQRTX to avoid C lib. conflicts. *>> 1996-04-27 SCOMQR Krogh Changes to use .C. and C%%. *>> 1996-03-30 SCOMQR Krogh Added external statement. *>> 1996-01-18 SCOMQR Krogh Added M77CON statements for conv. to C. *>> 1995-01-03 SCOMQR WVS Added EXTERNAL CQUO, ZSQRT so VAX won't gripe *>> 1992-03-13 SCOMQR FTK Removed implicit statements. *>> 1987-11-12 SCOMQR Lawson Initial code. *--S Replaces "?": ?COMQR *--C (Type) Replaces "?": ?QUO, ?SQRTX * ------------------------------------------------------------------ * Version of the Eispack subr, COMQR, for use in the JPL MATH77 * library. C. L. Lawson, JPL, 1987 Feb 17. * ------------------------------------------------------------------ */ /*++ Default NO_COMPLEX = .C. | (.N. == 'D') *++ Default COMPLEX = ~NO_COMPLEX *++ CODE for COMPLEX is inactive * COMPLEX Z(N), Z3 *++ CODE for NO_COMPLEX is active *++ Replace "SCABS" = .N.//.Y.//"ABS" */ /*++ END * * THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE * ALGOL PROCEDURE COMLR, NUM. MATH. 12, 369-376(1968) BY MARTIN * AND WILKINSON. * HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(1971). * THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS * (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM. * * THIS SUBROUTINE FINDS THE EIGENVALUES OF A COMPLEX * UPPER HESSENBERG MATRIX BY THE QR METHOD. * * ON INPUT- * * NM 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, * * LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING * SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, * SET LOW=1, IGH=N, * * HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, * RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. * THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN * INFORMATION ABOUT THE UNITARY TRANSFORMATIONS USED IN * THE REDUCTION BY CORTH, IF PERFORMED. * * ON OUTPUT- * * THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN * DESTROYED. THEREFORE, THEY MUST BE SAVED BEFORE * CALLING COMQR IF SUBSEQUENT CALCULATION OF * EIGENVECTORS IS TO BE PERFORMED, * * WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, * RESPECTIVELY, OF THE EIGENVALUES. IF AN ERROR * EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT * FOR INDICES IERR+1,...,N, * * IERR IS SET TO * ZERO FOR NORMAL RETURN, * J IF THE J-TH EIGENVALUE HAS NOT BEEN * DETERMINED AFTER 30 ITERATIONS. * * ARITHMETIC IS REAL EXCEPT FOR THE REPLACEMENT OF THE ALGOL * PROCEDURE CDIV BY COMPLEX DIVISION AND USE OF THE SUBROUTINES * SQRT AND CMPLX IN COMPUTING COMPLEX SQUARE ROOTS. * * QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, * APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY * * ------------------------------------------------------------------ * * ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING * THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC. * * ********** */ machep = FLT_EPSILON; *ierr = 0; if (low == igh) goto L_180; /* ********** CREATE REAL SUBDIAGONAL ELEMENTS ********** */ l = low + 1; for (i = l; i <= igh; i++) { ll = min( i + 1, igh ); if (HI(i - 2,i - 1) == 0.0) goto L_170; /*++ CODE for COMPLEX is inactive * NORM = ABS(CMPLX(HR(I,I-1),HI(I,I-1))) *++ CODE for NO_COMPLEX is active */ Z1[1] = HR(i - 2,i - 1); Z1[2] = HR(i - 2,i - 1); norm = scabs( z1 ); /*++ END */ yr = HR(i - 2,i - 1)/norm; yi = HI(i - 2,i - 1)/norm; HR(i - 2,i - 1) = norm; HI(i - 2,i - 1) = 0.0; for (j = i; j <= igh; j++) { si = yr*HI(j - 1,i - 1) - yi*HR(j - 1,i - 1); HR(j - 1,i - 1) = yr*HR(j - 1,i - 1) + yi*HI(j - 1,i - 1); HI(j - 1,i - 1) = si; } for (j = low; j <= ll; j++) { si = yr*HI(i - 1,j - 1) + yi*HR(i - 1,j - 1); HR(i - 1,j - 1) = yr*HR(i - 1,j - 1) - yi*HI(i - 1,j - 1); HI(i - 1,j - 1) = si; } L_170: ; } /* ********** STORE ROOTS ISOLATED BY CBAL ********** */ L_180: for (i = 1; i <= n; i++) { if (i >= low && i <= igh) goto L_200; /*++ CODE for COMPLEX is inactive * Z(I) = CMPLX(HR(I,I),HI(I,I)) *++ CODE for NO_COMPLEX is active */ Z[2*i - 1] = HR(i - 1,i - 1); Z[2*i] = HI(i - 1,i - 1); /*++ END */ L_200: ; } en = igh; tr = 0.0; ti = 0.0; /* ********** SEARCH FOR NEXT EIGENVALUE ********** */ L_220: if (en < low) goto L_1001; its = 0; enm1 = en - 1; /* ********** LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT * FOR L=EN STEP -1 UNTIL LOW -- ********** */ L_240: for (ll = low; ll <= en; ll++) { l = en + low - ll; if (l == low) goto L_300; if (fabsf( HR(l - 2,l - 1) ) <= machep*(fabsf( HR(l - 2,l - 2) ) + fabsf( HI(l - 2,l - 2) ) + fabsf( HR(l - 1,l - 1) ) + fabsf( HI(l - 1,l - 1) ))) goto L_300; } /* ********** FORM SHIFT ********** */ L_300: if (l == en) goto L_660; if (its == 30) goto L_1000; if (its == 10 || its == 20) goto L_320; sr = HR(en - 1,en - 1); si = HI(en - 1,en - 1); xr = HR(en - 1,enm1 - 1)*HR(enm1 - 1,en - 1); xi = HI(en - 1,enm1 - 1)*HR(enm1 - 1,en - 1); if (xr == 0.0 && xi == 0.0) goto L_340; yr = (HR(enm1 - 1,enm1 - 1) - sr)/2.0; yi = (HI(enm1 - 1,enm1 - 1) - si)/2.0; /*++ CODE for COMPLEX is inactive * Z3 = SQRT(CMPLX(YR**2-YI**2+XR,2.0*YR*YI+XI)) * ZZR = REAL(Z3) * ZZI = AIMAG(Z3) *++ CODE for NO_COMPLEX is active */ Zin[1] = SQ(yr) - SQ(yi) + xr; Zin[2] = 2.e0*yr*yi + xi; csqrtx( zin, zout ); zzr = Zout[1]; zzi = Zout[2]; /*++ END */ if (yr*zzr + yi*zzi >= 0.0) goto L_310; zzr = -zzr; zzi = -zzi; L_310: ; /*++ CODE for COMPLEX is inactive * Z3 = CMPLX(XR,XI) / CMPLX(YR+ZZR,YI+ZZI) * SR = SR - REAL(Z3) * SI = SI - AIMAG(Z3) *++ CODE for NO_COMPLEX is active */ Znum[1] = xr; Znum[2] = xi; Zdenom[1] = yr + zzr; Zdenom[2] = yi + zzi; /* CQUO computes quotient of two complex numbers. */ cquo( znum, zdenom, zquot ); sr -= Zquot[1]; si -= Zquot[2]; /*++ END */ goto L_340; /* ********** FORM EXCEPTIONAL SHIFT ********** */ L_320: sr = fabsf( HR(enm1 - 1,en - 1) ) + fabsf( HR(en - 3,enm1 - 1) ); si = 0.0; L_340: for (i = low; i <= en; i++) { HR(i - 1,i - 1) -= sr; HI(i - 1,i - 1) -= si; } tr += sr; ti += si; its += 1; /* ********** REDUCE TO TRIANGLE (ROWS) ********** */ lp1 = l + 1; for (i = lp1; i <= en; i++) { sr = HR(i - 2,i - 1); HR(i - 2,i - 1) = 0.0; norm = sqrtf( HR(i - 2,i - 2)*HR(i - 2,i - 2) + HI(i - 2,i - 2)* HI(i - 2,i - 2) + sr*sr ); xr = HR(i - 2,i - 2)/norm; xi = HI(i - 2,i - 2)/norm; /*++ CODE for COMPLEX is inactive * Z(I-1) = CMPLX(XR,XI) *++ CODE for NO_COMPLEX is active */ Z[2*i - 3] = xr; Z[2*i - 2] = xi; /*++ END */ HR(i - 2,i - 2) = norm; HI(i - 2,i - 2) = 0.0; HI(i - 2,i - 1) = sr/norm; for (j = i; j <= en; j++) { yr = HR(j - 1,i - 2); yi = HI(j - 1,i - 2); zzr = HR(j - 1,i - 1); zzi = HI(j - 1,i - 1); HR(j - 1,i - 2) = xr*yr + xi*yi + HI(i - 2,i - 1)*zzr; HI(j - 1,i - 2) = xr*yi - xi*yr + HI(i - 2,i - 1)*zzi; HR(j - 1,i - 1) = xr*zzr - xi*zzi - HI(i - 2,i - 1)*yr; HI(j - 1,i - 1) = xr*zzi + xi*zzr - HI(i - 2,i - 1)*yi; } } si = HI(en - 1,en - 1); if (si == 0.0) goto L_540; /*++ CODE for COMPLEX is inactive * NORM = ABS(CMPLX(HR(EN,EN),SI)) *++ CODE for NO_COMPLEX is active */ Z2[1] = HR(en - 1,en - 1); Z2[2] = si; norm = scabs( z2 ); /*++ END */ sr = HR(en - 1,en - 1)/norm; si /= norm; HR(en - 1,en - 1) = norm; HI(en - 1,en - 1) = 0.0; /* ********** INVERSE OPERATION (COLUMNS) ********** */ L_540: for (j = lp1; j <= en; j++) { /*++ CODE for COMPLEX is inactive * XR = REAL(Z(J-1)) * XI = AIMAG(Z(J-1)) *++ CODE for NO_COMPLEX is active */ xr = Z[2*j - 3]; xi = Z[2*j - 2]; /*++ END * */ for (i = l; i <= j; i++) { yr = HR(j - 2,i - 1); yi = 0.0; zzr = HR(j - 1,i - 1); zzi = HI(j - 1,i - 1); if (i == j) goto L_560; yi = HI(j - 2,i - 1); HI(j - 2,i - 1) = xr*yi + xi*yr + HI(j - 2,j - 1)*zzi; L_560: HR(j - 2,i - 1) = xr*yr - xi*yi + HI(j - 2,j - 1)*zzr; HR(j - 1,i - 1) = xr*zzr + xi*zzi - HI(j - 2,j - 1)*yr; HI(j - 1,i - 1) = xr*zzi - xi*zzr - HI(j - 2,j - 1)*yi; } } if (si == 0.0) goto L_240; for (i = l; i <= en; i++) { yr = HR(en - 1,i - 1); yi = HI(en - 1,i - 1); HR(en - 1,i - 1) = sr*yr - si*yi; HI(en - 1,i - 1) = sr*yi + si*yr; } goto L_240; /* ********** A ROOT FOUND ********** */ L_660: ; /*++ CODE for COMPLEX is inactive * Z(EN) = CMPLX(HR(EN,EN)+TR,HI(EN,EN)+TI) *++ CODE for NO_COMPLEX is active */ Z[2*en - 1] = HR(en - 1,en - 1) + tr; Z[2*en] = HI(en - 1,en - 1) + ti; /*++ END */ en = enm1; goto L_220; /* ********** SET ERROR -- NO CONVERGENCE TO AN * EIGENVALUE AFTER 30 ITERATIONS ********** */ L_1000: *ierr = en; L_1001: return; /* ********** LAST CARD OF COMQR ********** */ #undef HR #undef HI } /* end of function */