swiftnav: cla_gbrcond_x.c Source File

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