/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:04 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "devbh.h" #include #include void /*FUNCTION*/ devbh( double *a, long lda, long n, long *low, long *igh, long iflag[], double scale[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) LOGICAL32 noconv; long int _l0, i, iexc, j, k, kp1, l, la, m, mm1; double b2, base, c, f, g, r, s, x, y; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Iflag = &iflag[0] - 1; double *const Scale = &scale[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-05-25 DEVBH Krogh Minor change for making .f90 version. *>> 1996-03-30 DEVBH Krogh Added external statement. *>> 1994-10-20 DEVBH Krogh Changes to use M77CON *>> 1992-04-23 DEVBH CLL Declaring all variables. *>> 1992-04-22 DEVBH Krogh Removed unused label 80. *>> 1991-10-25 DEVBH Krogh Initial version, converted from EISPACK. * * This subroutine is a slight modification of EISPACK subroutines * BALANC and ELMHES (dated August 1983). These routines in turn are * translations of the ALGOL procedures BALANCE, Num. Math. 13, * 293-304(1969) by Parlett and Reinsch. Handbook for Auto. Comp., * Vol.ii-Linear Algebra, 315-326(1971), and * ELMHES, Num. Math. 12, 349-368(1968) by Martin and Wilkinson. * Handbook for Auto. Comp., Vol.ii-Linear Algebra, 339-358(1971). * BALANC balances a real matrix & isolates eigenvalues if possible. * ELMHES reduces a submatrix situated in rows and columns LOW * through IGH of a real general matrix to upper Hessenberg form by * stabilized elementary similarity transformations. * ------------------------------------------------------------------ *--D replaces "?": ?EVBH * Both versions use I1MACH * ------------------------------------------------------------------ */ /* On input * A contains the input matrix to be balanced and reduced to * upper Hessenburg form. * 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 contains the balanced Hessenburg matrix. The multipliers * which were used in the reduction are stored in the remaining * triangle under the Hessenberg matrix. * LOW and IGH are two integers such that A(i,j) * is equal to zero if * (1) i is greater than j and * (2) j=1,...,LOW-1 or i=IGH+1,...,N. * SCALE contains information determining the permutations and * scaling factors used. Suppose that the principal submatrix * in rows LOW through IGH has been balanced, that p(j) denotes * the index interchanged with j during the permutation step, * and that the elements of the diagonal matrix used are * denoted by d(i,j). Then * SCALE(j) = p(j), for j = 1,...,LOW-1 * = d(j,j), j = LOW,...,IGH * = p(j) j = IGH+1,...,N. * The order in which the interchanges are made is N to IGH+1, * then 1 to LOW-1. * Note that 1 is returned for IGH if IGH is zero formally. * IFLAG contains information on the rows and columns interchanged in * the reduction. Only elements LOW through IGH are used. * * The ALGOL procedure exc contained in balance appears in * BALANC inline. (Note that the ALGOL roles of identifiers * K,L have been reversed.) * * --------------------------- BALANC ------------------------------- * */ base = FLT_RADIX; b2 = base*base; k = 1; l = n; goto L_100; /* .......... in-line procedure for row and column exchange ......... */ L_20: Scale[m] = j; if (j != m) { for (i = 1; i <= l; i++) { f = A(j - 1,i - 1); A(j - 1,i - 1) = A(m - 1,i - 1); A(m - 1,i - 1) = f; } for (i = k; i <= n; i++) { f = A(i - 1,j - 1); A(i - 1,j - 1) = A(i - 1,m - 1); A(i - 1,m - 1) = f; } } if (iexc != 0) goto L_130; /* .......... search for rows isolating an eigenvalue * and push them down .......... */ if (l <= 1) goto L_280; l -= 1; L_100: for (j = l; j >= 1; j--) { for (i = 1; i <= l; i++) { if (i == j) goto L_110; if (A(i - 1,j - 1) != 0.0e0) goto L_120; L_110: ; } m = l; iexc = 0; goto L_20; L_120: ; } goto L_140; /* .......... search for columns isolating an eigenvalue * and push them left .......... */ L_130: k += 1; L_140: for (j = k; j <= l; j++) { for (i = k; i <= l; i++) { if (i == j) goto L_150; if (A(j - 1,i - 1) != 0.0e0) goto L_170; L_150: ; } m = k; iexc = 1; goto L_20; L_170: ; } /* .......... now balance the submatrix in rows K to L .......... */ for (i = k; i <= l; i++) { Scale[i] = 1.0e0; } /* .......... iterative loop for norm reduction .......... */ L_190: noconv = FALSE; for (i = k; i <= l; i++) { c = 0.0e0; r = 0.0e0; for (j = k; j <= l; j++) { if (j != i) { c += fabs( A(i - 1,j - 1) ); r += fabs( A(j - 1,i - 1) ); } } /* .......... guard against zero C or R due to underflow .......... */ if ((c != 0.0e0) && (r != 0.0e0)) { g = r/base; f = 1.0e0; s = c + r; L_210: if (c < g) { f *= base; c *= b2; goto L_210; } g = r*base; L_230: if (c >= g) { f /= base; c /= b2; goto L_230; } /* .......... now balance .......... */ if ((c + r)/f < 0.95e0*s) { g = 1.0e0/f; Scale[i] *= f; noconv = TRUE; for (j = k; j <= n; j++) { A(j - 1,i - 1) *= g; } for (j = 1; j <= l; j++) { A(i - 1,j - 1) *= f; } } } } if (noconv) goto L_190; L_280: *low = k; *igh = l; /* --------------------------- ELMHES ------------------------------- * */ la = *igh - 1; kp1 = *low + 1; if (la < kp1) return; for (m = kp1; m <= la; m++) { mm1 = m - 1; x = 0.0e0; i = m; for (j = m; j <= *igh; j++) { if (fabs( A(mm1 - 1,j - 1) ) <= fabs( x )) goto L_300; x = A(mm1 - 1,j - 1); i = j; L_300: ; } Iflag[m] = i; if (i == m) goto L_330; /* .......... interchange rows and columns of A .......... */ for (j = mm1; j <= n; j++) { y = A(j - 1,i - 1); A(j - 1,i - 1) = A(j - 1,m - 1); A(j - 1,m - 1) = y; } for (j = 1; j <= *igh; j++) { y = A(i - 1,j - 1); A(i - 1,j - 1) = A(m - 1,j - 1); A(m - 1,j - 1) = y; } /* .......... end interchange .......... */ L_330: if (x == 0.0e0) goto L_380; for (i = m + 1; i <= *igh; i++) { y = A(mm1 - 1,i - 1); if (y == 0.0e0) goto L_360; y /= x; A(mm1 - 1,i - 1) = y; for (j = m; j <= n; j++) { A(j - 1,i - 1) += -y*A(j - 1,m - 1); } for (j = 1; j <= *igh; j++) { A(m - 1,j - 1) += y*A(i - 1,j - 1); } L_360: ; } L_380: ; } return; #undef A } /* end of function */