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/cgees.c | |
parent | 01f64c1ecd0d4ffa9e3a74478335f1745f26cc75 (diff) | |
download | ydb-90d450f74722da7859d6f510a869f6c6908fd12f.tar.gz |
[] add metering mode to CLI
Diffstat (limited to 'contrib/libs/clapack/cgees.c')
-rw-r--r-- | contrib/libs/clapack/cgees.c | 404 |
1 files changed, 404 insertions, 0 deletions
diff --git a/contrib/libs/clapack/cgees.c b/contrib/libs/clapack/cgees.c new file mode 100644 index 0000000000..d113eb3749 --- /dev/null +++ b/contrib/libs/clapack/cgees.c @@ -0,0 +1,404 @@ +/* cgees.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 integer c__1 = 1; +static integer c__0 = 0; +static integer c_n1 = -1; + +/* Subroutine */ int cgees_(char *jobvs, char *sort, L_fp select, integer *n, + complex *a, integer *lda, integer *sdim, complex *w, complex *vs, + integer *ldvs, complex *work, integer *lwork, real *rwork, logical * + bwork, integer *info) +{ + /* System generated locals */ + integer a_dim1, a_offset, vs_dim1, vs_offset, i__1, i__2; + + /* Builtin functions */ + double sqrt(doublereal); + + /* Local variables */ + integer i__; + real s; + integer ihi, ilo; + real dum[1], eps, sep; + integer ibal; + real anrm; + integer ierr, itau, iwrk, icond, ieval; + extern logical lsame_(char *, char *); + extern /* Subroutine */ int ccopy_(integer *, complex *, integer *, + complex *, integer *), cgebak_(char *, char *, integer *, integer + *, integer *, real *, integer *, complex *, integer *, integer *), cgebal_(char *, integer *, complex *, integer *, + integer *, integer *, real *, integer *), slabad_(real *, + real *); + logical scalea; + extern doublereal clange_(char *, integer *, integer *, complex *, + integer *, real *); + real cscale; + extern /* Subroutine */ int cgehrd_(integer *, integer *, integer *, + complex *, integer *, complex *, complex *, integer *, integer *), + clascl_(char *, integer *, integer *, real *, real *, integer *, + integer *, complex *, integer *, integer *); + extern doublereal slamch_(char *); + extern /* Subroutine */ int clacpy_(char *, integer *, integer *, complex + *, integer *, complex *, integer *), xerbla_(char *, + integer *); + extern integer ilaenv_(integer *, char *, char *, integer *, integer *, + integer *, integer *); + real bignum; + extern /* Subroutine */ int chseqr_(char *, char *, integer *, integer *, + integer *, complex *, integer *, complex *, complex *, integer *, + complex *, integer *, integer *), cunghr_(integer + *, integer *, integer *, complex *, integer *, complex *, complex + *, integer *, integer *), ctrsen_(char *, char *, logical *, + integer *, complex *, integer *, complex *, integer *, complex *, + integer *, real *, real *, complex *, integer *, integer *); + integer minwrk, maxwrk; + real smlnum; + integer hswork; + logical wantst, lquery, wantvs; + + +/* -- LAPACK driver routine (version 3.2) -- */ +/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ +/* November 2006 */ + +/* .. Scalar Arguments .. */ +/* .. */ +/* .. Array Arguments .. */ +/* .. */ +/* .. Function Arguments .. */ +/* .. */ + +/* Purpose */ +/* ======= */ + +/* CGEES computes for an N-by-N complex nonsymmetric matrix A, the */ +/* eigenvalues, the Schur form T, and, optionally, the matrix of Schur */ +/* vectors Z. This gives the Schur factorization A = Z*T*(Z**H). */ + +/* Optionally, it also orders the eigenvalues on the diagonal of the */ +/* Schur form so that selected eigenvalues are at the top left. */ +/* The leading columns of Z then form an orthonormal basis for the */ +/* invariant subspace corresponding to the selected eigenvalues. */ +/* A complex matrix is in Schur form if it is upper triangular. */ + +/* Arguments */ +/* ========= */ + +/* JOBVS (input) CHARACTER*1 */ +/* = 'N': Schur vectors are not computed; */ +/* = 'V': Schur vectors are computed. */ + +/* SORT (input) CHARACTER*1 */ +/* Specifies whether or not to order the eigenvalues on the */ +/* diagonal of the Schur form. */ +/* = 'N': Eigenvalues are not ordered: */ +/* = 'S': Eigenvalues are ordered (see SELECT). */ + +/* SELECT (external procedure) LOGICAL FUNCTION of one COMPLEX argument */ +/* SELECT must be declared EXTERNAL in the calling subroutine. */ +/* If SORT = 'S', SELECT is used to select eigenvalues to order */ +/* to the top left of the Schur form. */ +/* IF SORT = 'N', SELECT is not referenced. */ +/* The eigenvalue W(j) is selected if SELECT(W(j)) is true. */ + +/* N (input) INTEGER */ +/* The order of the matrix A. N >= 0. */ + +/* A (input/output) COMPLEX array, dimension (LDA,N) */ +/* On entry, the N-by-N matrix A. */ +/* On exit, A has been overwritten by its Schur form T. */ + +/* LDA (input) INTEGER */ +/* The leading dimension of the array A. LDA >= max(1,N). */ + +/* SDIM (output) INTEGER */ +/* If SORT = 'N', SDIM = 0. */ +/* If SORT = 'S', SDIM = number of eigenvalues for which */ +/* SELECT is true. */ + +/* W (output) COMPLEX array, dimension (N) */ +/* W contains the computed eigenvalues, in the same order that */ +/* they appear on the diagonal of the output Schur form T. */ + +/* VS (output) COMPLEX array, dimension (LDVS,N) */ +/* If JOBVS = 'V', VS contains the unitary matrix Z of Schur */ +/* vectors. */ +/* If JOBVS = 'N', VS is not referenced. */ + +/* LDVS (input) INTEGER */ +/* The leading dimension of the array VS. LDVS >= 1; if */ +/* JOBVS = 'V', LDVS >= N. */ + +/* WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) */ +/* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */ + +/* LWORK (input) INTEGER */ +/* The dimension of the array WORK. LWORK >= max(1,2*N). */ +/* For good performance, LWORK must generally be larger. */ + +/* If LWORK = -1, then a workspace query is assumed; the routine */ +/* only calculates the optimal size of the WORK array, returns */ +/* this value as the first entry of the WORK array, and no error */ +/* message related to LWORK is issued by XERBLA. */ + +/* RWORK (workspace) REAL array, dimension (N) */ + +/* BWORK (workspace) LOGICAL array, dimension (N) */ +/* Not referenced if SORT = 'N'. */ + +/* INFO (output) INTEGER */ +/* = 0: successful exit */ +/* < 0: if INFO = -i, the i-th argument had an illegal value. */ +/* > 0: if INFO = i, and i is */ +/* <= N: the QR algorithm failed to compute all the */ +/* eigenvalues; elements 1:ILO-1 and i+1:N of W */ +/* contain those eigenvalues which have converged; */ +/* if JOBVS = 'V', VS contains the matrix which */ +/* reduces A to its partially converged Schur form. */ +/* = N+1: the eigenvalues could not be reordered because */ +/* some eigenvalues were too close to separate (the */ +/* problem is very ill-conditioned); */ +/* = N+2: after reordering, roundoff changed values of */ +/* some complex eigenvalues so that leading */ +/* eigenvalues in the Schur form no longer satisfy */ +/* SELECT = .TRUE.. This could also be caused by */ +/* underflow due to scaling. */ + +/* ===================================================================== */ + +/* .. Parameters .. */ +/* .. */ +/* .. Local Scalars .. */ +/* .. */ +/* .. Local Arrays .. */ +/* .. */ +/* .. External Subroutines .. */ +/* .. */ +/* .. External Functions .. */ +/* .. */ +/* .. Intrinsic Functions .. */ +/* .. */ +/* .. Executable Statements .. */ + +/* Test the input arguments */ + + /* Parameter adjustments */ + a_dim1 = *lda; + a_offset = 1 + a_dim1; + a -= a_offset; + --w; + vs_dim1 = *ldvs; + vs_offset = 1 + vs_dim1; + vs -= vs_offset; + --work; + --rwork; + --bwork; + + /* Function Body */ + *info = 0; + lquery = *lwork == -1; + wantvs = lsame_(jobvs, "V"); + wantst = lsame_(sort, "S"); + if (! wantvs && ! lsame_(jobvs, "N")) { + *info = -1; + } else if (! wantst && ! lsame_(sort, "N")) { + *info = -2; + } else if (*n < 0) { + *info = -4; + } else if (*lda < max(1,*n)) { + *info = -6; + } else if (*ldvs < 1 || wantvs && *ldvs < *n) { + *info = -10; + } + +/* Compute workspace */ +/* (Note: Comments in the code beginning "Workspace:" describe the */ +/* minimal amount of workspace needed at that point in the code, */ +/* as well as the preferred amount for good performance. */ +/* CWorkspace refers to complex workspace, and RWorkspace to real */ +/* workspace. NB refers to the optimal block size for the */ +/* immediately following subroutine, as returned by ILAENV. */ +/* HSWORK refers to the workspace preferred by CHSEQR, as */ +/* calculated below. HSWORK is computed assuming ILO=1 and IHI=N, */ +/* the worst case.) */ + + if (*info == 0) { + if (*n == 0) { + minwrk = 1; + maxwrk = 1; + } else { + maxwrk = *n + *n * ilaenv_(&c__1, "CGEHRD", " ", n, &c__1, n, & + c__0); + minwrk = *n << 1; + + chseqr_("S", jobvs, n, &c__1, n, &a[a_offset], lda, &w[1], &vs[ + vs_offset], ldvs, &work[1], &c_n1, &ieval); + hswork = work[1].r; + + if (! wantvs) { + maxwrk = max(maxwrk,hswork); + } else { +/* Computing MAX */ + i__1 = maxwrk, i__2 = *n + (*n - 1) * ilaenv_(&c__1, "CUNGHR", + " ", n, &c__1, n, &c_n1); + maxwrk = max(i__1,i__2); + maxwrk = max(maxwrk,hswork); + } + } + work[1].r = (real) maxwrk, work[1].i = 0.f; + + if (*lwork < minwrk && ! lquery) { + *info = -12; + } + } + + if (*info != 0) { + i__1 = -(*info); + xerbla_("CGEES ", &i__1); + return 0; + } else if (lquery) { + return 0; + } + +/* Quick return if possible */ + + if (*n == 0) { + *sdim = 0; + return 0; + } + +/* Get machine constants */ + + eps = slamch_("P"); + smlnum = slamch_("S"); + bignum = 1.f / smlnum; + slabad_(&smlnum, &bignum); + smlnum = sqrt(smlnum) / eps; + bignum = 1.f / smlnum; + +/* Scale A if max element outside range [SMLNUM,BIGNUM] */ + + anrm = clange_("M", n, n, &a[a_offset], lda, dum); + scalea = FALSE_; + if (anrm > 0.f && anrm < smlnum) { + scalea = TRUE_; + cscale = smlnum; + } else if (anrm > bignum) { + scalea = TRUE_; + cscale = bignum; + } + if (scalea) { + clascl_("G", &c__0, &c__0, &anrm, &cscale, n, n, &a[a_offset], lda, & + ierr); + } + +/* Permute the matrix to make it more nearly triangular */ +/* (CWorkspace: none) */ +/* (RWorkspace: need N) */ + + ibal = 1; + cgebal_("P", n, &a[a_offset], lda, &ilo, &ihi, &rwork[ibal], &ierr); + +/* Reduce to upper Hessenberg form */ +/* (CWorkspace: need 2*N, prefer N+N*NB) */ +/* (RWorkspace: none) */ + + itau = 1; + iwrk = *n + itau; + i__1 = *lwork - iwrk + 1; + cgehrd_(n, &ilo, &ihi, &a[a_offset], lda, &work[itau], &work[iwrk], &i__1, + &ierr); + + if (wantvs) { + +/* Copy Householder vectors to VS */ + + clacpy_("L", n, n, &a[a_offset], lda, &vs[vs_offset], ldvs) + ; + +/* Generate unitary matrix in VS */ +/* (CWorkspace: need 2*N-1, prefer N+(N-1)*NB) */ +/* (RWorkspace: none) */ + + i__1 = *lwork - iwrk + 1; + cunghr_(n, &ilo, &ihi, &vs[vs_offset], ldvs, &work[itau], &work[iwrk], + &i__1, &ierr); + } + + *sdim = 0; + +/* Perform QR iteration, accumulating Schur vectors in VS if desired */ +/* (CWorkspace: need 1, prefer HSWORK (see comments) ) */ +/* (RWorkspace: none) */ + + iwrk = itau; + i__1 = *lwork - iwrk + 1; + chseqr_("S", jobvs, n, &ilo, &ihi, &a[a_offset], lda, &w[1], &vs[ + vs_offset], ldvs, &work[iwrk], &i__1, &ieval); + if (ieval > 0) { + *info = ieval; + } + +/* Sort eigenvalues if desired */ + + if (wantst && *info == 0) { + if (scalea) { + clascl_("G", &c__0, &c__0, &cscale, &anrm, n, &c__1, &w[1], n, & + ierr); + } + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + bwork[i__] = (*select)(&w[i__]); +/* L10: */ + } + +/* Reorder eigenvalues and transform Schur vectors */ +/* (CWorkspace: none) */ +/* (RWorkspace: none) */ + + i__1 = *lwork - iwrk + 1; + ctrsen_("N", jobvs, &bwork[1], n, &a[a_offset], lda, &vs[vs_offset], + ldvs, &w[1], sdim, &s, &sep, &work[iwrk], &i__1, &icond); + } + + if (wantvs) { + +/* Undo balancing */ +/* (CWorkspace: none) */ +/* (RWorkspace: need N) */ + + cgebak_("P", "R", n, &ilo, &ihi, &rwork[ibal], n, &vs[vs_offset], + ldvs, &ierr); + } + + if (scalea) { + +/* Undo scaling for the Schur form of A */ + + clascl_("U", &c__0, &c__0, &cscale, &anrm, n, n, &a[a_offset], lda, & + ierr); + i__1 = *lda + 1; + ccopy_(n, &a[a_offset], &i__1, &w[1], &c__1); + } + + work[1].r = (real) maxwrk, work[1].i = 0.f; + return 0; + +/* End of CGEES */ + +} /* cgees_ */ |