#include "blaswrap.h" /* dget53.f -- translated by f2c (version 20061008). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib; on Linux or Unix systems, link with .../path/to/libf2c.a -lm or, if you install libf2c.a in a standard place, with -lf2c -lm -- in that order, at the end of the command line, as in cc *.o -lf2c -lm Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., http://www.netlib.org/f2c/libf2c.zip */ #include "f2c.h" /* Subroutine */ int dget53_(doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *scale, doublereal *wr, doublereal *wi, doublereal *result, integer *info) { /* System generated locals */ integer a_dim1, a_offset, b_dim1, b_offset; doublereal d__1, d__2, d__3, d__4, d__5, d__6; /* Builtin functions */ double sqrt(doublereal); /* Local variables */ static doublereal s1, ci11, ci12, ci22, cr11, cr12, cr21, cr22, ulp, wis, wrs, deti, absw, detr, temp, anorm, bnorm, cnorm; extern doublereal dlamch_(char *); static doublereal cscale, scales, safmin, sigmin; /* -- LAPACK test routine (version 3.1) -- Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. November 2006 Purpose ======= DGET53 checks the generalized eigenvalues computed by DLAG2. The basic test for an eigenvalue is: | det( s A - w B ) | RESULT = --------------------------------------------------- ulp max( s norm(A), |w| norm(B) )*norm( s A - w B ) Two "safety checks" are performed: (1) ulp*max( s*norm(A), |w|*norm(B) ) must be at least safe_minimum. This insures that the test performed is not essentially det(0*A + 0*B)=0. (2) s*norm(A) + |w|*norm(B) must be less than 1/safe_minimum. This insures that s*A - w*B will not overflow. If these tests are not passed, then s and w are scaled and tested anyway, if this is possible. Arguments ========= A (input) DOUBLE PRECISION array, dimension (LDA, 2) The 2x2 matrix A. LDA (input) INTEGER The leading dimension of A. It must be at least 2. B (input) DOUBLE PRECISION array, dimension (LDB, N) The 2x2 upper-triangular matrix B. LDB (input) INTEGER The leading dimension of B. It must be at least 2. SCALE (input) DOUBLE PRECISION The "scale factor" s in the formula s A - w B . It is assumed to be non-negative. WR (input) DOUBLE PRECISION The real part of the eigenvalue w in the formula s A - w B . WI (input) DOUBLE PRECISION The imaginary part of the eigenvalue w in the formula s A - w B . RESULT (output) DOUBLE PRECISION If INFO is 2 or less, the value computed by the test described above. If INFO=3, this will just be 1/ulp. INFO (output) INTEGER =0: The input data pass the "safety checks". =1: s*norm(A) + |w|*norm(B) > 1/safe_minimum. =2: ulp*max( s*norm(A), |w|*norm(B) ) < safe_minimum =3: same as INFO=2, but s and w could not be scaled so as to compute the test. ===================================================================== Initialize Parameter adjustments */ a_dim1 = *lda; a_offset = 1 + a_dim1; a -= a_offset; b_dim1 = *ldb; b_offset = 1 + b_dim1; b -= b_offset; /* Function Body */ *info = 0; *result = 0.; scales = *scale; wrs = *wr; wis = *wi; /* Machine constants and norms */ safmin = dlamch_("Safe minimum"); ulp = dlamch_("Epsilon") * dlamch_("Base"); absw = abs(wrs) + abs(wis); /* Computing MAX */ d__5 = (d__1 = a[a_dim1 + 1], abs(d__1)) + (d__2 = a[a_dim1 + 2], abs( d__2)), d__6 = (d__3 = a[(a_dim1 << 1) + 1], abs(d__3)) + (d__4 = a[(a_dim1 << 1) + 2], abs(d__4)), d__5 = max(d__5,d__6); anorm = max(d__5,safmin); /* Computing MAX */ d__4 = (d__3 = b[b_dim1 + 1], abs(d__3)), d__5 = (d__1 = b[(b_dim1 << 1) + 1], abs(d__1)) + (d__2 = b[(b_dim1 << 1) + 2], abs(d__2)), d__4 = max(d__4,d__5); bnorm = max(d__4,safmin); /* Check for possible overflow. */ temp = safmin * bnorm * absw + safmin * anorm * scales; if (temp >= 1.) { /* Scale down to avoid overflow */ *info = 1; temp = 1. / temp; scales *= temp; wrs *= temp; wis *= temp; absw = abs(wrs) + abs(wis); } /* Computing MAX Computing MAX */ d__3 = scales * anorm, d__4 = absw * bnorm; d__1 = ulp * max(d__3,d__4), d__2 = safmin * max(scales,absw); s1 = max(d__1,d__2); /* Check for W and SCALE essentially zero. */ if (s1 < safmin) { *info = 2; if (scales < safmin && absw < safmin) { *info = 3; *result = 1. / ulp; return 0; } /* Scale up to avoid underflow Computing MAX */ d__1 = scales * anorm + absw * bnorm; temp = 1. / max(d__1,safmin); scales *= temp; wrs *= temp; wis *= temp; absw = abs(wrs) + abs(wis); /* Computing MAX Computing MAX */ d__3 = scales * anorm, d__4 = absw * bnorm; d__1 = ulp * max(d__3,d__4), d__2 = safmin * max(scales,absw); s1 = max(d__1,d__2); if (s1 < safmin) { *info = 3; *result = 1. / ulp; return 0; } } /* Compute C = s A - w B */ cr11 = scales * a[a_dim1 + 1] - wrs * b[b_dim1 + 1]; ci11 = -wis * b[b_dim1 + 1]; cr21 = scales * a[a_dim1 + 2]; cr12 = scales * a[(a_dim1 << 1) + 1] - wrs * b[(b_dim1 << 1) + 1]; ci12 = -wis * b[(b_dim1 << 1) + 1]; cr22 = scales * a[(a_dim1 << 1) + 2] - wrs * b[(b_dim1 << 1) + 2]; ci22 = -wis * b[(b_dim1 << 1) + 2]; /* Compute the smallest singular value of s A - w B: |det( s A - w B )| sigma_min = ------------------ norm( s A - w B ) Computing MAX */ d__1 = abs(cr11) + abs(ci11) + abs(cr21), d__2 = abs(cr12) + abs(ci12) + abs(cr22) + abs(ci22), d__1 = max(d__1,d__2); cnorm = max(d__1,safmin); cscale = 1. / sqrt(cnorm); detr = cscale * cr11 * (cscale * cr22) - cscale * ci11 * (cscale * ci22) - cscale * cr12 * (cscale * cr21); deti = cscale * cr11 * (cscale * ci22) + cscale * ci11 * (cscale * cr22) - cscale * ci12 * (cscale * cr21); sigmin = abs(detr) + abs(deti); *result = sigmin / s1; return 0; /* End of DGET53 */ } /* dget53_ */