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/dgeqp3.c | |
parent | 01f64c1ecd0d4ffa9e3a74478335f1745f26cc75 (diff) | |
download | ydb-90d450f74722da7859d6f510a869f6c6908fd12f.tar.gz |
[] add metering mode to CLI
Diffstat (limited to 'contrib/libs/clapack/dgeqp3.c')
-rw-r--r-- | contrib/libs/clapack/dgeqp3.c | 358 |
1 files changed, 358 insertions, 0 deletions
diff --git a/contrib/libs/clapack/dgeqp3.c b/contrib/libs/clapack/dgeqp3.c new file mode 100644 index 0000000000..e3fda572c2 --- /dev/null +++ b/contrib/libs/clapack/dgeqp3.c @@ -0,0 +1,358 @@ +/* dgeqp3.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_n1 = -1; +static integer c__3 = 3; +static integer c__2 = 2; + +/* Subroutine */ int dgeqp3_(integer *m, integer *n, doublereal *a, integer * + lda, integer *jpvt, doublereal *tau, doublereal *work, integer *lwork, + integer *info) +{ + /* System generated locals */ + integer a_dim1, a_offset, i__1, i__2, i__3; + + /* Local variables */ + integer j, jb, na, nb, sm, sn, nx, fjb, iws, nfxd; + extern doublereal dnrm2_(integer *, doublereal *, integer *); + integer nbmin, minmn; + extern /* Subroutine */ int dswap_(integer *, doublereal *, integer *, + doublereal *, integer *); + integer minws; + extern /* Subroutine */ int dlaqp2_(integer *, integer *, integer *, + doublereal *, integer *, integer *, doublereal *, doublereal *, + doublereal *, doublereal *), dgeqrf_(integer *, integer *, + doublereal *, integer *, doublereal *, doublereal *, integer *, + integer *), xerbla_(char *, integer *); + extern integer ilaenv_(integer *, char *, char *, integer *, integer *, + integer *, integer *); + extern /* Subroutine */ int dlaqps_(integer *, integer *, integer *, + integer *, integer *, doublereal *, integer *, integer *, + doublereal *, doublereal *, doublereal *, doublereal *, + doublereal *, integer *); + integer topbmn, sminmn; + extern /* Subroutine */ int dormqr_(char *, char *, integer *, integer *, + integer *, doublereal *, integer *, doublereal *, doublereal *, + integer *, doublereal *, integer *, integer *); + integer lwkopt; + logical lquery; + + +/* -- LAPACK routine (version 3.2) -- */ +/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ +/* November 2006 */ + +/* .. Scalar Arguments .. */ +/* .. */ +/* .. Array Arguments .. */ +/* .. */ + +/* Purpose */ +/* ======= */ + +/* DGEQP3 computes a QR factorization with column pivoting of a */ +/* matrix A: A*P = Q*R using Level 3 BLAS. */ + +/* Arguments */ +/* ========= */ + +/* M (input) INTEGER */ +/* The number of rows of the matrix A. M >= 0. */ + +/* N (input) INTEGER */ +/* The number of columns of the matrix A. N >= 0. */ + +/* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */ +/* On entry, the M-by-N matrix A. */ +/* On exit, the upper triangle of the array contains the */ +/* min(M,N)-by-N upper trapezoidal matrix R; the elements below */ +/* the diagonal, together with the array TAU, represent the */ +/* orthogonal matrix Q as a product of min(M,N) elementary */ +/* reflectors. */ + +/* LDA (input) INTEGER */ +/* The leading dimension of the array A. LDA >= max(1,M). */ + +/* JPVT (input/output) INTEGER array, dimension (N) */ +/* On entry, if JPVT(J).ne.0, the J-th column of A is permuted */ +/* to the front of A*P (a leading column); if JPVT(J)=0, */ +/* the J-th column of A is a free column. */ +/* On exit, if JPVT(J)=K, then the J-th column of A*P was the */ +/* the K-th column of A. */ + +/* TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) */ +/* The scalar factors of the elementary reflectors. */ + +/* WORK (workspace/output) DOUBLE PRECISION 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 >= 3*N+1. */ +/* For optimal performance LWORK >= 2*N+( N+1 )*NB, where NB */ +/* is the optimal blocksize. */ + +/* 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. */ + +/* INFO (output) INTEGER */ +/* = 0: successful exit. */ +/* < 0: if INFO = -i, the i-th argument had an illegal value. */ + +/* Further Details */ +/* =============== */ + +/* The matrix Q is represented as a product of elementary reflectors */ + +/* Q = H(1) H(2) . . . H(k), where k = min(m,n). */ + +/* Each H(i) has the form */ + +/* H(i) = I - tau * v * v' */ + +/* where tau is a real/complex scalar, and v is a real/complex vector */ +/* with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in */ +/* A(i+1:m,i), and tau in TAU(i). */ + +/* Based on contributions by */ +/* G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain */ +/* X. Sun, Computer Science Dept., Duke University, USA */ + +/* ===================================================================== */ + +/* .. Parameters .. */ +/* .. */ +/* .. Local Scalars .. */ +/* .. */ +/* .. External Subroutines .. */ +/* .. */ +/* .. External Functions .. */ +/* .. */ +/* .. Intrinsic Functions .. */ +/* .. */ +/* .. Executable Statements .. */ + +/* Test input arguments */ +/* ==================== */ + + /* Parameter adjustments */ + a_dim1 = *lda; + a_offset = 1 + a_dim1; + a -= a_offset; + --jpvt; + --tau; + --work; + + /* Function Body */ + *info = 0; + lquery = *lwork == -1; + if (*m < 0) { + *info = -1; + } else if (*n < 0) { + *info = -2; + } else if (*lda < max(1,*m)) { + *info = -4; + } + + if (*info == 0) { + minmn = min(*m,*n); + if (minmn == 0) { + iws = 1; + lwkopt = 1; + } else { + iws = *n * 3 + 1; + nb = ilaenv_(&c__1, "DGEQRF", " ", m, n, &c_n1, &c_n1); + lwkopt = (*n << 1) + (*n + 1) * nb; + } + work[1] = (doublereal) lwkopt; + + if (*lwork < iws && ! lquery) { + *info = -8; + } + } + + if (*info != 0) { + i__1 = -(*info); + xerbla_("DGEQP3", &i__1); + return 0; + } else if (lquery) { + return 0; + } + +/* Quick return if possible. */ + + if (minmn == 0) { + return 0; + } + +/* Move initial columns up front. */ + + nfxd = 1; + i__1 = *n; + for (j = 1; j <= i__1; ++j) { + if (jpvt[j] != 0) { + if (j != nfxd) { + dswap_(m, &a[j * a_dim1 + 1], &c__1, &a[nfxd * a_dim1 + 1], & + c__1); + jpvt[j] = jpvt[nfxd]; + jpvt[nfxd] = j; + } else { + jpvt[j] = j; + } + ++nfxd; + } else { + jpvt[j] = j; + } +/* L10: */ + } + --nfxd; + +/* Factorize fixed columns */ +/* ======================= */ + +/* Compute the QR factorization of fixed columns and update */ +/* remaining columns. */ + + if (nfxd > 0) { + na = min(*m,nfxd); +/* CC CALL DGEQR2( M, NA, A, LDA, TAU, WORK, INFO ) */ + dgeqrf_(m, &na, &a[a_offset], lda, &tau[1], &work[1], lwork, info); +/* Computing MAX */ + i__1 = iws, i__2 = (integer) work[1]; + iws = max(i__1,i__2); + if (na < *n) { +/* CC CALL DORM2R( 'Left', 'Transpose', M, N-NA, NA, A, LDA, */ +/* CC $ TAU, A( 1, NA+1 ), LDA, WORK, INFO ) */ + i__1 = *n - na; + dormqr_("Left", "Transpose", m, &i__1, &na, &a[a_offset], lda, & + tau[1], &a[(na + 1) * a_dim1 + 1], lda, &work[1], lwork, + info); +/* Computing MAX */ + i__1 = iws, i__2 = (integer) work[1]; + iws = max(i__1,i__2); + } + } + +/* Factorize free columns */ +/* ====================== */ + + if (nfxd < minmn) { + + sm = *m - nfxd; + sn = *n - nfxd; + sminmn = minmn - nfxd; + +/* Determine the block size. */ + + nb = ilaenv_(&c__1, "DGEQRF", " ", &sm, &sn, &c_n1, &c_n1); + nbmin = 2; + nx = 0; + + if (nb > 1 && nb < sminmn) { + +/* Determine when to cross over from blocked to unblocked code. */ + +/* Computing MAX */ + i__1 = 0, i__2 = ilaenv_(&c__3, "DGEQRF", " ", &sm, &sn, &c_n1, & + c_n1); + nx = max(i__1,i__2); + + + if (nx < sminmn) { + +/* Determine if workspace is large enough for blocked code. */ + + minws = (sn << 1) + (sn + 1) * nb; + iws = max(iws,minws); + if (*lwork < minws) { + +/* Not enough workspace to use optimal NB: Reduce NB and */ +/* determine the minimum value of NB. */ + + nb = (*lwork - (sn << 1)) / (sn + 1); +/* Computing MAX */ + i__1 = 2, i__2 = ilaenv_(&c__2, "DGEQRF", " ", &sm, &sn, & + c_n1, &c_n1); + nbmin = max(i__1,i__2); + + + } + } + } + +/* Initialize partial column norms. The first N elements of work */ +/* store the exact column norms. */ + + i__1 = *n; + for (j = nfxd + 1; j <= i__1; ++j) { + work[j] = dnrm2_(&sm, &a[nfxd + 1 + j * a_dim1], &c__1); + work[*n + j] = work[j]; +/* L20: */ + } + + if (nb >= nbmin && nb < sminmn && nx < sminmn) { + +/* Use blocked code initially. */ + + j = nfxd + 1; + +/* Compute factorization: while loop. */ + + + topbmn = minmn - nx; +L30: + if (j <= topbmn) { +/* Computing MIN */ + i__1 = nb, i__2 = topbmn - j + 1; + jb = min(i__1,i__2); + +/* Factorize JB columns among columns J:N. */ + + i__1 = *n - j + 1; + i__2 = j - 1; + i__3 = *n - j + 1; + dlaqps_(m, &i__1, &i__2, &jb, &fjb, &a[j * a_dim1 + 1], lda, & + jpvt[j], &tau[j], &work[j], &work[*n + j], &work[(*n + << 1) + 1], &work[(*n << 1) + jb + 1], &i__3); + + j += fjb; + goto L30; + } + } else { + j = nfxd + 1; + } + +/* Use unblocked code to factor the last or only block. */ + + + if (j <= minmn) { + i__1 = *n - j + 1; + i__2 = j - 1; + dlaqp2_(m, &i__1, &i__2, &a[j * a_dim1 + 1], lda, &jpvt[j], &tau[ + j], &work[j], &work[*n + j], &work[(*n << 1) + 1]); + } + + } + + work[1] = (doublereal) iws; + return 0; + +/* End of DGEQP3 */ + +} /* dgeqp3_ */ |