swiftnav: zla_gbrcond_x.c Source File

Go to the documentation of this file.
00001 /* zla_gbrcond_x.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 doublereal zla_gbrcond_x__(char *trans, integer *n, integer *kl, integer *ku, 
00021         doublecomplex *ab, integer *ldab, doublecomplex *afb, integer *ldafb, 
00022         integer *ipiv, doublecomplex *x, integer *info, doublecomplex *work, 
00023         doublereal *rwork, ftnlen trans_len)
00024 {
00025     /* System generated locals */
00026     integer ab_dim1, ab_offset, afb_dim1, afb_offset, i__1, i__2, i__3, i__4;
00027     doublereal ret_val, d__1, d__2;
00028     doublecomplex z__1, z__2;
00029 
00030     /* Builtin functions */
00031     double d_imag(doublecomplex *);
00032     void z_div(doublecomplex *, doublecomplex *, doublecomplex *);
00033 
00034     /* Local variables */
00035     integer i__, j, kd, ke;
00036     doublereal tmp;
00037     integer kase;
00038     extern logical lsame_(char *, char *);
00039     integer isave[3];
00040     doublereal anorm;
00041     extern /* Subroutine */ int zlacn2_(integer *, doublecomplex *, 
00042             doublecomplex *, doublereal *, integer *, integer *), xerbla_(
00043             char *, integer *);
00044     doublereal ainvnm;
00045     extern /* Subroutine */ int zgbtrs_(char *, integer *, integer *, integer 
00046             *, integer *, doublecomplex *, integer *, integer *, 
00047             doublecomplex *, integer *, integer *);
00048     logical notrans;
00049 
00050 
00051 /*     -- LAPACK routine (version 3.2.1)                               -- */
00052 /*     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- */
00053 /*     -- Jason Riedy of Univ. of California Berkeley.                 -- */
00054 /*     -- April 2009                                                   -- */
00055 
00056 /*     -- LAPACK is a software package provided by Univ. of Tennessee, -- */
00057 /*     -- Univ. of California Berkeley and NAG Ltd.                    -- */
00058 
00059 /*     .. */
00060 /*     .. Scalar Arguments .. */
00061 /*     .. */
00062 /*     .. Array Arguments .. */
00063 
00064 
00065 /*  Purpose */
00066 /*  ======= */
00067 
00068 /*     ZLA_GBRCOND_X Computes the infinity norm condition number of */
00069 /*     op(A) * diag(X) where X is a COMPLEX*16 vector. */
00070 
00071 /*  Arguments */
00072 /*  ========= */
00073 
00074 /*     TRANS   (input) CHARACTER*1 */
00075 /*     Specifies the form of the system of equations: */
00076 /*       = 'N':  A * X = B     (No transpose) */
00077 /*       = 'T':  A**T * X = B  (Transpose) */
00078 /*       = 'C':  A**H * X = B  (Conjugate Transpose = Transpose) */
00079 
00080 /*     N       (input) INTEGER */
00081 /*     The number of linear equations, i.e., the order of the */
00082 /*     matrix A.  N >= 0. */
00083 
00084 /*     KL      (input) INTEGER */
00085 /*     The number of subdiagonals within the band of A.  KL >= 0. */
00086 
00087 /*     KU      (input) INTEGER */
00088 /*     The number of superdiagonals within the band of A.  KU >= 0. */
00089 
00090 /*     AB      (input) COMPLEX*16 array, dimension (LDAB,N) */
00091 /*     On entry, the matrix A in band storage, in rows 1 to KL+KU+1. */
00092 /*     The j-th column of A is stored in the j-th column of the */
00093 /*     array AB as follows: */
00094 /*     AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl) */
00095 
00096 /*     LDAB    (input) INTEGER */
00097 /*     The leading dimension of the array AB.  LDAB >= KL+KU+1. */
00098 
00099 /*     AFB     (input) COMPLEX*16 array, dimension (LDAFB,N) */
00100 /*     Details of the LU factorization of the band matrix A, as */
00101 /*     computed by ZGBTRF.  U is stored as an upper triangular */
00102 /*     band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, */
00103 /*     and the multipliers used during the factorization are stored */
00104 /*     in rows KL+KU+2 to 2*KL+KU+1. */
00105 
00106 /*     LDAFB   (input) INTEGER */
00107 /*     The leading dimension of the array AFB.  LDAFB >= 2*KL+KU+1. */
00108 
00109 /*     IPIV    (input) INTEGER array, dimension (N) */
00110 /*     The pivot indices from the factorization A = P*L*U */
00111 /*     as computed by ZGBTRF; row i of the matrix was interchanged */
00112 /*     with row IPIV(i). */
00113 
00114 /*     X       (input) COMPLEX*16 array, dimension (N) */
00115 /*     The vector X in the formula op(A) * diag(X). */
00116 
00117 /*     INFO    (output) INTEGER */
00118 /*       = 0:  Successful exit. */
00119 /*     i > 0:  The ith argument is invalid. */
00120 
00121 /*     WORK    (input) COMPLEX*16 array, dimension (2*N). */
00122 /*     Workspace. */
00123 
00124 /*     RWORK   (input) DOUBLE PRECISION array, dimension (N). */
00125 /*     Workspace. */
00126 
00127 /*  ===================================================================== */
00128 
00129 /*     .. Local Scalars .. */
00130 /*     .. */
00131 /*     .. Local Arrays .. */
00132 /*     .. */
00133 /*     .. External Functions .. */
00134 /*     .. */
00135 /*     .. External Subroutines .. */
00136 /*     .. */
00137 /*     .. Intrinsic Functions .. */
00138 /*     .. */
00139 /*     .. Statement Functions .. */
00140 /*     .. */
00141 /*     .. Statement Function Definitions .. */
00142 /*     .. */
00143 /*     .. Executable Statements .. */
00144 
00145     /* Parameter adjustments */
00146     ab_dim1 = *ldab;
00147     ab_offset = 1 + ab_dim1;
00148     ab -= ab_offset;
00149     afb_dim1 = *ldafb;
00150     afb_offset = 1 + afb_dim1;
00151     afb -= afb_offset;
00152     --ipiv;
00153     --x;
00154     --work;
00155     --rwork;
00156 
00157     /* Function Body */
00158     ret_val = 0.;
00159 
00160     *info = 0;
00161     notrans = lsame_(trans, "N");
00162     if (! notrans && ! lsame_(trans, "T") && ! lsame_(
00163             trans, "C")) {
00164         *info = -1;
00165     } else if (*n < 0) {
00166         *info = -2;
00167     } else if (*kl < 0 || *kl > *n - 1) {
00168         *info = -3;
00169     } else if (*ku < 0 || *ku > *n - 1) {
00170         *info = -4;
00171     } else if (*ldab < *kl + *ku + 1) {
00172         *info = -6;
00173     } else if (*ldafb < (*kl << 1) + *ku + 1) {
00174         *info = -8;
00175     }
00176     if (*info != 0) {
00177         i__1 = -(*info);
00178         xerbla_("ZLA_GBRCOND_X", &i__1);
00179         return ret_val;
00180     }
00181 
00182 /*     Compute norm of op(A)*op2(C). */
00183 
00184     kd = *ku + 1;
00185     ke = *kl + 1;
00186     anorm = 0.;
00187     if (notrans) {
00188         i__1 = *n;
00189         for (i__ = 1; i__ <= i__1; ++i__) {
00190             tmp = 0.;
00191 /* Computing MAX */
00192             i__2 = i__ - *kl;
00193 /* Computing MIN */
00194             i__4 = i__ + *ku;
00195             i__3 = min(i__4,*n);
00196             for (j = max(i__2,1); j <= i__3; ++j) {
00197                 i__2 = kd + i__ - j + j * ab_dim1;
00198                 i__4 = j;
00199                 z__2.r = ab[i__2].r * x[i__4].r - ab[i__2].i * x[i__4].i, 
00200                         z__2.i = ab[i__2].r * x[i__4].i + ab[i__2].i * x[i__4]
00201                         .r;
00202                 z__1.r = z__2.r, z__1.i = z__2.i;
00203                 tmp += (d__1 = z__1.r, abs(d__1)) + (d__2 = d_imag(&z__1), 
00204                         abs(d__2));
00205             }
00206             rwork[i__] = tmp;
00207             anorm = max(anorm,tmp);
00208         }
00209     } else {
00210         i__1 = *n;
00211         for (i__ = 1; i__ <= i__1; ++i__) {
00212             tmp = 0.;
00213 /* Computing MAX */
00214             i__3 = i__ - *kl;
00215 /* Computing MIN */
00216             i__4 = i__ + *ku;
00217             i__2 = min(i__4,*n);
00218             for (j = max(i__3,1); j <= i__2; ++j) {
00219                 i__3 = ke - i__ + j + i__ * ab_dim1;
00220                 i__4 = j;
00221                 z__2.r = ab[i__3].r * x[i__4].r - ab[i__3].i * x[i__4].i, 
00222                         z__2.i = ab[i__3].r * x[i__4].i + ab[i__3].i * x[i__4]
00223                         .r;
00224                 z__1.r = z__2.r, z__1.i = z__2.i;
00225                 tmp += (d__1 = z__1.r, abs(d__1)) + (d__2 = d_imag(&z__1), 
00226                         abs(d__2));
00227             }
00228             rwork[i__] = tmp;
00229             anorm = max(anorm,tmp);
00230         }
00231     }
00232 
00233 /*     Quick return if possible. */
00234 
00235     if (*n == 0) {
00236         ret_val = 1.;
00237         return ret_val;
00238     } else if (anorm == 0.) {
00239         return ret_val;
00240     }
00241 
00242 /*     Estimate the norm of inv(op(A)). */
00243 
00244     ainvnm = 0.;
00245 
00246     kase = 0;
00247 L10:
00248     zlacn2_(n, &work[*n + 1], &work[1], &ainvnm, &kase, isave);
00249     if (kase != 0) {
00250         if (kase == 2) {
00251 
00252 /*           Multiply by R. */
00253 
00254             i__1 = *n;
00255             for (i__ = 1; i__ <= i__1; ++i__) {
00256                 i__2 = i__;
00257                 i__3 = i__;
00258                 i__4 = i__;
00259                 z__1.r = rwork[i__4] * work[i__3].r, z__1.i = rwork[i__4] * 
00260                         work[i__3].i;
00261                 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00262             }
00263 
00264             if (notrans) {
00265                 zgbtrs_("No transpose", n, kl, ku, &c__1, &afb[afb_offset], 
00266                         ldafb, &ipiv[1], &work[1], n, info);
00267             } else {
00268                 zgbtrs_("Conjugate transpose", n, kl, ku, &c__1, &afb[
00269                         afb_offset], ldafb, &ipiv[1], &work[1], n, info);
00270             }
00271 
00272 /*           Multiply by inv(X). */
00273 
00274             i__1 = *n;
00275             for (i__ = 1; i__ <= i__1; ++i__) {
00276                 i__2 = i__;
00277                 z_div(&z__1, &work[i__], &x[i__]);
00278                 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00279             }
00280         } else {
00281 
00282 /*           Multiply by inv(X'). */
00283 
00284             i__1 = *n;
00285             for (i__ = 1; i__ <= i__1; ++i__) {
00286                 i__2 = i__;
00287                 z_div(&z__1, &work[i__], &x[i__]);
00288                 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00289             }
00290 
00291             if (notrans) {
00292                 zgbtrs_("Conjugate transpose", n, kl, ku, &c__1, &afb[
00293                         afb_offset], ldafb, &ipiv[1], &work[1], n, info);
00294             } else {
00295                 zgbtrs_("No transpose", n, kl, ku, &c__1, &afb[afb_offset], 
00296                         ldafb, &ipiv[1], &work[1], n, info);
00297             }
00298 
00299 /*           Multiply by R. */
00300 
00301             i__1 = *n;
00302             for (i__ = 1; i__ <= i__1; ++i__) {
00303                 i__2 = i__;
00304                 i__3 = i__;
00305                 i__4 = i__;
00306                 z__1.r = rwork[i__4] * work[i__3].r, z__1.i = rwork[i__4] * 
00307                         work[i__3].i;
00308                 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00309             }
00310         }
00311         goto L10;
00312     }
00313 
00314 /*     Compute the estimate of the reciprocal condition number. */
00315 
00316     if (ainvnm != 0.) {
00317         ret_val = 1. / ainvnm;
00318     }
00319 
00320     return ret_val;
00321 
00322 } /* zla_gbrcond_x__ */