swiftnav: dlagv2.c Source File

Go to the documentation of this file.
00001 /* dlagv2.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Table of constant values */
00017 
00018 static integer c__2 = 2;
00019 static integer c__1 = 1;
00020 
00021 /* Subroutine */ int dlagv2_(doublereal *a, integer *lda, doublereal *b, 
00022         integer *ldb, doublereal *alphar, doublereal *alphai, doublereal *
00023         beta, doublereal *csl, doublereal *snl, doublereal *csr, doublereal *
00024         snr)
00025 {
00026     /* System generated locals */
00027     integer a_dim1, a_offset, b_dim1, b_offset;
00028     doublereal d__1, d__2, d__3, d__4, d__5, d__6;
00029 
00030     /* Local variables */
00031     doublereal r__, t, h1, h2, h3, wi, qq, rr, wr1, wr2, ulp;
00032     extern /* Subroutine */ int drot_(integer *, doublereal *, integer *, 
00033             doublereal *, integer *, doublereal *, doublereal *), dlag2_(
00034             doublereal *, integer *, doublereal *, integer *, doublereal *, 
00035             doublereal *, doublereal *, doublereal *, doublereal *, 
00036             doublereal *);
00037     doublereal anorm, bnorm, scale1, scale2;
00038     extern /* Subroutine */ int dlasv2_(doublereal *, doublereal *, 
00039             doublereal *, doublereal *, doublereal *, doublereal *, 
00040             doublereal *, doublereal *, doublereal *);
00041     extern doublereal dlapy2_(doublereal *, doublereal *);
00042     doublereal ascale, bscale;
00043     extern doublereal dlamch_(char *);
00044     doublereal safmin;
00045     extern /* Subroutine */ int dlartg_(doublereal *, doublereal *, 
00046             doublereal *, doublereal *, doublereal *);
00047 
00048 
00049 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00050 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00051 /*     November 2006 */
00052 
00053 /*     .. Scalar Arguments .. */
00054 /*     .. */
00055 /*     .. Array Arguments .. */
00056 /*     .. */
00057 
00058 /*  Purpose */
00059 /*  ======= */
00060 
00061 /*  DLAGV2 computes the Generalized Schur factorization of a real 2-by-2 */
00062 /*  matrix pencil (A,B) where B is upper triangular. This routine */
00063 /*  computes orthogonal (rotation) matrices given by CSL, SNL and CSR, */
00064 /*  SNR such that */
00065 
00066 /*  1) if the pencil (A,B) has two real eigenvalues (include 0/0 or 1/0 */
00067 /*     types), then */
00068 
00069 /*     [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ] */
00070 /*     [  0  a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ] */
00071 
00072 /*     [ b11 b12 ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ] */
00073 /*     [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ], */
00074 
00075 /*  2) if the pencil (A,B) has a pair of complex conjugate eigenvalues, */
00076 /*     then */
00077 
00078 /*     [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ] */
00079 /*     [ a21 a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ] */
00080 
00081 /*     [ b11  0  ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ] */
00082 /*     [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ] */
00083 
00084 /*     where b11 >= b22 > 0. */
00085 
00086 
00087 /*  Arguments */
00088 /*  ========= */
00089 
00090 /*  A       (input/output) DOUBLE PRECISION array, dimension (LDA, 2) */
00091 /*          On entry, the 2 x 2 matrix A. */
00092 /*          On exit, A is overwritten by the ``A-part'' of the */
00093 /*          generalized Schur form. */
00094 
00095 /*  LDA     (input) INTEGER */
00096 /*          THe leading dimension of the array A.  LDA >= 2. */
00097 
00098 /*  B       (input/output) DOUBLE PRECISION array, dimension (LDB, 2) */
00099 /*          On entry, the upper triangular 2 x 2 matrix B. */
00100 /*          On exit, B is overwritten by the ``B-part'' of the */
00101 /*          generalized Schur form. */
00102 
00103 /*  LDB     (input) INTEGER */
00104 /*          THe leading dimension of the array B.  LDB >= 2. */
00105 
00106 /*  ALPHAR  (output) DOUBLE PRECISION array, dimension (2) */
00107 /*  ALPHAI  (output) DOUBLE PRECISION array, dimension (2) */
00108 /*  BETA    (output) DOUBLE PRECISION array, dimension (2) */
00109 /*          (ALPHAR(k)+i*ALPHAI(k))/BETA(k) are the eigenvalues of the */
00110 /*          pencil (A,B), k=1,2, i = sqrt(-1).  Note that BETA(k) may */
00111 /*          be zero. */
00112 
00113 /*  CSL     (output) DOUBLE PRECISION */
00114 /*          The cosine of the left rotation matrix. */
00115 
00116 /*  SNL     (output) DOUBLE PRECISION */
00117 /*          The sine of the left rotation matrix. */
00118 
00119 /*  CSR     (output) DOUBLE PRECISION */
00120 /*          The cosine of the right rotation matrix. */
00121 
00122 /*  SNR     (output) DOUBLE PRECISION */
00123 /*          The sine of the right rotation matrix. */
00124 
00125 /*  Further Details */
00126 /*  =============== */
00127 
00128 /*  Based on contributions by */
00129 /*     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA */
00130 
00131 /*  ===================================================================== */
00132 
00133 /*     .. Parameters .. */
00134 /*     .. */
00135 /*     .. Local Scalars .. */
00136 /*     .. */
00137 /*     .. External Subroutines .. */
00138 /*     .. */
00139 /*     .. External Functions .. */
00140 /*     .. */
00141 /*     .. Intrinsic Functions .. */
00142 /*     .. */
00143 /*     .. Executable Statements .. */
00144 
00145     /* Parameter adjustments */
00146     a_dim1 = *lda;
00147     a_offset = 1 + a_dim1;
00148     a -= a_offset;
00149     b_dim1 = *ldb;
00150     b_offset = 1 + b_dim1;
00151     b -= b_offset;
00152     --alphar;
00153     --alphai;
00154     --beta;
00155 
00156     /* Function Body */
00157     safmin = dlamch_("S");
00158     ulp = dlamch_("P");
00159 
00160 /*     Scale A */
00161 
00162 /* Computing MAX */
00163     d__5 = (d__1 = a[a_dim1 + 1], abs(d__1)) + (d__2 = a[a_dim1 + 2], abs(
00164             d__2)), d__6 = (d__3 = a[(a_dim1 << 1) + 1], abs(d__3)) + (d__4 = 
00165             a[(a_dim1 << 1) + 2], abs(d__4)), d__5 = max(d__5,d__6);
00166     anorm = max(d__5,safmin);
00167     ascale = 1. / anorm;
00168     a[a_dim1 + 1] = ascale * a[a_dim1 + 1];
00169     a[(a_dim1 << 1) + 1] = ascale * a[(a_dim1 << 1) + 1];
00170     a[a_dim1 + 2] = ascale * a[a_dim1 + 2];
00171     a[(a_dim1 << 1) + 2] = ascale * a[(a_dim1 << 1) + 2];
00172 
00173 /*     Scale B */
00174 
00175 /* Computing MAX */
00176     d__4 = (d__3 = b[b_dim1 + 1], abs(d__3)), d__5 = (d__1 = b[(b_dim1 << 1) 
00177             + 1], abs(d__1)) + (d__2 = b[(b_dim1 << 1) + 2], abs(d__2)), d__4 
00178             = max(d__4,d__5);
00179     bnorm = max(d__4,safmin);
00180     bscale = 1. / bnorm;
00181     b[b_dim1 + 1] = bscale * b[b_dim1 + 1];
00182     b[(b_dim1 << 1) + 1] = bscale * b[(b_dim1 << 1) + 1];
00183     b[(b_dim1 << 1) + 2] = bscale * b[(b_dim1 << 1) + 2];
00184 
00185 /*     Check if A can be deflated */
00186 
00187     if ((d__1 = a[a_dim1 + 2], abs(d__1)) <= ulp) {
00188         *csl = 1.;
00189         *snl = 0.;
00190         *csr = 1.;
00191         *snr = 0.;
00192         a[a_dim1 + 2] = 0.;
00193         b[b_dim1 + 2] = 0.;
00194 
00195 /*     Check if B is singular */
00196 
00197     } else if ((d__1 = b[b_dim1 + 1], abs(d__1)) <= ulp) {
00198         dlartg_(&a[a_dim1 + 1], &a[a_dim1 + 2], csl, snl, &r__);
00199         *csr = 1.;
00200         *snr = 0.;
00201         drot_(&c__2, &a[a_dim1 + 1], lda, &a[a_dim1 + 2], lda, csl, snl);
00202         drot_(&c__2, &b[b_dim1 + 1], ldb, &b[b_dim1 + 2], ldb, csl, snl);
00203         a[a_dim1 + 2] = 0.;
00204         b[b_dim1 + 1] = 0.;
00205         b[b_dim1 + 2] = 0.;
00206 
00207     } else if ((d__1 = b[(b_dim1 << 1) + 2], abs(d__1)) <= ulp) {
00208         dlartg_(&a[(a_dim1 << 1) + 2], &a[a_dim1 + 2], csr, snr, &t);
00209         *snr = -(*snr);
00210         drot_(&c__2, &a[a_dim1 + 1], &c__1, &a[(a_dim1 << 1) + 1], &c__1, csr, 
00211                  snr);
00212         drot_(&c__2, &b[b_dim1 + 1], &c__1, &b[(b_dim1 << 1) + 1], &c__1, csr, 
00213                  snr);
00214         *csl = 1.;
00215         *snl = 0.;
00216         a[a_dim1 + 2] = 0.;
00217         b[b_dim1 + 2] = 0.;
00218         b[(b_dim1 << 1) + 2] = 0.;
00219 
00220     } else {
00221 
00222 /*        B is nonsingular, first compute the eigenvalues of (A,B) */
00223 
00224         dlag2_(&a[a_offset], lda, &b[b_offset], ldb, &safmin, &scale1, &
00225                 scale2, &wr1, &wr2, &wi);
00226 
00227         if (wi == 0.) {
00228 
00229 /*           two real eigenvalues, compute s*A-w*B */
00230 
00231             h1 = scale1 * a[a_dim1 + 1] - wr1 * b[b_dim1 + 1];
00232             h2 = scale1 * a[(a_dim1 << 1) + 1] - wr1 * b[(b_dim1 << 1) + 1];
00233             h3 = scale1 * a[(a_dim1 << 1) + 2] - wr1 * b[(b_dim1 << 1) + 2];
00234 
00235             rr = dlapy2_(&h1, &h2);
00236             d__1 = scale1 * a[a_dim1 + 2];
00237             qq = dlapy2_(&d__1, &h3);
00238 
00239             if (rr > qq) {
00240 
00241 /*              find right rotation matrix to zero 1,1 element of */
00242 /*              (sA - wB) */
00243 
00244                 dlartg_(&h2, &h1, csr, snr, &t);
00245 
00246             } else {
00247 
00248 /*              find right rotation matrix to zero 2,1 element of */
00249 /*              (sA - wB) */
00250 
00251                 d__1 = scale1 * a[a_dim1 + 2];
00252                 dlartg_(&h3, &d__1, csr, snr, &t);
00253 
00254             }
00255 
00256             *snr = -(*snr);
00257             drot_(&c__2, &a[a_dim1 + 1], &c__1, &a[(a_dim1 << 1) + 1], &c__1, 
00258                     csr, snr);
00259             drot_(&c__2, &b[b_dim1 + 1], &c__1, &b[(b_dim1 << 1) + 1], &c__1, 
00260                     csr, snr);
00261 
00262 /*           compute inf norms of A and B */
00263 
00264 /* Computing MAX */
00265             d__5 = (d__1 = a[a_dim1 + 1], abs(d__1)) + (d__2 = a[(a_dim1 << 1)
00266                      + 1], abs(d__2)), d__6 = (d__3 = a[a_dim1 + 2], abs(d__3)
00267                     ) + (d__4 = a[(a_dim1 << 1) + 2], abs(d__4));
00268             h1 = max(d__5,d__6);
00269 /* Computing MAX */
00270             d__5 = (d__1 = b[b_dim1 + 1], abs(d__1)) + (d__2 = b[(b_dim1 << 1)
00271                      + 1], abs(d__2)), d__6 = (d__3 = b[b_dim1 + 2], abs(d__3)
00272                     ) + (d__4 = b[(b_dim1 << 1) + 2], abs(d__4));
00273             h2 = max(d__5,d__6);
00274 
00275             if (scale1 * h1 >= abs(wr1) * h2) {
00276 
00277 /*              find left rotation matrix Q to zero out B(2,1) */
00278 
00279                 dlartg_(&b[b_dim1 + 1], &b[b_dim1 + 2], csl, snl, &r__);
00280 
00281             } else {
00282 
00283 /*              find left rotation matrix Q to zero out A(2,1) */
00284 
00285                 dlartg_(&a[a_dim1 + 1], &a[a_dim1 + 2], csl, snl, &r__);
00286 
00287             }
00288 
00289             drot_(&c__2, &a[a_dim1 + 1], lda, &a[a_dim1 + 2], lda, csl, snl);
00290             drot_(&c__2, &b[b_dim1 + 1], ldb, &b[b_dim1 + 2], ldb, csl, snl);
00291 
00292             a[a_dim1 + 2] = 0.;
00293             b[b_dim1 + 2] = 0.;
00294 
00295         } else {
00296 
00297 /*           a pair of complex conjugate eigenvalues */
00298 /*           first compute the SVD of the matrix B */
00299 
00300             dlasv2_(&b[b_dim1 + 1], &b[(b_dim1 << 1) + 1], &b[(b_dim1 << 1) + 
00301                     2], &r__, &t, snr, csr, snl, csl);
00302 
00303 /*           Form (A,B) := Q(A,B)Z' where Q is left rotation matrix and */
00304 /*           Z is right rotation matrix computed from DLASV2 */
00305 
00306             drot_(&c__2, &a[a_dim1 + 1], lda, &a[a_dim1 + 2], lda, csl, snl);
00307             drot_(&c__2, &b[b_dim1 + 1], ldb, &b[b_dim1 + 2], ldb, csl, snl);
00308             drot_(&c__2, &a[a_dim1 + 1], &c__1, &a[(a_dim1 << 1) + 1], &c__1, 
00309                     csr, snr);
00310             drot_(&c__2, &b[b_dim1 + 1], &c__1, &b[(b_dim1 << 1) + 1], &c__1, 
00311                     csr, snr);
00312 
00313             b[b_dim1 + 2] = 0.;
00314             b[(b_dim1 << 1) + 1] = 0.;
00315 
00316         }
00317 
00318     }
00319 
00320 /*     Unscaling */
00321 
00322     a[a_dim1 + 1] = anorm * a[a_dim1 + 1];
00323     a[a_dim1 + 2] = anorm * a[a_dim1 + 2];
00324     a[(a_dim1 << 1) + 1] = anorm * a[(a_dim1 << 1) + 1];
00325     a[(a_dim1 << 1) + 2] = anorm * a[(a_dim1 << 1) + 2];
00326     b[b_dim1 + 1] = bnorm * b[b_dim1 + 1];
00327     b[b_dim1 + 2] = bnorm * b[b_dim1 + 2];
00328     b[(b_dim1 << 1) + 1] = bnorm * b[(b_dim1 << 1) + 1];
00329     b[(b_dim1 << 1) + 2] = bnorm * b[(b_dim1 << 1) + 2];
00330 
00331     if (wi == 0.) {
00332         alphar[1] = a[a_dim1 + 1];
00333         alphar[2] = a[(a_dim1 << 1) + 2];
00334         alphai[1] = 0.;
00335         alphai[2] = 0.;
00336         beta[1] = b[b_dim1 + 1];
00337         beta[2] = b[(b_dim1 << 1) + 2];
00338     } else {
00339         alphar[1] = anorm * wr1 / scale1 / bnorm;
00340         alphai[1] = anorm * wi / scale1 / bnorm;
00341         alphar[2] = alphar[1];
00342         alphai[2] = -alphai[1];
00343         beta[1] = 1.;
00344         beta[2] = 1.;
00345     }
00346 
00347     return 0;
00348 
00349 /*     End of DLAGV2 */
00350 
00351 } /* dlagv2_ */