diff options
author | shmel1k <shmel1k@ydb.tech> | 2022-09-02 12:44:59 +0300 |
---|---|---|
committer | shmel1k <shmel1k@ydb.tech> | 2022-09-02 12:44:59 +0300 |
commit | 90d450f74722da7859d6f510a869f6c6908fd12f (patch) | |
tree | 538c718dedc76cdfe37ad6d01ff250dd930d9278 /contrib/libs/clapack/zhsein.c | |
parent | 01f64c1ecd0d4ffa9e3a74478335f1745f26cc75 (diff) | |
download | ydb-90d450f74722da7859d6f510a869f6c6908fd12f.tar.gz |
[] add metering mode to CLI
Diffstat (limited to 'contrib/libs/clapack/zhsein.c')
-rw-r--r-- | contrib/libs/clapack/zhsein.c | 433 |
1 files changed, 433 insertions, 0 deletions
diff --git a/contrib/libs/clapack/zhsein.c b/contrib/libs/clapack/zhsein.c new file mode 100644 index 0000000000..0e27b9bea9 --- /dev/null +++ b/contrib/libs/clapack/zhsein.c @@ -0,0 +1,433 @@ +/* zhsein.f -- translated by f2c (version 20061008). + You must link the resulting object file with libf2c: + on Microsoft Windows system, link with libf2c.lib; + on Linux or Unix systems, link with .../path/to/libf2c.a -lm + or, if you install libf2c.a in a standard place, with -lf2c -lm + -- in that order, at the end of the command line, as in + cc *.o -lf2c -lm + Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., + + http://www.netlib.org/f2c/libf2c.zip +*/ + +#include "f2c.h" +#include "blaswrap.h" + +/* Table of constant values */ + +static logical c_false = FALSE_; +static logical c_true = TRUE_; + +/* Subroutine */ int zhsein_(char *side, char *eigsrc, char *initv, logical * + select, integer *n, doublecomplex *h__, integer *ldh, doublecomplex * + w, doublecomplex *vl, integer *ldvl, doublecomplex *vr, integer *ldvr, + integer *mm, integer *m, doublecomplex *work, doublereal *rwork, + integer *ifaill, integer *ifailr, integer *info) +{ + /* System generated locals */ + integer h_dim1, h_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, i__1, + i__2, i__3; + doublereal d__1, d__2; + doublecomplex z__1, z__2; + + /* Builtin functions */ + double d_imag(doublecomplex *); + + /* Local variables */ + integer i__, k, kl, kr, ks; + doublecomplex wk; + integer kln; + doublereal ulp, eps3, unfl; + extern logical lsame_(char *, char *); + integer iinfo; + logical leftv, bothv; + doublereal hnorm; + extern doublereal dlamch_(char *); + extern /* Subroutine */ int xerbla_(char *, integer *), zlaein_( + logical *, logical *, integer *, doublecomplex *, integer *, + doublecomplex *, doublecomplex *, doublecomplex *, integer *, + doublereal *, doublereal *, doublereal *, integer *); + extern doublereal zlanhs_(char *, integer *, doublecomplex *, integer *, + doublereal *); + logical noinit; + integer ldwork; + logical rightv, fromqr; + doublereal smlnum; + + +/* -- LAPACK routine (version 3.2) -- */ +/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ +/* November 2006 */ + +/* .. Scalar Arguments .. */ +/* .. */ +/* .. Array Arguments .. */ +/* .. */ + +/* Purpose */ +/* ======= */ + +/* ZHSEIN uses inverse iteration to find specified right and/or left */ +/* eigenvectors of a complex upper Hessenberg matrix H. */ + +/* The right eigenvector x and the left eigenvector y of the matrix H */ +/* corresponding to an eigenvalue w are defined by: */ + +/* H * x = w * x, y**h * H = w * y**h */ + +/* where y**h denotes the conjugate transpose of the vector y. */ + +/* Arguments */ +/* ========= */ + +/* SIDE (input) CHARACTER*1 */ +/* = 'R': compute right eigenvectors only; */ +/* = 'L': compute left eigenvectors only; */ +/* = 'B': compute both right and left eigenvectors. */ + +/* EIGSRC (input) CHARACTER*1 */ +/* Specifies the source of eigenvalues supplied in W: */ +/* = 'Q': the eigenvalues were found using ZHSEQR; thus, if */ +/* H has zero subdiagonal elements, and so is */ +/* block-triangular, then the j-th eigenvalue can be */ +/* assumed to be an eigenvalue of the block containing */ +/* the j-th row/column. This property allows ZHSEIN to */ +/* perform inverse iteration on just one diagonal block. */ +/* = 'N': no assumptions are made on the correspondence */ +/* between eigenvalues and diagonal blocks. In this */ +/* case, ZHSEIN must always perform inverse iteration */ +/* using the whole matrix H. */ + +/* INITV (input) CHARACTER*1 */ +/* = 'N': no initial vectors are supplied; */ +/* = 'U': user-supplied initial vectors are stored in the arrays */ +/* VL and/or VR. */ + +/* SELECT (input) LOGICAL array, dimension (N) */ +/* Specifies the eigenvectors to be computed. To select the */ +/* eigenvector corresponding to the eigenvalue W(j), */ +/* SELECT(j) must be set to .TRUE.. */ + +/* N (input) INTEGER */ +/* The order of the matrix H. N >= 0. */ + +/* H (input) COMPLEX*16 array, dimension (LDH,N) */ +/* The upper Hessenberg matrix H. */ + +/* LDH (input) INTEGER */ +/* The leading dimension of the array H. LDH >= max(1,N). */ + +/* W (input/output) COMPLEX*16 array, dimension (N) */ +/* On entry, the eigenvalues of H. */ +/* On exit, the real parts of W may have been altered since */ +/* close eigenvalues are perturbed slightly in searching for */ +/* independent eigenvectors. */ + +/* VL (input/output) COMPLEX*16 array, dimension (LDVL,MM) */ +/* On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must */ +/* contain starting vectors for the inverse iteration for the */ +/* left eigenvectors; the starting vector for each eigenvector */ +/* must be in the same column in which the eigenvector will be */ +/* stored. */ +/* On exit, if SIDE = 'L' or 'B', the left eigenvectors */ +/* specified by SELECT will be stored consecutively in the */ +/* columns of VL, in the same order as their eigenvalues. */ +/* If SIDE = 'R', VL is not referenced. */ + +/* LDVL (input) INTEGER */ +/* The leading dimension of the array VL. */ +/* LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise. */ + +/* VR (input/output) COMPLEX*16 array, dimension (LDVR,MM) */ +/* On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must */ +/* contain starting vectors for the inverse iteration for the */ +/* right eigenvectors; the starting vector for each eigenvector */ +/* must be in the same column in which the eigenvector will be */ +/* stored. */ +/* On exit, if SIDE = 'R' or 'B', the right eigenvectors */ +/* specified by SELECT will be stored consecutively in the */ +/* columns of VR, in the same order as their eigenvalues. */ +/* If SIDE = 'L', VR is not referenced. */ + +/* LDVR (input) INTEGER */ +/* The leading dimension of the array VR. */ +/* LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise. */ + +/* MM (input) INTEGER */ +/* The number of columns in the arrays VL and/or VR. MM >= M. */ + +/* M (output) INTEGER */ +/* The number of columns in the arrays VL and/or VR required to */ +/* store the eigenvectors (= the number of .TRUE. elements in */ +/* SELECT). */ + +/* WORK (workspace) COMPLEX*16 array, dimension (N*N) */ + +/* RWORK (workspace) DOUBLE PRECISION array, dimension (N) */ + +/* IFAILL (output) INTEGER array, dimension (MM) */ +/* If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left */ +/* eigenvector in the i-th column of VL (corresponding to the */ +/* eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the */ +/* eigenvector converged satisfactorily. */ +/* If SIDE = 'R', IFAILL is not referenced. */ + +/* IFAILR (output) INTEGER array, dimension (MM) */ +/* If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right */ +/* eigenvector in the i-th column of VR (corresponding to the */ +/* eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the */ +/* eigenvector converged satisfactorily. */ +/* If SIDE = 'L', IFAILR is not referenced. */ + +/* INFO (output) INTEGER */ +/* = 0: successful exit */ +/* < 0: if INFO = -i, the i-th argument had an illegal value */ +/* > 0: if INFO = i, i is the number of eigenvectors which */ +/* failed to converge; see IFAILL and IFAILR for further */ +/* details. */ + +/* Further Details */ +/* =============== */ + +/* Each eigenvector is normalized so that the element of largest */ +/* magnitude has magnitude 1; here the magnitude of a complex number */ +/* (x,y) is taken to be |x|+|y|. */ + +/* ===================================================================== */ + +/* .. Parameters .. */ +/* .. */ +/* .. Local Scalars .. */ +/* .. */ +/* .. External Functions .. */ +/* .. */ +/* .. External Subroutines .. */ +/* .. */ +/* .. Intrinsic Functions .. */ +/* .. */ +/* .. Statement Functions .. */ +/* .. */ +/* .. Statement Function definitions .. */ +/* .. */ +/* .. Executable Statements .. */ + +/* Decode and test the input parameters. */ + + /* Parameter adjustments */ + --select; + h_dim1 = *ldh; + h_offset = 1 + h_dim1; + h__ -= h_offset; + --w; + vl_dim1 = *ldvl; + vl_offset = 1 + vl_dim1; + vl -= vl_offset; + vr_dim1 = *ldvr; + vr_offset = 1 + vr_dim1; + vr -= vr_offset; + --work; + --rwork; + --ifaill; + --ifailr; + + /* Function Body */ + bothv = lsame_(side, "B"); + rightv = lsame_(side, "R") || bothv; + leftv = lsame_(side, "L") || bothv; + + fromqr = lsame_(eigsrc, "Q"); + + noinit = lsame_(initv, "N"); + +/* Set M to the number of columns required to store the selected */ +/* eigenvectors. */ + + *m = 0; + i__1 = *n; + for (k = 1; k <= i__1; ++k) { + if (select[k]) { + ++(*m); + } +/* L10: */ + } + + *info = 0; + if (! rightv && ! leftv) { + *info = -1; + } else if (! fromqr && ! lsame_(eigsrc, "N")) { + *info = -2; + } else if (! noinit && ! lsame_(initv, "U")) { + *info = -3; + } else if (*n < 0) { + *info = -5; + } else if (*ldh < max(1,*n)) { + *info = -7; + } else if (*ldvl < 1 || leftv && *ldvl < *n) { + *info = -10; + } else if (*ldvr < 1 || rightv && *ldvr < *n) { + *info = -12; + } else if (*mm < *m) { + *info = -13; + } + if (*info != 0) { + i__1 = -(*info); + xerbla_("ZHSEIN", &i__1); + return 0; + } + +/* Quick return if possible. */ + + if (*n == 0) { + return 0; + } + +/* Set machine-dependent constants. */ + + unfl = dlamch_("Safe minimum"); + ulp = dlamch_("Precision"); + smlnum = unfl * (*n / ulp); + + ldwork = *n; + + kl = 1; + kln = 0; + if (fromqr) { + kr = 0; + } else { + kr = *n; + } + ks = 1; + + i__1 = *n; + for (k = 1; k <= i__1; ++k) { + if (select[k]) { + +/* Compute eigenvector(s) corresponding to W(K). */ + + if (fromqr) { + +/* If affiliation of eigenvalues is known, check whether */ +/* the matrix splits. */ + +/* Determine KL and KR such that 1 <= KL <= K <= KR <= N */ +/* and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or */ +/* KR = N). */ + +/* Then inverse iteration can be performed with the */ +/* submatrix H(KL:N,KL:N) for a left eigenvector, and with */ +/* the submatrix H(1:KR,1:KR) for a right eigenvector. */ + + i__2 = kl + 1; + for (i__ = k; i__ >= i__2; --i__) { + i__3 = i__ + (i__ - 1) * h_dim1; + if (h__[i__3].r == 0. && h__[i__3].i == 0.) { + goto L30; + } +/* L20: */ + } +L30: + kl = i__; + if (k > kr) { + i__2 = *n - 1; + for (i__ = k; i__ <= i__2; ++i__) { + i__3 = i__ + 1 + i__ * h_dim1; + if (h__[i__3].r == 0. && h__[i__3].i == 0.) { + goto L50; + } +/* L40: */ + } +L50: + kr = i__; + } + } + + if (kl != kln) { + kln = kl; + +/* Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it */ +/* has not ben computed before. */ + + i__2 = kr - kl + 1; + hnorm = zlanhs_("I", &i__2, &h__[kl + kl * h_dim1], ldh, & + rwork[1]); + if (hnorm > 0.) { + eps3 = hnorm * ulp; + } else { + eps3 = smlnum; + } + } + +/* Perturb eigenvalue if it is close to any previous */ +/* selected eigenvalues affiliated to the submatrix */ +/* H(KL:KR,KL:KR). Close roots are modified by EPS3. */ + + i__2 = k; + wk.r = w[i__2].r, wk.i = w[i__2].i; +L60: + i__2 = kl; + for (i__ = k - 1; i__ >= i__2; --i__) { + i__3 = i__; + z__2.r = w[i__3].r - wk.r, z__2.i = w[i__3].i - wk.i; + z__1.r = z__2.r, z__1.i = z__2.i; + if (select[i__] && (d__1 = z__1.r, abs(d__1)) + (d__2 = + d_imag(&z__1), abs(d__2)) < eps3) { + z__1.r = wk.r + eps3, z__1.i = wk.i; + wk.r = z__1.r, wk.i = z__1.i; + goto L60; + } +/* L70: */ + } + i__2 = k; + w[i__2].r = wk.r, w[i__2].i = wk.i; + + if (leftv) { + +/* Compute left eigenvector. */ + + i__2 = *n - kl + 1; + zlaein_(&c_false, &noinit, &i__2, &h__[kl + kl * h_dim1], ldh, + &wk, &vl[kl + ks * vl_dim1], &work[1], &ldwork, & + rwork[1], &eps3, &smlnum, &iinfo); + if (iinfo > 0) { + ++(*info); + ifaill[ks] = k; + } else { + ifaill[ks] = 0; + } + i__2 = kl - 1; + for (i__ = 1; i__ <= i__2; ++i__) { + i__3 = i__ + ks * vl_dim1; + vl[i__3].r = 0., vl[i__3].i = 0.; +/* L80: */ + } + } + if (rightv) { + +/* Compute right eigenvector. */ + + zlaein_(&c_true, &noinit, &kr, &h__[h_offset], ldh, &wk, &vr[ + ks * vr_dim1 + 1], &work[1], &ldwork, &rwork[1], & + eps3, &smlnum, &iinfo); + if (iinfo > 0) { + ++(*info); + ifailr[ks] = k; + } else { + ifailr[ks] = 0; + } + i__2 = *n; + for (i__ = kr + 1; i__ <= i__2; ++i__) { + i__3 = i__ + ks * vr_dim1; + vr[i__3].r = 0., vr[i__3].i = 0.; +/* L90: */ + } + } + ++ks; + } +/* L100: */ + } + + return 0; + +/* End of ZHSEIN */ + +} /* zhsein_ */ |