zla_herpvgrw.c
Go to the documentation of this file.
00001 /* zla_herpvgrw.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 doublereal zla_herpvgrw__(char *uplo, integer *n, integer *info, 
00017         doublecomplex *a, integer *lda, doublecomplex *af, integer *ldaf, 
00018         integer *ipiv, doublereal *work, ftnlen uplo_len)
00019 {
00020     /* System generated locals */
00021     integer a_dim1, a_offset, af_dim1, af_offset, i__1, i__2, i__3;
00022     doublereal ret_val, d__1, d__2, d__3, d__4;
00023 
00024     /* Builtin functions */
00025     double d_imag(doublecomplex *);
00026 
00027     /* Local variables */
00028     integer i__, j, k, kp;
00029     doublereal tmp, amax, umax;
00030     extern logical lsame_(char *, char *);
00031     integer ncols;
00032     logical upper;
00033     doublereal rpvgrw;
00034 
00035 
00036 /*     -- LAPACK routine (version 3.2.1)                                 -- */
00037 /*     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- */
00038 /*     -- Jason Riedy of Univ. of California Berkeley.                 -- */
00039 /*     -- April 2009                                                   -- */
00040 
00041 /*     -- LAPACK is a software package provided by Univ. of Tennessee, -- */
00042 /*     -- Univ. of California Berkeley and NAG Ltd.                    -- */
00043 
00044 /*     .. */
00045 /*     .. Scalar Arguments .. */
00046 /*     .. */
00047 /*     .. Array Arguments .. */
00048 /*     .. */
00049 
00050 /*  Purpose */
00051 /*  ======= */
00052 
00053 /*  ZLA_HERPVGRW computes the reciprocal pivot growth factor */
00054 /*  norm(A)/norm(U). The "max absolute element" norm is used. If this is */
00055 /*  much less than 1, the stability of the LU factorization of the */
00056 /*  (equilibrated) matrix A could be poor. This also means that the */
00057 /*  solution X, estimated condition numbers, and error bounds could be */
00058 /*  unreliable. */
00059 
00060 /*  Arguments */
00061 /*  ========= */
00062 
00063 /*     UPLO    (input) CHARACTER*1 */
00064 /*       = 'U':  Upper triangle of A is stored; */
00065 /*       = 'L':  Lower triangle of A is stored. */
00066 
00067 /*     N       (input) INTEGER */
00068 /*     The number of linear equations, i.e., the order of the */
00069 /*     matrix A.  N >= 0. */
00070 
00071 /*     INFO    (input) INTEGER */
00072 /*     The value of INFO returned from ZHETRF, .i.e., the pivot in */
00073 /*     column INFO is exactly 0. */
00074 
00075 /*     NCOLS   (input) INTEGER */
00076 /*     The number of columns of the matrix A. NCOLS >= 0. */
00077 
00078 /*     A       (input) COMPLEX*16 array, dimension (LDA,N) */
00079 /*     On entry, the N-by-N matrix A. */
00080 
00081 /*     LDA     (input) INTEGER */
00082 /*     The leading dimension of the array A.  LDA >= max(1,N). */
00083 
00084 /*     AF      (input) COMPLEX*16 array, dimension (LDAF,N) */
00085 /*     The block diagonal matrix D and the multipliers used to */
00086 /*     obtain the factor U or L as computed by ZHETRF. */
00087 
00088 /*     LDAF    (input) INTEGER */
00089 /*     The leading dimension of the array AF.  LDAF >= max(1,N). */
00090 
00091 /*     IPIV    (input) INTEGER array, dimension (N) */
00092 /*     Details of the interchanges and the block structure of D */
00093 /*     as determined by ZHETRF. */
00094 
00095 /*     WORK    (input) COMPLEX*16 array, dimension (2*N) */
00096 
00097 /*  ===================================================================== */
00098 
00099 /*     .. Local Scalars .. */
00100 /*     .. */
00101 /*     .. External Functions .. */
00102 /*     .. */
00103 /*     .. Intrinsic Functions .. */
00104 /*     .. */
00105 /*     .. Statement Functions .. */
00106 /*     .. */
00107 /*     .. Statement Function Definitions .. */
00108 /*     .. */
00109 /*     .. Executable Statements .. */
00110 
00111     /* Parameter adjustments */
00112     a_dim1 = *lda;
00113     a_offset = 1 + a_dim1;
00114     a -= a_offset;
00115     af_dim1 = *ldaf;
00116     af_offset = 1 + af_dim1;
00117     af -= af_offset;
00118     --ipiv;
00119     --work;
00120 
00121     /* Function Body */
00122     upper = lsame_("Upper", uplo);
00123     if (*info == 0) {
00124         if (upper) {
00125             ncols = 1;
00126         } else {
00127             ncols = *n;
00128         }
00129     } else {
00130         ncols = *info;
00131     }
00132     rpvgrw = 1.;
00133     i__1 = *n << 1;
00134     for (i__ = 1; i__ <= i__1; ++i__) {
00135         work[i__] = 0.;
00136     }
00137 
00138 /*     Find the max magnitude entry of each column of A.  Compute the max */
00139 /*     for all N columns so we can apply the pivot permutation while */
00140 /*     looping below.  Assume a full factorization is the common case. */
00141 
00142     if (upper) {
00143         i__1 = *n;
00144         for (j = 1; j <= i__1; ++j) {
00145             i__2 = j;
00146             for (i__ = 1; i__ <= i__2; ++i__) {
00147 /* Computing MAX */
00148                 i__3 = i__ + j * a_dim1;
00149                 d__3 = (d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[i__ 
00150                         + j * a_dim1]), abs(d__2)), d__4 = work[*n + i__];
00151                 work[*n + i__] = max(d__3,d__4);
00152 /* Computing MAX */
00153                 i__3 = i__ + j * a_dim1;
00154                 d__3 = (d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[i__ 
00155                         + j * a_dim1]), abs(d__2)), d__4 = work[*n + j];
00156                 work[*n + j] = max(d__3,d__4);
00157             }
00158         }
00159     } else {
00160         i__1 = *n;
00161         for (j = 1; j <= i__1; ++j) {
00162             i__2 = *n;
00163             for (i__ = j; i__ <= i__2; ++i__) {
00164 /* Computing MAX */
00165                 i__3 = i__ + j * a_dim1;
00166                 d__3 = (d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[i__ 
00167                         + j * a_dim1]), abs(d__2)), d__4 = work[*n + i__];
00168                 work[*n + i__] = max(d__3,d__4);
00169 /* Computing MAX */
00170                 i__3 = i__ + j * a_dim1;
00171                 d__3 = (d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a[i__ 
00172                         + j * a_dim1]), abs(d__2)), d__4 = work[*n + j];
00173                 work[*n + j] = max(d__3,d__4);
00174             }
00175         }
00176     }
00177 
00178 /*     Now find the max magnitude entry of each column of U or L.  Also */
00179 /*     permute the magnitudes of A above so they're in the same order as */
00180 /*     the factor. */
00181 
00182 /*     The iteration orders and permutations were copied from zsytrs. */
00183 /*     Calls to SSWAP would be severe overkill. */
00184 
00185     if (upper) {
00186         k = *n;
00187         while(k < ncols && k > 0) {
00188             if (ipiv[k] > 0) {
00189 /*              1x1 pivot */
00190                 kp = ipiv[k];
00191                 if (kp != k) {
00192                     tmp = work[*n + k];
00193                     work[*n + k] = work[*n + kp];
00194                     work[*n + kp] = tmp;
00195                 }
00196                 i__1 = k;
00197                 for (i__ = 1; i__ <= i__1; ++i__) {
00198 /* Computing MAX */
00199                     i__2 = i__ + k * af_dim1;
00200                     d__3 = (d__1 = af[i__2].r, abs(d__1)) + (d__2 = d_imag(&
00201                             af[i__ + k * af_dim1]), abs(d__2)), d__4 = work[k]
00202                             ;
00203                     work[k] = max(d__3,d__4);
00204                 }
00205                 --k;
00206             } else {
00207 /*              2x2 pivot */
00208                 kp = -ipiv[k];
00209                 tmp = work[*n + k - 1];
00210                 work[*n + k - 1] = work[*n + kp];
00211                 work[*n + kp] = tmp;
00212                 i__1 = k - 1;
00213                 for (i__ = 1; i__ <= i__1; ++i__) {
00214 /* Computing MAX */
00215                     i__2 = i__ + k * af_dim1;
00216                     d__3 = (d__1 = af[i__2].r, abs(d__1)) + (d__2 = d_imag(&
00217                             af[i__ + k * af_dim1]), abs(d__2)), d__4 = work[k]
00218                             ;
00219                     work[k] = max(d__3,d__4);
00220 /* Computing MAX */
00221                     i__2 = i__ + (k - 1) * af_dim1;
00222                     d__3 = (d__1 = af[i__2].r, abs(d__1)) + (d__2 = d_imag(&
00223                             af[i__ + (k - 1) * af_dim1]), abs(d__2)), d__4 = 
00224                             work[k - 1];
00225                     work[k - 1] = max(d__3,d__4);
00226                 }
00227 /* Computing MAX */
00228                 i__1 = k + k * af_dim1;
00229                 d__3 = (d__1 = af[i__1].r, abs(d__1)) + (d__2 = d_imag(&af[k 
00230                         + k * af_dim1]), abs(d__2)), d__4 = work[k];
00231                 work[k] = max(d__3,d__4);
00232                 k += -2;
00233             }
00234         }
00235         k = ncols;
00236         while(k <= *n) {
00237             if (ipiv[k] > 0) {
00238                 kp = ipiv[k];
00239                 if (kp != k) {
00240                     tmp = work[*n + k];
00241                     work[*n + k] = work[*n + kp];
00242                     work[*n + kp] = tmp;
00243                 }
00244                 ++k;
00245             } else {
00246                 kp = -ipiv[k];
00247                 tmp = work[*n + k];
00248                 work[*n + k] = work[*n + kp];
00249                 work[*n + kp] = tmp;
00250                 k += 2;
00251             }
00252         }
00253     } else {
00254         k = 1;
00255         while(k <= ncols) {
00256             if (ipiv[k] > 0) {
00257 /*              1x1 pivot */
00258                 kp = ipiv[k];
00259                 if (kp != k) {
00260                     tmp = work[*n + k];
00261                     work[*n + k] = work[*n + kp];
00262                     work[*n + kp] = tmp;
00263                 }
00264                 i__1 = *n;
00265                 for (i__ = k; i__ <= i__1; ++i__) {
00266 /* Computing MAX */
00267                     i__2 = i__ + k * af_dim1;
00268                     d__3 = (d__1 = af[i__2].r, abs(d__1)) + (d__2 = d_imag(&
00269                             af[i__ + k * af_dim1]), abs(d__2)), d__4 = work[k]
00270                             ;
00271                     work[k] = max(d__3,d__4);
00272                 }
00273                 ++k;
00274             } else {
00275 /*              2x2 pivot */
00276                 kp = -ipiv[k];
00277                 tmp = work[*n + k + 1];
00278                 work[*n + k + 1] = work[*n + kp];
00279                 work[*n + kp] = tmp;
00280                 i__1 = *n;
00281                 for (i__ = k + 1; i__ <= i__1; ++i__) {
00282 /* Computing MAX */
00283                     i__2 = i__ + k * af_dim1;
00284                     d__3 = (d__1 = af[i__2].r, abs(d__1)) + (d__2 = d_imag(&
00285                             af[i__ + k * af_dim1]), abs(d__2)), d__4 = work[k]
00286                             ;
00287                     work[k] = max(d__3,d__4);
00288 /* Computing MAX */
00289                     i__2 = i__ + (k + 1) * af_dim1;
00290                     d__3 = (d__1 = af[i__2].r, abs(d__1)) + (d__2 = d_imag(&
00291                             af[i__ + (k + 1) * af_dim1]), abs(d__2)), d__4 = 
00292                             work[k + 1];
00293                     work[k + 1] = max(d__3,d__4);
00294                 }
00295 /* Computing MAX */
00296                 i__1 = k + k * af_dim1;
00297                 d__3 = (d__1 = af[i__1].r, abs(d__1)) + (d__2 = d_imag(&af[k 
00298                         + k * af_dim1]), abs(d__2)), d__4 = work[k];
00299                 work[k] = max(d__3,d__4);
00300                 k += 2;
00301             }
00302         }
00303         k = ncols;
00304         while(k >= 1) {
00305             if (ipiv[k] > 0) {
00306                 kp = ipiv[k];
00307                 if (kp != k) {
00308                     tmp = work[*n + k];
00309                     work[*n + k] = work[*n + kp];
00310                     work[*n + kp] = tmp;
00311                 }
00312                 --k;
00313             } else {
00314                 kp = -ipiv[k];
00315                 tmp = work[*n + k];
00316                 work[*n + k] = work[*n + kp];
00317                 work[*n + kp] = tmp;
00318                 k += -2;
00319             }
00320         }
00321     }
00322 
00323 /*     Compute the *inverse* of the max element growth factor.  Dividing */
00324 /*     by zero would imply the largest entry of the factor's column is */
00325 /*     zero.  Than can happen when either the column of A is zero or */
00326 /*     massive pivots made the factor underflow to zero.  Neither counts */
00327 /*     as growth in itself, so simply ignore terms with zero */
00328 /*     denominators. */
00329 
00330     if (upper) {
00331         i__1 = *n;
00332         for (i__ = ncols; i__ <= i__1; ++i__) {
00333             umax = work[i__];
00334             amax = work[*n + i__];
00335             if (umax != 0.) {
00336 /* Computing MIN */
00337                 d__1 = amax / umax;
00338                 rpvgrw = min(d__1,rpvgrw);
00339             }
00340         }
00341     } else {
00342         i__1 = ncols;
00343         for (i__ = 1; i__ <= i__1; ++i__) {
00344             umax = work[i__];
00345             amax = work[*n + i__];
00346             if (umax != 0.) {
00347 /* Computing MIN */
00348                 d__1 = amax / umax;
00349                 rpvgrw = min(d__1,rpvgrw);
00350             }
00351         }
00352     }
00353     ret_val = rpvgrw;
00354     return ret_val;
00355 } /* zla_herpvgrw__ */


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