zla_syrcond_c.c
Go to the documentation of this file.
00001 /* zla_syrcond_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 zla_syrcond_c__(char *uplo, integer *n, doublecomplex *a, integer *
00021         lda, doublecomplex *af, integer *ldaf, integer *ipiv, doublereal *c__,
00022          logical *capply, integer *info, doublecomplex *work, doublereal *
00023         rwork, ftnlen uplo_len)
00024 {
00025     /* System generated locals */
00026     integer a_dim1, a_offset, af_dim1, af_offset, i__1, i__2, i__3, i__4;
00027     doublereal ret_val, d__1, d__2;
00028     doublecomplex z__1;
00029 
00030     /* Builtin functions */
00031     double d_imag(doublecomplex *);
00032 
00033     /* Local variables */
00034     integer i__, j;
00035     logical up;
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 zsytrs_(char *, integer *, integer *, 
00046             doublecomplex *, integer *, integer *, doublecomplex *, integer *, 
00047              integer *);
00048 
00049 
00050 /*     -- LAPACK routine (version 3.2.1)                                 -- */
00051 /*     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- */
00052 /*     -- Jason Riedy of Univ. of California Berkeley.                 -- */
00053 /*     -- April 2009                                                   -- */
00054 
00055 /*     -- LAPACK is a software package provided by Univ. of Tennessee, -- */
00056 /*     -- Univ. of California Berkeley and NAG Ltd.                    -- */
00057 
00058 /*     .. */
00059 /*     .. Scalar Arguments .. */
00060 /*     .. */
00061 /*     .. Array Arguments .. */
00062 /*     .. */
00063 
00064 /*  Purpose */
00065 /*  ======= */
00066 
00067 /*     ZLA_SYRCOND_C Computes the infinity norm condition number of */
00068 /*     op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector. */
00069 
00070 /*  Arguments */
00071 /*  ========= */
00072 
00073 /*     UPLO    (input) CHARACTER*1 */
00074 /*       = 'U':  Upper triangle of A is stored; */
00075 /*       = 'L':  Lower triangle of A is stored. */
00076 
00077 /*     N       (input) INTEGER */
00078 /*     The number of linear equations, i.e., the order of the */
00079 /*     matrix A.  N >= 0. */
00080 
00081 /*     A       (input) COMPLEX*16 array, dimension (LDA,N) */
00082 /*     On entry, the N-by-N matrix A */
00083 
00084 /*     LDA     (input) INTEGER */
00085 /*     The leading dimension of the array A.  LDA >= max(1,N). */
00086 
00087 /*     AF      (input) COMPLEX*16 array, dimension (LDAF,N) */
00088 /*     The block diagonal matrix D and the multipliers used to */
00089 /*     obtain the factor U or L as computed by ZSYTRF. */
00090 
00091 /*     LDAF    (input) INTEGER */
00092 /*     The leading dimension of the array AF.  LDAF >= max(1,N). */
00093 
00094 /*     IPIV    (input) INTEGER array, dimension (N) */
00095 /*     Details of the interchanges and the block structure of D */
00096 /*     as determined by ZSYTRF. */
00097 
00098 /*     C       (input) DOUBLE PRECISION array, dimension (N) */
00099 /*     The vector C in the formula op(A) * inv(diag(C)). */
00100 
00101 /*     CAPPLY  (input) LOGICAL */
00102 /*     If .TRUE. then access the vector C in the formula above. */
00103 
00104 /*     INFO    (output) INTEGER */
00105 /*       = 0:  Successful exit. */
00106 /*     i > 0:  The ith argument is invalid. */
00107 
00108 /*     WORK    (input) COMPLEX*16 array, dimension (2*N). */
00109 /*     Workspace. */
00110 
00111 /*     RWORK   (input) DOUBLE PRECISION array, dimension (N). */
00112 /*     Workspace. */
00113 
00114 /*  ===================================================================== */
00115 
00116 /*     .. Local Scalars .. */
00117 /*     .. */
00118 /*     .. Local Arrays .. */
00119 /*     .. */
00120 /*     .. External Functions .. */
00121 /*     .. */
00122 /*     .. External Subroutines .. */
00123 /*     .. */
00124 /*     .. Intrinsic Functions .. */
00125 /*     .. */
00126 /*     .. Statement Functions .. */
00127 /*     .. */
00128 /*     .. Statement Function Definitions .. */
00129 /*     .. */
00130 /*     .. Executable Statements .. */
00131 
00132     /* Parameter adjustments */
00133     a_dim1 = *lda;
00134     a_offset = 1 + a_dim1;
00135     a -= a_offset;
00136     af_dim1 = *ldaf;
00137     af_offset = 1 + af_dim1;
00138     af -= af_offset;
00139     --ipiv;
00140     --c__;
00141     --work;
00142     --rwork;
00143 
00144     /* Function Body */
00145     ret_val = 0.;
00146 
00147     *info = 0;
00148     if (*n < 0) {
00149         *info = -2;
00150     }
00151     if (*info != 0) {
00152         i__1 = -(*info);
00153         xerbla_("ZLA_SYRCOND_C", &i__1);
00154         return ret_val;
00155     }
00156     up = FALSE_;
00157     if (lsame_(uplo, "U")) {
00158         up = TRUE_;
00159     }
00160 
00161 /*     Compute norm of op(A)*op2(C). */
00162 
00163     anorm = 0.;
00164     if (up) {
00165         i__1 = *n;
00166         for (i__ = 1; i__ <= i__1; ++i__) {
00167             tmp = 0.;
00168             if (*capply) {
00169                 i__2 = i__;
00170                 for (j = 1; j <= i__2; ++j) {
00171                     i__3 = j + i__ * a_dim1;
00172                     tmp += ((d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[
00173                             j + i__ * a_dim1]), abs(d__2))) / c__[j];
00174                 }
00175                 i__2 = *n;
00176                 for (j = i__ + 1; j <= i__2; ++j) {
00177                     i__3 = i__ + j * a_dim1;
00178                     tmp += ((d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[
00179                             i__ + j * a_dim1]), abs(d__2))) / c__[j];
00180                 }
00181             } else {
00182                 i__2 = i__;
00183                 for (j = 1; j <= i__2; ++j) {
00184                     i__3 = j + i__ * a_dim1;
00185                     tmp += (d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[
00186                             j + i__ * a_dim1]), abs(d__2));
00187                 }
00188                 i__2 = *n;
00189                 for (j = i__ + 1; j <= i__2; ++j) {
00190                     i__3 = i__ + j * a_dim1;
00191                     tmp += (d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[
00192                             i__ + j * a_dim1]), abs(d__2));
00193                 }
00194             }
00195             rwork[i__] = tmp;
00196             anorm = max(anorm,tmp);
00197         }
00198     } else {
00199         i__1 = *n;
00200         for (i__ = 1; i__ <= i__1; ++i__) {
00201             tmp = 0.;
00202             if (*capply) {
00203                 i__2 = i__;
00204                 for (j = 1; j <= i__2; ++j) {
00205                     i__3 = i__ + j * a_dim1;
00206                     tmp += ((d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[
00207                             i__ + j * a_dim1]), abs(d__2))) / c__[j];
00208                 }
00209                 i__2 = *n;
00210                 for (j = i__ + 1; j <= i__2; ++j) {
00211                     i__3 = j + i__ * a_dim1;
00212                     tmp += ((d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[
00213                             j + i__ * a_dim1]), abs(d__2))) / c__[j];
00214                 }
00215             } else {
00216                 i__2 = i__;
00217                 for (j = 1; j <= i__2; ++j) {
00218                     i__3 = i__ + j * a_dim1;
00219                     tmp += (d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[
00220                             i__ + j * a_dim1]), abs(d__2));
00221                 }
00222                 i__2 = *n;
00223                 for (j = i__ + 1; j <= i__2; ++j) {
00224                     i__3 = j + i__ * a_dim1;
00225                     tmp += (d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[
00226                             j + i__ * a_dim1]), abs(d__2));
00227                 }
00228             }
00229             rwork[i__] = tmp;
00230             anorm = max(anorm,tmp);
00231         }
00232     }
00233 
00234 /*     Quick return if possible. */
00235 
00236     if (*n == 0) {
00237         ret_val = 1.;
00238         return ret_val;
00239     } else if (anorm == 0.) {
00240         return ret_val;
00241     }
00242 
00243 /*     Estimate the norm of inv(op(A)). */
00244 
00245     ainvnm = 0.;
00246 
00247     kase = 0;
00248 L10:
00249     zlacn2_(n, &work[*n + 1], &work[1], &ainvnm, &kase, isave);
00250     if (kase != 0) {
00251         if (kase == 2) {
00252 
00253 /*           Multiply by R. */
00254 
00255             i__1 = *n;
00256             for (i__ = 1; i__ <= i__1; ++i__) {
00257                 i__2 = i__;
00258                 i__3 = i__;
00259                 i__4 = i__;
00260                 z__1.r = rwork[i__4] * work[i__3].r, z__1.i = rwork[i__4] * 
00261                         work[i__3].i;
00262                 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00263             }
00264 
00265             if (up) {
00266                 zsytrs_("U", n, &c__1, &af[af_offset], ldaf, &ipiv[1], &work[
00267                         1], n, info);
00268             } else {
00269                 zsytrs_("L", n, &c__1, &af[af_offset], ldaf, &ipiv[1], &work[
00270                         1], n, info);
00271             }
00272 
00273 /*           Multiply by inv(C). */
00274 
00275             if (*capply) {
00276                 i__1 = *n;
00277                 for (i__ = 1; i__ <= i__1; ++i__) {
00278                     i__2 = i__;
00279                     i__3 = i__;
00280                     i__4 = i__;
00281                     z__1.r = c__[i__4] * work[i__3].r, z__1.i = c__[i__4] * 
00282                             work[i__3].i;
00283                     work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00284                 }
00285             }
00286         } else {
00287 
00288 /*           Multiply by inv(C'). */
00289 
00290             if (*capply) {
00291                 i__1 = *n;
00292                 for (i__ = 1; i__ <= i__1; ++i__) {
00293                     i__2 = i__;
00294                     i__3 = i__;
00295                     i__4 = i__;
00296                     z__1.r = c__[i__4] * work[i__3].r, z__1.i = c__[i__4] * 
00297                             work[i__3].i;
00298                     work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00299                 }
00300             }
00301 
00302             if (up) {
00303                 zsytrs_("U", n, &c__1, &af[af_offset], ldaf, &ipiv[1], &work[
00304                         1], n, info);
00305             } else {
00306                 zsytrs_("L", n, &c__1, &af[af_offset], ldaf, &ipiv[1], &work[
00307                         1], n, info);
00308             }
00309 
00310 /*           Multiply by R. */
00311 
00312             i__1 = *n;
00313             for (i__ = 1; i__ <= i__1; ++i__) {
00314                 i__2 = i__;
00315                 i__3 = i__;
00316                 i__4 = i__;
00317                 z__1.r = rwork[i__4] * work[i__3].r, z__1.i = rwork[i__4] * 
00318                         work[i__3].i;
00319                 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00320             }
00321         }
00322         goto L10;
00323     }
00324 
00325 /*     Compute the estimate of the reciprocal condition number. */
00326 
00327     if (ainvnm != 0.) {
00328         ret_val = 1. / ainvnm;
00329     }
00330 
00331     return ret_val;
00332 
00333 } /* zla_syrcond_c__ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:40