zla_porcond_c.c
Go to the documentation of this file.
00001 /* zla_porcond_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_porcond_c__(char *uplo, integer *n, doublecomplex *a, integer *
00021         lda, doublecomplex *af, integer *ldaf, doublereal *c__, logical *
00022         capply, integer *info, doublecomplex *work, doublereal *rwork, ftnlen 
00023         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 zpotrs_(char *, integer *, integer *, 
00046             doublecomplex *, integer *, doublecomplex *, integer *, integer *);
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 /*     ZLA_PORCOND_C Computes the infinity norm condition number of */
00067 /*     op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector */
00068 
00069 /*  Arguments */
00070 /*  ========= */
00071 
00072 /*     UPLO    (input) CHARACTER*1 */
00073 /*       = 'U':  Upper triangle of A is stored; */
00074 /*       = 'L':  Lower triangle of A is stored. */
00075 
00076 /*     N       (input) INTEGER */
00077 /*     The number of linear equations, i.e., the order of the */
00078 /*     matrix A.  N >= 0. */
00079 
00080 /*     A       (input) COMPLEX*16 array, dimension (LDA,N) */
00081 /*     On entry, the N-by-N matrix A */
00082 
00083 /*     LDA     (input) INTEGER */
00084 /*     The leading dimension of the array A.  LDA >= max(1,N). */
00085 
00086 /*     AF      (input) COMPLEX*16 array, dimension (LDAF,N) */
00087 /*     The triangular factor U or L from the Cholesky factorization */
00088 /*     A = U**T*U or A = L*L**T, as computed by ZPOTRF. */
00089 
00090 /*     LDAF    (input) INTEGER */
00091 /*     The leading dimension of the array AF.  LDAF >= max(1,N). */
00092 
00093 /*     C       (input) DOUBLE PRECISION array, dimension (N) */
00094 /*     The vector C in the formula op(A) * inv(diag(C)). */
00095 
00096 /*     CAPPLY  (input) LOGICAL */
00097 /*     If .TRUE. then access the vector C in the formula above. */
00098 
00099 /*     INFO    (output) INTEGER */
00100 /*       = 0:  Successful exit. */
00101 /*     i > 0:  The ith argument is invalid. */
00102 
00103 /*     WORK    (input) COMPLEX*16 array, dimension (2*N). */
00104 /*     Workspace. */
00105 
00106 /*     RWORK   (input) DOUBLE PRECISION 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     --c__;
00135     --work;
00136     --rwork;
00137 
00138     /* Function Body */
00139     ret_val = 0.;
00140 
00141     *info = 0;
00142     if (*n < 0) {
00143         *info = -2;
00144     }
00145     if (*info != 0) {
00146         i__1 = -(*info);
00147         xerbla_("ZLA_PORCOND_C", &i__1);
00148         return ret_val;
00149     }
00150     up = FALSE_;
00151     if (lsame_(uplo, "U")) {
00152         up = TRUE_;
00153     }
00154 
00155 /*     Compute norm of op(A)*op2(C). */
00156 
00157     anorm = 0.;
00158     if (up) {
00159         i__1 = *n;
00160         for (i__ = 1; i__ <= i__1; ++i__) {
00161             tmp = 0.;
00162             if (*capply) {
00163                 i__2 = i__;
00164                 for (j = 1; j <= i__2; ++j) {
00165                     i__3 = j + i__ * a_dim1;
00166                     tmp += ((d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[
00167                             j + i__ * a_dim1]), abs(d__2))) / c__[j];
00168                 }
00169                 i__2 = *n;
00170                 for (j = i__ + 1; j <= i__2; ++j) {
00171                     i__3 = i__ + j * a_dim1;
00172                     tmp += ((d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[
00173                             i__ + j * a_dim1]), abs(d__2))) / c__[j];
00174                 }
00175             } else {
00176                 i__2 = i__;
00177                 for (j = 1; j <= i__2; ++j) {
00178                     i__3 = j + i__ * a_dim1;
00179                     tmp += (d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[
00180                             j + i__ * a_dim1]), abs(d__2));
00181                 }
00182                 i__2 = *n;
00183                 for (j = i__ + 1; j <= i__2; ++j) {
00184                     i__3 = i__ + j * a_dim1;
00185                     tmp += (d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[
00186                             i__ + j * a_dim1]), abs(d__2));
00187                 }
00188             }
00189             rwork[i__] = tmp;
00190             anorm = max(anorm,tmp);
00191         }
00192     } else {
00193         i__1 = *n;
00194         for (i__ = 1; i__ <= i__1; ++i__) {
00195             tmp = 0.;
00196             if (*capply) {
00197                 i__2 = i__;
00198                 for (j = 1; j <= i__2; ++j) {
00199                     i__3 = i__ + j * a_dim1;
00200                     tmp += ((d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[
00201                             i__ + j * a_dim1]), abs(d__2))) / c__[j];
00202                 }
00203                 i__2 = *n;
00204                 for (j = i__ + 1; j <= i__2; ++j) {
00205                     i__3 = j + i__ * a_dim1;
00206                     tmp += ((d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[
00207                             j + i__ * a_dim1]), abs(d__2))) / c__[j];
00208                 }
00209             } else {
00210                 i__2 = i__;
00211                 for (j = 1; j <= i__2; ++j) {
00212                     i__3 = i__ + j * a_dim1;
00213                     tmp += (d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[
00214                             i__ + j * a_dim1]), abs(d__2));
00215                 }
00216                 i__2 = *n;
00217                 for (j = i__ + 1; j <= i__2; ++j) {
00218                     i__3 = j + i__ * a_dim1;
00219                     tmp += (d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[
00220                             j + i__ * a_dim1]), abs(d__2));
00221                 }
00222             }
00223             rwork[i__] = tmp;
00224             anorm = max(anorm,tmp);
00225         }
00226     }
00227 
00228 /*     Quick return if possible. */
00229 
00230     if (*n == 0) {
00231         ret_val = 1.;
00232         return ret_val;
00233     } else if (anorm == 0.) {
00234         return ret_val;
00235     }
00236 
00237 /*     Estimate the norm of inv(op(A)). */
00238 
00239     ainvnm = 0.;
00240 
00241     kase = 0;
00242 L10:
00243     zlacn2_(n, &work[*n + 1], &work[1], &ainvnm, &kase, isave);
00244     if (kase != 0) {
00245         if (kase == 2) {
00246 
00247 /*           Multiply by R. */
00248 
00249             i__1 = *n;
00250             for (i__ = 1; i__ <= i__1; ++i__) {
00251                 i__2 = i__;
00252                 i__3 = i__;
00253                 i__4 = i__;
00254                 z__1.r = rwork[i__4] * work[i__3].r, z__1.i = rwork[i__4] * 
00255                         work[i__3].i;
00256                 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00257             }
00258 
00259             if (up) {
00260                 zpotrs_("U", n, &c__1, &af[af_offset], ldaf, &work[1], n, 
00261                         info);
00262             } else {
00263                 zpotrs_("L", n, &c__1, &af[af_offset], ldaf, &work[1], n, 
00264                         info);
00265             }
00266 
00267 /*           Multiply by inv(C). */
00268 
00269             if (*capply) {
00270                 i__1 = *n;
00271                 for (i__ = 1; i__ <= i__1; ++i__) {
00272                     i__2 = i__;
00273                     i__3 = i__;
00274                     i__4 = i__;
00275                     z__1.r = c__[i__4] * work[i__3].r, z__1.i = c__[i__4] * 
00276                             work[i__3].i;
00277                     work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00278                 }
00279             }
00280         } else {
00281 
00282 /*           Multiply by inv(C'). */
00283 
00284             if (*capply) {
00285                 i__1 = *n;
00286                 for (i__ = 1; i__ <= i__1; ++i__) {
00287                     i__2 = i__;
00288                     i__3 = i__;
00289                     i__4 = i__;
00290                     z__1.r = c__[i__4] * work[i__3].r, z__1.i = c__[i__4] * 
00291                             work[i__3].i;
00292                     work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00293                 }
00294             }
00295 
00296             if (up) {
00297                 zpotrs_("U", n, &c__1, &af[af_offset], ldaf, &work[1], n, 
00298                         info);
00299             } else {
00300                 zpotrs_("L", n, &c__1, &af[af_offset], ldaf, &work[1], n, 
00301                         info);
00302             }
00303 
00304 /*           Multiply by R. */
00305 
00306             i__1 = *n;
00307             for (i__ = 1; i__ <= i__1; ++i__) {
00308                 i__2 = i__;
00309                 i__3 = i__;
00310                 i__4 = i__;
00311                 z__1.r = rwork[i__4] * work[i__3].r, z__1.i = rwork[i__4] * 
00312                         work[i__3].i;
00313                 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00314             }
00315         }
00316         goto L10;
00317     }
00318 
00319 /*     Compute the estimate of the reciprocal condition number. */
00320 
00321     if (ainvnm != 0.) {
00322         ret_val = 1. / ainvnm;
00323     }
00324 
00325     return ret_val;
00326 
00327 } /* zla_porcond_c__ */

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