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

swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:29