cla_syrcond_x.c
Go to the documentation of this file.
00001 /* cla_syrcond_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_syrcond_x__(char *uplo, integer *n, complex *a, integer *lda, 
00021         complex *af, integer *ldaf, integer *ipiv, complex *x, integer *info, 
00022         complex *work, real *rwork, ftnlen uplo_len)
00023 {
00024     /* System generated locals */
00025     integer a_dim1, a_offset, af_dim1, af_offset, i__1, i__2, i__3, i__4;
00026     real ret_val, r__1, r__2;
00027     complex q__1, q__2;
00028 
00029     /* Builtin functions */
00030     double r_imag(complex *);
00031     void c_div(complex *, complex *, complex *);
00032 
00033     /* Local variables */
00034     integer i__, j;
00035     logical up;
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     real ainvnm;
00044     extern /* Subroutine */ int csytrs_(char *, integer *, integer *, complex 
00045             *, integer *, integer *, complex *, integer *, integer *);
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_SYRCOND_X Computes the infinity norm condition number of */
00066 /*     op(A) * diag(X) where X is a COMPLEX vector. */
00067 
00068 /*  Arguments */
00069 /*  ========= */
00070 
00071 /*     UPLO    (input) CHARACTER*1 */
00072 /*       = 'U':  Upper triangle of A is stored; */
00073 /*       = 'L':  Lower triangle of A is stored. */
00074 
00075 /*     N       (input) INTEGER */
00076 /*     The number of linear equations, i.e., the order of the */
00077 /*     matrix A.  N >= 0. */
00078 
00079 /*     A       (input) COMPLEX array, dimension (LDA,N) */
00080 /*     On entry, the N-by-N matrix A. */
00081 
00082 /*     LDA     (input) INTEGER */
00083 /*     The leading dimension of the array A.  LDA >= max(1,N). */
00084 
00085 /*     AF      (input) COMPLEX array, dimension (LDAF,N) */
00086 /*     The block diagonal matrix D and the multipliers used to */
00087 /*     obtain the factor U or L as computed by CSYTRF. */
00088 
00089 /*     LDAF    (input) INTEGER */
00090 /*     The leading dimension of the array AF.  LDAF >= max(1,N). */
00091 
00092 /*     IPIV    (input) INTEGER array, dimension (N) */
00093 /*     Details of the interchanges and the block structure of D */
00094 /*     as determined by CSYTRF. */
00095 
00096 /*     X       (input) COMPLEX array, dimension (N) */
00097 /*     The vector X in the formula op(A) * diag(X). */
00098 
00099 /*     INFO    (output) INTEGER */
00100 /*       = 0:  Successful exit. */
00101 /*     i > 0:  The ith argument is invalid. */
00102 
00103 /*     WORK    (input) COMPLEX array, dimension (2*N). */
00104 /*     Workspace. */
00105 
00106 /*     RWORK   (input) REAL array, dimension (N). */
00107 /*     Workspace. */
00108 
00109 /*  ===================================================================== */
00110 
00111 /*     .. Local Scalars .. */
00112 /*     .. */
00113 /*     .. Local Arrays .. */
00114 /*     .. */
00115 /*     .. External Functions .. */
00116 /*     .. */
00117 /*     .. External Subroutines .. */
00118 /*     .. */
00119 /*     .. Intrinsic Functions .. */
00120 /*     .. */
00121 /*     .. Statement Functions .. */
00122 /*     .. */
00123 /*     .. Statement Function Definitions .. */
00124 /*     .. */
00125 /*     .. Executable Statements .. */
00126 
00127     /* Parameter adjustments */
00128     a_dim1 = *lda;
00129     a_offset = 1 + a_dim1;
00130     a -= a_offset;
00131     af_dim1 = *ldaf;
00132     af_offset = 1 + af_dim1;
00133     af -= af_offset;
00134     --ipiv;
00135     --x;
00136     --work;
00137     --rwork;
00138 
00139     /* Function Body */
00140     ret_val = 0.f;
00141 
00142     *info = 0;
00143     if (*n < 0) {
00144         *info = -2;
00145     }
00146     if (*info != 0) {
00147         i__1 = -(*info);
00148         xerbla_("CLA_SYRCOND_X", &i__1);
00149         return ret_val;
00150     }
00151     up = FALSE_;
00152     if (lsame_(uplo, "U")) {
00153         up = TRUE_;
00154     }
00155 
00156 /*     Compute norm of op(A)*op2(C). */
00157 
00158     anorm = 0.f;
00159     if (up) {
00160         i__1 = *n;
00161         for (i__ = 1; i__ <= i__1; ++i__) {
00162             tmp = 0.f;
00163             i__2 = i__;
00164             for (j = 1; j <= i__2; ++j) {
00165                 i__3 = j + i__ * a_dim1;
00166                 i__4 = j;
00167                 q__2.r = a[i__3].r * x[i__4].r - a[i__3].i * x[i__4].i, 
00168                         q__2.i = a[i__3].r * x[i__4].i + a[i__3].i * x[i__4]
00169                         .r;
00170                 q__1.r = q__2.r, q__1.i = q__2.i;
00171                 tmp += (r__1 = q__1.r, dabs(r__1)) + (r__2 = r_imag(&q__1), 
00172                         dabs(r__2));
00173             }
00174             i__2 = *n;
00175             for (j = i__ + 1; j <= i__2; ++j) {
00176                 i__3 = i__ + j * a_dim1;
00177                 i__4 = j;
00178                 q__2.r = a[i__3].r * x[i__4].r - a[i__3].i * x[i__4].i, 
00179                         q__2.i = a[i__3].r * x[i__4].i + a[i__3].i * x[i__4]
00180                         .r;
00181                 q__1.r = q__2.r, q__1.i = q__2.i;
00182                 tmp += (r__1 = q__1.r, dabs(r__1)) + (r__2 = r_imag(&q__1), 
00183                         dabs(r__2));
00184             }
00185             rwork[i__] = tmp;
00186             anorm = dmax(anorm,tmp);
00187         }
00188     } else {
00189         i__1 = *n;
00190         for (i__ = 1; i__ <= i__1; ++i__) {
00191             tmp = 0.f;
00192             i__2 = i__;
00193             for (j = 1; j <= i__2; ++j) {
00194                 i__3 = i__ + j * a_dim1;
00195                 i__4 = j;
00196                 q__2.r = a[i__3].r * x[i__4].r - a[i__3].i * x[i__4].i, 
00197                         q__2.i = a[i__3].r * x[i__4].i + a[i__3].i * x[i__4]
00198                         .r;
00199                 q__1.r = q__2.r, q__1.i = q__2.i;
00200                 tmp += (r__1 = q__1.r, dabs(r__1)) + (r__2 = r_imag(&q__1), 
00201                         dabs(r__2));
00202             }
00203             i__2 = *n;
00204             for (j = i__ + 1; j <= i__2; ++j) {
00205                 i__3 = j + i__ * a_dim1;
00206                 i__4 = j;
00207                 q__2.r = a[i__3].r * x[i__4].r - a[i__3].i * x[i__4].i, 
00208                         q__2.i = a[i__3].r * x[i__4].i + a[i__3].i * x[i__4]
00209                         .r;
00210                 q__1.r = q__2.r, q__1.i = q__2.i;
00211                 tmp += (r__1 = q__1.r, dabs(r__1)) + (r__2 = r_imag(&q__1), 
00212                         dabs(r__2));
00213             }
00214             rwork[i__] = tmp;
00215             anorm = dmax(anorm,tmp);
00216         }
00217     }
00218 
00219 /*     Quick return if possible. */
00220 
00221     if (*n == 0) {
00222         ret_val = 1.f;
00223         return ret_val;
00224     } else if (anorm == 0.f) {
00225         return ret_val;
00226     }
00227 
00228 /*     Estimate the norm of inv(op(A)). */
00229 
00230     ainvnm = 0.f;
00231 
00232     kase = 0;
00233 L10:
00234     clacn2_(n, &work[*n + 1], &work[1], &ainvnm, &kase, isave);
00235     if (kase != 0) {
00236         if (kase == 2) {
00237 
00238 /*           Multiply by R. */
00239 
00240             i__1 = *n;
00241             for (i__ = 1; i__ <= i__1; ++i__) {
00242                 i__2 = i__;
00243                 i__3 = i__;
00244                 i__4 = i__;
00245                 q__1.r = rwork[i__4] * work[i__3].r, q__1.i = rwork[i__4] * 
00246                         work[i__3].i;
00247                 work[i__2].r = q__1.r, work[i__2].i = q__1.i;
00248             }
00249 
00250             if (up) {
00251                 csytrs_("U", n, &c__1, &af[af_offset], ldaf, &ipiv[1], &work[
00252                         1], n, info);
00253             } else {
00254                 csytrs_("L", n, &c__1, &af[af_offset], ldaf, &ipiv[1], &work[
00255                         1], n, info);
00256             }
00257 
00258 /*           Multiply by inv(X). */
00259 
00260             i__1 = *n;
00261             for (i__ = 1; i__ <= i__1; ++i__) {
00262                 i__2 = i__;
00263                 c_div(&q__1, &work[i__], &x[i__]);
00264                 work[i__2].r = q__1.r, work[i__2].i = q__1.i;
00265             }
00266         } else {
00267 
00268 /*           Multiply by inv(X'). */
00269 
00270             i__1 = *n;
00271             for (i__ = 1; i__ <= i__1; ++i__) {
00272                 i__2 = i__;
00273                 c_div(&q__1, &work[i__], &x[i__]);
00274                 work[i__2].r = q__1.r, work[i__2].i = q__1.i;
00275             }
00276 
00277             if (up) {
00278                 csytrs_("U", n, &c__1, &af[af_offset], ldaf, &ipiv[1], &work[
00279                         1], n, info);
00280             } else {
00281                 csytrs_("L", n, &c__1, &af[af_offset], ldaf, &ipiv[1], &work[
00282                         1], n, info);
00283             }
00284 
00285 /*           Multiply by R. */
00286 
00287             i__1 = *n;
00288             for (i__ = 1; i__ <= i__1; ++i__) {
00289                 i__2 = i__;
00290                 i__3 = i__;
00291                 i__4 = i__;
00292                 q__1.r = rwork[i__4] * work[i__3].r, q__1.i = rwork[i__4] * 
00293                         work[i__3].i;
00294                 work[i__2].r = q__1.r, work[i__2].i = q__1.i;
00295             }
00296         }
00297         goto L10;
00298     }
00299 
00300 /*     Compute the estimate of the reciprocal condition number. */
00301 
00302     if (ainvnm != 0.f) {
00303         ret_val = 1.f / ainvnm;
00304     }
00305 
00306     return ret_val;
00307 
00308 } /* cla_syrcond_x__ */

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