swiftnav: ctrsyl.c Source File

Go to the documentation of this file.
00001 /* ctrsyl.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 
00020 /* Subroutine */ int ctrsyl_(char *trana, char *tranb, integer *isgn, integer 
00021         *m, integer *n, complex *a, integer *lda, complex *b, integer *ldb, 
00022         complex *c__, integer *ldc, real *scale, integer *info)
00023 {
00024     /* System generated locals */
00025     integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2, 
00026             i__3, i__4;
00027     real r__1, r__2;
00028     complex q__1, q__2, q__3, q__4;
00029 
00030     /* Builtin functions */
00031     double r_imag(complex *);
00032     void r_cnjg(complex *, complex *);
00033 
00034     /* Local variables */
00035     integer j, k, l;
00036     complex a11;
00037     real db;
00038     complex x11;
00039     real da11;
00040     complex vec;
00041     real dum[1], eps, sgn, smin;
00042     complex suml, sumr;
00043     extern /* Complex */ VOID cdotc_(complex *, integer *, complex *, integer 
00044             *, complex *, integer *);
00045     extern logical lsame_(char *, char *);
00046     extern /* Complex */ VOID cdotu_(complex *, integer *, complex *, integer 
00047             *, complex *, integer *);
00048     extern /* Subroutine */ int slabad_(real *, real *);
00049     extern doublereal clange_(char *, integer *, integer *, complex *, 
00050             integer *, real *);
00051     extern /* Complex */ VOID cladiv_(complex *, complex *, complex *);
00052     real scaloc;
00053     extern doublereal slamch_(char *);
00054     extern /* Subroutine */ int csscal_(integer *, real *, complex *, integer 
00055             *), xerbla_(char *, integer *);
00056     real bignum;
00057     logical notrna, notrnb;
00058     real smlnum;
00059 
00060 
00061 /*  -- LAPACK routine (version 3.2) -- */
00062 /*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */
00063 /*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
00064 /*     November 2006 */
00065 
00066 /*     .. Scalar Arguments .. */
00067 /*     .. */
00068 /*     .. Array Arguments .. */
00069 /*     .. */
00070 
00071 /*  Purpose */
00072 /*  ======= */
00073 
00074 /*  CTRSYL solves the complex Sylvester matrix equation: */
00075 
00076 /*     op(A)*X + X*op(B) = scale*C or */
00077 /*     op(A)*X - X*op(B) = scale*C, */
00078 
00079 /*  where op(A) = A or A**H, and A and B are both upper triangular. A is */
00080 /*  M-by-M and B is N-by-N; the right hand side C and the solution X are */
00081 /*  M-by-N; and scale is an output scale factor, set <= 1 to avoid */
00082 /*  overflow in X. */
00083 
00084 /*  Arguments */
00085 /*  ========= */
00086 
00087 /*  TRANA   (input) CHARACTER*1 */
00088 /*          Specifies the option op(A): */
00089 /*          = 'N': op(A) = A    (No transpose) */
00090 /*          = 'C': op(A) = A**H (Conjugate transpose) */
00091 
00092 /*  TRANB   (input) CHARACTER*1 */
00093 /*          Specifies the option op(B): */
00094 /*          = 'N': op(B) = B    (No transpose) */
00095 /*          = 'C': op(B) = B**H (Conjugate transpose) */
00096 
00097 /*  ISGN    (input) INTEGER */
00098 /*          Specifies the sign in the equation: */
00099 /*          = +1: solve op(A)*X + X*op(B) = scale*C */
00100 /*          = -1: solve op(A)*X - X*op(B) = scale*C */
00101 
00102 /*  M       (input) INTEGER */
00103 /*          The order of the matrix A, and the number of rows in the */
00104 /*          matrices X and C. M >= 0. */
00105 
00106 /*  N       (input) INTEGER */
00107 /*          The order of the matrix B, and the number of columns in the */
00108 /*          matrices X and C. N >= 0. */
00109 
00110 /*  A       (input) COMPLEX array, dimension (LDA,M) */
00111 /*          The upper triangular matrix A. */
00112 
00113 /*  LDA     (input) INTEGER */
00114 /*          The leading dimension of the array A. LDA >= max(1,M). */
00115 
00116 /*  B       (input) COMPLEX array, dimension (LDB,N) */
00117 /*          The upper triangular matrix B. */
00118 
00119 /*  LDB     (input) INTEGER */
00120 /*          The leading dimension of the array B. LDB >= max(1,N). */
00121 
00122 /*  C       (input/output) COMPLEX array, dimension (LDC,N) */
00123 /*          On entry, the M-by-N right hand side matrix C. */
00124 /*          On exit, C is overwritten by the solution matrix X. */
00125 
00126 /*  LDC     (input) INTEGER */
00127 /*          The leading dimension of the array C. LDC >= max(1,M) */
00128 
00129 /*  SCALE   (output) REAL */
00130 /*          The scale factor, scale, set <= 1 to avoid overflow in X. */
00131 
00132 /*  INFO    (output) INTEGER */
00133 /*          = 0: successful exit */
00134 /*          < 0: if INFO = -i, the i-th argument had an illegal value */
00135 /*          = 1: A and B have common or very close eigenvalues; perturbed */
00136 /*               values were used to solve the equation (but the matrices */
00137 /*               A and B are unchanged). */
00138 
00139 /*  ===================================================================== */
00140 
00141 /*     .. Parameters .. */
00142 /*     .. */
00143 /*     .. Local Scalars .. */
00144 /*     .. */
00145 /*     .. Local Arrays .. */
00146 /*     .. */
00147 /*     .. External Functions .. */
00148 /*     .. */
00149 /*     .. External Subroutines .. */
00150 /*     .. */
00151 /*     .. Intrinsic Functions .. */
00152 /*     .. */
00153 /*     .. Executable Statements .. */
00154 
00155 /*     Decode and Test input parameters */
00156 
00157     /* Parameter adjustments */
00158     a_dim1 = *lda;
00159     a_offset = 1 + a_dim1;
00160     a -= a_offset;
00161     b_dim1 = *ldb;
00162     b_offset = 1 + b_dim1;
00163     b -= b_offset;
00164     c_dim1 = *ldc;
00165     c_offset = 1 + c_dim1;
00166     c__ -= c_offset;
00167 
00168     /* Function Body */
00169     notrna = lsame_(trana, "N");
00170     notrnb = lsame_(tranb, "N");
00171 
00172     *info = 0;
00173     if (! notrna && ! lsame_(trana, "C")) {
00174         *info = -1;
00175     } else if (! notrnb && ! lsame_(tranb, "C")) {
00176         *info = -2;
00177     } else if (*isgn != 1 && *isgn != -1) {
00178         *info = -3;
00179     } else if (*m < 0) {
00180         *info = -4;
00181     } else if (*n < 0) {
00182         *info = -5;
00183     } else if (*lda < max(1,*m)) {
00184         *info = -7;
00185     } else if (*ldb < max(1,*n)) {
00186         *info = -9;
00187     } else if (*ldc < max(1,*m)) {
00188         *info = -11;
00189     }
00190     if (*info != 0) {
00191         i__1 = -(*info);
00192         xerbla_("CTRSYL", &i__1);
00193         return 0;
00194     }
00195 
00196 /*     Quick return if possible */
00197 
00198     *scale = 1.f;
00199     if (*m == 0 || *n == 0) {
00200         return 0;
00201     }
00202 
00203 /*     Set constants to control overflow */
00204 
00205     eps = slamch_("P");
00206     smlnum = slamch_("S");
00207     bignum = 1.f / smlnum;
00208     slabad_(&smlnum, &bignum);
00209     smlnum = smlnum * (real) (*m * *n) / eps;
00210     bignum = 1.f / smlnum;
00211 /* Computing MAX */
00212     r__1 = smlnum, r__2 = eps * clange_("M", m, m, &a[a_offset], lda, dum), r__1 = max(r__1,r__2), r__2 = eps * clange_("M", n, n, 
00213             &b[b_offset], ldb, dum);
00214     smin = dmax(r__1,r__2);
00215     sgn = (real) (*isgn);
00216 
00217     if (notrna && notrnb) {
00218 
00219 /*        Solve    A*X + ISGN*X*B = scale*C. */
00220 
00221 /*        The (K,L)th block of X is determined starting from */
00222 /*        bottom-left corner column by column by */
00223 
00224 /*            A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) */
00225 
00226 /*        Where */
00227 /*                    M                        L-1 */
00228 /*          R(K,L) = SUM [A(K,I)*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]. */
00229 /*                  I=K+1                      J=1 */
00230 
00231         i__1 = *n;
00232         for (l = 1; l <= i__1; ++l) {
00233             for (k = *m; k >= 1; --k) {
00234 
00235                 i__2 = *m - k;
00236 /* Computing MIN */
00237                 i__3 = k + 1;
00238 /* Computing MIN */
00239                 i__4 = k + 1;
00240                 cdotu_(&q__1, &i__2, &a[k + min(i__3, *m)* a_dim1], lda, &c__[
00241                         min(i__4, *m)+ l * c_dim1], &c__1);
00242                 suml.r = q__1.r, suml.i = q__1.i;
00243                 i__2 = l - 1;
00244                 cdotu_(&q__1, &i__2, &c__[k + c_dim1], ldc, &b[l * b_dim1 + 1]
00245 , &c__1);
00246                 sumr.r = q__1.r, sumr.i = q__1.i;
00247                 i__2 = k + l * c_dim1;
00248                 q__3.r = sgn * sumr.r, q__3.i = sgn * sumr.i;
00249                 q__2.r = suml.r + q__3.r, q__2.i = suml.i + q__3.i;
00250                 q__1.r = c__[i__2].r - q__2.r, q__1.i = c__[i__2].i - q__2.i;
00251                 vec.r = q__1.r, vec.i = q__1.i;
00252 
00253                 scaloc = 1.f;
00254                 i__2 = k + k * a_dim1;
00255                 i__3 = l + l * b_dim1;
00256                 q__2.r = sgn * b[i__3].r, q__2.i = sgn * b[i__3].i;
00257                 q__1.r = a[i__2].r + q__2.r, q__1.i = a[i__2].i + q__2.i;
00258                 a11.r = q__1.r, a11.i = q__1.i;
00259                 da11 = (r__1 = a11.r, dabs(r__1)) + (r__2 = r_imag(&a11), 
00260                         dabs(r__2));
00261                 if (da11 <= smin) {
00262                     a11.r = smin, a11.i = 0.f;
00263                     da11 = smin;
00264                     *info = 1;
00265                 }
00266                 db = (r__1 = vec.r, dabs(r__1)) + (r__2 = r_imag(&vec), dabs(
00267                         r__2));
00268                 if (da11 < 1.f && db > 1.f) {
00269                     if (db > bignum * da11) {
00270                         scaloc = 1.f / db;
00271                     }
00272                 }
00273                 q__3.r = scaloc, q__3.i = 0.f;
00274                 q__2.r = vec.r * q__3.r - vec.i * q__3.i, q__2.i = vec.r * 
00275                         q__3.i + vec.i * q__3.r;
00276                 cladiv_(&q__1, &q__2, &a11);
00277                 x11.r = q__1.r, x11.i = q__1.i;
00278 
00279                 if (scaloc != 1.f) {
00280                     i__2 = *n;
00281                     for (j = 1; j <= i__2; ++j) {
00282                         csscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
00283 /* L10: */
00284                     }
00285                     *scale *= scaloc;
00286                 }
00287                 i__2 = k + l * c_dim1;
00288                 c__[i__2].r = x11.r, c__[i__2].i = x11.i;
00289 
00290 /* L20: */
00291             }
00292 /* L30: */
00293         }
00294 
00295     } else if (! notrna && notrnb) {
00296 
00297 /*        Solve    A' *X + ISGN*X*B = scale*C. */
00298 
00299 /*        The (K,L)th block of X is determined starting from */
00300 /*        upper-left corner column by column by */
00301 
00302 /*            A'(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) */
00303 
00304 /*        Where */
00305 /*                   K-1                         L-1 */
00306 /*          R(K,L) = SUM [A'(I,K)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)] */
00307 /*                   I=1                         J=1 */
00308 
00309         i__1 = *n;
00310         for (l = 1; l <= i__1; ++l) {
00311             i__2 = *m;
00312             for (k = 1; k <= i__2; ++k) {
00313 
00314                 i__3 = k - 1;
00315                 cdotc_(&q__1, &i__3, &a[k * a_dim1 + 1], &c__1, &c__[l * 
00316                         c_dim1 + 1], &c__1);
00317                 suml.r = q__1.r, suml.i = q__1.i;
00318                 i__3 = l - 1;
00319                 cdotu_(&q__1, &i__3, &c__[k + c_dim1], ldc, &b[l * b_dim1 + 1]
00320 , &c__1);
00321                 sumr.r = q__1.r, sumr.i = q__1.i;
00322                 i__3 = k + l * c_dim1;
00323                 q__3.r = sgn * sumr.r, q__3.i = sgn * sumr.i;
00324                 q__2.r = suml.r + q__3.r, q__2.i = suml.i + q__3.i;
00325                 q__1.r = c__[i__3].r - q__2.r, q__1.i = c__[i__3].i - q__2.i;
00326                 vec.r = q__1.r, vec.i = q__1.i;
00327 
00328                 scaloc = 1.f;
00329                 r_cnjg(&q__2, &a[k + k * a_dim1]);
00330                 i__3 = l + l * b_dim1;
00331                 q__3.r = sgn * b[i__3].r, q__3.i = sgn * b[i__3].i;
00332                 q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
00333                 a11.r = q__1.r, a11.i = q__1.i;
00334                 da11 = (r__1 = a11.r, dabs(r__1)) + (r__2 = r_imag(&a11), 
00335                         dabs(r__2));
00336                 if (da11 <= smin) {
00337                     a11.r = smin, a11.i = 0.f;
00338                     da11 = smin;
00339                     *info = 1;
00340                 }
00341                 db = (r__1 = vec.r, dabs(r__1)) + (r__2 = r_imag(&vec), dabs(
00342                         r__2));
00343                 if (da11 < 1.f && db > 1.f) {
00344                     if (db > bignum * da11) {
00345                         scaloc = 1.f / db;
00346                     }
00347                 }
00348 
00349                 q__3.r = scaloc, q__3.i = 0.f;
00350                 q__2.r = vec.r * q__3.r - vec.i * q__3.i, q__2.i = vec.r * 
00351                         q__3.i + vec.i * q__3.r;
00352                 cladiv_(&q__1, &q__2, &a11);
00353                 x11.r = q__1.r, x11.i = q__1.i;
00354 
00355                 if (scaloc != 1.f) {
00356                     i__3 = *n;
00357                     for (j = 1; j <= i__3; ++j) {
00358                         csscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
00359 /* L40: */
00360                     }
00361                     *scale *= scaloc;
00362                 }
00363                 i__3 = k + l * c_dim1;
00364                 c__[i__3].r = x11.r, c__[i__3].i = x11.i;
00365 
00366 /* L50: */
00367             }
00368 /* L60: */
00369         }
00370 
00371     } else if (! notrna && ! notrnb) {
00372 
00373 /*        Solve    A'*X + ISGN*X*B' = C. */
00374 
00375 /*        The (K,L)th block of X is determined starting from */
00376 /*        upper-right corner column by column by */
00377 
00378 /*            A'(K,K)*X(K,L) + ISGN*X(K,L)*B'(L,L) = C(K,L) - R(K,L) */
00379 
00380 /*        Where */
00381 /*                    K-1 */
00382 /*           R(K,L) = SUM [A'(I,K)*X(I,L)] + */
00383 /*                    I=1 */
00384 /*                           N */
00385 /*                     ISGN*SUM [X(K,J)*B'(L,J)]. */
00386 /*                          J=L+1 */
00387 
00388         for (l = *n; l >= 1; --l) {
00389             i__1 = *m;
00390             for (k = 1; k <= i__1; ++k) {
00391 
00392                 i__2 = k - 1;
00393                 cdotc_(&q__1, &i__2, &a[k * a_dim1 + 1], &c__1, &c__[l * 
00394                         c_dim1 + 1], &c__1);
00395                 suml.r = q__1.r, suml.i = q__1.i;
00396                 i__2 = *n - l;
00397 /* Computing MIN */
00398                 i__3 = l + 1;
00399 /* Computing MIN */
00400                 i__4 = l + 1;
00401                 cdotc_(&q__1, &i__2, &c__[k + min(i__3, *n)* c_dim1], ldc, &b[
00402                         l + min(i__4, *n)* b_dim1], ldb);
00403                 sumr.r = q__1.r, sumr.i = q__1.i;
00404                 i__2 = k + l * c_dim1;
00405                 r_cnjg(&q__4, &sumr);
00406                 q__3.r = sgn * q__4.r, q__3.i = sgn * q__4.i;
00407                 q__2.r = suml.r + q__3.r, q__2.i = suml.i + q__3.i;
00408                 q__1.r = c__[i__2].r - q__2.r, q__1.i = c__[i__2].i - q__2.i;
00409                 vec.r = q__1.r, vec.i = q__1.i;
00410 
00411                 scaloc = 1.f;
00412                 i__2 = k + k * a_dim1;
00413                 i__3 = l + l * b_dim1;
00414                 q__3.r = sgn * b[i__3].r, q__3.i = sgn * b[i__3].i;
00415                 q__2.r = a[i__2].r + q__3.r, q__2.i = a[i__2].i + q__3.i;
00416                 r_cnjg(&q__1, &q__2);
00417                 a11.r = q__1.r, a11.i = q__1.i;
00418                 da11 = (r__1 = a11.r, dabs(r__1)) + (r__2 = r_imag(&a11), 
00419                         dabs(r__2));
00420                 if (da11 <= smin) {
00421                     a11.r = smin, a11.i = 0.f;
00422                     da11 = smin;
00423                     *info = 1;
00424                 }
00425                 db = (r__1 = vec.r, dabs(r__1)) + (r__2 = r_imag(&vec), dabs(
00426                         r__2));
00427                 if (da11 < 1.f && db > 1.f) {
00428                     if (db > bignum * da11) {
00429                         scaloc = 1.f / db;
00430                     }
00431                 }
00432 
00433                 q__3.r = scaloc, q__3.i = 0.f;
00434                 q__2.r = vec.r * q__3.r - vec.i * q__3.i, q__2.i = vec.r * 
00435                         q__3.i + vec.i * q__3.r;
00436                 cladiv_(&q__1, &q__2, &a11);
00437                 x11.r = q__1.r, x11.i = q__1.i;
00438 
00439                 if (scaloc != 1.f) {
00440                     i__2 = *n;
00441                     for (j = 1; j <= i__2; ++j) {
00442                         csscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
00443 /* L70: */
00444                     }
00445                     *scale *= scaloc;
00446                 }
00447                 i__2 = k + l * c_dim1;
00448                 c__[i__2].r = x11.r, c__[i__2].i = x11.i;
00449 
00450 /* L80: */
00451             }
00452 /* L90: */
00453         }
00454 
00455     } else if (notrna && ! notrnb) {
00456 
00457 /*        Solve    A*X + ISGN*X*B' = C. */
00458 
00459 /*        The (K,L)th block of X is determined starting from */
00460 /*        bottom-left corner column by column by */
00461 
00462 /*           A(K,K)*X(K,L) + ISGN*X(K,L)*B'(L,L) = C(K,L) - R(K,L) */
00463 
00464 /*        Where */
00465 /*                    M                          N */
00466 /*          R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B'(L,J)] */
00467 /*                  I=K+1                      J=L+1 */
00468 
00469         for (l = *n; l >= 1; --l) {
00470             for (k = *m; k >= 1; --k) {
00471 
00472                 i__1 = *m - k;
00473 /* Computing MIN */
00474                 i__2 = k + 1;
00475 /* Computing MIN */
00476                 i__3 = k + 1;
00477                 cdotu_(&q__1, &i__1, &a[k + min(i__2, *m)* a_dim1], lda, &c__[
00478                         min(i__3, *m)+ l * c_dim1], &c__1);
00479                 suml.r = q__1.r, suml.i = q__1.i;
00480                 i__1 = *n - l;
00481 /* Computing MIN */
00482                 i__2 = l + 1;
00483 /* Computing MIN */
00484                 i__3 = l + 1;
00485                 cdotc_(&q__1, &i__1, &c__[k + min(i__2, *n)* c_dim1], ldc, &b[
00486                         l + min(i__3, *n)* b_dim1], ldb);
00487                 sumr.r = q__1.r, sumr.i = q__1.i;
00488                 i__1 = k + l * c_dim1;
00489                 r_cnjg(&q__4, &sumr);
00490                 q__3.r = sgn * q__4.r, q__3.i = sgn * q__4.i;
00491                 q__2.r = suml.r + q__3.r, q__2.i = suml.i + q__3.i;
00492                 q__1.r = c__[i__1].r - q__2.r, q__1.i = c__[i__1].i - q__2.i;
00493                 vec.r = q__1.r, vec.i = q__1.i;
00494 
00495                 scaloc = 1.f;
00496                 i__1 = k + k * a_dim1;
00497                 r_cnjg(&q__3, &b[l + l * b_dim1]);
00498                 q__2.r = sgn * q__3.r, q__2.i = sgn * q__3.i;
00499                 q__1.r = a[i__1].r + q__2.r, q__1.i = a[i__1].i + q__2.i;
00500                 a11.r = q__1.r, a11.i = q__1.i;
00501                 da11 = (r__1 = a11.r, dabs(r__1)) + (r__2 = r_imag(&a11), 
00502                         dabs(r__2));
00503                 if (da11 <= smin) {
00504                     a11.r = smin, a11.i = 0.f;
00505                     da11 = smin;
00506                     *info = 1;
00507                 }
00508                 db = (r__1 = vec.r, dabs(r__1)) + (r__2 = r_imag(&vec), dabs(
00509                         r__2));
00510                 if (da11 < 1.f && db > 1.f) {
00511                     if (db > bignum * da11) {
00512                         scaloc = 1.f / db;
00513                     }
00514                 }
00515 
00516                 q__3.r = scaloc, q__3.i = 0.f;
00517                 q__2.r = vec.r * q__3.r - vec.i * q__3.i, q__2.i = vec.r * 
00518                         q__3.i + vec.i * q__3.r;
00519                 cladiv_(&q__1, &q__2, &a11);
00520                 x11.r = q__1.r, x11.i = q__1.i;
00521 
00522                 if (scaloc != 1.f) {
00523                     i__1 = *n;
00524                     for (j = 1; j <= i__1; ++j) {
00525                         csscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
00526 /* L100: */
00527                     }
00528                     *scale *= scaloc;
00529                 }
00530                 i__1 = k + l * c_dim1;
00531                 c__[i__1].r = x11.r, c__[i__1].i = x11.i;
00532 
00533 /* L110: */
00534             }
00535 /* L120: */
00536         }
00537 
00538     }
00539 
00540     return 0;
00541 
00542 /*     End of CTRSYL */
00543 
00544 } /* ctrsyl_ */