stgexc.c
Go to the documentation of this file.
00001 /* stgexc.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__1 = 1;
00019 static integer c__2 = 2;
00020 
00021 /* Subroutine */ int stgexc_(logical *wantq, logical *wantz, integer *n, real 
00022         *a, integer *lda, real *b, integer *ldb, real *q, integer *ldq, real *
00023         z__, integer *ldz, integer *ifst, integer *ilst, real *work, integer *
00024         lwork, integer *info)
00025 {
00026     /* System generated locals */
00027     integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1, 
00028             z_offset, i__1;
00029 
00030     /* Local variables */
00031     integer nbf, nbl, here, lwmin;
00032     extern /* Subroutine */ int stgex2_(logical *, logical *, integer *, real 
00033             *, integer *, real *, integer *, real *, integer *, real *, 
00034             integer *, integer *, integer *, integer *, real *, integer *, 
00035             integer *), xerbla_(char *, integer *);
00036     integer nbnext;
00037     logical lquery;
00038 
00039 
00040 /*  -- LAPACK routine (version 3.2) -- */
00041 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00042 /*     November 2006 */
00043 
00044 /*     .. Scalar Arguments .. */
00045 /*     .. */
00046 /*     .. Array Arguments .. */
00047 /*     .. */
00048 
00049 /*  Purpose */
00050 /*  ======= */
00051 
00052 /*  STGEXC reorders the generalized real Schur decomposition of a real */
00053 /*  matrix pair (A,B) using an orthogonal equivalence transformation */
00054 
00055 /*                 (A, B) = Q * (A, B) * Z', */
00056 
00057 /*  so that the diagonal block of (A, B) with row index IFST is moved */
00058 /*  to row ILST. */
00059 
00060 /*  (A, B) must be in generalized real Schur canonical form (as returned */
00061 /*  by SGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2 */
00062 /*  diagonal blocks. B is upper triangular. */
00063 
00064 /*  Optionally, the matrices Q and Z of generalized Schur vectors are */
00065 /*  updated. */
00066 
00067 /*         Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)' */
00068 /*         Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)' */
00069 
00070 
00071 /*  Arguments */
00072 /*  ========= */
00073 
00074 /*  WANTQ   (input) LOGICAL */
00075 /*          .TRUE. : update the left transformation matrix Q; */
00076 /*          .FALSE.: do not update Q. */
00077 
00078 /*  WANTZ   (input) LOGICAL */
00079 /*          .TRUE. : update the right transformation matrix Z; */
00080 /*          .FALSE.: do not update Z. */
00081 
00082 /*  N       (input) INTEGER */
00083 /*          The order of the matrices A and B. N >= 0. */
00084 
00085 /*  A       (input/output) REAL array, dimension (LDA,N) */
00086 /*          On entry, the matrix A in generalized real Schur canonical */
00087 /*          form. */
00088 /*          On exit, the updated matrix A, again in generalized */
00089 /*          real Schur canonical form. */
00090 
00091 /*  LDA     (input)  INTEGER */
00092 /*          The leading dimension of the array A. LDA >= max(1,N). */
00093 
00094 /*  B       (input/output) REAL array, dimension (LDB,N) */
00095 /*          On entry, the matrix B in generalized real Schur canonical */
00096 /*          form (A,B). */
00097 /*          On exit, the updated matrix B, again in generalized */
00098 /*          real Schur canonical form (A,B). */
00099 
00100 /*  LDB     (input)  INTEGER */
00101 /*          The leading dimension of the array B. LDB >= max(1,N). */
00102 
00103 /*  Q       (input/output) REAL array, dimension (LDZ,N) */
00104 /*          On entry, if WANTQ = .TRUE., the orthogonal matrix Q. */
00105 /*          On exit, the updated matrix Q. */
00106 /*          If WANTQ = .FALSE., Q is not referenced. */
00107 
00108 /*  LDQ     (input) INTEGER */
00109 /*          The leading dimension of the array Q. LDQ >= 1. */
00110 /*          If WANTQ = .TRUE., LDQ >= N. */
00111 
00112 /*  Z       (input/output) REAL array, dimension (LDZ,N) */
00113 /*          On entry, if WANTZ = .TRUE., the orthogonal matrix Z. */
00114 /*          On exit, the updated matrix Z. */
00115 /*          If WANTZ = .FALSE., Z is not referenced. */
00116 
00117 /*  LDZ     (input) INTEGER */
00118 /*          The leading dimension of the array Z. LDZ >= 1. */
00119 /*          If WANTZ = .TRUE., LDZ >= N. */
00120 
00121 /*  IFST    (input/output) INTEGER */
00122 /*  ILST    (input/output) INTEGER */
00123 /*          Specify the reordering of the diagonal blocks of (A, B). */
00124 /*          The block with row index IFST is moved to row ILST, by a */
00125 /*          sequence of swapping between adjacent blocks. */
00126 /*          On exit, if IFST pointed on entry to the second row of */
00127 /*          a 2-by-2 block, it is changed to point to the first row; */
00128 /*          ILST always points to the first row of the block in its */
00129 /*          final position (which may differ from its input value by */
00130 /*          +1 or -1). 1 <= IFST, ILST <= N. */
00131 
00132 /*  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK)) */
00133 /*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
00134 
00135 /*  LWORK   (input) INTEGER */
00136 /*          The dimension of the array WORK. */
00137 /*          LWORK >= 1 when N <= 1, otherwise LWORK >= 4*N + 16. */
00138 
00139 /*          If LWORK = -1, then a workspace query is assumed; the routine */
00140 /*          only calculates the optimal size of the WORK array, returns */
00141 /*          this value as the first entry of the WORK array, and no error */
00142 /*          message related to LWORK is issued by XERBLA. */
00143 
00144 /*  INFO    (output) INTEGER */
00145 /*           =0:  successful exit. */
00146 /*           <0:  if INFO = -i, the i-th argument had an illegal value. */
00147 /*           =1:  The transformed matrix pair (A, B) would be too far */
00148 /*                from generalized Schur form; the problem is ill- */
00149 /*                conditioned. (A, B) may have been partially reordered, */
00150 /*                and ILST points to the first row of the current */
00151 /*                position of the block being moved. */
00152 
00153 /*  Further Details */
00154 /*  =============== */
00155 
00156 /*  Based on contributions by */
00157 /*     Bo Kagstrom and Peter Poromaa, Department of Computing Science, */
00158 /*     Umea University, S-901 87 Umea, Sweden. */
00159 
00160 /*  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the */
00161 /*      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in */
00162 /*      M.S. Moonen et al (eds), Linear Algebra for Large Scale and */
00163 /*      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. */
00164 
00165 /*  ===================================================================== */
00166 
00167 /*     .. Parameters .. */
00168 /*     .. */
00169 /*     .. Local Scalars .. */
00170 /*     .. */
00171 /*     .. External Subroutines .. */
00172 /*     .. */
00173 /*     .. Intrinsic Functions .. */
00174 /*     .. */
00175 /*     .. Executable Statements .. */
00176 
00177 /*     Decode and test input arguments. */
00178 
00179     /* Parameter adjustments */
00180     a_dim1 = *lda;
00181     a_offset = 1 + a_dim1;
00182     a -= a_offset;
00183     b_dim1 = *ldb;
00184     b_offset = 1 + b_dim1;
00185     b -= b_offset;
00186     q_dim1 = *ldq;
00187     q_offset = 1 + q_dim1;
00188     q -= q_offset;
00189     z_dim1 = *ldz;
00190     z_offset = 1 + z_dim1;
00191     z__ -= z_offset;
00192     --work;
00193 
00194     /* Function Body */
00195     *info = 0;
00196     lquery = *lwork == -1;
00197     if (*n < 0) {
00198         *info = -3;
00199     } else if (*lda < max(1,*n)) {
00200         *info = -5;
00201     } else if (*ldb < max(1,*n)) {
00202         *info = -7;
00203     } else if (*ldq < 1 || *wantq && *ldq < max(1,*n)) {
00204         *info = -9;
00205     } else if (*ldz < 1 || *wantz && *ldz < max(1,*n)) {
00206         *info = -11;
00207     } else if (*ifst < 1 || *ifst > *n) {
00208         *info = -12;
00209     } else if (*ilst < 1 || *ilst > *n) {
00210         *info = -13;
00211     }
00212 
00213     if (*info == 0) {
00214         if (*n <= 1) {
00215             lwmin = 1;
00216         } else {
00217             lwmin = (*n << 2) + 16;
00218         }
00219         work[1] = (real) lwmin;
00220 
00221         if (*lwork < lwmin && ! lquery) {
00222             *info = -15;
00223         }
00224     }
00225 
00226     if (*info != 0) {
00227         i__1 = -(*info);
00228         xerbla_("STGEXC", &i__1);
00229         return 0;
00230     } else if (lquery) {
00231         return 0;
00232     }
00233 
00234 /*     Quick return if possible */
00235 
00236     if (*n <= 1) {
00237         return 0;
00238     }
00239 
00240 /*     Determine the first row of the specified block and find out */
00241 /*     if it is 1-by-1 or 2-by-2. */
00242 
00243     if (*ifst > 1) {
00244         if (a[*ifst + (*ifst - 1) * a_dim1] != 0.f) {
00245             --(*ifst);
00246         }
00247     }
00248     nbf = 1;
00249     if (*ifst < *n) {
00250         if (a[*ifst + 1 + *ifst * a_dim1] != 0.f) {
00251             nbf = 2;
00252         }
00253     }
00254 
00255 /*     Determine the first row of the final block */
00256 /*     and find out if it is 1-by-1 or 2-by-2. */
00257 
00258     if (*ilst > 1) {
00259         if (a[*ilst + (*ilst - 1) * a_dim1] != 0.f) {
00260             --(*ilst);
00261         }
00262     }
00263     nbl = 1;
00264     if (*ilst < *n) {
00265         if (a[*ilst + 1 + *ilst * a_dim1] != 0.f) {
00266             nbl = 2;
00267         }
00268     }
00269     if (*ifst == *ilst) {
00270         return 0;
00271     }
00272 
00273     if (*ifst < *ilst) {
00274 
00275 /*        Update ILST. */
00276 
00277         if (nbf == 2 && nbl == 1) {
00278             --(*ilst);
00279         }
00280         if (nbf == 1 && nbl == 2) {
00281             ++(*ilst);
00282         }
00283 
00284         here = *ifst;
00285 
00286 L10:
00287 
00288 /*        Swap with next one below. */
00289 
00290         if (nbf == 1 || nbf == 2) {
00291 
00292 /*           Current block either 1-by-1 or 2-by-2. */
00293 
00294             nbnext = 1;
00295             if (here + nbf + 1 <= *n) {
00296                 if (a[here + nbf + 1 + (here + nbf) * a_dim1] != 0.f) {
00297                     nbnext = 2;
00298                 }
00299             }
00300             stgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[
00301                     q_offset], ldq, &z__[z_offset], ldz, &here, &nbf, &nbnext, 
00302                      &work[1], lwork, info);
00303             if (*info != 0) {
00304                 *ilst = here;
00305                 return 0;
00306             }
00307             here += nbnext;
00308 
00309 /*           Test if 2-by-2 block breaks into two 1-by-1 blocks. */
00310 
00311             if (nbf == 2) {
00312                 if (a[here + 1 + here * a_dim1] == 0.f) {
00313                     nbf = 3;
00314                 }
00315             }
00316 
00317         } else {
00318 
00319 /*           Current block consists of two 1-by-1 blocks, each of which */
00320 /*           must be swapped individually. */
00321 
00322             nbnext = 1;
00323             if (here + 3 <= *n) {
00324                 if (a[here + 3 + (here + 2) * a_dim1] != 0.f) {
00325                     nbnext = 2;
00326                 }
00327             }
00328             i__1 = here + 1;
00329             stgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[
00330                     q_offset], ldq, &z__[z_offset], ldz, &i__1, &c__1, &
00331                     nbnext, &work[1], lwork, info);
00332             if (*info != 0) {
00333                 *ilst = here;
00334                 return 0;
00335             }
00336             if (nbnext == 1) {
00337 
00338 /*              Swap two 1-by-1 blocks. */
00339 
00340                 stgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, 
00341                          &q[q_offset], ldq, &z__[z_offset], ldz, &here, &c__1, 
00342                          &c__1, &work[1], lwork, info);
00343                 if (*info != 0) {
00344                     *ilst = here;
00345                     return 0;
00346                 }
00347                 ++here;
00348 
00349             } else {
00350 
00351 /*              Recompute NBNEXT in case of 2-by-2 split. */
00352 
00353                 if (a[here + 2 + (here + 1) * a_dim1] == 0.f) {
00354                     nbnext = 1;
00355                 }
00356                 if (nbnext == 2) {
00357 
00358 /*                 2-by-2 block did not split. */
00359 
00360                     stgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], 
00361                             ldb, &q[q_offset], ldq, &z__[z_offset], ldz, &
00362                             here, &c__1, &nbnext, &work[1], lwork, info);
00363                     if (*info != 0) {
00364                         *ilst = here;
00365                         return 0;
00366                     }
00367                     here += 2;
00368                 } else {
00369 
00370 /*                 2-by-2 block did split. */
00371 
00372                     stgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], 
00373                             ldb, &q[q_offset], ldq, &z__[z_offset], ldz, &
00374                             here, &c__1, &c__1, &work[1], lwork, info);
00375                     if (*info != 0) {
00376                         *ilst = here;
00377                         return 0;
00378                     }
00379                     ++here;
00380                     stgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], 
00381                             ldb, &q[q_offset], ldq, &z__[z_offset], ldz, &
00382                             here, &c__1, &c__1, &work[1], lwork, info);
00383                     if (*info != 0) {
00384                         *ilst = here;
00385                         return 0;
00386                     }
00387                     ++here;
00388                 }
00389 
00390             }
00391         }
00392         if (here < *ilst) {
00393             goto L10;
00394         }
00395     } else {
00396         here = *ifst;
00397 
00398 L20:
00399 
00400 /*        Swap with next one below. */
00401 
00402         if (nbf == 1 || nbf == 2) {
00403 
00404 /*           Current block either 1-by-1 or 2-by-2. */
00405 
00406             nbnext = 1;
00407             if (here >= 3) {
00408                 if (a[here - 1 + (here - 2) * a_dim1] != 0.f) {
00409                     nbnext = 2;
00410                 }
00411             }
00412             i__1 = here - nbnext;
00413             stgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[
00414                     q_offset], ldq, &z__[z_offset], ldz, &i__1, &nbnext, &nbf, 
00415                      &work[1], lwork, info);
00416             if (*info != 0) {
00417                 *ilst = here;
00418                 return 0;
00419             }
00420             here -= nbnext;
00421 
00422 /*           Test if 2-by-2 block breaks into two 1-by-1 blocks. */
00423 
00424             if (nbf == 2) {
00425                 if (a[here + 1 + here * a_dim1] == 0.f) {
00426                     nbf = 3;
00427                 }
00428             }
00429 
00430         } else {
00431 
00432 /*           Current block consists of two 1-by-1 blocks, each of which */
00433 /*           must be swapped individually. */
00434 
00435             nbnext = 1;
00436             if (here >= 3) {
00437                 if (a[here - 1 + (here - 2) * a_dim1] != 0.f) {
00438                     nbnext = 2;
00439                 }
00440             }
00441             i__1 = here - nbnext;
00442             stgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[
00443                     q_offset], ldq, &z__[z_offset], ldz, &i__1, &nbnext, &
00444                     c__1, &work[1], lwork, info);
00445             if (*info != 0) {
00446                 *ilst = here;
00447                 return 0;
00448             }
00449             if (nbnext == 1) {
00450 
00451 /*              Swap two 1-by-1 blocks. */
00452 
00453                 stgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, 
00454                          &q[q_offset], ldq, &z__[z_offset], ldz, &here, &
00455                         nbnext, &c__1, &work[1], lwork, info);
00456                 if (*info != 0) {
00457                     *ilst = here;
00458                     return 0;
00459                 }
00460                 --here;
00461             } else {
00462 
00463 /*             Recompute NBNEXT in case of 2-by-2 split. */
00464 
00465                 if (a[here + (here - 1) * a_dim1] == 0.f) {
00466                     nbnext = 1;
00467                 }
00468                 if (nbnext == 2) {
00469 
00470 /*                 2-by-2 block did not split. */
00471 
00472                     i__1 = here - 1;
00473                     stgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], 
00474                             ldb, &q[q_offset], ldq, &z__[z_offset], ldz, &
00475                             i__1, &c__2, &c__1, &work[1], lwork, info);
00476                     if (*info != 0) {
00477                         *ilst = here;
00478                         return 0;
00479                     }
00480                     here += -2;
00481                 } else {
00482 
00483 /*                 2-by-2 block did split. */
00484 
00485                     stgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], 
00486                             ldb, &q[q_offset], ldq, &z__[z_offset], ldz, &
00487                             here, &c__1, &c__1, &work[1], lwork, info);
00488                     if (*info != 0) {
00489                         *ilst = here;
00490                         return 0;
00491                     }
00492                     --here;
00493                     stgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], 
00494                             ldb, &q[q_offset], ldq, &z__[z_offset], ldz, &
00495                             here, &c__1, &c__1, &work[1], lwork, info);
00496                     if (*info != 0) {
00497                         *ilst = here;
00498                         return 0;
00499                     }
00500                     --here;
00501                 }
00502             }
00503         }
00504         if (here > *ilst) {
00505             goto L20;
00506         }
00507     }
00508     *ilst = here;
00509     work[1] = (real) lwmin;
00510     return 0;
00511 
00512 /*     End of STGEXC */
00513 
00514 } /* stgexc_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:14