.TH "SRC/zgeqp3rk.f" 3 "Version 3.12.0" "LAPACK" \" -*- nroff -*-
.ad l
.nh
.SH NAME
SRC/zgeqp3rk.f
.SH SYNOPSIS
.br
.PP
.SS "Functions/Subroutines"

.in +1c
.ti -1c
.RI "subroutine \fBzgeqp3rk\fP (m, n, nrhs, kmax, abstol, reltol, a, lda, k, maxc2nrmk, relmaxc2nrmk, jpiv, tau, work, lwork, rwork, iwork, info)"
.br
.RI "\fBZGEQP3RK\fP computes a truncated Householder QR factorization with column pivoting of a complex m-by-n matrix A by using Level 3 BLAS and overwrites m-by-nrhs matrix B with Q**H * B\&. "
.in -1c
.SH "Function/Subroutine Documentation"
.PP 
.SS "subroutine zgeqp3rk (integer m, integer n, integer nrhs, integer kmax, double precision abstol, double precision reltol, complex*16, dimension( lda, * ) a, integer lda, integer k, double precision maxc2nrmk, double precision relmaxc2nrmk, integer, dimension( * ) jpiv, complex*16, dimension( * ) tau, complex*16, dimension( * ) work, integer lwork, double precision, dimension( * ) rwork, integer, dimension( * ) iwork, integer info)"

.PP
\fBZGEQP3RK\fP computes a truncated Householder QR factorization with column pivoting of a complex m-by-n matrix A by using Level 3 BLAS and overwrites m-by-nrhs matrix B with Q**H * B\&.  
.PP
\fBPurpose:\fP
.RS 4

.PP
.nf
!>
!> ZGEQP3RK performs two tasks simultaneously:
!>
!> Task 1: The routine computes a truncated (rank K) or full rank
!> Householder QR factorization with column pivoting of a complex
!> M-by-N matrix A using Level 3 BLAS\&. K is the number of columns
!> that were factorized, i\&.e\&. factorization rank of the
!> factor R, K <= min(M,N)\&.
!>
!>  A * P(K) = Q(K) * R(K)  =
!>
!>        = Q(K) * ( R11(K) R12(K) ) = Q(K) * (   R(K)_approx    )
!>                 ( 0      R22(K) )          ( 0  R(K)_residual ),
!>
!> where:
!>
!>  P(K)            is an N-by-N permutation matrix;
!>  Q(K)            is an M-by-M unitary matrix;
!>  R(K)_approx   = ( R11(K), R12(K) ) is a rank K approximation of the
!>                    full rank factor R with K-by-K upper-triangular
!>                    R11(K) and K-by-N rectangular R12(K)\&. The diagonal
!>                    entries of R11(K) appear in non-increasing order
!>                    of absolute value, and absolute values of all of
!>                    them exceed the maximum column 2-norm of R22(K)
!>                    up to roundoff error\&.
!>  R(K)_residual = R22(K) is the residual of a rank K approximation
!>                    of the full rank factor R\&. It is a
!>                    an (M-K)-by-(N-K) rectangular matrix;
!>  0               is a an (M-K)-by-K zero matrix\&.
!>
!> Task 2: At the same time, the routine overwrites a complex M-by-NRHS
!> matrix B with  Q(K)**H * B  using Level 3 BLAS\&.
!>
!> =====================================================================
!>
!> The matrices A and B are stored on input in the array A as
!> the left and right blocks A(1:M,1:N) and A(1:M, N+1:N+NRHS)
!> respectively\&.
!>
!>                                  N     NRHS
!>             array_A   =   M  [ mat_A, mat_B ]
!>
!> The truncation criteria (i\&.e\&. when to stop the factorization)
!> can be any of the following:
!>
!>   1) The input parameter KMAX, the maximum number of columns
!>      KMAX to factorize, i\&.e\&. the factorization rank is limited
!>      to KMAX\&. If KMAX >= min(M,N), the criterion is not used\&.
!>
!>   2) The input parameter ABSTOL, the absolute tolerance for
!>      the maximum column 2-norm of the residual matrix R22(K)\&. This
!>      means that the factorization stops if this norm is less or
!>      equal to ABSTOL\&. If ABSTOL < 0\&.0, the criterion is not used\&.
!>
!>   3) The input parameter RELTOL, the tolerance for the maximum
!>      column 2-norm matrix of the residual matrix R22(K) divided
!>      by the maximum column 2-norm of the original matrix A, which
!>      is equal to abs(R(1,1))\&. This means that the factorization stops
!>      when the ratio of the maximum column 2-norm of R22(K) to
!>      the maximum column 2-norm of A is less than or equal to RELTOL\&.
!>      If RELTOL < 0\&.0, the criterion is not used\&.
!>
!>   4) In case both stopping criteria ABSTOL or RELTOL are not used,
!>      and when the residual matrix R22(K) is a zero matrix in some
!>      factorization step K\&. ( This stopping criterion is implicit\&. )
!>
!>  The algorithm stops when any of these conditions is first
!>  satisfied, otherwise the whole matrix A is factorized\&.
!>
!>  To factorize the whole matrix A, use the values
!>  KMAX >= min(M,N), ABSTOL < 0\&.0 and RELTOL < 0\&.0\&.
!>
!>  The routine returns:
!>     a) Q(K), R(K)_approx = ( R11(K), R12(K) ),
!>        R(K)_residual = R22(K), P(K), i\&.e\&. the resulting matrices
!>        of the factorization; P(K) is represented by JPIV,
!>        ( if K = min(M,N), R(K)_approx is the full factor R,
!>        and there is no residual matrix R(K)_residual);
!>     b) K, the number of columns that were factorized,
!>        i\&.e\&. factorization rank;
!>     c) MAXC2NRMK, the maximum column 2-norm of the residual
!>        matrix R(K)_residual = R22(K),
!>        ( if K = min(M,N), MAXC2NRMK = 0\&.0 );
!>     d) RELMAXC2NRMK equals MAXC2NRMK divided by MAXC2NRM, the maximum
!>        column 2-norm of the original matrix A, which is equal
!>        to abs(R(1,1)), ( if K = min(M,N), RELMAXC2NRMK = 0\&.0 );
!>     e) Q(K)**H * B, the matrix B with the unitary
!>        transformation Q(K)**H applied on the left\&.
!>
!> The N-by-N permutation matrix P(K) is stored in a compact form in
!> the integer array JPIV\&. For 1 <= j <= N, column j
!> of the matrix A was interchanged with column JPIV(j)\&.
!>
!> The M-by-M unitary matrix Q is represented as a product
!> of elementary Householder reflectors
!>
!>     Q(K) = H(1) *  H(2) * \&. \&. \&. * H(K),
!>
!> where K is the number of columns that were factorized\&.
!>
!> Each H(j) has the form
!>
!>     H(j) = I - tau * v * v**H,
!>
!> where 1 <= j <= K and
!>   I    is an M-by-M identity matrix,
!>   tau  is a complex scalar,
!>   v    is a complex vector with v(1:j-1) = 0 and v(j) = 1\&.
!>
!> v(j+1:M) is stored on exit in A(j+1:M,j) and tau in TAU(j)\&.
!>
!> See the Further Details section for more information\&.
!> 
.fi
.PP
 
.RE
.PP
\fBParameters\fP
.RS 4
\fIM\fP 
.PP
.nf
!>          M is INTEGER
!>          The number of rows of the matrix A\&. M >= 0\&.
!> 
.fi
.PP
.br
\fIN\fP 
.PP
.nf
!>          N is INTEGER
!>          The number of columns of the matrix A\&. N >= 0\&.
!> 
.fi
.PP
.br
\fINRHS\fP 
.PP
.nf
!>          NRHS is INTEGER
!>          The number of right hand sides, i\&.e\&. the number of
!>          columns of the matrix B\&. NRHS >= 0\&.
!> 
.fi
.PP
.br
\fIKMAX\fP 
.PP
.nf
!>          KMAX is INTEGER
!>
!>          The first factorization stopping criterion\&. KMAX >= 0\&.
!>
!>          The maximum number of columns of the matrix A to factorize,
!>          i\&.e\&. the maximum factorization rank\&.
!>
!>          a) If KMAX >= min(M,N), then this stopping criterion
!>                is not used, the routine factorizes columns
!>                depending on ABSTOL and RELTOL\&.
!>
!>          b) If KMAX = 0, then this stopping criterion is
!>                satisfied on input and the routine exits immediately\&.
!>                This means that the factorization is not performed,
!>                the matrices A and B are not modified, and
!>                the matrix A is itself the residual\&.
!> 
.fi
.PP
.br
\fIABSTOL\fP 
.PP
.nf
!>          ABSTOL is DOUBLE PRECISION
!>
!>          The second factorization stopping criterion, cannot be NaN\&.
!>
!>          The absolute tolerance (stopping threshold) for
!>          maximum column 2-norm of the residual matrix R22(K)\&.
!>          The algorithm converges (stops the factorization) when
!>          the maximum column 2-norm of the residual matrix R22(K)
!>          is less than or equal to ABSTOL\&. Let SAFMIN = DLAMCH('S')\&.
!>
!>          a) If ABSTOL is NaN, then no computation is performed
!>                and an error message ( INFO = -5 ) is issued
!>                by XERBLA\&.
!>
!>          b) If ABSTOL < 0\&.0, then this stopping criterion is not
!>                used, the routine factorizes columns depending
!>                on KMAX and RELTOL\&.
!>                This includes the case ABSTOL = -Inf\&.
!>
!>          c) If 0\&.0 <= ABSTOL < 2*SAFMIN, then ABSTOL = 2*SAFMIN
!>                is used\&. This includes the case ABSTOL = -0\&.0\&.
!>
!>          d) If 2*SAFMIN <= ABSTOL then the input value
!>                of ABSTOL is used\&.
!>
!>          Let MAXC2NRM be the maximum column 2-norm of the
!>          whole original matrix A\&.
!>          If ABSTOL chosen above is >= MAXC2NRM, then this
!>          stopping criterion is satisfied on input and routine exits
!>          immediately after MAXC2NRM is computed\&. The routine
!>          returns MAXC2NRM in MAXC2NORMK,
!>          and 1\&.0 in RELMAXC2NORMK\&.
!>          This includes the case ABSTOL = +Inf\&. This means that the
!>          factorization is not performed, the matrices A and B are not
!>          modified, and the matrix A is itself the residual\&.
!> 
.fi
.PP
.br
\fIRELTOL\fP 
.PP
.nf
!>          RELTOL is DOUBLE PRECISION
!>
!>          The third factorization stopping criterion, cannot be NaN\&.
!>
!>          The tolerance (stopping threshold) for the ratio
!>          abs(R(K+1,K+1))/abs(R(1,1)) of the maximum column 2-norm of
!>          the residual matrix R22(K) to the maximum column 2-norm of
!>          the original matrix A\&. The algorithm converges (stops the
!>          factorization), when abs(R(K+1,K+1))/abs(R(1,1)) A is less
!>          than or equal to RELTOL\&. Let EPS = DLAMCH('E')\&.
!>
!>          a) If RELTOL is NaN, then no computation is performed
!>                and an error message ( INFO = -6 ) is issued
!>                by XERBLA\&.
!>
!>          b) If RELTOL < 0\&.0, then this stopping criterion is not
!>                used, the routine factorizes columns depending
!>                on KMAX and ABSTOL\&.
!>                This includes the case RELTOL = -Inf\&.
!>
!>          c) If 0\&.0 <= RELTOL < EPS, then RELTOL = EPS is used\&.
!>                This includes the case RELTOL = -0\&.0\&.
!>
!>          d) If EPS <= RELTOL then the input value of RELTOL
!>                is used\&.
!>
!>          Let MAXC2NRM be the maximum column 2-norm of the
!>          whole original matrix A\&.
!>          If RELTOL chosen above is >= 1\&.0, then this stopping
!>          criterion is satisfied on input and routine exits
!>          immediately after MAXC2NRM is computed\&.
!>          The routine returns MAXC2NRM in MAXC2NORMK,
!>          and 1\&.0 in RELMAXC2NORMK\&.
!>          This includes the case RELTOL = +Inf\&. This means that the
!>          factorization is not performed, the matrices A and B are not
!>          modified, and the matrix A is itself the residual\&.
!>
!>          NOTE: We recommend that RELTOL satisfy
!>                min( 10*max(M,N)*EPS, sqrt(EPS) ) <= RELTOL
!> 
.fi
.PP
.br
\fIA\fP 
.PP
.nf
!>          A is COMPLEX*16 array, dimension (LDA,N+NRHS)
!>
!>          On entry:
!>
!>          a) The subarray A(1:M,1:N) contains the M-by-N matrix A\&.
!>          b) The subarray A(1:M,N+1:N+NRHS) contains the M-by-NRHS
!>             matrix B\&.
!>
!>                                  N     NRHS
!>              array_A   =   M  [ mat_A, mat_B ]
!>
!>          On exit:
!>
!>          a) The subarray A(1:M,1:N) contains parts of the factors
!>             of the matrix A:
!>
!>            1) If K = 0, A(1:M,1:N) contains the original matrix A\&.
!>            2) If K > 0, A(1:M,1:N) contains parts of the
!>            factors:
!>
!>              1\&. The elements below the diagonal of the subarray
!>                 A(1:M,1:K) together with TAU(1:K) represent the
!>                 unitary matrix Q(K) as a product of K Householder
!>                 elementary reflectors\&.
!>
!>              2\&. The elements on and above the diagonal of
!>                 the subarray A(1:K,1:N) contain K-by-N
!>                 upper-trapezoidal matrix
!>                 R(K)_approx = ( R11(K), R12(K) )\&.
!>                 NOTE: If K=min(M,N), i\&.e\&. full rank factorization,
!>                       then R_approx(K) is the full factor R which
!>                       is upper-trapezoidal\&. If, in addition, M>=N,
!>                       then R is upper-triangular\&.
!>
!>              3\&. The subarray A(K+1:M,K+1:N) contains (M-K)-by-(N-K)
!>                 rectangular matrix R(K)_residual = R22(K)\&.
!>
!>          b) If NRHS > 0, the subarray A(1:M,N+1:N+NRHS) contains
!>             the M-by-NRHS product Q(K)**H * B\&.
!> 
.fi
.PP
.br
\fILDA\fP 
.PP
.nf
!>          LDA is INTEGER
!>          The leading dimension of the array A\&. LDA >= max(1,M)\&.
!>          This is the leading dimension for both matrices, A and B\&.
!> 
.fi
.PP
.br
\fIK\fP 
.PP
.nf
!>          K is INTEGER
!>          Factorization rank of the matrix A, i\&.e\&. the rank of
!>          the factor R, which is the same as the number of non-zero
!>          rows of the factor R\&. 0 <= K <= min(M,KMAX,N)\&.
!>
!>          K also represents the number of non-zero Householder
!>          vectors\&.
!>
!>          NOTE: If K = 0, a) the arrays A and B are not modified;
!>                          b) the array TAU(1:min(M,N)) is set to ZERO,
!>                             if the matrix A does not contain NaN,
!>                             otherwise the elements TAU(1:min(M,N))
!>                             are undefined;
!>                          c) the elements of the array JPIV are set
!>                             as follows: for j = 1:N, JPIV(j) = j\&.
!> 
.fi
.PP
.br
\fIMAXC2NRMK\fP 
.PP
.nf
!>          MAXC2NRMK is DOUBLE PRECISION
!>          The maximum column 2-norm of the residual matrix R22(K),
!>          when the factorization stopped at rank K\&. MAXC2NRMK >= 0\&.
!>
!>          a) If K = 0, i\&.e\&. the factorization was not performed,
!>             the matrix A was not modified and is itself a residual
!>             matrix, then MAXC2NRMK equals the maximum column 2-norm
!>             of the original matrix A\&.
!>
!>          b) If 0 < K < min(M,N), then MAXC2NRMK is returned\&.
!>
!>          c) If K = min(M,N), i\&.e\&. the whole matrix A was
!>             factorized and there is no residual matrix,
!>             then MAXC2NRMK = 0\&.0\&.
!>
!>          NOTE: MAXC2NRMK in the factorization step K would equal
!>                R(K+1,K+1) in the next factorization step K+1\&.
!> 
.fi
.PP
.br
\fIRELMAXC2NRMK\fP 
.PP
.nf
!>          RELMAXC2NRMK is DOUBLE PRECISION
!>          The ratio MAXC2NRMK / MAXC2NRM of the maximum column
!>          2-norm of the residual matrix R22(K) (when the factorization
!>          stopped at rank K) to the maximum column 2-norm of the
!>          whole original matrix A\&. RELMAXC2NRMK >= 0\&.
!>
!>          a) If K = 0, i\&.e\&. the factorization was not performed,
!>             the matrix A was not modified and is itself a residual
!>             matrix, then RELMAXC2NRMK = 1\&.0\&.
!>
!>          b) If 0 < K < min(M,N), then
!>                RELMAXC2NRMK = MAXC2NRMK / MAXC2NRM is returned\&.
!>
!>          c) If K = min(M,N), i\&.e\&. the whole matrix A was
!>             factorized and there is no residual matrix,
!>             then RELMAXC2NRMK = 0\&.0\&.
!>
!>         NOTE: RELMAXC2NRMK in the factorization step K would equal
!>               abs(R(K+1,K+1))/abs(R(1,1)) in the next factorization
!>               step K+1\&.
!> 
.fi
.PP
.br
\fIJPIV\fP 
.PP
.nf
!>          JPIV is INTEGER array, dimension (N)
!>          Column pivot indices\&. For 1 <= j <= N, column j
!>          of the matrix A was interchanged with column JPIV(j)\&.
!>
!>          The elements of the array JPIV(1:N) are always set
!>          by the routine, for example, even  when no columns
!>          were factorized, i\&.e\&. when K = 0, the elements are
!>          set as JPIV(j) = j for j = 1:N\&.
!> 
.fi
.PP
.br
\fITAU\fP 
.PP
.nf
!>          TAU is COMPLEX*16 array, dimension (min(M,N))
!>          The scalar factors of the elementary reflectors\&.
!>
!>          If 0 < K <= min(M,N), only the elements TAU(1:K) of
!>          the array TAU are modified by the factorization\&.
!>          After the factorization computed, if no NaN was found
!>          during the factorization, the remaining elements
!>          TAU(K+1:min(M,N)) are set to zero, otherwise the
!>          elements TAU(K+1:min(M,N)) are not set and therefore
!>          undefined\&.
!>          ( If K = 0, all elements of TAU are set to zero, if
!>          the matrix A does not contain NaN\&. )
!> 
.fi
.PP
.br
\fIWORK\fP 
.PP
.nf
!>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK\&.
!> 
.fi
.PP
.br
\fILWORK\fP 
.PP
.nf
!>          LWORK is INTEGER
!>          The dimension of the array WORK\&.
!>          LWORK >= 1, if MIN(M,N) = 0, and
!>          LWORK >= N+NRHS-1, otherwise\&.
!>          For optimal performance LWORK >= NB*( N+NRHS+1 ),
!>          where NB is the optimal block size for ZGEQP3RK returned
!>          by ILAENV\&. Minimal block size MINNB=2\&.
!>
!>          NOTE: The decision, whether to use unblocked BLAS 2
!>          or blocked BLAS 3 code is based not only on the dimension
!>          LWORK of the availbale workspace WORK, but also also on the
!>          matrix A dimension N via crossover point NX returned
!>          by ILAENV\&. (For N less than NX, unblocked code should be
!>          used\&.)
!>
!>          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\&.
!> 
.fi
.PP
.br
\fIRWORK\fP 
.PP
.nf
!>          RWORK is DOUBLE PRECISION array, dimension (2*N)
!> 
.fi
.PP
.br
\fIIWORK\fP 
.PP
.nf
!>          IWORK is INTEGER array, dimension (N-1)\&.
!>          Is a work array\&. ( IWORK is used to store indices
!>          of  columns for norm downdating in the residual
!>          matrix in the blocked step auxiliary subroutine ZLAQP3RK )\&.
!> 
.fi
.PP
.br
\fIINFO\fP 
.PP
.nf
!>          INFO is INTEGER
!>          1) INFO = 0: successful exit\&.
!>          2) INFO < 0: if INFO = -i, the i-th argument had an
!>                       illegal value\&.
!>          3) If INFO = j_1, where 1 <= j_1 <= N, then NaN was
!>             detected and the routine stops the computation\&.
!>             The j_1-th column of the matrix A or the j_1-th
!>             element of array TAU contains the first occurrence
!>             of NaN in the factorization step K+1 ( when K columns
!>             have been factorized )\&.
!>
!>             On exit:
!>             K                  is set to the number of
!>                                   factorized columns without
!>                                   exception\&.
!>             MAXC2NRMK          is set to NaN\&.
!>             RELMAXC2NRMK       is set to NaN\&.
!>             TAU(K+1:min(M,N))  is not set and contains undefined
!>                                   elements\&. If j_1=K+1, TAU(K+1)
!>                                   may contain NaN\&.
!>          4) If INFO = j_2, where N+1 <= j_2 <= 2*N, then no NaN
!>             was detected, but +Inf (or -Inf) was detected and
!>             the routine continues the computation until completion\&.
!>             The (j_2-N)-th column of the matrix A contains the first
!>             occurrence of +Inf (or -Inf) in the factorization
!>             step K+1 ( when K columns have been factorized )\&.
!> 
.fi
.PP
 
.RE
.PP
\fBAuthor\fP
.RS 4
Univ\&. of Tennessee 

.PP
Univ\&. of California Berkeley 

.PP
Univ\&. of Colorado Denver 

.PP
NAG Ltd\&. 
.RE
.PP
\fBFurther Details:\fP
.RS 4

.PP
.nf
!> ZGEQP3RK is based on the same BLAS3 Householder QR factorization
!> algorithm with column pivoting as in ZGEQP3 routine which uses
!> ZLARFG routine to generate Householder reflectors
!> for QR factorization\&.
!>
!> We can also write:
!>
!>   A = A_approx(K) + A_residual(K)
!>
!> The low rank approximation matrix A(K)_approx from
!> the truncated QR factorization of rank K of the matrix A is:
!>
!>   A(K)_approx = Q(K) * ( R(K)_approx ) * P(K)**T
!>                        (     0     0 )
!>
!>               = Q(K) * ( R11(K) R12(K) ) * P(K)**T
!>                        (      0      0 )
!>
!> The residual A_residual(K) of the matrix A is:
!>
!>   A_residual(K) = Q(K) * ( 0              0 ) * P(K)**T =
!>                          ( 0  R(K)_residual )
!>
!>                 = Q(K) * ( 0        0 ) * P(K)**T
!>                          ( 0   R22(K) )
!>
!> The truncated (rank K) factorization guarantees that
!> the maximum column 2-norm of A_residual(K) is less than
!> or equal to MAXC2NRMK up to roundoff error\&.
!>
!> NOTE: An approximation of the null vectors
!>       of A can be easily computed from R11(K)
!>       and R12(K):
!>
!>       Null( A(K) )_approx = P * ( inv(R11(K)) * R12(K) )
!>                                 (         -I           )
!>
!> 
.fi
.PP
 
.RE
.PP
\fBReferences:\fP
.RS 4
[1] A Level 3 BLAS QR factorization algorithm with column pivoting developed in 1996\&. G\&. Quintana-Orti, Depto\&. de Informatica, Universidad Jaime I, Spain\&. X\&. Sun, Computer Science Dept\&., Duke University, USA\&. C\&. H\&. Bischof, Math\&. and Comp\&. Sci\&. Div\&., Argonne National Lab, USA\&. A BLAS-3 version of the QR factorization with column pivoting\&. LAPACK Working Note 114  and in SIAM J\&. Sci\&. Comput\&., 19(5):1486-1494, Sept\&. 1998\&. 
.RE
.PP
[2] A partial column norm updating strategy developed in 2006\&. Z\&. Drmac and Z\&. Bujanovic, Dept\&. of Math\&., University of Zagreb, Croatia\&. On the failure of rank revealing QR factorization software – a case study\&. LAPACK Working Note 176\&.  and in ACM Trans\&. Math\&. Softw\&. 35, 2, Article 12 (July 2008), 28 pages\&.  
.PP
\fBContributors:\fP
.RS 4

.PP
.nf
!>
!>  November  2023, Igor Kozachenko, James Demmel,
!>                  EECS Department,
!>                  University of California, Berkeley, USA\&.
!>
!> 
.fi
.PP
 
.RE
.PP

.PP
Definition at line \fB589\fP of file \fBzgeqp3rk\&.f\fP\&.
.SH "Author"
.PP 
Generated automatically by Doxygen for LAPACK from the source code\&.