- dgebak
- dgebal
- dgebd2
- dgebrd
- dgecon
- dgeequ
- dgeequb
- dgees
- dgeesx
- dgeev
- dgeevx
- dgegs
- dgegv
- dgehd2
- dgehrd
- dgejsv
- dgelq2
- dgelqf
- dgels
- dgelsd
- dgelss
- dgelsx
- dgelsy
- dgeql2
- dgeqlf
- dgeqp3
- dgeqpf
- dgeqr2
- dgeqr2p
- dgeqrf
- dgeqrfp
- dgerfs
- dgerfsx
- dgerq2
- dgerqf
- dgesc2
- dgesdd
- dgesv
- dgesvd
- dgesvj
- dgesvx
- dgesvxx
- dgetc2
- dgetf2
- dgetrf
- dgetri
- dgetrs

USAGE: info, v = NumRu::Lapack.dgebak( job, side, ilo, ihi, scale, v, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEBAK( JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO ) * Purpose * ======= * * DGEBAK forms the right or left eigenvectors of a real general matrix * by backward transformation on the computed eigenvectors of the * balanced matrix output by DGEBAL. * * Arguments * ========= * * JOB (input) CHARACTER*1 * Specifies the type of backward transformation required: * = 'N', do nothing, return immediately; * = 'P', do backward transformation for permutation only; * = 'S', do backward transformation for scaling only; * = 'B', do backward transformations for both permutation and * scaling. * JOB must be the same as the argument JOB supplied to DGEBAL. * * SIDE (input) CHARACTER*1 * = 'R': V contains right eigenvectors; * = 'L': V contains left eigenvectors. * * N (input) INTEGER * The number of rows of the matrix V. N >= 0. * * ILO (input) INTEGER * IHI (input) INTEGER * The integers ILO and IHI determined by DGEBAL. * 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. * * SCALE (input) DOUBLE PRECISION array, dimension (N) * Details of the permutation and scaling factors, as returned * by DGEBAL. * * M (input) INTEGER * The number of columns of the matrix V. M >= 0. * * V (input/output) DOUBLE PRECISION array, dimension (LDV,M) * On entry, the matrix of right or left eigenvectors to be * transformed, as returned by DHSEIN or DTREVC. * On exit, V is overwritten by the transformed eigenvectors. * * LDV (input) INTEGER * The leading dimension of the array V. LDV >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value. * * ===================================================================== *go to the page top

USAGE: ilo, ihi, scale, info, a = NumRu::Lapack.dgebal( job, a, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO ) * Purpose * ======= * * DGEBAL balances a general real matrix A. This involves, first, * permuting A by a similarity transformation to isolate eigenvalues * in the first 1 to ILO-1 and last IHI+1 to N elements on the * diagonal; and second, applying a diagonal similarity transformation * to rows and columns ILO to IHI to make the rows and columns as * close in norm as possible. Both steps are optional. * * Balancing may reduce the 1-norm of the matrix, and improve the * accuracy of the computed eigenvalues and/or eigenvectors. * * Arguments * ========= * * JOB (input) CHARACTER*1 * Specifies the operations to be performed on A: * = 'N': none: simply set ILO = 1, IHI = N, SCALE(I) = 1.0 * for i = 1,...,N; * = 'P': permute only; * = 'S': scale only; * = 'B': both permute and scale. * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the input matrix A. * On exit, A is overwritten by the balanced matrix. * If JOB = 'N', A is not referenced. * See Further Details. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * ILO (output) INTEGER * IHI (output) INTEGER * ILO and IHI are set to integers such that on exit * A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N. * If JOB = 'N' or 'S', ILO = 1 and IHI = N. * * SCALE (output) DOUBLE PRECISION array, dimension (N) * Details of the permutations and scaling factors applied to * A. If P(j) is the index of the row and column interchanged * with row and column j and D(j) is the scaling factor * applied to row and column j, then * SCALE(j) = P(j) for j = 1,...,ILO-1 * = D(j) for j = ILO,...,IHI * = P(j) for j = IHI+1,...,N. * The order in which the interchanges are made is N to IHI+1, * then 1 to ILO-1. * * INFO (output) INTEGER * = 0: successful exit. * < 0: if INFO = -i, the i-th argument had an illegal value. * * Further Details * =============== * * The permutations consist of row and column interchanges which put * the matrix in the form * * ( T1 X Y ) * P A P = ( 0 B Z ) * ( 0 0 T2 ) * * where T1 and T2 are upper triangular matrices whose eigenvalues lie * along the diagonal. The column indices ILO and IHI mark the starting * and ending columns of the submatrix B. Balancing consists of applying * a diagonal similarity transformation inv(D) * B * D to make the * 1-norms of each row of B and its corresponding column nearly equal. * The output matrix is * * ( T1 X*D Y ) * ( 0 inv(D)*B*D inv(D)*Z ). * ( 0 0 T2 ) * * Information about the permutations P and the diagonal matrix D is * returned in the vector SCALE. * * This subroutine is based on the EISPACK routine BALANC. * * Modified by Tzu-Yi Chen, Computer Science Division, University of * California at Berkeley, USA * * ===================================================================== *go to the page top

USAGE: d, e, tauq, taup, info, a = NumRu::Lapack.dgebd2( m, a, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEBD2( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO ) * Purpose * ======= * * DGEBD2 reduces a real general m by n matrix A to upper or lower * bidiagonal form B by an orthogonal transformation: Q' * A * P = B. * * If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal. * * Arguments * ========= * * M (input) INTEGER * The number of rows in the matrix A. M >= 0. * * N (input) INTEGER * The number of columns in the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the m by n general matrix to be reduced. * On exit, * if m >= n, the diagonal and the first superdiagonal are * overwritten with the upper bidiagonal matrix B; the * elements below the diagonal, with the array TAUQ, represent * the orthogonal matrix Q as a product of elementary * reflectors, and the elements above the first superdiagonal, * with the array TAUP, represent the orthogonal matrix P as * a product of elementary reflectors; * if m < n, the diagonal and the first subdiagonal are * overwritten with the lower bidiagonal matrix B; the * elements below the first subdiagonal, with the array TAUQ, * represent the orthogonal matrix Q as a product of * elementary reflectors, and the elements above the diagonal, * with the array TAUP, represent the orthogonal matrix P as * a product of elementary reflectors. * See Further Details. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * D (output) DOUBLE PRECISION array, dimension (min(M,N)) * The diagonal elements of the bidiagonal matrix B: * D(i) = A(i,i). * * E (output) DOUBLE PRECISION array, dimension (min(M,N)-1) * The off-diagonal elements of the bidiagonal matrix B: * if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1; * if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1. * * TAUQ (output) DOUBLE PRECISION array dimension (min(M,N)) * The scalar factors of the elementary reflectors which * represent the orthogonal matrix Q. See Further Details. * * TAUP (output) DOUBLE PRECISION array, dimension (min(M,N)) * The scalar factors of the elementary reflectors which * represent the orthogonal matrix P. See Further Details. * * WORK (workspace) DOUBLE PRECISION array, dimension (max(M,N)) * * INFO (output) INTEGER * = 0: successful exit. * < 0: if INFO = -i, the i-th argument had an illegal value. * * Further Details * =============== * * The matrices Q and P are represented as products of elementary * reflectors: * * If m >= n, * * Q = H(1) H(2) . . . H(n) and P = G(1) G(2) . . . G(n-1) * * Each H(i) and G(i) has the form: * * H(i) = I - tauq * v * v' and G(i) = I - taup * u * u' * * where tauq and taup are real scalars, and v and u are real vectors; * v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in A(i+1:m,i); * u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in A(i,i+2:n); * tauq is stored in TAUQ(i) and taup in TAUP(i). * * If m < n, * * Q = H(1) H(2) . . . H(m-1) and P = G(1) G(2) . . . G(m) * * Each H(i) and G(i) has the form: * * H(i) = I - tauq * v * v' and G(i) = I - taup * u * u' * * where tauq and taup are real scalars, and v and u are real vectors; * v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in A(i+2:m,i); * u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n); * tauq is stored in TAUQ(i) and taup in TAUP(i). * * The contents of A on exit are illustrated by the following examples: * * m = 6 and n = 5 (m > n): m = 5 and n = 6 (m < n): * * ( d e u1 u1 u1 ) ( d u1 u1 u1 u1 u1 ) * ( v1 d e u2 u2 ) ( e d u2 u2 u2 u2 ) * ( v1 v2 d e u3 ) ( v1 e d u3 u3 u3 ) * ( v1 v2 v3 d e ) ( v1 v2 e d u4 u4 ) * ( v1 v2 v3 v4 d ) ( v1 v2 v3 e d u5 ) * ( v1 v2 v3 v4 v5 ) * * where d and e denote diagonal and off-diagonal elements of B, vi * denotes an element of the vector defining H(i), and ui an element of * the vector defining G(i). * * ===================================================================== *go to the page top

USAGE: d, e, tauq, taup, work, info, a = NumRu::Lapack.dgebrd( m, a, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO ) * Purpose * ======= * * DGEBRD reduces a general real M-by-N matrix A to upper or lower * bidiagonal form B by an orthogonal transformation: Q**T * A * P = B. * * If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal. * * Arguments * ========= * * M (input) INTEGER * The number of rows in the matrix A. M >= 0. * * N (input) INTEGER * The number of columns in the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the M-by-N general matrix to be reduced. * On exit, * if m >= n, the diagonal and the first superdiagonal are * overwritten with the upper bidiagonal matrix B; the * elements below the diagonal, with the array TAUQ, represent * the orthogonal matrix Q as a product of elementary * reflectors, and the elements above the first superdiagonal, * with the array TAUP, represent the orthogonal matrix P as * a product of elementary reflectors; * if m < n, the diagonal and the first subdiagonal are * overwritten with the lower bidiagonal matrix B; the * elements below the first subdiagonal, with the array TAUQ, * represent the orthogonal matrix Q as a product of * elementary reflectors, and the elements above the diagonal, * with the array TAUP, represent the orthogonal matrix P as * a product of elementary reflectors. * See Further Details. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * D (output) DOUBLE PRECISION array, dimension (min(M,N)) * The diagonal elements of the bidiagonal matrix B: * D(i) = A(i,i). * * E (output) DOUBLE PRECISION array, dimension (min(M,N)-1) * The off-diagonal elements of the bidiagonal matrix B: * if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1; * if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1. * * TAUQ (output) DOUBLE PRECISION array dimension (min(M,N)) * The scalar factors of the elementary reflectors which * represent the orthogonal matrix Q. See Further Details. * * TAUP (output) DOUBLE PRECISION array, dimension (min(M,N)) * The scalar factors of the elementary reflectors which * represent the orthogonal matrix P. See Further Details. * * 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 length of the array WORK. LWORK >= max(1,M,N). * For optimum performance LWORK >= (M+N)*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 matrices Q and P are represented as products of elementary * reflectors: * * If m >= n, * * Q = H(1) H(2) . . . H(n) and P = G(1) G(2) . . . G(n-1) * * Each H(i) and G(i) has the form: * * H(i) = I - tauq * v * v' and G(i) = I - taup * u * u' * * where tauq and taup are real scalars, and v and u are real vectors; * v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in A(i+1:m,i); * u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in A(i,i+2:n); * tauq is stored in TAUQ(i) and taup in TAUP(i). * * If m < n, * * Q = H(1) H(2) . . . H(m-1) and P = G(1) G(2) . . . G(m) * * Each H(i) and G(i) has the form: * * H(i) = I - tauq * v * v' and G(i) = I - taup * u * u' * * where tauq and taup are real scalars, and v and u are real vectors; * v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in A(i+2:m,i); * u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n); * tauq is stored in TAUQ(i) and taup in TAUP(i). * * The contents of A on exit are illustrated by the following examples: * * m = 6 and n = 5 (m > n): m = 5 and n = 6 (m < n): * * ( d e u1 u1 u1 ) ( d u1 u1 u1 u1 u1 ) * ( v1 d e u2 u2 ) ( e d u2 u2 u2 u2 ) * ( v1 v2 d e u3 ) ( v1 e d u3 u3 u3 ) * ( v1 v2 v3 d e ) ( v1 v2 e d u4 u4 ) * ( v1 v2 v3 v4 d ) ( v1 v2 v3 e d u5 ) * ( v1 v2 v3 v4 v5 ) * * where d and e denote diagonal and off-diagonal elements of B, vi * denotes an element of the vector defining H(i), and ui an element of * the vector defining G(i). * * ===================================================================== *go to the page top

USAGE: rcond, info = NumRu::Lapack.dgecon( norm, a, anorm, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO ) * Purpose * ======= * * DGECON estimates the reciprocal of the condition number of a general * real matrix A, in either the 1-norm or the infinity-norm, using * the LU factorization computed by DGETRF. * * An estimate is obtained for norm(inv(A)), and the reciprocal of the * condition number is computed as * RCOND = 1 / ( norm(A) * norm(inv(A)) ). * * Arguments * ========= * * NORM (input) CHARACTER*1 * Specifies whether the 1-norm condition number or the * infinity-norm condition number is required: * = '1' or 'O': 1-norm; * = 'I': Infinity-norm. * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input) DOUBLE PRECISION array, dimension (LDA,N) * The factors L and U from the factorization A = P*L*U * as computed by DGETRF. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * ANORM (input) DOUBLE PRECISION * If NORM = '1' or 'O', the 1-norm of the original matrix A. * If NORM = 'I', the infinity-norm of the original matrix A. * * RCOND (output) DOUBLE PRECISION * The reciprocal of the condition number of the matrix A, * computed as RCOND = 1/(norm(A) * norm(inv(A))). * * WORK (workspace) DOUBLE PRECISION array, dimension (4*N) * * IWORK (workspace) INTEGER array, dimension (N) * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * * ===================================================================== *go to the page top

USAGE: r, c, rowcnd, colcnd, amax, info = NumRu::Lapack.dgeequ( a, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEEQU( M, N, A, LDA, R, C, ROWCND, COLCND, AMAX, INFO ) * Purpose * ======= * * DGEEQU computes row and column scalings intended to equilibrate an * M-by-N matrix A and reduce its condition number. R returns the row * scale factors and C the column scale factors, chosen to try to make * the largest element in each row and column of the matrix B with * elements B(i,j)=R(i)*A(i,j)*C(j) have absolute value 1. * * R(i) and C(j) are restricted to be between SMLNUM = smallest safe * number and BIGNUM = largest safe number. Use of these scaling * factors is not guaranteed to reduce the condition number of A but * works well in practice. * * 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) DOUBLE PRECISION array, dimension (LDA,N) * The M-by-N matrix whose equilibration factors are * to be computed. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * R (output) DOUBLE PRECISION array, dimension (M) * If INFO = 0 or INFO > M, R contains the row scale factors * for A. * * C (output) DOUBLE PRECISION array, dimension (N) * If INFO = 0, C contains the column scale factors for A. * * ROWCND (output) DOUBLE PRECISION * If INFO = 0 or INFO > M, ROWCND contains the ratio of the * smallest R(i) to the largest R(i). If ROWCND >= 0.1 and * AMAX is neither too large nor too small, it is not worth * scaling by R. * * COLCND (output) DOUBLE PRECISION * If INFO = 0, COLCND contains the ratio of the smallest * C(i) to the largest C(i). If COLCND >= 0.1, it is not * worth scaling by C. * * AMAX (output) DOUBLE PRECISION * Absolute value of largest matrix element. If AMAX is very * close to overflow or very close to underflow, the matrix * should be scaled. * * 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 * <= M: the i-th row of A is exactly zero * > M: the (i-M)-th column of A is exactly zero * * ===================================================================== *go to the page top

USAGE: r, c, rowcnd, colcnd, amax, info = NumRu::Lapack.dgeequb( a, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEEQUB( M, N, A, LDA, R, C, ROWCND, COLCND, AMAX, INFO ) * Purpose * ======= * * DGEEQUB computes row and column scalings intended to equilibrate an * M-by-N matrix A and reduce its condition number. R returns the row * scale factors and C the column scale factors, chosen to try to make * the largest element in each row and column of the matrix B with * elements B(i,j)=R(i)*A(i,j)*C(j) have an absolute value of at most * the radix. * * R(i) and C(j) are restricted to be a power of the radix between * SMLNUM = smallest safe number and BIGNUM = largest safe number. Use * of these scaling factors is not guaranteed to reduce the condition * number of A but works well in practice. * * This routine differs from DGEEQU by restricting the scaling factors * to a power of the radix. Baring over- and underflow, scaling by * these factors introduces no additional rounding errors. However, the * scaled entries' magnitured are no longer approximately 1 but lie * between sqrt(radix) and 1/sqrt(radix). * * 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) DOUBLE PRECISION array, dimension (LDA,N) * The M-by-N matrix whose equilibration factors are * to be computed. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * R (output) DOUBLE PRECISION array, dimension (M) * If INFO = 0 or INFO > M, R contains the row scale factors * for A. * * C (output) DOUBLE PRECISION array, dimension (N) * If INFO = 0, C contains the column scale factors for A. * * ROWCND (output) DOUBLE PRECISION * If INFO = 0 or INFO > M, ROWCND contains the ratio of the * smallest R(i) to the largest R(i). If ROWCND >= 0.1 and * AMAX is neither too large nor too small, it is not worth * scaling by R. * * COLCND (output) DOUBLE PRECISION * If INFO = 0, COLCND contains the ratio of the smallest * C(i) to the largest C(i). If COLCND >= 0.1, it is not * worth scaling by C. * * AMAX (output) DOUBLE PRECISION * Absolute value of largest matrix element. If AMAX is very * close to overflow or very close to underflow, the matrix * should be scaled. * * 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 * <= M: the i-th row of A is exactly zero * > M: the (i-M)-th column of A is exactly zero * * ===================================================================== *go to the page top

USAGE: sdim, wr, wi, vs, work, info, a = NumRu::Lapack.dgees( jobvs, sort, a, [:lwork => lwork, :usage => usage, :help => help]){|a,b| ... } FORTRAN MANUAL SUBROUTINE DGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, WR, WI, VS, LDVS, WORK, LWORK, BWORK, INFO ) * Purpose * ======= * * DGEES computes for an N-by-N real nonsymmetric matrix A, the * eigenvalues, the real Schur form T, and, optionally, the matrix of * Schur vectors Z. This gives the Schur factorization A = Z*T*(Z**T). * * Optionally, it also orders the eigenvalues on the diagonal of the * real 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 matrix is in real Schur form if it is upper quasi-triangular with * 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in the * form * [ a b ] * [ c a ] * * where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc). * * 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 two DOUBLE PRECISION arguments * SELECT must be declared EXTERNAL in the calling subroutine. * If SORT = 'S', SELECT is used to select eigenvalues to sort * to the top left of the Schur form. * If SORT = 'N', SELECT is not referenced. * An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if * SELECT(WR(j),WI(j)) is true; i.e., if either one of a complex * conjugate pair of eigenvalues is selected, then both complex * eigenvalues are selected. * Note that a selected complex eigenvalue may no longer * satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since * ordering may change the value of complex eigenvalues * (especially if the eigenvalue is ill-conditioned); in this * case INFO is set to N+2 (see INFO below). * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the N-by-N matrix A. * On exit, A has been overwritten by its real 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 (after sorting) * for which SELECT is true. (Complex conjugate * pairs for which SELECT is true for either * eigenvalue count as 2.) * * WR (output) DOUBLE PRECISION array, dimension (N) * WI (output) DOUBLE PRECISION array, dimension (N) * WR and WI contain the real and imaginary parts, * respectively, of the computed eigenvalues in the same order * that they appear on the diagonal of the output Schur form T. * Complex conjugate pairs of eigenvalues will appear * consecutively with the eigenvalue having the positive * imaginary part first. * * VS (output) DOUBLE PRECISION array, dimension (LDVS,N) * If JOBVS = 'V', VS contains the orthogonal 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) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) * On exit, if INFO = 0, WORK(1) contains the optimal LWORK. * * LWORK (input) INTEGER * The dimension of the array WORK. LWORK >= max(1,3*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. * * 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 WR and WI * 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. * * ===================================================================== *go to the page top

USAGE: sdim, wr, wi, vs, rconde, rcondv, work, iwork, info, a = NumRu::Lapack.dgeesx( jobvs, sort, sense, a, liwork, [:lwork => lwork, :usage => usage, :help => help]){|a,b| ... } FORTRAN MANUAL SUBROUTINE DGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, WR, WI, VS, LDVS, RCONDE, RCONDV, WORK, LWORK, IWORK, LIWORK, BWORK, INFO ) * Purpose * ======= * * DGEESX computes for an N-by-N real nonsymmetric matrix A, the * eigenvalues, the real Schur form T, and, optionally, the matrix of * Schur vectors Z. This gives the Schur factorization A = Z*T*(Z**T). * * Optionally, it also orders the eigenvalues on the diagonal of the * real Schur form so that selected eigenvalues are at the top left; * computes a reciprocal condition number for the average of the * selected eigenvalues (RCONDE); and computes a reciprocal condition * number for the right invariant subspace corresponding to the * selected eigenvalues (RCONDV). The leading columns of Z form an * orthonormal basis for this invariant subspace. * * For further explanation of the reciprocal condition numbers RCONDE * and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where * these quantities are called s and sep respectively). * * A real matrix is in real Schur form if it is upper quasi-triangular * with 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in * the form * [ a b ] * [ c a ] * * where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc). * * 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 two DOUBLE PRECISION arguments * SELECT must be declared EXTERNAL in the calling subroutine. * If SORT = 'S', SELECT is used to select eigenvalues to sort * to the top left of the Schur form. * If SORT = 'N', SELECT is not referenced. * An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if * SELECT(WR(j),WI(j)) is true; i.e., if either one of a * complex conjugate pair of eigenvalues is selected, then both * are. Note that a selected complex eigenvalue may no longer * satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since * ordering may change the value of complex eigenvalues * (especially if the eigenvalue is ill-conditioned); in this * case INFO may be set to N+3 (see INFO below). * * SENSE (input) CHARACTER*1 * Determines which reciprocal condition numbers are computed. * = 'N': None are computed; * = 'E': Computed for average of selected eigenvalues only; * = 'V': Computed for selected right invariant subspace only; * = 'B': Computed for both. * If SENSE = 'E', 'V' or 'B', SORT must equal 'S'. * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA, N) * On entry, the N-by-N matrix A. * On exit, A is overwritten by its real 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 (after sorting) * for which SELECT is true. (Complex conjugate * pairs for which SELECT is true for either * eigenvalue count as 2.) * * WR (output) DOUBLE PRECISION array, dimension (N) * WI (output) DOUBLE PRECISION array, dimension (N) * WR and WI contain the real and imaginary parts, respectively, * of the computed eigenvalues, in the same order that they * appear on the diagonal of the output Schur form T. Complex * conjugate pairs of eigenvalues appear consecutively with the * eigenvalue having the positive imaginary part first. * * VS (output) DOUBLE PRECISION array, dimension (LDVS,N) * If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur * vectors. * If JOBVS = 'N', VS is not referenced. * * LDVS (input) INTEGER * The leading dimension of the array VS. LDVS >= 1, and if * JOBVS = 'V', LDVS >= N. * * RCONDE (output) DOUBLE PRECISION * If SENSE = 'E' or 'B', RCONDE contains the reciprocal * condition number for the average of the selected eigenvalues. * Not referenced if SENSE = 'N' or 'V'. * * RCONDV (output) DOUBLE PRECISION * If SENSE = 'V' or 'B', RCONDV contains the reciprocal * condition number for the selected right invariant subspace. * Not referenced if SENSE = 'N' or 'E'. * * 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 >= max(1,3*N). * Also, if SENSE = 'E' or 'V' or 'B', * LWORK >= N+2*SDIM*(N-SDIM), where SDIM is the number of * selected eigenvalues computed by this routine. Note that * N+2*SDIM*(N-SDIM) <= N+N*N/2. Note also that an error is only * returned if LWORK < max(1,3*N), but if SENSE = 'E' or 'V' or * 'B' this may not be large enough. * For good performance, LWORK must generally be larger. * * If LWORK = -1, then a workspace query is assumed; the routine * only calculates upper bounds on the optimal sizes of the * arrays WORK and IWORK, returns these values as the first * entries of the WORK and IWORK arrays, and no error messages * related to LWORK or LIWORK are issued by XERBLA. * * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. * * LIWORK (input) INTEGER * The dimension of the array IWORK. * LIWORK >= 1; if SENSE = 'V' or 'B', LIWORK >= SDIM*(N-SDIM). * Note that SDIM*(N-SDIM) <= N*N/4. Note also that an error is * only returned if LIWORK < 1, but if SENSE = 'V' or 'B' this * may not be large enough. * * If LIWORK = -1, then a workspace query is assumed; the * routine only calculates upper bounds on the optimal sizes of * the arrays WORK and IWORK, returns these values as the first * entries of the WORK and IWORK arrays, and no error messages * related to LWORK or LIWORK are issued by XERBLA. * * 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 WR and WI * contain those eigenvalues which have converged; if * JOBVS = 'V', VS contains the transformation 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. * * ===================================================================== *go to the page top

USAGE: wr, wi, vl, vr, work, info, a = NumRu::Lapack.dgeev( jobvl, jobvr, a, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR, LDVR, WORK, LWORK, INFO ) * Purpose * ======= * * DGEEV computes for an N-by-N real nonsymmetric matrix A, the * eigenvalues and, optionally, the left and/or right eigenvectors. * * The right eigenvector v(j) of A satisfies * A * v(j) = lambda(j) * v(j) * where lambda(j) is its eigenvalue. * The left eigenvector u(j) of A satisfies * u(j)**H * A = lambda(j) * u(j)**H * where u(j)**H denotes the conjugate transpose of u(j). * * The computed eigenvectors are normalized to have Euclidean norm * equal to 1 and largest component real. * * Arguments * ========= * * JOBVL (input) CHARACTER*1 * = 'N': left eigenvectors of A are not computed; * = 'V': left eigenvectors of A are computed. * * JOBVR (input) CHARACTER*1 * = 'N': right eigenvectors of A are not computed; * = 'V': right eigenvectors of A are computed. * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the N-by-N matrix A. * On exit, A has been overwritten. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * WR (output) DOUBLE PRECISION array, dimension (N) * WI (output) DOUBLE PRECISION array, dimension (N) * WR and WI contain the real and imaginary parts, * respectively, of the computed eigenvalues. Complex * conjugate pairs of eigenvalues appear consecutively * with the eigenvalue having the positive imaginary part * first. * * VL (output) DOUBLE PRECISION array, dimension (LDVL,N) * If JOBVL = 'V', the left eigenvectors u(j) are stored one * after another in the columns of VL, in the same order * as their eigenvalues. * If JOBVL = 'N', VL is not referenced. * If the j-th eigenvalue is real, then u(j) = VL(:,j), * the j-th column of VL. * If the j-th and (j+1)-st eigenvalues form a complex * conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and * u(j+1) = VL(:,j) - i*VL(:,j+1). * * LDVL (input) INTEGER * The leading dimension of the array VL. LDVL >= 1; if * JOBVL = 'V', LDVL >= N. * * VR (output) DOUBLE PRECISION array, dimension (LDVR,N) * If JOBVR = 'V', the right eigenvectors v(j) are stored one * after another in the columns of VR, in the same order * as their eigenvalues. * If JOBVR = 'N', VR is not referenced. * If the j-th eigenvalue is real, then v(j) = VR(:,j), * the j-th column of VR. * If the j-th and (j+1)-st eigenvalues form a complex * conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and * v(j+1) = VR(:,j) - i*VR(:,j+1). * * LDVR (input) INTEGER * The leading dimension of the array VR. LDVR >= 1; if * JOBVR = 'V', LDVR >= N. * * 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 >= max(1,3*N), and * if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*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. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value. * > 0: if INFO = i, the QR algorithm failed to compute all the * eigenvalues, and no eigenvectors have been computed; * elements i+1:N of WR and WI contain eigenvalues which * have converged. * * ===================================================================== *go to the page top

USAGE: wr, wi, vl, vr, ilo, ihi, scale, abnrm, rconde, rcondv, work, info, a = NumRu::Lapack.dgeevx( balanc, jobvl, jobvr, sense, a, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI, VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE, RCONDV, WORK, LWORK, IWORK, INFO ) * Purpose * ======= * * DGEEVX computes for an N-by-N real nonsymmetric matrix A, the * eigenvalues and, optionally, the left and/or right eigenvectors. * * Optionally also, it computes a balancing transformation to improve * the conditioning of the eigenvalues and eigenvectors (ILO, IHI, * SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues * (RCONDE), and reciprocal condition numbers for the right * eigenvectors (RCONDV). * * The right eigenvector v(j) of A satisfies * A * v(j) = lambda(j) * v(j) * where lambda(j) is its eigenvalue. * The left eigenvector u(j) of A satisfies * u(j)**H * A = lambda(j) * u(j)**H * where u(j)**H denotes the conjugate transpose of u(j). * * The computed eigenvectors are normalized to have Euclidean norm * equal to 1 and largest component real. * * Balancing a matrix means permuting the rows and columns to make it * more nearly upper triangular, and applying a diagonal similarity * transformation D * A * D**(-1), where D is a diagonal matrix, to * make its rows and columns closer in norm and the condition numbers * of its eigenvalues and eigenvectors smaller. The computed * reciprocal condition numbers correspond to the balanced matrix. * Permuting rows and columns will not change the condition numbers * (in exact arithmetic) but diagonal scaling will. For further * explanation of balancing, see section 4.10.2 of the LAPACK * Users' Guide. * * Arguments * ========= * * BALANC (input) CHARACTER*1 * Indicates how the input matrix should be diagonally scaled * and/or permuted to improve the conditioning of its * eigenvalues. * = 'N': Do not diagonally scale or permute; * = 'P': Perform permutations to make the matrix more nearly * upper triangular. Do not diagonally scale; * = 'S': Diagonally scale the matrix, i.e. replace A by * D*A*D**(-1), where D is a diagonal matrix chosen * to make the rows and columns of A more equal in * norm. Do not permute; * = 'B': Both diagonally scale and permute A. * * Computed reciprocal condition numbers will be for the matrix * after balancing and/or permuting. Permuting does not change * condition numbers (in exact arithmetic), but balancing does. * * JOBVL (input) CHARACTER*1 * = 'N': left eigenvectors of A are not computed; * = 'V': left eigenvectors of A are computed. * If SENSE = 'E' or 'B', JOBVL must = 'V'. * * JOBVR (input) CHARACTER*1 * = 'N': right eigenvectors of A are not computed; * = 'V': right eigenvectors of A are computed. * If SENSE = 'E' or 'B', JOBVR must = 'V'. * * SENSE (input) CHARACTER*1 * Determines which reciprocal condition numbers are computed. * = 'N': None are computed; * = 'E': Computed for eigenvalues only; * = 'V': Computed for right eigenvectors only; * = 'B': Computed for eigenvalues and right eigenvectors. * * If SENSE = 'E' or 'B', both left and right eigenvectors * must also be computed (JOBVL = 'V' and JOBVR = 'V'). * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the N-by-N matrix A. * On exit, A has been overwritten. If JOBVL = 'V' or * JOBVR = 'V', A contains the real Schur form of the balanced * version of the input matrix A. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * WR (output) DOUBLE PRECISION array, dimension (N) * WI (output) DOUBLE PRECISION array, dimension (N) * WR and WI contain the real and imaginary parts, * respectively, of the computed eigenvalues. Complex * conjugate pairs of eigenvalues will appear consecutively * with the eigenvalue having the positive imaginary part * first. * * VL (output) DOUBLE PRECISION array, dimension (LDVL,N) * If JOBVL = 'V', the left eigenvectors u(j) are stored one * after another in the columns of VL, in the same order * as their eigenvalues. * If JOBVL = 'N', VL is not referenced. * If the j-th eigenvalue is real, then u(j) = VL(:,j), * the j-th column of VL. * If the j-th and (j+1)-st eigenvalues form a complex * conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and * u(j+1) = VL(:,j) - i*VL(:,j+1). * * LDVL (input) INTEGER * The leading dimension of the array VL. LDVL >= 1; if * JOBVL = 'V', LDVL >= N. * * VR (output) DOUBLE PRECISION array, dimension (LDVR,N) * If JOBVR = 'V', the right eigenvectors v(j) are stored one * after another in the columns of VR, in the same order * as their eigenvalues. * If JOBVR = 'N', VR is not referenced. * If the j-th eigenvalue is real, then v(j) = VR(:,j), * the j-th column of VR. * If the j-th and (j+1)-st eigenvalues form a complex * conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and * v(j+1) = VR(:,j) - i*VR(:,j+1). * * LDVR (input) INTEGER * The leading dimension of the array VR. LDVR >= 1, and if * JOBVR = 'V', LDVR >= N. * * ILO (output) INTEGER * IHI (output) INTEGER * ILO and IHI are integer values determined when A was * balanced. The balanced A(i,j) = 0 if I > J and * J = 1,...,ILO-1 or I = IHI+1,...,N. * * SCALE (output) DOUBLE PRECISION array, dimension (N) * Details of the permutations and scaling factors applied * when balancing A. If P(j) is the index of the row and column * interchanged with row and column j, and D(j) is the scaling * factor applied to row and column j, then * SCALE(J) = P(J), for J = 1,...,ILO-1 * = D(J), for J = ILO,...,IHI * = P(J) for J = IHI+1,...,N. * The order in which the interchanges are made is N to IHI+1, * then 1 to ILO-1. * * ABNRM (output) DOUBLE PRECISION * The one-norm of the balanced matrix (the maximum * of the sum of absolute values of elements of any column). * * RCONDE (output) DOUBLE PRECISION array, dimension (N) * RCONDE(j) is the reciprocal condition number of the j-th * eigenvalue. * * RCONDV (output) DOUBLE PRECISION array, dimension (N) * RCONDV(j) is the reciprocal condition number of the j-th * right eigenvector. * * 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. If SENSE = 'N' or 'E', * LWORK >= max(1,2*N), and if JOBVL = 'V' or JOBVR = 'V', * LWORK >= 3*N. If SENSE = 'V' or 'B', LWORK >= N*(N+6). * 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. * * IWORK (workspace) INTEGER array, dimension (2*N-2) * If SENSE = 'N' or 'E', not referenced. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value. * > 0: if INFO = i, the QR algorithm failed to compute all the * eigenvalues, and no eigenvectors or condition numbers * have been computed; elements 1:ILO-1 and i+1:N of WR * and WI contain eigenvalues which have converged. * * ===================================================================== *go to the page top

USAGE: alphar, alphai, beta, vsl, vsr, work, info, a, b = NumRu::Lapack.dgegs( jobvsl, jobvsr, a, b, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK, LWORK, INFO ) * Purpose * ======= * * This routine is deprecated and has been replaced by routine DGGES. * * DGEGS computes the eigenvalues, real Schur form, and, optionally, * left and or/right Schur vectors of a real matrix pair (A,B). * Given two square matrices A and B, the generalized real Schur * factorization has the form * * A = Q*S*Z**T, B = Q*T*Z**T * * where Q and Z are orthogonal matrices, T is upper triangular, and S * is an upper quasi-triangular matrix with 1-by-1 and 2-by-2 diagonal * blocks, the 2-by-2 blocks corresponding to complex conjugate pairs * of eigenvalues of (A,B). The columns of Q are the left Schur vectors * and the columns of Z are the right Schur vectors. * * If only the eigenvalues of (A,B) are needed, the driver routine * DGEGV should be used instead. See DGEGV for a description of the * eigenvalues of the generalized nonsymmetric eigenvalue problem * (GNEP). * * Arguments * ========= * * JOBVSL (input) CHARACTER*1 * = 'N': do not compute the left Schur vectors; * = 'V': compute the left Schur vectors (returned in VSL). * * JOBVSR (input) CHARACTER*1 * = 'N': do not compute the right Schur vectors; * = 'V': compute the right Schur vectors (returned in VSR). * * N (input) INTEGER * The order of the matrices A, B, VSL, and VSR. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA, N) * On entry, the matrix A. * On exit, the upper quasi-triangular matrix S from the * generalized real Schur factorization. * * LDA (input) INTEGER * The leading dimension of A. LDA >= max(1,N). * * B (input/output) DOUBLE PRECISION array, dimension (LDB, N) * On entry, the matrix B. * On exit, the upper triangular matrix T from the generalized * real Schur factorization. * * LDB (input) INTEGER * The leading dimension of B. LDB >= max(1,N). * * ALPHAR (output) DOUBLE PRECISION array, dimension (N) * The real parts of each scalar alpha defining an eigenvalue * of GNEP. * * ALPHAI (output) DOUBLE PRECISION array, dimension (N) * The imaginary parts of each scalar alpha defining an * eigenvalue of GNEP. If ALPHAI(j) is zero, then the j-th * eigenvalue is real; if positive, then the j-th and (j+1)-st * eigenvalues are a complex conjugate pair, with * ALPHAI(j+1) = -ALPHAI(j). * * BETA (output) DOUBLE PRECISION array, dimension (N) * The scalars beta that define the eigenvalues of GNEP. * Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and * beta = BETA(j) represent the j-th eigenvalue of the matrix * pair (A,B), in one of the forms lambda = alpha/beta or * mu = beta/alpha. Since either lambda or mu may overflow, * they should not, in general, be computed. * * VSL (output) DOUBLE PRECISION array, dimension (LDVSL,N) * If JOBVSL = 'V', the matrix of left Schur vectors Q. * Not referenced if JOBVSL = 'N'. * * LDVSL (input) INTEGER * The leading dimension of the matrix VSL. LDVSL >=1, and * if JOBVSL = 'V', LDVSL >= N. * * VSR (output) DOUBLE PRECISION array, dimension (LDVSR,N) * If JOBVSR = 'V', the matrix of right Schur vectors Z. * Not referenced if JOBVSR = 'N'. * * LDVSR (input) INTEGER * The leading dimension of the matrix VSR. LDVSR >= 1, and * if JOBVSR = 'V', LDVSR >= N. * * 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 >= max(1,4*N). * For good performance, LWORK must generally be larger. * To compute the optimal value of LWORK, call ILAENV to get * blocksizes (for DGEQRF, DORMQR, and DORGQR.) Then compute: * NB -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR * The optimal LWORK is 2*N + N*(NB+1). * * 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. * = 1,...,N: * The QZ iteration failed. (A,B) are not in Schur * form, but ALPHAR(j), ALPHAI(j), and BETA(j) should * be correct for j=INFO+1,...,N. * > N: errors that usually indicate LAPACK problems: * =N+1: error return from DGGBAL * =N+2: error return from DGEQRF * =N+3: error return from DORMQR * =N+4: error return from DORGQR * =N+5: error return from DGGHRD * =N+6: error return from DHGEQZ (other than failed * iteration) * =N+7: error return from DGGBAK (computing VSL) * =N+8: error return from DGGBAK (computing VSR) * =N+9: error return from DLASCL (various places) * * ===================================================================== *go to the page top

USAGE: alphar, alphai, beta, vl, vr, work, info, a, b = NumRu::Lapack.dgegv( jobvl, jobvr, a, b, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO ) * Purpose * ======= * * This routine is deprecated and has been replaced by routine DGGEV. * * DGEGV computes the eigenvalues and, optionally, the left and/or right * eigenvectors of a real matrix pair (A,B). * Given two square matrices A and B, * the generalized nonsymmetric eigenvalue problem (GNEP) is to find the * eigenvalues lambda and corresponding (non-zero) eigenvectors x such * that * * A*x = lambda*B*x. * * An alternate form is to find the eigenvalues mu and corresponding * eigenvectors y such that * * mu*A*y = B*y. * * These two forms are equivalent with mu = 1/lambda and x = y if * neither lambda nor mu is zero. In order to deal with the case that * lambda or mu is zero or small, two values alpha and beta are returned * for each eigenvalue, such that lambda = alpha/beta and * mu = beta/alpha. * * The vectors x and y in the above equations are right eigenvectors of * the matrix pair (A,B). Vectors u and v satisfying * * u**H*A = lambda*u**H*B or mu*v**H*A = v**H*B * * are left eigenvectors of (A,B). * * Note: this routine performs "full balancing" on A and B -- see * "Further Details", below. * * Arguments * ========= * * JOBVL (input) CHARACTER*1 * = 'N': do not compute the left generalized eigenvectors; * = 'V': compute the left generalized eigenvectors (returned * in VL). * * JOBVR (input) CHARACTER*1 * = 'N': do not compute the right generalized eigenvectors; * = 'V': compute the right generalized eigenvectors (returned * in VR). * * N (input) INTEGER * The order of the matrices A, B, VL, and VR. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA, N) * On entry, the matrix A. * If JOBVL = 'V' or JOBVR = 'V', then on exit A * contains the real Schur form of A from the generalized Schur * factorization of the pair (A,B) after balancing. * If no eigenvectors were computed, then only the diagonal * blocks from the Schur form will be correct. See DGGHRD and * DHGEQZ for details. * * LDA (input) INTEGER * The leading dimension of A. LDA >= max(1,N). * * B (input/output) DOUBLE PRECISION array, dimension (LDB, N) * On entry, the matrix B. * If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the * upper triangular matrix obtained from B in the generalized * Schur factorization of the pair (A,B) after balancing. * If no eigenvectors were computed, then only those elements of * B corresponding to the diagonal blocks from the Schur form of * A will be correct. See DGGHRD and DHGEQZ for details. * * LDB (input) INTEGER * The leading dimension of B. LDB >= max(1,N). * * ALPHAR (output) DOUBLE PRECISION array, dimension (N) * The real parts of each scalar alpha defining an eigenvalue of * GNEP. * * ALPHAI (output) DOUBLE PRECISION array, dimension (N) * The imaginary parts of each scalar alpha defining an * eigenvalue of GNEP. If ALPHAI(j) is zero, then the j-th * eigenvalue is real; if positive, then the j-th and * (j+1)-st eigenvalues are a complex conjugate pair, with * ALPHAI(j+1) = -ALPHAI(j). * * BETA (output) DOUBLE PRECISION array, dimension (N) * The scalars beta that define the eigenvalues of GNEP. * * Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and * beta = BETA(j) represent the j-th eigenvalue of the matrix * pair (A,B), in one of the forms lambda = alpha/beta or * mu = beta/alpha. Since either lambda or mu may overflow, * they should not, in general, be computed. * * VL (output) DOUBLE PRECISION array, dimension (LDVL,N) * If JOBVL = 'V', the left eigenvectors u(j) are stored * in the columns of VL, in the same order as their eigenvalues. * If the j-th eigenvalue is real, then u(j) = VL(:,j). * If the j-th and (j+1)-st eigenvalues form a complex conjugate * pair, then * u(j) = VL(:,j) + i*VL(:,j+1) * and * u(j+1) = VL(:,j) - i*VL(:,j+1). * * Each eigenvector is scaled so that its largest component has * abs(real part) + abs(imag. part) = 1, except for eigenvectors * corresponding to an eigenvalue with alpha = beta = 0, which * are set to zero. * Not referenced if JOBVL = 'N'. * * LDVL (input) INTEGER * The leading dimension of the matrix VL. LDVL >= 1, and * if JOBVL = 'V', LDVL >= N. * * VR (output) DOUBLE PRECISION array, dimension (LDVR,N) * If JOBVR = 'V', the right eigenvectors x(j) are stored * in the columns of VR, in the same order as their eigenvalues. * If the j-th eigenvalue is real, then x(j) = VR(:,j). * If the j-th and (j+1)-st eigenvalues form a complex conjugate * pair, then * x(j) = VR(:,j) + i*VR(:,j+1) * and * x(j+1) = VR(:,j) - i*VR(:,j+1). * * Each eigenvector is scaled so that its largest component has * abs(real part) + abs(imag. part) = 1, except for eigenvalues * corresponding to an eigenvalue with alpha = beta = 0, which * are set to zero. * Not referenced if JOBVR = 'N'. * * LDVR (input) INTEGER * The leading dimension of the matrix VR. LDVR >= 1, and * if JOBVR = 'V', LDVR >= N. * * 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 >= max(1,8*N). * For good performance, LWORK must generally be larger. * To compute the optimal value of LWORK, call ILAENV to get * blocksizes (for DGEQRF, DORMQR, and DORGQR.) Then compute: * NB -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR; * The optimal LWORK is: * 2*N + MAX( 6*N, N*(NB+1) ). * * 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. * = 1,...,N: * The QZ iteration failed. No eigenvectors have been * calculated, but ALPHAR(j), ALPHAI(j), and BETA(j) * should be correct for j=INFO+1,...,N. * > N: errors that usually indicate LAPACK problems: * =N+1: error return from DGGBAL * =N+2: error return from DGEQRF * =N+3: error return from DORMQR * =N+4: error return from DORGQR * =N+5: error return from DGGHRD * =N+6: error return from DHGEQZ (other than failed * iteration) * =N+7: error return from DTGEVC * =N+8: error return from DGGBAK (computing VL) * =N+9: error return from DGGBAK (computing VR) * =N+10: error return from DLASCL (various calls) * * Further Details * =============== * * Balancing * --------- * * This driver calls DGGBAL to both permute and scale rows and columns * of A and B. The permutations PL and PR are chosen so that PL*A*PR * and PL*B*R will be upper triangular except for the diagonal blocks * A(i:j,i:j) and B(i:j,i:j), with i and j as close together as * possible. The diagonal scaling matrices DL and DR are chosen so * that the pair DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to * one (except for the elements that start out zero.) * * After the eigenvalues and eigenvectors of the balanced matrices * have been computed, DGGBAK transforms the eigenvectors back to what * they would have been (in perfect arithmetic) if they had not been * balanced. * * Contents of A and B on Exit * -------- -- - --- - -- ---- * * If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or * both), then on exit the arrays A and B will contain the real Schur * form[*] of the "balanced" versions of A and B. If no eigenvectors * are computed, then only the diagonal blocks will be correct. * * [*] See DHGEQZ, DGEGS, or read the book "Matrix Computations", * by Golub & van Loan, pub. by Johns Hopkins U. Press. * * ===================================================================== *go to the page top

USAGE: tau, info, a = NumRu::Lapack.dgehd2( ilo, ihi, a, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEHD2( N, ILO, IHI, A, LDA, TAU, WORK, INFO ) * Purpose * ======= * * DGEHD2 reduces a real general matrix A to upper Hessenberg form H by * an orthogonal similarity transformation: Q' * A * Q = H . * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. N >= 0. * * ILO (input) INTEGER * IHI (input) INTEGER * It is assumed that A is already upper triangular in rows * and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally * set by a previous call to DGEBAL; otherwise they should be * set to 1 and N respectively. See Further Details. * 1 <= ILO <= IHI <= max(1,N). * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the n by n general matrix to be reduced. * On exit, the upper triangle and the first subdiagonal of A * are overwritten with the upper Hessenberg matrix H, and the * elements below the first subdiagonal, with the array TAU, * represent the orthogonal matrix Q as a product of elementary * reflectors. See Further Details. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * TAU (output) DOUBLE PRECISION array, dimension (N-1) * The scalar factors of the elementary reflectors (see Further * Details). * * WORK (workspace) DOUBLE PRECISION array, dimension (N) * * 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 (ihi-ilo) elementary * reflectors * * Q = H(ilo) H(ilo+1) . . . H(ihi-1). * * Each H(i) has the form * * H(i) = I - tau * v * v' * * where tau is a real scalar, and v is a real vector with * v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on * exit in A(i+2:ihi,i), and tau in TAU(i). * * The contents of A are illustrated by the following example, with * n = 7, ilo = 2 and ihi = 6: * * on entry, on exit, * * ( a a a a a a a ) ( a a h h h h a ) * ( a a a a a a ) ( a h h h h a ) * ( a a a a a a ) ( h h h h h h ) * ( a a a a a a ) ( v2 h h h h h ) * ( a a a a a a ) ( v2 v3 h h h h ) * ( a a a a a a ) ( v2 v3 v4 h h h ) * ( a ) ( a ) * * where a denotes an element of the original matrix A, h denotes a * modified element of the upper Hessenberg matrix H, and vi denotes an * element of the vector defining H(i). * * ===================================================================== *go to the page top

USAGE: tau, work, info, a = NumRu::Lapack.dgehrd( ilo, ihi, a, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO ) * Purpose * ======= * * DGEHRD reduces a real general matrix A to upper Hessenberg form H by * an orthogonal similarity transformation: Q' * A * Q = H . * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. N >= 0. * * ILO (input) INTEGER * IHI (input) INTEGER * It is assumed that A is already upper triangular in rows * and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally * set by a previous call to DGEBAL; otherwise they should be * set to 1 and N respectively. See Further Details. * 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the N-by-N general matrix to be reduced. * On exit, the upper triangle and the first subdiagonal of A * are overwritten with the upper Hessenberg matrix H, and the * elements below the first subdiagonal, with the array TAU, * represent the orthogonal matrix Q as a product of elementary * reflectors. See Further Details. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * TAU (output) DOUBLE PRECISION array, dimension (N-1) * The scalar factors of the elementary reflectors (see Further * Details). Elements 1:ILO-1 and IHI:N-1 of TAU are set to * zero. * * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. * * LWORK (input) INTEGER * The length of the array WORK. LWORK >= max(1,N). * For optimum performance LWORK >= N*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 (ihi-ilo) elementary * reflectors * * Q = H(ilo) H(ilo+1) . . . H(ihi-1). * * Each H(i) has the form * * H(i) = I - tau * v * v' * * where tau is a real scalar, and v is a real vector with * v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on * exit in A(i+2:ihi,i), and tau in TAU(i). * * The contents of A are illustrated by the following example, with * n = 7, ilo = 2 and ihi = 6: * * on entry, on exit, * * ( a a a a a a a ) ( a a h h h h a ) * ( a a a a a a ) ( a h h h h a ) * ( a a a a a a ) ( h h h h h h ) * ( a a a a a a ) ( v2 h h h h h ) * ( a a a a a a ) ( v2 v3 h h h h ) * ( a a a a a a ) ( v2 v3 v4 h h h ) * ( a ) ( a ) * * where a denotes an element of the original matrix A, h denotes a * modified element of the upper Hessenberg matrix H, and vi denotes an * element of the vector defining H(i). * * This file is a slight modification of LAPACK-3.0's DGEHRD * subroutine incorporating improvements proposed by Quintana-Orti and * Van de Geijn (2006). (See DLAHR2.) * * ===================================================================== *go to the page top

USAGE: sva, u, v, iwork, info, work = NumRu::Lapack.dgejsv( joba, jobu, jobv, jobr, jobt, jobp, m, a, work, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP, M, N, A, LDA, SVA, U, LDU, V, LDV, WORK, LWORK, IWORK, INFO ) * Purpose * ======= * * DGEJSV computes the singular value decomposition (SVD) of a real M-by-N * matrix [A], where M >= N. The SVD of [A] is written as * * [A] = [U] * [SIGMA] * [V]^t, * * where [SIGMA] is an N-by-N (M-by-N) matrix which is zero except for its N * diagonal elements, [U] is an M-by-N (or M-by-M) orthonormal matrix, and * [V] is an N-by-N orthogonal matrix. The diagonal elements of [SIGMA] are * the singular values of [A]. The columns of [U] and [V] are the left and * the right singular vectors of [A], respectively. The matrices [U] and [V] * are computed and stored in the arrays U and V, respectively. The diagonal * of [SIGMA] is computed and stored in the array SVA. * * Arguments * ========= * * JOBA (input) CHARACTER*1 * Specifies the level of accuracy: * = 'C': This option works well (high relative accuracy) if A = B * D, * with well-conditioned B and arbitrary diagonal matrix D. * The accuracy cannot be spoiled by COLUMN scaling. The * accuracy of the computed output depends on the condition of * B, and the procedure aims at the best theoretical accuracy. * The relative error max_{i=1:N}|d sigma_i| / sigma_i is * bounded by f(M,N)*epsilon* cond(B), independent of D. * The input matrix is preprocessed with the QRF with column * pivoting. This initial preprocessing and preconditioning by * a rank revealing QR factorization is common for all values of * JOBA. Additional actions are specified as follows: * = 'E': Computation as with 'C' with an additional estimate of the * condition number of B. It provides a realistic error bound. * = 'F': If A = D1 * C * D2 with ill-conditioned diagonal scalings * D1, D2, and well-conditioned matrix C, this option gives * higher accuracy than the 'C' option. If the structure of the * input matrix is not known, and relative accuracy is * desirable, then this option is advisable. The input matrix A * is preprocessed with QR factorization with FULL (row and * column) pivoting. * = 'G' Computation as with 'F' with an additional estimate of the * condition number of B, where A=D*B. If A has heavily weighted * rows, then using this condition number gives too pessimistic * error bound. * = 'A': Small singular values are the noise and the matrix is treated * as numerically rank defficient. The error in the computed * singular values is bounded by f(m,n)*epsilon*||A||. * The computed SVD A = U * S * V^t restores A up to * f(m,n)*epsilon*||A||. * This gives the procedure the licence to discard (set to zero) * all singular values below N*epsilon*||A||. * = 'R': Similar as in 'A'. Rank revealing property of the initial * QR factorization is used do reveal (using triangular factor) * a gap sigma_{r+1} < epsilon * sigma_r in which case the * numerical RANK is declared to be r. The SVD is computed with * absolute error bounds, but more accurately than with 'A'. * * JOBU (input) CHARACTER*1 * Specifies whether to compute the columns of U: * = 'U': N columns of U are returned in the array U. * = 'F': full set of M left sing. vectors is returned in the array U. * = 'W': U may be used as workspace of length M*N. See the description * of U. * = 'N': U is not computed. * * JOBV (input) CHARACTER*1 * Specifies whether to compute the matrix V: * = 'V': N columns of V are returned in the array V; Jacobi rotations * are not explicitly accumulated. * = 'J': N columns of V are returned in the array V, but they are * computed as the product of Jacobi rotations. This option is * allowed only if JOBU .NE. 'N', i.e. in computing the full SVD. * = 'W': V may be used as workspace of length N*N. See the description * of V. * = 'N': V is not computed. * * JOBR (input) CHARACTER*1 * Specifies the RANGE for the singular values. Issues the licence to * set to zero small positive singular values if they are outside * specified range. If A .NE. 0 is scaled so that the largest singular * value of c*A is around DSQRT(BIG), BIG=SLAMCH('O'), then JOBR issues * the licence to kill columns of A whose norm in c*A is less than * DSQRT(SFMIN) (for JOBR.EQ.'R'), or less than SMALL=SFMIN/EPSLN, * where SFMIN=SLAMCH('S'), EPSLN=SLAMCH('E'). * = 'N': Do not kill small columns of c*A. This option assumes that * BLAS and QR factorizations and triangular solvers are * implemented to work in that range. If the condition of A * is greater than BIG, use DGESVJ. * = 'R': RESTRICTED range for sigma(c*A) is [DSQRT(SFMIN), DSQRT(BIG)] * (roughly, as described above). This option is recommended. * ~~~~~~~~~~~~~~~~~~~~~~~~~~~ * For computing the singular values in the FULL range [SFMIN,BIG] * use DGESVJ. * * JOBT (input) CHARACTER*1 * If the matrix is square then the procedure may determine to use * transposed A if A^t seems to be better with respect to convergence. * If the matrix is not square, JOBT is ignored. This is subject to * changes in the future. * The decision is based on two values of entropy over the adjoint * orbit of A^t * A. See the descriptions of WORK(6) and WORK(7). * = 'T': transpose if entropy test indicates possibly faster * convergence of Jacobi process if A^t is taken as input. If A is * replaced with A^t, then the row pivoting is included automatically. * = 'N': do not speculate. * This option can be used to compute only the singular values, or the * full SVD (U, SIGMA and V). For only one set of singular vectors * (U or V), the caller should provide both U and V, as one of the * matrices is used as workspace if the matrix A is transposed. * The implementer can easily remove this constraint and make the * code more complicated. See the descriptions of U and V. * * JOBP (input) CHARACTER*1 * Issues the licence to introduce structured perturbations to drown * denormalized numbers. This licence should be active if the * denormals are poorly implemented, causing slow computation, * especially in cases of fast convergence (!). For details see [1,2]. * For the sake of simplicity, this perturbations are included only * when the full SVD or only the singular values are requested. The * implementer/user can easily add the perturbation for the cases of * computing one set of singular vectors. * = 'P': introduce perturbation * = 'N': do not perturb * * M (input) INTEGER * The number of rows of the input matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the input matrix A. M >= N >= 0. * * A (input/workspace) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the M-by-N matrix A. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * SVA (workspace/output) DOUBLE PRECISION array, dimension (N) * On exit, * - For WORK(1)/WORK(2) = ONE: The singular values of A. During the * computation SVA contains Euclidean column norms of the * iterated matrices in the array A. * - For WORK(1) .NE. WORK(2): The singular values of A are * (WORK(1)/WORK(2)) * SVA(1:N). This factored form is used if * sigma_max(A) overflows or if small singular values have been * saved from underflow by scaling the input matrix A. * - If JOBR='R' then some of the singular values may be returned * as exact zeros obtained by "set to zero" because they are * below the numerical rank threshold or are denormalized numbers. * * U (workspace/output) DOUBLE PRECISION array, dimension ( LDU, N ) * If JOBU = 'U', then U contains on exit the M-by-N matrix of * the left singular vectors. * If JOBU = 'F', then U contains on exit the M-by-M matrix of * the left singular vectors, including an ONB * of the orthogonal complement of the Range(A). * If JOBU = 'W' .AND. (JOBV.EQ.'V' .AND. JOBT.EQ.'T' .AND. M.EQ.N), * then U is used as workspace if the procedure * replaces A with A^t. In that case, [V] is computed * in U as left singular vectors of A^t and then * copied back to the V array. This 'W' option is just * a reminder to the caller that in this case U is * reserved as workspace of length N*N. * If JOBU = 'N' U is not referenced. * * LDU (input) INTEGER * The leading dimension of the array U, LDU >= 1. * IF JOBU = 'U' or 'F' or 'W', then LDU >= M. * * V (workspace/output) DOUBLE PRECISION array, dimension ( LDV, N ) * If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of * the right singular vectors; * If JOBV = 'W', AND (JOBU.EQ.'U' AND JOBT.EQ.'T' AND M.EQ.N), * then V is used as workspace if the pprocedure * replaces A with A^t. In that case, [U] is computed * in V as right singular vectors of A^t and then * copied back to the U array. This 'W' option is just * a reminder to the caller that in this case V is * reserved as workspace of length N*N. * If JOBV = 'N' V is not referenced. * * LDV (input) INTEGER * The leading dimension of the array V, LDV >= 1. * If JOBV = 'V' or 'J' or 'W', then LDV >= N. * * WORK (workspace/output) DOUBLE PRECISION array, dimension at least LWORK. * On exit, * WORK(1) = SCALE = WORK(2) / WORK(1) is the scaling factor such * that SCALE*SVA(1:N) are the computed singular values * of A. (See the description of SVA().) * WORK(2) = See the description of WORK(1). * WORK(3) = SCONDA is an estimate for the condition number of * column equilibrated A. (If JOBA .EQ. 'E' or 'G') * SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1). * It is computed using DPOCON. It holds * N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA * where R is the triangular factor from the QRF of A. * However, if R is truncated and the numerical rank is * determined to be strictly smaller than N, SCONDA is * returned as -1, thus indicating that the smallest * singular values might be lost. * * If full SVD is needed, the following two condition numbers are * useful for the analysis of the algorithm. They are provied for * a developer/implementer who is familiar with the details of * the method. * * WORK(4) = an estimate of the scaled condition number of the * triangular factor in the first QR factorization. * WORK(5) = an estimate of the scaled condition number of the * triangular factor in the second QR factorization. * The following two parameters are computed if JOBT .EQ. 'T'. * They are provided for a developer/implementer who is familiar * with the details of the method. * * WORK(6) = the entropy of A^t*A :: this is the Shannon entropy * of diag(A^t*A) / Trace(A^t*A) taken as point in the * probability simplex. * WORK(7) = the entropy of A*A^t. * * LWORK (input) INTEGER * Length of WORK to confirm proper allocation of work space. * LWORK depends on the job: * * If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and * -> .. no scaled condition estimate required ( JOBE.EQ.'N'): * LWORK >= max(2*M+N,4*N+1,7). This is the minimal requirement. * For optimal performance (blocked code) the optimal value * is LWORK >= max(2*M+N,3*N+(N+1)*NB,7). Here NB is the optimal * block size for xGEQP3/xGEQRF. * -> .. an estimate of the scaled condition number of A is * required (JOBA='E', 'G'). In this case, LWORK is the maximum * of the above and N*N+4*N, i.e. LWORK >= max(2*M+N,N*N+4N,7). * * If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'), * -> the minimal requirement is LWORK >= max(2*N+M,7). * -> For optimal performance, LWORK >= max(2*N+M,2*N+N*NB,7), * where NB is the optimal block size. * * If SIGMA and the left singular vectors are needed * -> the minimal requirement is LWORK >= max(2*N+M,7). * -> For optimal performance, LWORK >= max(2*N+M,2*N+N*NB,7), * where NB is the optimal block size. * * If full SVD is needed ( JOBU.EQ.'U' or 'F', JOBV.EQ.'V' ) and * -> .. the singular vectors are computed without explicit * accumulation of the Jacobi rotations, LWORK >= 6*N+2*N*N * -> .. in the iterative part, the Jacobi rotations are * explicitly accumulated (option, see the description of JOBV), * then the minimal requirement is LWORK >= max(M+3*N+N*N,7). * For better performance, if NB is the optimal block size, * LWORK >= max(3*N+N*N+M,3*N+N*N+N*NB,7). * * IWORK (workspace/output) INTEGER array, dimension M+3*N. * On exit, * IWORK(1) = the numerical rank determined after the initial * QR factorization with pivoting. See the descriptions * of JOBA and JOBR. * IWORK(2) = the number of the computed nonzero singular values * IWORK(3) = if nonzero, a warning message: * If IWORK(3).EQ.1 then some of the column norms of A * were denormalized floats. The requested high accuracy * is not warranted by the data. * * INFO (output) INTEGER * < 0 : if INFO = -i, then the i-th argument had an illegal value. * = 0 : successfull exit; * > 0 : DGEJSV did not converge in the maximal allowed number * of sweeps. The computed values may be inaccurate. * * Further Details * =============== * * DGEJSV implements a preconditioned Jacobi SVD algorithm. It uses SGEQP3, * SGEQRF, and SGELQF as preprocessors and preconditioners. Optionally, an * additional row pivoting can be used as a preprocessor, which in some * cases results in much higher accuracy. An example is matrix A with the * structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned * diagonal matrices and C is well-conditioned matrix. In that case, complete * pivoting in the first QR factorizations provides accuracy dependent on the * condition number of C, and independent of D1, D2. Such higher accuracy is * not completely understood theoretically, but it works well in practice. * Further, if A can be written as A = B*D, with well-conditioned B and some * diagonal D, then the high accuracy is guaranteed, both theoretically and * in software, independent of D. For more details see [1], [2]. * The computational range for the singular values can be the full range * ( UNDERFLOW,OVERFLOW ), provided that the machine arithmetic and the BLAS * & LAPACK routines called by DGEJSV are implemented to work in that range. * If that is not the case, then the restriction for safe computation with * the singular values in the range of normalized IEEE numbers is that the * spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not * overflow. This code (DGEJSV) is best used in this restricted range, * meaning that singular values of magnitude below ||A||_2 / SLAMCH('O') are * returned as zeros. See JOBR for details on this. * Further, this implementation is somewhat slower than the one described * in [1,2] due to replacement of some non-LAPACK components, and because * the choice of some tuning parameters in the iterative part (DGESVJ) is * left to the implementer on a particular machine. * The rank revealing QR factorization (in this code: SGEQP3) should be * implemented as in [3]. We have a new version of SGEQP3 under development * that is more robust than the current one in LAPACK, with a cleaner cut in * rank defficient cases. It will be available in the SIGMA library [4]. * If M is much larger than N, it is obvious that the inital QRF with * column pivoting can be preprocessed by the QRF without pivoting. That * well known trick is not used in DGEJSV because in some cases heavy row * weighting can be treated with complete pivoting. The overhead in cases * M much larger than N is then only due to pivoting, but the benefits in * terms of accuracy have prevailed. The implementer/user can incorporate * this extra QRF step easily. The implementer can also improve data movement * (matrix transpose, matrix copy, matrix transposed copy) - this * implementation of DGEJSV uses only the simplest, naive data movement. * * Contributors * * Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany) * * References * * [1] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I. * SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342. * LAPACK Working note 169. * [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II. * SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362. * LAPACK Working note 170. * [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR * factorization software - a case study. * ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28. * LAPACK Working note 176. * [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, * QSVD, (H,K)-SVD computations. * Department of Mathematics, University of Zagreb, 2008. * * Bugs, examples and comments * * Please report all bugs and send interesting examples and/or comments to * drmac@math.hr. Thank you. * * ========================================================================== * * .. Local Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) * .. * .. Local Scalars .. DOUBLE PRECISION AAPP, AAQQ, AATMAX, AATMIN, BIG, BIG1, COND_OK, & CONDR1, CONDR2, ENTRA, ENTRAT, EPSLN, MAXPRJ, SCALEM, & SCONDA, SFMIN, SMALL, TEMP1, USCAL1, USCAL2, XSC INTEGER IERR, N1, NR, NUMRANK, p, q, WARNING LOGICAL ALMORT, DEFR, ERREST, GOSCAL, JRACC, KILL, LSVEC, & L2ABER, L2KILL, L2PERT, L2RANK, L2TRAN, & NOSCAL, ROWPIV, RSVEC, TRANSP * .. * .. Intrinsic Functions .. INTRINSIC DABS, DLOG, DMAX1, DMIN1, DBLE, & MAX0, MIN0, IDNINT, DSIGN, DSQRT * .. * .. External Functions .. DOUBLE PRECISION DLAMCH, DNRM2 INTEGER IDAMAX LOGICAL LSAME EXTERNAL IDAMAX, LSAME, DLAMCH, DNRM2 * .. * .. External Subroutines .. EXTERNAL DCOPY, DGELQF, DGEQP3, DGEQRF, DLACPY, DLASCL, & DLASET, DLASSQ, DLASWP, DORGQR, DORMLQ, & DORMQR, DPOCON, DSCAL, DSWAP, DTRSM, XERBLA * EXTERNAL DGESVJ * .. * * Test the input arguments * LSVEC = LSAME( JOBU, 'U' ) .OR. LSAME( JOBU, 'F' ) JRACC = LSAME( JOBV, 'J' ) RSVEC = LSAME( JOBV, 'V' ) .OR. JRACC ROWPIV = LSAME( JOBA, 'F' ) .OR. LSAME( JOBA, 'G' ) L2RANK = LSAME( JOBA, 'R' ) L2ABER = LSAME( JOBA, 'A' ) ERREST = LSAME( JOBA, 'E' ) .OR. LSAME( JOBA, 'G' ) L2TRAN = LSAME( JOBT, 'T' ) L2KILL = LSAME( JOBR, 'R' ) DEFR = LSAME( JOBR, 'N' ) L2PERT = LSAME( JOBP, 'P' ) * IF ( .NOT.(ROWPIV .OR. L2RANK .OR. L2ABER .OR. & ERREST .OR. LSAME( JOBA, 'C' ) )) THEN INFO = - 1 ELSE IF ( .NOT.( LSVEC .OR. LSAME( JOBU, 'N' ) .OR. & LSAME( JOBU, 'W' )) ) THEN INFO = - 2 ELSE IF ( .NOT.( RSVEC .OR. LSAME( JOBV, 'N' ) .OR. & LSAME( JOBV, 'W' )) .OR. ( JRACC .AND. (.NOT.LSVEC) ) ) THEN INFO = - 3 ELSE IF ( .NOT. ( L2KILL .OR. DEFR ) ) THEN INFO = - 4 ELSE IF ( .NOT. ( L2TRAN .OR. LSAME( JOBT, 'N' ) ) ) THEN INFO = - 5 ELSE IF ( .NOT. ( L2PERT .OR. LSAME( JOBP, 'N' ) ) ) THEN INFO = - 6 ELSE IF ( M .LT. 0 ) THEN INFO = - 7 ELSE IF ( ( N .LT. 0 ) .OR. ( N .GT. M ) ) THEN INFO = - 8 ELSE IF ( LDA .LT. M ) THEN INFO = - 10 ELSE IF ( LSVEC .AND. ( LDU .LT. M ) ) THEN INFO = - 13 ELSE IF ( RSVEC .AND. ( LDV .LT. N ) ) THEN INFO = - 14 ELSE IF ( (.NOT.(LSVEC .OR. RSVEC .OR. ERREST).AND. & (LWORK .LT. MAX0(7,4*N+1,2*M+N))) .OR. & (.NOT.(LSVEC .OR. LSVEC) .AND. ERREST .AND. & (LWORK .LT. MAX0(7,4*N+N*N,2*M+N))) .OR. & (LSVEC .AND. (.NOT.RSVEC) .AND. (LWORK .LT. MAX0(7,2*N+M))) .OR. & (RSVEC .AND. (.NOT.LSVEC) .AND. (LWORK .LT. MAX0(7,2*N+M))) .OR. & (LSVEC .AND. RSVEC .AND. .NOT.JRACC .AND. (LWORK.LT.6*N+2*N*N)) & .OR. (LSVEC.AND.RSVEC.AND.JRACC.AND.LWORK.LT.MAX0(7,M+3*N+N*N))) & THEN INFO = - 17 ELSE * #:) INFO = 0 END IF * IF ( INFO .NE. 0 ) THEN * #:( CALL XERBLA( 'DGEJSV', - INFO ) END IF * * Quick return for void matrix (Y3K safe) * #:) IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) RETURN * * Determine whether the matrix U should be M x N or M x M * IF ( LSVEC ) THEN N1 = N IF ( LSAME( JOBU, 'F' ) ) N1 = M END IF * * Set numerical parameters * *! NOTE: Make sure DLAMCH() does not fail on the target architecture. * EPSLN = DLAMCH('Epsilon') SFMIN = DLAMCH('SafeMinimum') SMALL = SFMIN / EPSLN BIG = DLAMCH('O') * BIG = ONE / SFMIN * * Initialize SVA(1:N) = diag( ||A e_i||_2 )_1^N * *(!) If necessary, scale SVA() to protect the largest norm from * overflow. It is possible that this scaling pushes the smallest * column norm left from the underflow threshold (extreme case). * SCALEM = ONE / DSQRT(DBLE(M)*DBLE(N)) NOSCAL = .TRUE. GOSCAL = .TRUE. DO 1874 p = 1, N AAPP = ZERO AAQQ = ONE CALL DLASSQ( M, A(1,p), 1, AAPP, AAQQ ) IF ( AAPP .GT. BIG ) THEN INFO = - 9 CALL XERBLA( 'DGEJSV', -INFO ) RETURN END IF AAQQ = DSQRT(AAQQ) IF ( ( AAPP .LT. (BIG / AAQQ) ) .AND. NOSCAL ) THEN SVA(p) = AAPP * AAQQ ELSE NOSCAL = .FALSE. SVA(p) = AAPP * ( AAQQ * SCALEM ) IF ( GOSCAL ) THEN GOSCAL = .FALSE. CALL DSCAL( p-1, SCALEM, SVA, 1 ) END IF END IF 1874 CONTINUE * IF ( NOSCAL ) SCALEM = ONE * AAPP = ZERO AAQQ = BIG DO 4781 p = 1, N AAPP = DMAX1( AAPP, SVA(p) ) IF ( SVA(p) .NE. ZERO ) AAQQ = DMIN1( AAQQ, SVA(p) ) 4781 CONTINUE * * Quick return for zero M x N matrix * #:) IF ( AAPP .EQ. ZERO ) THEN IF ( LSVEC ) CALL DLASET( 'G', M, N1, ZERO, ONE, U, LDU ) IF ( RSVEC ) CALL DLASET( 'G', N, N, ZERO, ONE, V, LDV ) WORK(1) = ONE WORK(2) = ONE IF ( ERREST ) WORK(3) = ONE IF ( LSVEC .AND. RSVEC ) THEN WORK(4) = ONE WORK(5) = ONE END IF IF ( L2TRAN ) THEN WORK(6) = ZERO WORK(7) = ZERO END IF IWORK(1) = 0 IWORK(2) = 0 RETURN END IF * * Issue warning if denormalized column norms detected. Override the * high relative accuracy request. Issue licence to kill columns * (set them to zero) whose norm is less than sigma_max / BIG (roughly). * #:( WARNING = 0 IF ( AAQQ .LE. SFMIN ) THEN L2RANK = .TRUE. L2KILL = .TRUE. WARNING = 1 END IF * * Quick return for one-column matrix * #:) IF ( N .EQ. 1 ) THEN * IF ( LSVEC ) THEN CALL DLASCL( 'G',0,0,SVA(1),SCALEM, M,1,A(1,1),LDA,IERR ) CALL DLACPY( 'A', M, 1, A, LDA, U, LDU ) * computing all M left singular vectors of the M x 1 matrix IF ( N1 .NE. N ) THEN CALL DGEQRF( M, N, U,LDU, WORK, WORK(N+1),LWORK-N,IERR ) CALL DORGQR( M,N1,1, U,LDU,WORK,WORK(N+1),LWORK-N,IERR ) CALL DCOPY( M, A(1,1), 1, U(1,1), 1 ) END IF END IF IF ( RSVEC ) THEN V(1,1) = ONE END IF IF ( SVA(1) .LT. (BIG*SCALEM) ) THEN SVA(1) = SVA(1) / SCALEM SCALEM = ONE END IF WORK(1) = ONE / SCALEM WORK(2) = ONE IF ( SVA(1) .NE. ZERO ) THEN IWORK(1) = 1 IF ( ( SVA(1) / SCALEM) .GE. SFMIN ) THEN IWORK(2) = 1 ELSE IWORK(2) = 0 END IF ELSE IWORK(1) = 0 IWORK(2) = 0 END IF IF ( ERREST ) WORK(3) = ONE IF ( LSVEC .AND. RSVEC ) THEN WORK(4) = ONE WORK(5) = ONE END IF IF ( L2TRAN ) THEN WORK(6) = ZERO WORK(7) = ZERO END IF RETURN * END IF * TRANSP = .FALSE. L2TRAN = L2TRAN .AND. ( M .EQ. N ) * AATMAX = -ONE AATMIN = BIG IF ( ROWPIV .OR. L2TRAN ) THEN * * Compute the row norms, needed to determine row pivoting sequence * (in the case of heavily row weighted A, row pivoting is strongly * advised) and to collect information needed to compare the * structures of A * A^t and A^t * A (in the case L2TRAN.EQ..TRUE.). * IF ( L2TRAN ) THEN DO 1950 p = 1, M XSC = ZERO TEMP1 = ONE CALL DLASSQ( N, A(p,1), LDA, XSC, TEMP1 ) * DLASSQ gets both the ell_2 and the ell_infinity norm * in one pass through the vector WORK(M+N+p) = XSC * SCALEM WORK(N+p) = XSC * (SCALEM*DSQRT(TEMP1)) AATMAX = DMAX1( AATMAX, WORK(N+p) ) IF (WORK(N+p) .NE. ZERO) AATMIN = DMIN1(AATMIN,WORK(N+p)) 1950 CONTINUE ELSE DO 1904 p = 1, M WORK(M+N+p) = SCALEM*DABS( A(p,IDAMAX(N,A(p,1),LDA)) ) AATMAX = DMAX1( AATMAX, WORK(M+N+p) ) AATMIN = DMIN1( AATMIN, WORK(M+N+p) ) 1904 CONTINUE END IF * END IF * * For square matrix A try to determine whether A^t would be better * input for the preconditioned Jacobi SVD, with faster convergence. * The decision is based on an O(N) function of the vector of column * and row norms of A, based on the Shannon entropy. This should give * the right choice in most cases when the difference actually matters. * It may fail and pick the slower converging side. * ENTRA = ZERO ENTRAT = ZERO IF ( L2TRAN ) THEN * XSC = ZERO TEMP1 = ONE CALL DLASSQ( N, SVA, 1, XSC, TEMP1 ) TEMP1 = ONE / TEMP1 * ENTRA = ZERO DO 1113 p = 1, N BIG1 = ( ( SVA(p) / XSC )**2 ) * TEMP1 IF ( BIG1 .NE. ZERO ) ENTRA = ENTRA + BIG1 * DLOG(BIG1) 1113 CONTINUE ENTRA = - ENTRA / DLOG(DBLE(N)) * * Now, SVA().^2/Trace(A^t * A) is a point in the probability simplex. * It is derived from the diagonal of A^t * A. Do the same with the * diagonal of A * A^t, compute the entropy of the corresponding * probability distribution. Note that A * A^t and A^t * A have the * same trace. * ENTRAT = ZERO DO 1114 p = N+1, N+M BIG1 = ( ( WORK(p) / XSC )**2 ) * TEMP1 IF ( BIG1 .NE. ZERO ) ENTRAT = ENTRAT + BIG1 * DLOG(BIG1) 1114 CONTINUE ENTRAT = - ENTRAT / DLOG(DBLE(M)) * * Analyze the entropies and decide A or A^t. Smaller entropy * usually means better input for the algorithm. * TRANSP = ( ENTRAT .LT. ENTRA ) * * If A^t is better than A, transpose A. * IF ( TRANSP ) THEN * In an optimal implementation, this trivial transpose * should be replaced with faster transpose. DO 1115 p = 1, N - 1 DO 1116 q = p + 1, N TEMP1 = A(q,p) A(q,p) = A(p,q) A(p,q) = TEMP1 1116 CONTINUE 1115 CONTINUE DO 1117 p = 1, N WORK(M+N+p) = SVA(p) SVA(p) = WORK(N+p) 1117 CONTINUE TEMP1 = AAPP AAPP = AATMAX AATMAX = TEMP1 TEMP1 = AAQQ AAQQ = AATMIN AATMIN = TEMP1 KILL = LSVEC LSVEC = RSVEC RSVEC = KILL IF ( LSVEC ) N1 = N * ROWPIV = .TRUE. END IF * END IF * END IF L2TRAN * * Scale the matrix so that its maximal singular value remains less * than DSQRT(BIG) -- the matrix is scaled so that its maximal column * has Euclidean norm equal to DSQRT(BIG/N). The only reason to keep * DSQRT(BIG) instead of BIG is the fact that DGEJSV uses LAPACK and * BLAS routines that, in some implementations, are not capable of * working in the full interval [SFMIN,BIG] and that they may provoke * overflows in the intermediate results. If the singular values spread * from SFMIN to BIG, then DGESVJ will compute them. So, in that case, * one should use DGESVJ instead of DGEJSV. * BIG1 = DSQRT( BIG ) TEMP1 = DSQRT( BIG / DBLE(N) ) * CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, N, 1, SVA, N, IERR ) IF ( AAQQ .GT. (AAPP * SFMIN) ) THEN AAQQ = ( AAQQ / AAPP ) * TEMP1 ELSE AAQQ = ( AAQQ * TEMP1 ) / AAPP END IF TEMP1 = TEMP1 * SCALEM CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, M, N, A, LDA, IERR ) * * To undo scaling at the end of this procedure, multiply the * computed singular values with USCAL2 / USCAL1. * USCAL1 = TEMP1 USCAL2 = AAPP * IF ( L2KILL ) THEN * L2KILL enforces computation of nonzero singular values in * the restricted range of condition number of the initial A, * sigma_max(A) / sigma_min(A) approx. DSQRT(BIG)/DSQRT(SFMIN). XSC = DSQRT( SFMIN ) ELSE XSC = SMALL * * Now, if the condition number of A is too big, * sigma_max(A) / sigma_min(A) .GT. DSQRT(BIG/N) * EPSLN / SFMIN, * as a precaution measure, the full SVD is computed using DGESVJ * with accumulated Jacobi rotations. This provides numerically * more robust computation, at the cost of slightly increased run * time. Depending on the concrete implementation of BLAS and LAPACK * (i.e. how they behave in presence of extreme ill-conditioning) the * implementor may decide to remove this switch. IF ( ( AAQQ.LT.DSQRT(SFMIN) ) .AND. LSVEC .AND. RSVEC ) THEN JRACC = .TRUE. END IF * END IF IF ( AAQQ .LT. XSC ) THEN DO 700 p = 1, N IF ( SVA(p) .LT. XSC ) THEN CALL DLASET( 'A', M, 1, ZERO, ZERO, A(1,p), LDA ) SVA(p) = ZERO END IF 700 CONTINUE END IF * * Preconditioning using QR factorization with pivoting * IF ( ROWPIV ) THEN * Optional row permutation (Bjoerck row pivoting): * A result by Cox and Higham shows that the Bjoerck's * row pivoting combined with standard column pivoting * has similar effect as Powell-Reid complete pivoting. * The ell-infinity norms of A are made nonincreasing. DO 1952 p = 1, M - 1 q = IDAMAX( M-p+1, WORK(M+N+p), 1 ) + p - 1 IWORK(2*N+p) = q IF ( p .NE. q ) THEN TEMP1 = WORK(M+N+p) WORK(M+N+p) = WORK(M+N+q) WORK(M+N+q) = TEMP1 END IF 1952 CONTINUE CALL DLASWP( N, A, LDA, 1, M-1, IWORK(2*N+1), 1 ) END IF * * End of the preparation phase (scaling, optional sorting and * transposing, optional flushing of small columns). * * Preconditioning * * If the full SVD is needed, the right singular vectors are computed * from a matrix equation, and for that we need theoretical analysis * of the Businger-Golub pivoting. So we use DGEQP3 as the first RR QRF. * In all other cases the first RR QRF can be chosen by other criteria * (eg speed by replacing global with restricted window pivoting, such * as in SGEQPX from TOMS # 782). Good results will be obtained using * SGEQPX with properly (!) chosen numerical parameters. * Any improvement of DGEQP3 improves overal performance of DGEJSV. * * A * P1 = Q1 * [ R1^t 0]^t: DO 1963 p = 1, N * .. all columns are free columns IWORK(p) = 0 1963 CONTINUE CALL DGEQP3( M,N,A,LDA, IWORK,WORK, WORK(N+1),LWORK-N, IERR ) * * The upper triangular matrix R1 from the first QRF is inspected for * rank deficiency and possibilities for deflation, or possible * ill-conditioning. Depending on the user specified flag L2RANK, * the procedure explores possibilities to reduce the numerical * rank by inspecting the computed upper triangular factor. If * L2RANK or L2ABER are up, then DGEJSV will compute the SVD of * A + dA, where ||dA|| <= f(M,N)*EPSLN. * NR = 1 IF ( L2ABER ) THEN * Standard absolute error bound suffices. All sigma_i with * sigma_i < N*EPSLN*||A|| are flushed to zero. This is an * agressive enforcement of lower numerical rank by introducing a * backward error of the order of N*EPSLN*||A||. TEMP1 = DSQRT(DBLE(N))*EPSLN DO 3001 p = 2, N IF ( DABS(A(p,p)) .GE. (TEMP1*DABS(A(1,1))) ) THEN NR = NR + 1 ELSE GO TO 3002 END IF 3001 CONTINUE 3002 CONTINUE ELSE IF ( L2RANK ) THEN * .. similarly as above, only slightly more gentle (less agressive). * Sudden drop on the diagonal of R1 is used as the criterion for * close-to-rank-defficient. TEMP1 = DSQRT(SFMIN) DO 3401 p = 2, N IF ( ( DABS(A(p,p)) .LT. (EPSLN*DABS(A(p-1,p-1))) ) .OR. & ( DABS(A(p,p)) .LT. SMALL ) .OR. & ( L2KILL .AND. (DABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3402 NR = NR + 1 3401 CONTINUE 3402 CONTINUE * ELSE * The goal is high relative accuracy. However, if the matrix * has high scaled condition number the relative accuracy is in * general not feasible. Later on, a condition number estimator * will be deployed to estimate the scaled condition number. * Here we just remove the underflowed part of the triangular * factor. This prevents the situation in which the code is * working hard to get the accuracy not warranted by the data. TEMP1 = DSQRT(SFMIN) DO 3301 p = 2, N IF ( ( DABS(A(p,p)) .LT. SMALL ) .OR. & ( L2KILL .AND. (DABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3302 NR = NR + 1 3301 CONTINUE 3302 CONTINUE * END IF * ALMORT = .FALSE. IF ( NR .EQ. N ) THEN MAXPRJ = ONE DO 3051 p = 2, N TEMP1 = DABS(A(p,p)) / SVA(IWORK(p)) MAXPRJ = DMIN1( MAXPRJ, TEMP1 ) 3051 CONTINUE IF ( MAXPRJ**2 .GE. ONE - DBLE(N)*EPSLN ) ALMORT = .TRUE. END IF * * SCONDA = - ONE CONDR1 = - ONE CONDR2 = - ONE * IF ( ERREST ) THEN IF ( N .EQ. NR ) THEN IF ( RSVEC ) THEN * .. V is available as workspace CALL DLACPY( 'U', N, N, A, LDA, V, LDV ) DO 3053 p = 1, N TEMP1 = SVA(IWORK(p)) CALL DSCAL( p, ONE/TEMP1, V(1,p), 1 ) 3053 CONTINUE CALL DPOCON( 'U', N, V, LDV, ONE, TEMP1, & WORK(N+1), IWORK(2*N+M+1), IERR ) ELSE IF ( LSVEC ) THEN * .. U is available as workspace CALL DLACPY( 'U', N, N, A, LDA, U, LDU ) DO 3054 p = 1, N TEMP1 = SVA(IWORK(p)) CALL DSCAL( p, ONE/TEMP1, U(1,p), 1 ) 3054 CONTINUE CALL DPOCON( 'U', N, U, LDU, ONE, TEMP1, & WORK(N+1), IWORK(2*N+M+1), IERR ) ELSE CALL DLACPY( 'U', N, N, A, LDA, WORK(N+1), N ) DO 3052 p = 1, N TEMP1 = SVA(IWORK(p)) CALL DSCAL( p, ONE/TEMP1, WORK(N+(p-1)*N+1), 1 ) 3052 CONTINUE * .. the columns of R are scaled to have unit Euclidean lengths. CALL DPOCON( 'U', N, WORK(N+1), N, ONE, TEMP1, & WORK(N+N*N+1), IWORK(2*N+M+1), IERR ) END IF SCONDA = ONE / DSQRT(TEMP1) * SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1). * N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA ELSE SCONDA = - ONE END IF END IF * L2PERT = L2PERT .AND. ( DABS( A(1,1)/A(NR,NR) ) .GT. DSQRT(BIG1) ) * If there is no violent scaling, artificial perturbation is not needed. * * Phase 3: * IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN * * Singular Values only * * .. transpose A(1:NR,1:N) DO 1946 p = 1, MIN0( N-1, NR ) CALL DCOPY( N-p, A(p,p+1), LDA, A(p+1,p), 1 ) 1946 CONTINUE * * The following two DO-loops introduce small relative perturbation * into the strict upper triangle of the lower triangular matrix. * Small entries below the main diagonal are also changed. * This modification is useful if the computing environment does not * provide/allow FLUSH TO ZERO underflow, for it prevents many * annoying denormalized numbers in case of strongly scaled matrices. * The perturbation is structured so that it does not introduce any * new perturbation of the singular values, and it does not destroy * the job done by the preconditioner. * The licence for this perturbation is in the variable L2PERT, which * should be .FALSE. if FLUSH TO ZERO underflow is active. * IF ( .NOT. ALMORT ) THEN * IF ( L2PERT ) THEN * XSC = DSQRT(SMALL) XSC = EPSLN / DBLE(N) DO 4947 q = 1, NR TEMP1 = XSC*DABS(A(q,q)) DO 4949 p = 1, N IF ( ( (p.GT.q) .AND. (DABS(A(p,q)).LE.TEMP1) ) & .OR. ( p .LT. q ) ) & A(p,q) = DSIGN( TEMP1, A(p,q) ) 4949 CONTINUE 4947 CONTINUE ELSE CALL DLASET( 'U', NR-1,NR-1, ZERO,ZERO, A(1,2),LDA ) END IF * * .. second preconditioning using the QR factorization * CALL DGEQRF( N,NR, A,LDA, WORK, WORK(N+1),LWORK-N, IERR ) * * .. and transpose upper to lower triangular DO 1948 p = 1, NR - 1 CALL DCOPY( NR-p, A(p,p+1), LDA, A(p+1,p), 1 ) 1948 CONTINUE * END IF * * Row-cyclic Jacobi SVD algorithm with column pivoting * * .. again some perturbation (a "background noise") is added * to drown denormals IF ( L2PERT ) THEN * XSC = DSQRT(SMALL) XSC = EPSLN / DBLE(N) DO 1947 q = 1, NR TEMP1 = XSC*DABS(A(q,q)) DO 1949 p = 1, NR IF ( ( (p.GT.q) .AND. (DABS(A(p,q)).LE.TEMP1) ) & .OR. ( p .LT. q ) ) & A(p,q) = DSIGN( TEMP1, A(p,q) ) 1949 CONTINUE 1947 CONTINUE ELSE CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, A(1,2), LDA ) END IF * * .. and one-sided Jacobi rotations are started on a lower * triangular matrix (plus perturbation which is ignored in * the part which destroys triangular form (confusing?!)) * CALL DGESVJ( 'L', 'NoU', 'NoV', NR, NR, A, LDA, SVA, & N, V, LDV, WORK, LWORK, INFO ) * SCALEM = WORK(1) NUMRANK = IDNINT(WORK(2)) * * ELSE IF ( RSVEC .AND. ( .NOT. LSVEC ) ) THEN * * -> Singular Values and Right Singular Vectors <- * IF ( ALMORT ) THEN * * .. in this case NR equals N DO 1998 p = 1, NR CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 ) 1998 CONTINUE CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV ) * CALL DGESVJ( 'L','U','N', N, NR, V,LDV, SVA, NR, A,LDA, & WORK, LWORK, INFO ) SCALEM = WORK(1) NUMRANK = IDNINT(WORK(2)) ELSE * * .. two more QR factorizations ( one QRF is not enough, two require * accumulated product of Jacobi rotations, three are perfect ) * CALL DLASET( 'Lower', NR-1, NR-1, ZERO, ZERO, A(2,1), LDA ) CALL DGELQF( NR, N, A, LDA, WORK, WORK(N+1), LWORK-N, IERR) CALL DLACPY( 'Lower', NR, NR, A, LDA, V, LDV ) CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV ) CALL DGEQRF( NR, NR, V, LDV, WORK(N+1), WORK(2*N+1), & LWORK-2*N, IERR ) DO 8998 p = 1, NR CALL DCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 ) 8998 CONTINUE CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV ) * CALL DGESVJ( 'Lower', 'U','N', NR, NR, V,LDV, SVA, NR, U, & LDU, WORK(N+1), LWORK, INFO ) SCALEM = WORK(N+1) NUMRANK = IDNINT(WORK(N+2)) IF ( NR .LT. N ) THEN CALL DLASET( 'A',N-NR, NR, ZERO,ZERO, V(NR+1,1), LDV ) CALL DLASET( 'A',NR, N-NR, ZERO,ZERO, V(1,NR+1), LDV ) CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE, V(NR+1,NR+1), LDV ) END IF * CALL DORMLQ( 'Left', 'Transpose', N, N, NR, A, LDA, WORK, & V, LDV, WORK(N+1), LWORK-N, IERR ) * END IF * DO 8991 p = 1, N CALL DCOPY( N, V(p,1), LDV, A(IWORK(p),1), LDA ) 8991 CONTINUE CALL DLACPY( 'All', N, N, A, LDA, V, LDV ) * IF ( TRANSP ) THEN CALL DLACPY( 'All', N, N, V, LDV, U, LDU ) END IF * ELSE IF ( LSVEC .AND. ( .NOT. RSVEC ) ) THEN * * .. Singular Values and Left Singular Vectors .. * * .. second preconditioning step to avoid need to accumulate * Jacobi rotations in the Jacobi iterations. DO 1965 p = 1, NR CALL DCOPY( N-p+1, A(p,p), LDA, U(p,p), 1 ) 1965 CONTINUE CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU ) * CALL DGEQRF( N, NR, U, LDU, WORK(N+1), WORK(2*N+1), & LWORK-2*N, IERR ) * DO 1967 p = 1, NR - 1 CALL DCOPY( NR-p, U(p,p+1), LDU, U(p+1,p), 1 ) 1967 CONTINUE CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU ) * CALL DGESVJ( 'Lower', 'U', 'N', NR,NR, U, LDU, SVA, NR, A, & LDA, WORK(N+1), LWORK-N, INFO ) SCALEM = WORK(N+1) NUMRANK = IDNINT(WORK(N+2)) * IF ( NR .LT. M ) THEN CALL DLASET( 'A', M-NR, NR,ZERO, ZERO, U(NR+1,1), LDU ) IF ( NR .LT. N1 ) THEN CALL DLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1), LDU ) CALL DLASET( 'A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1), LDU ) END IF END IF * CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U, & LDU, WORK(N+1), LWORK-N, IERR ) * IF ( ROWPIV ) & CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) * DO 1974 p = 1, N1 XSC = ONE / DNRM2( M, U(1,p), 1 ) CALL DSCAL( M, XSC, U(1,p), 1 ) 1974 CONTINUE * IF ( TRANSP ) THEN CALL DLACPY( 'All', N, N, U, LDU, V, LDV ) END IF * ELSE * * .. Full SVD .. * IF ( .NOT. JRACC ) THEN * IF ( .NOT. ALMORT ) THEN * * Second Preconditioning Step (QRF [with pivoting]) * Note that the composition of TRANSPOSE, QRF and TRANSPOSE is * equivalent to an LQF CALL. Since in many libraries the QRF * seems to be better optimized than the LQF, we do explicit * transpose and use the QRF. This is subject to changes in an * optimized implementation of DGEJSV. * DO 1968 p = 1, NR CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 ) 1968 CONTINUE * * .. the following two loops perturb small entries to avoid * denormals in the second QR factorization, where they are * as good as zeros. This is done to avoid painfully slow * computation with denormals. The relative size of the perturbation * is a parameter that can be changed by the implementer. * This perturbation device will be obsolete on machines with * properly implemented arithmetic. * To switch it off, set L2PERT=.FALSE. To remove it from the * code, remove the action under L2PERT=.TRUE., leave the ELSE part. * The following two loops should be blocked and fused with the * transposed copy above. * IF ( L2PERT ) THEN XSC = DSQRT(SMALL) DO 2969 q = 1, NR TEMP1 = XSC*DABS( V(q,q) ) DO 2968 p = 1, N IF ( ( p .GT. q ) .AND. ( DABS(V(p,q)) .LE. TEMP1 ) & .OR. ( p .LT. q ) ) & V(p,q) = DSIGN( TEMP1, V(p,q) ) IF ( p. LT. q ) V(p,q) = - V(p,q) 2968 CONTINUE 2969 CONTINUE ELSE CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV ) END IF * * Estimate the row scaled condition number of R1 * (If R1 is rectangular, N > NR, then the condition number * of the leading NR x NR submatrix is estimated.) * CALL DLACPY( 'L', NR, NR, V, LDV, WORK(2*N+1), NR ) DO 3950 p = 1, NR TEMP1 = DNRM2(NR-p+1,WORK(2*N+(p-1)*NR+p),1) CALL DSCAL(NR-p+1,ONE/TEMP1,WORK(2*N+(p-1)*NR+p),1) 3950 CONTINUE CALL DPOCON('Lower',NR,WORK(2*N+1),NR,ONE,TEMP1, & WORK(2*N+NR*NR+1),IWORK(M+2*N+1),IERR) CONDR1 = ONE / DSQRT(TEMP1) * .. here need a second oppinion on the condition number * .. then assume worst case scenario * R1 is OK for inverse <=> CONDR1 .LT. DBLE(N) * more conservative <=> CONDR1 .LT. DSQRT(DBLE(N)) * COND_OK = DSQRT(DBLE(NR)) *[TP] COND_OK is a tuning parameter. IF ( CONDR1 .LT. COND_OK ) THEN * .. the second QRF without pivoting. Note: in an optimized * implementation, this QRF should be implemented as the QRF * of a lower triangular matrix. * R1^t = Q2 * R2 CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1), & LWORK-2*N, IERR ) * IF ( L2PERT ) THEN XSC = DSQRT(SMALL)/EPSLN DO 3959 p = 2, NR DO 3958 q = 1, p - 1 TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q))) IF ( DABS(V(q,p)) .LE. TEMP1 ) & V(q,p) = DSIGN( TEMP1, V(q,p) ) 3958 CONTINUE 3959 CONTINUE END IF * IF ( NR .NE. N ) * .. save ... & CALL DLACPY( 'A', N, NR, V, LDV, WORK(2*N+1), N ) * * .. this transposed copy should be better than naive DO 1969 p = 1, NR - 1 CALL DCOPY( NR-p, V(p,p+1), LDV, V(p+1,p), 1 ) 1969 CONTINUE * CONDR2 = CONDR1 * ELSE * * .. ill-conditioned case: second QRF with pivoting * Note that windowed pivoting would be equaly good * numerically, and more run-time efficient. So, in * an optimal implementation, the next call to DGEQP3 * should be replaced with eg. CALL SGEQPX (ACM TOMS #782) * with properly (carefully) chosen parameters. * * R1^t * P2 = Q2 * R2 DO 3003 p = 1, NR IWORK(N+p) = 0 3003 CONTINUE CALL DGEQP3( N, NR, V, LDV, IWORK(N+1), WORK(N+1), & WORK(2*N+1), LWORK-2*N, IERR ) ** CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1), ** & LWORK-2*N, IERR ) IF ( L2PERT ) THEN XSC = DSQRT(SMALL) DO 3969 p = 2, NR DO 3968 q = 1, p - 1 TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q))) IF ( DABS(V(q,p)) .LE. TEMP1 ) & V(q,p) = DSIGN( TEMP1, V(q,p) ) 3968 CONTINUE 3969 CONTINUE END IF * CALL DLACPY( 'A', N, NR, V, LDV, WORK(2*N+1), N ) * IF ( L2PERT ) THEN XSC = DSQRT(SMALL) DO 8970 p = 2, NR DO 8971 q = 1, p - 1 TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q))) V(p,q) = - DSIGN( TEMP1, V(q,p) ) 8971 CONTINUE 8970 CONTINUE ELSE CALL DLASET( 'L',NR-1,NR-1,ZERO,ZERO,V(2,1),LDV ) END IF * Now, compute R2 = L3 * Q3, the LQ factorization. CALL DGELQF( NR, NR, V, LDV, WORK(2*N+N*NR+1), & WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, IERR ) * .. and estimate the condition number CALL DLACPY( 'L',NR,NR,V,LDV,WORK(2*N+N*NR+NR+1),NR ) DO 4950 p = 1, NR TEMP1 = DNRM2( p, WORK(2*N+N*NR+NR+p), NR ) CALL DSCAL( p, ONE/TEMP1, WORK(2*N+N*NR+NR+p), NR ) 4950 CONTINUE CALL DPOCON( 'L',NR,WORK(2*N+N*NR+NR+1),NR,ONE,TEMP1, & WORK(2*N+N*NR+NR+NR*NR+1),IWORK(M+2*N+1),IERR ) CONDR2 = ONE / DSQRT(TEMP1) * IF ( CONDR2 .GE. COND_OK ) THEN * .. save the Householder vectors used for Q3 * (this overwrittes the copy of R2, as it will not be * needed in this branch, but it does not overwritte the * Huseholder vectors of Q2.). CALL DLACPY( 'U', NR, NR, V, LDV, WORK(2*N+1), N ) * .. and the rest of the information on Q3 is in * WORK(2*N+N*NR+1:2*N+N*NR+N) END IF * END IF * IF ( L2PERT ) THEN XSC = DSQRT(SMALL) DO 4968 q = 2, NR TEMP1 = XSC * V(q,q) DO 4969 p = 1, q - 1 * V(p,q) = - DSIGN( TEMP1, V(q,p) ) V(p,q) = - DSIGN( TEMP1, V(p,q) ) 4969 CONTINUE 4968 CONTINUE ELSE CALL DLASET( 'U', NR-1,NR-1, ZERO,ZERO, V(1,2), LDV ) END IF * * Second preconditioning finished; continue with Jacobi SVD * The input matrix is lower trinagular. * * Recover the right singular vectors as solution of a well * conditioned triangular matrix equation. * IF ( CONDR1 .LT. COND_OK ) THEN * CALL DGESVJ( 'L','U','N',NR,NR,V,LDV,SVA,NR,U, & LDU,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,INFO ) SCALEM = WORK(2*N+N*NR+NR+1) NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2)) DO 3970 p = 1, NR CALL DCOPY( NR, V(1,p), 1, U(1,p), 1 ) CALL DSCAL( NR, SVA(p), V(1,p), 1 ) 3970 CONTINUE * .. pick the right matrix equation and solve it * IF ( NR. EQ. N ) THEN * :)) .. best case, R1 is inverted. The solution of this matrix * equation is Q2*V2 = the product of the Jacobi rotations * used in DGESVJ, premultiplied with the orthogonal matrix * from the second QR factorization. CALL DTRSM( 'L','U','N','N', NR,NR,ONE, A,LDA, V,LDV ) ELSE * .. R1 is well conditioned, but non-square. Transpose(R2) * is inverted to get the product of the Jacobi rotations * used in DGESVJ. The Q-factor from the second QR * factorization is then built in explicitly. CALL DTRSM('L','U','T','N',NR,NR,ONE,WORK(2*N+1), & N,V,LDV) IF ( NR .LT. N ) THEN CALL DLASET('A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV) CALL DLASET('A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV) CALL DLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV) END IF CALL DORMQR('L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1), & V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR) END IF * ELSE IF ( CONDR2 .LT. COND_OK ) THEN * * :) .. the input matrix A is very likely a relative of * the Kahan matrix :) * The matrix R2 is inverted. The solution of the matrix equation * is Q3^T*V3 = the product of the Jacobi rotations (appplied to * the lower triangular L3 from the LQ factorization of * R2=L3*Q3), pre-multiplied with the transposed Q3. CALL DGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR, U, & LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO ) SCALEM = WORK(2*N+N*NR+NR+1) NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2)) DO 3870 p = 1, NR CALL DCOPY( NR, V(1,p), 1, U(1,p), 1 ) CALL DSCAL( NR, SVA(p), U(1,p), 1 ) 3870 CONTINUE CALL DTRSM('L','U','N','N',NR,NR,ONE,WORK(2*N+1),N,U,LDU) * .. apply the permutation from the second QR factorization DO 873 q = 1, NR DO 872 p = 1, NR WORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q) 872 CONTINUE DO 874 p = 1, NR U(p,q) = WORK(2*N+N*NR+NR+p) 874 CONTINUE 873 CONTINUE IF ( NR .LT. N ) THEN CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV ) CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV ) CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV ) END IF CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1), & V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR ) ELSE * Last line of defense. * #:( This is a rather pathological case: no scaled condition * improvement after two pivoted QR factorizations. Other * possibility is that the rank revealing QR factorization * or the condition estimator has failed, or the COND_OK * is set very close to ONE (which is unnecessary). Normally, * this branch should never be executed, but in rare cases of * failure of the RRQR or condition estimator, the last line of * defense ensures that DGEJSV completes the task. * Compute the full SVD of L3 using DGESVJ with explicit * accumulation of Jacobi rotations. CALL DGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR, U, & LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO ) SCALEM = WORK(2*N+N*NR+NR+1) NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2)) IF ( NR .LT. N ) THEN CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV ) CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV ) CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV ) END IF CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1), & V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR ) * CALL DORMLQ( 'L', 'T', NR, NR, NR, WORK(2*N+1), N, & WORK(2*N+N*NR+1), U, LDU, WORK(2*N+N*NR+NR+1), & LWORK-2*N-N*NR-NR, IERR ) DO 773 q = 1, NR DO 772 p = 1, NR WORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q) 772 CONTINUE DO 774 p = 1, NR U(p,q) = WORK(2*N+N*NR+NR+p) 774 CONTINUE 773 CONTINUE * END IF * * Permute the rows of V using the (column) permutation from the * first QRF. Also, scale the columns to make them unit in * Euclidean norm. This applies to all cases. * TEMP1 = DSQRT(DBLE(N)) * EPSLN DO 1972 q = 1, N DO 972 p = 1, N WORK(2*N+N*NR+NR+IWORK(p)) = V(p,q) 972 CONTINUE DO 973 p = 1, N V(p,q) = WORK(2*N+N*NR+NR+p) 973 CONTINUE XSC = ONE / DNRM2( N, V(1,q), 1 ) IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) & CALL DSCAL( N, XSC, V(1,q), 1 ) 1972 CONTINUE * At this moment, V contains the right singular vectors of A. * Next, assemble the left singular vector matrix U (M x N). IF ( NR .LT. M ) THEN CALL DLASET( 'A', M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU ) IF ( NR .LT. N1 ) THEN CALL DLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU) CALL DLASET('A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1),LDU) END IF END IF * * The Q matrix from the first QRF is built into the left singular * matrix U. This applies to all cases. * CALL DORMQR( 'Left', 'No_Tr', M, N1, N, A, LDA, WORK, U, & LDU, WORK(N+1), LWORK-N, IERR ) * The columns of U are normalized. The cost is O(M*N) flops. TEMP1 = DSQRT(DBLE(M)) * EPSLN DO 1973 p = 1, NR XSC = ONE / DNRM2( M, U(1,p), 1 ) IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) & CALL DSCAL( M, XSC, U(1,p), 1 ) 1973 CONTINUE * * If the initial QRF is computed with row pivoting, the left * singular vectors must be adjusted. * IF ( ROWPIV ) & CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) * ELSE * * .. the initial matrix A has almost orthogonal columns and * the second QRF is not needed * CALL DLACPY( 'Upper', N, N, A, LDA, WORK(N+1), N ) IF ( L2PERT ) THEN XSC = DSQRT(SMALL) DO 5970 p = 2, N TEMP1 = XSC * WORK( N + (p-1)*N + p ) DO 5971 q = 1, p - 1 WORK(N+(q-1)*N+p)=-DSIGN(TEMP1,WORK(N+(p-1)*N+q)) 5971 CONTINUE 5970 CONTINUE ELSE CALL DLASET( 'Lower',N-1,N-1,ZERO,ZERO,WORK(N+2),N ) END IF * CALL DGESVJ( 'Upper', 'U', 'N', N, N, WORK(N+1), N, SVA, & N, U, LDU, WORK(N+N*N+1), LWORK-N-N*N, INFO ) * SCALEM = WORK(N+N*N+1) NUMRANK = IDNINT(WORK(N+N*N+2)) DO 6970 p = 1, N CALL DCOPY( N, WORK(N+(p-1)*N+1), 1, U(1,p), 1 ) CALL DSCAL( N, SVA(p), WORK(N+(p-1)*N+1), 1 ) 6970 CONTINUE * CALL DTRSM( 'Left', 'Upper', 'NoTrans', 'No UD', N, N, & ONE, A, LDA, WORK(N+1), N ) DO 6972 p = 1, N CALL DCOPY( N, WORK(N+p), N, V(IWORK(p),1), LDV ) 6972 CONTINUE TEMP1 = DSQRT(DBLE(N))*EPSLN DO 6971 p = 1, N XSC = ONE / DNRM2( N, V(1,p), 1 ) IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) & CALL DSCAL( N, XSC, V(1,p), 1 ) 6971 CONTINUE * * Assemble the left singular vector matrix U (M x N). * IF ( N .LT. M ) THEN CALL DLASET( 'A', M-N, N, ZERO, ZERO, U(N+1,1), LDU ) IF ( N .LT. N1 ) THEN CALL DLASET( 'A',N, N1-N, ZERO, ZERO, U(1,N+1),LDU ) CALL DLASET( 'A',M-N,N1-N, ZERO, ONE,U(N+1,N+1),LDU ) END IF END IF CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U, & LDU, WORK(N+1), LWORK-N, IERR ) TEMP1 = DSQRT(DBLE(M))*EPSLN DO 6973 p = 1, N1 XSC = ONE / DNRM2( M, U(1,p), 1 ) IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) & CALL DSCAL( M, XSC, U(1,p), 1 ) 6973 CONTINUE * IF ( ROWPIV ) & CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) * END IF * * end of the >> almost orthogonal case << in the full SVD * ELSE * * This branch deploys a preconditioned Jacobi SVD with explicitly * accumulated rotations. It is included as optional, mainly for * experimental purposes. It does perfom well, and can also be used. * In this implementation, this branch will be automatically activated * if the condition number sigma_max(A) / sigma_min(A) is predicted * to be greater than the overflow threshold. This is because the * a posteriori computation of the singular vectors assumes robust * implementation of BLAS and some LAPACK procedures, capable of working * in presence of extreme values. Since that is not always the case, ... * DO 7968 p = 1, NR CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 ) 7968 CONTINUE * IF ( L2PERT ) THEN XSC = DSQRT(SMALL/EPSLN) DO 5969 q = 1, NR TEMP1 = XSC*DABS( V(q,q) ) DO 5968 p = 1, N IF ( ( p .GT. q ) .AND. ( DABS(V(p,q)) .LE. TEMP1 ) & .OR. ( p .LT. q ) ) & V(p,q) = DSIGN( TEMP1, V(p,q) ) IF ( p. LT. q ) V(p,q) = - V(p,q) 5968 CONTINUE 5969 CONTINUE ELSE CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV ) END IF CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1), & LWORK-2*N, IERR ) CALL DLACPY( 'L', N, NR, V, LDV, WORK(2*N+1), N ) * DO 7969 p = 1, NR CALL DCOPY( NR-p+1, V(p,p), LDV, U(p,p), 1 ) 7969 CONTINUE IF ( L2PERT ) THEN XSC = DSQRT(SMALL/EPSLN) DO 9970 q = 2, NR DO 9971 p = 1, q - 1 TEMP1 = XSC * DMIN1(DABS(U(p,p)),DABS(U(q,q))) U(p,q) = - DSIGN( TEMP1, U(q,p) ) 9971 CONTINUE 9970 CONTINUE ELSE CALL DLASET('U', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU ) END IF CALL DGESVJ( 'G', 'U', 'V', NR, NR, U, LDU, SVA, & N, V, LDV, WORK(2*N+N*NR+1), LWORK-2*N-N*NR, INFO ) SCALEM = WORK(2*N+N*NR+1) NUMRANK = IDNINT(WORK(2*N+N*NR+2)) IF ( NR .LT. N ) THEN CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV ) CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV ) CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV ) END IF CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1), & V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR ) * * Permute the rows of V using the (column) permutation from the * first QRF. Also, scale the columns to make them unit in * Euclidean norm. This applies to all cases. * TEMP1 = DSQRT(DBLE(N)) * EPSLN DO 7972 q = 1, N DO 8972 p = 1, N WORK(2*N+N*NR+NR+IWORK(p)) = V(p,q) 8972 CONTINUE DO 8973 p = 1, N V(p,q) = WORK(2*N+N*NR+NR+p) 8973 CONTINUE XSC = ONE / DNRM2( N, V(1,q), 1 ) IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) & CALL DSCAL( N, XSC, V(1,q), 1 ) 7972 CONTINUE * * At this moment, V contains the right singular vectors of A. * Next, assemble the left singular vector matrix U (M x N). * IF ( NR .LT. M ) THEN CALL DLASET( 'A', M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU ) IF ( NR .LT. N1 ) THEN CALL DLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1),LDU ) CALL DLASET( 'A',M-NR,N1-NR, ZERO, ONE,U(NR+1,NR+1),LDU ) END IF END IF * CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U, & LDU, WORK(N+1), LWORK-N, IERR ) * IF ( ROWPIV ) & CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) * * END IF IF ( TRANSP ) THEN * .. swap U and V because the procedure worked on A^t DO 6974 p = 1, N CALL DSWAP( N, U(1,p), 1, V(1,p), 1 ) 6974 CONTINUE END IF * END IF * end of the full SVD * * Undo scaling, if necessary (and possible) * IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN CALL DLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR ) USCAL1 = ONE USCAL2 = ONE END IF * IF ( NR .LT. N ) THEN DO 3004 p = NR+1, N SVA(p) = ZERO 3004 CONTINUE END IF * WORK(1) = USCAL2 * SCALEM WORK(2) = USCAL1 IF ( ERREST ) WORK(3) = SCONDA IF ( LSVEC .AND. RSVEC ) THEN WORK(4) = CONDR1 WORK(5) = CONDR2 END IF IF ( L2TRAN ) THEN WORK(6) = ENTRA WORK(7) = ENTRAT END IF * IWORK(1) = NR IWORK(2) = NUMRANK IWORK(3) = WARNING * RETURN * .. * .. END OF DGEJSV * .. END *go to the page top

USAGE: tau, info, a = NumRu::Lapack.dgelq2( a, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGELQ2( M, N, A, LDA, TAU, WORK, INFO ) * Purpose * ======= * * DGELQ2 computes an LQ factorization of a real m by n matrix A: * A = L * Q. * * 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 elements on and below the diagonal of the array * contain the m by min(m,n) lower trapezoidal matrix L (L is * lower triangular if m <= n); the elements above the diagonal, * with the array TAU, represent the orthogonal matrix Q as a * product of elementary reflectors (see Further Details). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) * The scalar factors of the elementary reflectors (see Further * Details). * * WORK (workspace) DOUBLE PRECISION array, dimension (M) * * 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(k) . . . H(2) H(1), where k = min(m,n). * * Each H(i) has the form * * H(i) = I - tau * v * v' * * where tau is a real scalar, and v is a real vector with * v(1:i-1) = 0 and v(i) = 1; v(i+1:n) is stored on exit in A(i,i+1:n), * and tau in TAU(i). * * ===================================================================== *go to the page top

USAGE: tau, work, info, a = NumRu::Lapack.dgelqf( m, a, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGELQF( M, N, A, LDA, TAU, WORK, LWORK, INFO ) * Purpose * ======= * * DGELQF computes an LQ factorization of a real M-by-N matrix A: * A = L * Q. * * 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 elements on and below the diagonal of the array * contain the m-by-min(m,n) lower trapezoidal matrix L (L is * lower triangular if m <= n); the elements above the diagonal, * with the array TAU, represent the orthogonal matrix Q as a * product of elementary reflectors (see Further Details). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) * The scalar factors of the elementary reflectors (see Further * Details). * * 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 >= max(1,M). * For optimum performance LWORK >= M*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(k) . . . H(2) H(1), where k = min(m,n). * * Each H(i) has the form * * H(i) = I - tau * v * v' * * where tau is a real scalar, and v is a real vector with * v(1:i-1) = 0 and v(i) = 1; v(i+1:n) is stored on exit in A(i,i+1:n), * and tau in TAU(i). * * ===================================================================== * * .. Local Scalars .. LOGICAL LQUERY INTEGER I, IB, IINFO, IWS, K, LDWORK, LWKOPT, NB, $ NBMIN, NX * .. * .. External Subroutines .. EXTERNAL DGELQ2, DLARFB, DLARFT, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. External Functions .. INTEGER ILAENV EXTERNAL ILAENV * ..go to the page top

USAGE: work, info, a, b = NumRu::Lapack.dgels( trans, a, b, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO ) * Purpose * ======= * * DGELS solves overdetermined or underdetermined real linear systems * involving an M-by-N matrix A, or its transpose, using a QR or LQ * factorization of A. It is assumed that A has full rank. * * The following options are provided: * * 1. If TRANS = 'N' and m >= n: find the least squares solution of * an overdetermined system, i.e., solve the least squares problem * minimize || B - A*X ||. * * 2. If TRANS = 'N' and m < n: find the minimum norm solution of * an underdetermined system A * X = B. * * 3. If TRANS = 'T' and m >= n: find the minimum norm solution of * an undetermined system A**T * X = B. * * 4. If TRANS = 'T' and m < n: find the least squares solution of * an overdetermined system, i.e., solve the least squares problem * minimize || B - A**T * X ||. * * Several right hand side vectors b and solution vectors x can be * handled in a single call; they are stored as the columns of the * M-by-NRHS right hand side matrix B and the N-by-NRHS solution * matrix X. * * Arguments * ========= * * TRANS (input) CHARACTER*1 * = 'N': the linear system involves A; * = 'T': the linear system involves A**T. * * 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. * * NRHS (input) INTEGER * The number of right hand sides, i.e., the number of * columns of the matrices B and X. NRHS >=0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the M-by-N matrix A. * On exit, * if M >= N, A is overwritten by details of its QR * factorization as returned by DGEQRF; * if M < N, A is overwritten by details of its LQ * factorization as returned by DGELQF. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) * On entry, the matrix B of right hand side vectors, stored * columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS * if TRANS = 'T'. * On exit, if INFO = 0, B is overwritten by the solution * vectors, stored columnwise: * if TRANS = 'N' and m >= n, rows 1 to n of B contain the least * squares solution vectors; the residual sum of squares for the * solution in each column is given by the sum of squares of * elements N+1 to M in that column; * if TRANS = 'N' and m < n, rows 1 to N of B contain the * minimum norm solution vectors; * if TRANS = 'T' and m >= n, rows 1 to M of B contain the * minimum norm solution vectors; * if TRANS = 'T' and m < n, rows 1 to M of B contain the * least squares solution vectors; the residual sum of squares * for the solution in each column is given by the sum of * squares of elements M+1 to N in that column. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= MAX(1,M,N). * * 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 >= max( 1, MN + max( MN, NRHS ) ). * For optimal performance, * LWORK >= max( 1, MN + max( MN, NRHS )*NB ). * where MN = min(M,N) and NB is the optimum block size. * * 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 * > 0: if INFO = i, the i-th diagonal element of the * triangular factor of A is zero, so that A does not have * full rank; the least squares solution could not be * computed. * * ===================================================================== *go to the page top

USAGE: s, rank, work, info, b = NumRu::Lapack.dgelsd( a, b, rcond, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, IWORK, INFO ) * Purpose * ======= * * DGELSD computes the minimum-norm solution to a real linear least * squares problem: * minimize 2-norm(| b - A*x |) * using the singular value decomposition (SVD) of A. A is an M-by-N * matrix which may be rank-deficient. * * Several right hand side vectors b and solution vectors x can be * handled in a single call; they are stored as the columns of the * M-by-NRHS right hand side matrix B and the N-by-NRHS solution * matrix X. * * The problem is solved in three steps: * (1) Reduce the coefficient matrix A to bidiagonal form with * Householder transformations, reducing the original problem * into a "bidiagonal least squares problem" (BLS) * (2) Solve the BLS using a divide and conquer approach. * (3) Apply back all the Householder tranformations to solve * the original least squares problem. * * The effective rank of A is determined by treating as zero those * singular values which are less than RCOND times the largest singular * value. * * The divide and conquer algorithm makes very mild assumptions about * floating point arithmetic. It will work on machines with a guard * digit in add/subtract, or on those binary machines without guard * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or * Cray-2. It could conceivably fail on hexadecimal or decimal machines * without guard digits, but we know of none. * * Arguments * ========= * * M (input) INTEGER * The number of rows of A. M >= 0. * * N (input) INTEGER * The number of columns of A. N >= 0. * * NRHS (input) INTEGER * The number of right hand sides, i.e., the number of columns * of the matrices B and X. NRHS >= 0. * * A (input) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the M-by-N matrix A. * On exit, A has been destroyed. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) * On entry, the M-by-NRHS right hand side matrix B. * On exit, B is overwritten by the N-by-NRHS solution * matrix X. If m >= n and RANK = n, the residual * sum-of-squares for the solution in the i-th column is given * by the sum of squares of elements n+1:m in that column. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,max(M,N)). * * S (output) DOUBLE PRECISION array, dimension (min(M,N)) * The singular values of A in decreasing order. * The condition number of A in the 2-norm = S(1)/S(min(m,n)). * * RCOND (input) DOUBLE PRECISION * RCOND is used to determine the effective rank of A. * Singular values S(i) <= RCOND*S(1) are treated as zero. * If RCOND < 0, machine precision is used instead. * * RANK (output) INTEGER * The effective rank of A, i.e., the number of singular values * which are greater than RCOND*S(1). * * 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 must be at least 1. * The exact minimum amount of workspace needed depends on M, * N and NRHS. As long as LWORK is at least * 12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2, * if M is greater than or equal to N or * 12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + (SMLSIZ+1)**2, * if M is less than N, the code will execute correctly. * SMLSIZ is returned by ILAENV and is equal to the maximum * size of the subproblems at the bottom of the computation * tree (usually about 25), and * NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 ) * For good performance, LWORK should 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. * * IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK)) * LIWORK >= max(1, 3 * MINMN * NLVL + 11 * MINMN), * where MINMN = MIN( M,N ). * On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value. * > 0: the algorithm for computing the SVD failed to converge; * if INFO = i, i off-diagonal elements of an intermediate * bidiagonal form did not converge to zero. * * Further Details * =============== * * Based on contributions by * Ming Gu and Ren-Cang Li, Computer Science Division, University of * California at Berkeley, USA * Osni Marques, LBNL/NERSC, USA * * ===================================================================== *go to the page top

USAGE: s, rank, work, info, a, b = NumRu::Lapack.dgelss( a, b, rcond, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO ) * Purpose * ======= * * DGELSS computes the minimum norm solution to a real linear least * squares problem: * * Minimize 2-norm(| b - A*x |). * * using the singular value decomposition (SVD) of A. A is an M-by-N * matrix which may be rank-deficient. * * Several right hand side vectors b and solution vectors x can be * handled in a single call; they are stored as the columns of the * M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix * X. * * The effective rank of A is determined by treating as zero those * singular values which are less than RCOND times the largest singular * value. * * 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. * * NRHS (input) INTEGER * The number of right hand sides, i.e., the number of columns * of the matrices B and X. NRHS >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the M-by-N matrix A. * On exit, the first min(m,n) rows of A are overwritten with * its right singular vectors, stored rowwise. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) * On entry, the M-by-NRHS right hand side matrix B. * On exit, B is overwritten by the N-by-NRHS solution * matrix X. If m >= n and RANK = n, the residual * sum-of-squares for the solution in the i-th column is given * by the sum of squares of elements n+1:m in that column. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,max(M,N)). * * S (output) DOUBLE PRECISION array, dimension (min(M,N)) * The singular values of A in decreasing order. * The condition number of A in the 2-norm = S(1)/S(min(m,n)). * * RCOND (input) DOUBLE PRECISION * RCOND is used to determine the effective rank of A. * Singular values S(i) <= RCOND*S(1) are treated as zero. * If RCOND < 0, machine precision is used instead. * * RANK (output) INTEGER * The effective rank of A, i.e., the number of singular values * which are greater than RCOND*S(1). * * 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 >= 1, and also: * LWORK >= 3*min(M,N) + max( 2*min(M,N), max(M,N), NRHS ) * For good performance, LWORK should 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. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value. * > 0: the algorithm for computing the SVD failed to converge; * if INFO = i, i off-diagonal elements of an intermediate * bidiagonal form did not converge to zero. * * ===================================================================== *go to the page top

USAGE: rank, info, a, b, jpvt = NumRu::Lapack.dgelsx( m, a, b, jpvt, rcond, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, WORK, INFO ) * Purpose * ======= * * This routine is deprecated and has been replaced by routine DGELSY. * * DGELSX computes the minimum-norm solution to a real linear least * squares problem: * minimize || A * X - B || * using a complete orthogonal factorization of A. A is an M-by-N * matrix which may be rank-deficient. * * Several right hand side vectors b and solution vectors x can be * handled in a single call; they are stored as the columns of the * M-by-NRHS right hand side matrix B and the N-by-NRHS solution * matrix X. * * The routine first computes a QR factorization with column pivoting: * A * P = Q * [ R11 R12 ] * [ 0 R22 ] * with R11 defined as the largest leading submatrix whose estimated * condition number is less than 1/RCOND. The order of R11, RANK, * is the effective rank of A. * * Then, R22 is considered to be negligible, and R12 is annihilated * by orthogonal transformations from the right, arriving at the * complete orthogonal factorization: * A * P = Q * [ T11 0 ] * Z * [ 0 0 ] * The minimum-norm solution is then * X = P * Z' [ inv(T11)*Q1'*B ] * [ 0 ] * where Q1 consists of the first RANK columns of Q. * * 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. * * NRHS (input) INTEGER * The number of right hand sides, i.e., the number of * columns of matrices B and X. NRHS >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the M-by-N matrix A. * On exit, A has been overwritten by details of its * complete orthogonal factorization. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) * On entry, the M-by-NRHS right hand side matrix B. * On exit, the N-by-NRHS solution matrix X. * If m >= n and RANK = n, the residual sum-of-squares for * the solution in the i-th column is given by the sum of * squares of elements N+1:M in that column. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,M,N). * * JPVT (input/output) INTEGER array, dimension (N) * On entry, if JPVT(i) .ne. 0, the i-th column of A is an * initial column, otherwise it is a free column. Before * the QR factorization of A, all initial columns are * permuted to the leading positions; only the remaining * free columns are moved as a result of column pivoting * during the factorization. * On exit, if JPVT(i) = k, then the i-th column of A*P * was the k-th column of A. * * RCOND (input) DOUBLE PRECISION * RCOND is used to determine the effective rank of A, which * is defined as the order of the largest leading triangular * submatrix R11 in the QR factorization with pivoting of A, * whose estimated condition number < 1/RCOND. * * RANK (output) INTEGER * The effective rank of A, i.e., the order of the submatrix * R11. This is the same as the order of the submatrix T11 * in the complete orthogonal factorization of A. * * WORK (workspace) DOUBLE PRECISION array, dimension * (max( min(M,N)+3*N, 2*min(M,N)+NRHS )), * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * * ===================================================================== *go to the page top

USAGE: rank, work, info, a, b, jpvt = NumRu::Lapack.dgelsy( a, b, jpvt, rcond, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, WORK, LWORK, INFO ) * Purpose * ======= * * DGELSY computes the minimum-norm solution to a real linear least * squares problem: * minimize || A * X - B || * using a complete orthogonal factorization of A. A is an M-by-N * matrix which may be rank-deficient. * * Several right hand side vectors b and solution vectors x can be * handled in a single call; they are stored as the columns of the * M-by-NRHS right hand side matrix B and the N-by-NRHS solution * matrix X. * * The routine first computes a QR factorization with column pivoting: * A * P = Q * [ R11 R12 ] * [ 0 R22 ] * with R11 defined as the largest leading submatrix whose estimated * condition number is less than 1/RCOND. The order of R11, RANK, * is the effective rank of A. * * Then, R22 is considered to be negligible, and R12 is annihilated * by orthogonal transformations from the right, arriving at the * complete orthogonal factorization: * A * P = Q * [ T11 0 ] * Z * [ 0 0 ] * The minimum-norm solution is then * X = P * Z' [ inv(T11)*Q1'*B ] * [ 0 ] * where Q1 consists of the first RANK columns of Q. * * This routine is basically identical to the original xGELSX except * three differences: * o The call to the subroutine xGEQPF has been substituted by the * the call to the subroutine xGEQP3. This subroutine is a Blas-3 * version of the QR factorization with column pivoting. * o Matrix B (the right hand side) is updated with Blas-3. * o The permutation of matrix B (the right hand side) is faster and * more simple. * * 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. * * NRHS (input) INTEGER * The number of right hand sides, i.e., the number of * columns of matrices B and X. NRHS >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the M-by-N matrix A. * On exit, A has been overwritten by details of its * complete orthogonal factorization. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) * On entry, the M-by-NRHS right hand side matrix B. * On exit, the N-by-NRHS solution matrix X. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,M,N). * * JPVT (input/output) INTEGER array, dimension (N) * On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted * to the front of AP, otherwise column i is a free column. * On exit, if JPVT(i) = k, then the i-th column of AP * was the k-th column of A. * * RCOND (input) DOUBLE PRECISION * RCOND is used to determine the effective rank of A, which * is defined as the order of the largest leading triangular * submatrix R11 in the QR factorization with pivoting of A, * whose estimated condition number < 1/RCOND. * * RANK (output) INTEGER * The effective rank of A, i.e., the order of the submatrix * R11. This is the same as the order of the submatrix T11 * in the complete orthogonal factorization of A. * * 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. * The unblocked strategy requires that: * LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ), * where MN = min( M, N ). * The block algorithm requires that: * LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ), * where NB is an upper bound on the blocksize returned * by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR, * and DORMRZ. * * 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 * =============== * * Based on contributions by * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA * E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain * * ===================================================================== *go to the page top

USAGE: tau, info, a = NumRu::Lapack.dgeql2( m, a, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEQL2( M, N, A, LDA, TAU, WORK, INFO ) * Purpose * ======= * * DGEQL2 computes a QL factorization of a real m by n matrix A: * A = Q * L. * * 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, if m >= n, the lower triangle of the subarray * A(m-n+1:m,1:n) contains the n by n lower triangular matrix L; * if m <= n, the elements on and below the (n-m)-th * superdiagonal contain the m by n lower trapezoidal matrix L; * the remaining elements, with the array TAU, represent the * orthogonal matrix Q as a product of elementary reflectors * (see Further Details). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) * The scalar factors of the elementary reflectors (see Further * Details). * * WORK (workspace) DOUBLE PRECISION array, dimension (N) * * 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(k) . . . H(2) H(1), where k = min(m,n). * * Each H(i) has the form * * H(i) = I - tau * v * v' * * where tau is a real scalar, and v is a real vector with * v(m-k+i+1:m) = 0 and v(m-k+i) = 1; v(1:m-k+i-1) is stored on exit in * A(1:m-k+i-1,n-k+i), and tau in TAU(i). * * ===================================================================== *go to the page top

USAGE: tau, work, info, a = NumRu::Lapack.dgeqlf( m, a, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEQLF( M, N, A, LDA, TAU, WORK, LWORK, INFO ) * Purpose * ======= * * DGEQLF computes a QL factorization of a real M-by-N matrix A: * A = Q * L. * * 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, * if m >= n, the lower triangle of the subarray * A(m-n+1:m,1:n) contains the N-by-N lower triangular matrix L; * if m <= n, the elements on and below the (n-m)-th * superdiagonal contain the M-by-N lower trapezoidal matrix L; * the remaining elements, with the array TAU, represent the * orthogonal matrix Q as a product of elementary reflectors * (see Further Details). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) * The scalar factors of the elementary reflectors (see Further * Details). * * 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 >= max(1,N). * For optimum performance LWORK >= N*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(k) . . . H(2) H(1), where k = min(m,n). * * Each H(i) has the form * * H(i) = I - tau * v * v' * * where tau is a real scalar, and v is a real vector with * v(m-k+i+1:m) = 0 and v(m-k+i) = 1; v(1:m-k+i-1) is stored on exit in * A(1:m-k+i-1,n-k+i), and tau in TAU(i). * * ===================================================================== * * .. Local Scalars .. LOGICAL LQUERY INTEGER I, IB, IINFO, IWS, K, KI, KK, LDWORK, LWKOPT, $ MU, NB, NBMIN, NU, NX * .. * .. External Subroutines .. EXTERNAL DGEQL2, DLARFB, DLARFT, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. External Functions .. INTEGER ILAENV EXTERNAL ILAENV * ..go to the page top

USAGE: tau, work, info, a, jpvt = NumRu::Lapack.dgeqp3( m, a, jpvt, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO ) * 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 * * ===================================================================== *go to the page top

USAGE: tau, info, a, jpvt = NumRu::Lapack.dgeqpf( m, a, jpvt, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEQPF( M, N, A, LDA, JPVT, TAU, WORK, INFO ) * Purpose * ======= * * This routine is deprecated and has been replaced by routine DGEQP3. * * DGEQPF computes a QR factorization with column pivoting of a * real M-by-N matrix A: A*P = Q*R. * * 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 triangular 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(i) .ne. 0, the i-th column of A is permuted * to the front of A*P (a leading column); if JPVT(i) = 0, * the i-th column of A is a free column. * On exit, if JPVT(i) = k, then the i-th column of A*P * was the k-th column of A. * * TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) * The scalar factors of the elementary reflectors. * * WORK (workspace) DOUBLE PRECISION array, dimension (3*N) * * 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(n) * * Each H(i) has the form * * H = I - tau * v * v' * * where tau is a real scalar, and v is a real 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). * * The matrix P is represented in jpvt as follows: If * jpvt(j) = i * then the jth column of P is the ith canonical unit vector. * * Partial column norm updating strategy modified by * Z. Drmac and Z. Bujanovic, Dept. of Mathematics, * University of Zagreb, Croatia. * June 2010 * For more details see LAPACK Working Note 176. * * ===================================================================== *go to the page top

USAGE: tau, info, a = NumRu::Lapack.dgeqr2( m, a, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEQR2( M, N, A, LDA, TAU, WORK, INFO ) * Purpose * ======= * * DGEQR2 computes a QR factorization of a real m by n matrix A: * A = Q * R. * * 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 elements on and above the diagonal of the array * contain the min(m,n) by n upper trapezoidal matrix R (R is * upper triangular if m >= n); the elements below the diagonal, * with the array TAU, represent the orthogonal matrix Q as a * product of elementary reflectors (see Further Details). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) * The scalar factors of the elementary reflectors (see Further * Details). * * WORK (workspace) DOUBLE PRECISION array, dimension (N) * * 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 scalar, and v is a real 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). * * ===================================================================== *go to the page top

USAGE: tau, info, a = NumRu::Lapack.dgeqr2p( m, a, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEQR2P( M, N, A, LDA, TAU, WORK, INFO ) * Purpose * ======= * * DGEQR2 computes a QR factorization of a real m by n matrix A: * A = Q * R. * * 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 elements on and above the diagonal of the array * contain the min(m,n) by n upper trapezoidal matrix R (R is * upper triangular if m >= n); the elements below the diagonal, * with the array TAU, represent the orthogonal matrix Q as a * product of elementary reflectors (see Further Details). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) * The scalar factors of the elementary reflectors (see Further * Details). * * WORK (workspace) DOUBLE PRECISION array, dimension (N) * * 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 scalar, and v is a real 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). * * ===================================================================== *go to the page top

USAGE: tau, work, info, a = NumRu::Lapack.dgeqrf( m, a, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEQRF( M, N, A, LDA, TAU, WORK, LWORK, INFO ) * Purpose * ======= * * DGEQRF computes a QR factorization of a real M-by-N matrix A: * A = Q * R. * * 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 elements on and above the diagonal of the array * contain the min(M,N)-by-N upper trapezoidal matrix R (R is * upper triangular if m >= n); the elements below the diagonal, * with the array TAU, represent the orthogonal matrix Q as a * product of min(m,n) elementary reflectors (see Further * Details). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) * The scalar factors of the elementary reflectors (see Further * Details). * * 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 >= max(1,N). * For optimum performance LWORK >= N*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 scalar, and v is a real 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). * * ===================================================================== * * .. Local Scalars .. LOGICAL LQUERY INTEGER I, IB, IINFO, IWS, K, LDWORK, LWKOPT, NB, $ NBMIN, NX * .. * .. External Subroutines .. EXTERNAL DGEQR2, DLARFB, DLARFT, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. External Functions .. INTEGER ILAENV EXTERNAL ILAENV * ..go to the page top

USAGE: tau, work, info, a = NumRu::Lapack.dgeqrfp( m, a, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGEQRFP( M, N, A, LDA, TAU, WORK, LWORK, INFO ) * Purpose * ======= * * DGEQRFP computes a QR factorization of a real M-by-N matrix A: * A = Q * R. * * 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 elements on and above the diagonal of the array * contain the min(M,N)-by-N upper trapezoidal matrix R (R is * upper triangular if m >= n); the elements below the diagonal, * with the array TAU, represent the orthogonal matrix Q as a * product of min(m,n) elementary reflectors (see Further * Details). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) * The scalar factors of the elementary reflectors (see Further * Details). * * 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 >= max(1,N). * For optimum performance LWORK >= N*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 scalar, and v is a real 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). * * ===================================================================== * * .. Local Scalars .. LOGICAL LQUERY INTEGER I, IB, IINFO, IWS, K, LDWORK, LWKOPT, NB, $ NBMIN, NX * .. * .. External Subroutines .. EXTERNAL DGEQR2P, DLARFB, DLARFT, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. External Functions .. INTEGER ILAENV EXTERNAL ILAENV * ..go to the page top

USAGE: ferr, berr, info, x = NumRu::Lapack.dgerfs( trans, a, af, ipiv, b, x, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGERFS( TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO ) * Purpose * ======= * * DGERFS improves the computed solution to a system of linear * equations and provides error bounds and backward error estimates for * the solution. * * Arguments * ========= * * TRANS (input) CHARACTER*1 * Specifies the form of the system of equations: * = 'N': A * X = B (No transpose) * = 'T': A**T * X = B (Transpose) * = 'C': A**H * X = B (Conjugate transpose = Transpose) * * N (input) INTEGER * The order of the matrix A. N >= 0. * * NRHS (input) INTEGER * The number of right hand sides, i.e., the number of columns * of the matrices B and X. NRHS >= 0. * * A (input) DOUBLE PRECISION array, dimension (LDA,N) * The original N-by-N matrix A. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * AF (input) DOUBLE PRECISION array, dimension (LDAF,N) * The factors L and U from the factorization A = P*L*U * as computed by DGETRF. * * LDAF (input) INTEGER * The leading dimension of the array AF. LDAF >= max(1,N). * * IPIV (input) INTEGER array, dimension (N) * The pivot indices from DGETRF; for 1<=i<=N, row i of the * matrix was interchanged with row IPIV(i). * * B (input) DOUBLE PRECISION array, dimension (LDB,NRHS) * The right hand side matrix B. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,N). * * X (input/output) DOUBLE PRECISION array, dimension (LDX,NRHS) * On entry, the solution matrix X, as computed by DGETRS. * On exit, the improved solution matrix X. * * LDX (input) INTEGER * The leading dimension of the array X. LDX >= max(1,N). * * FERR (output) DOUBLE PRECISION array, dimension (NRHS) * The estimated forward error bound for each solution vector * X(j) (the j-th column of the solution matrix X). * If XTRUE is the true solution corresponding to X(j), FERR(j) * is an estimated upper bound for the magnitude of the largest * element in (X(j) - XTRUE) divided by the magnitude of the * largest element in X(j). The estimate is as reliable as * the estimate for RCOND, and is almost always a slight * overestimate of the true error. * * BERR (output) DOUBLE PRECISION array, dimension (NRHS) * The componentwise relative backward error of each solution * vector X(j) (i.e., the smallest relative change in * any element of A or B that makes X(j) an exact solution). * * WORK (workspace) DOUBLE PRECISION array, dimension (3*N) * * IWORK (workspace) INTEGER array, dimension (N) * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * * Internal Parameters * =================== * * ITMAX is the maximum number of steps of iterative refinement. * * ===================================================================== *go to the page top

USAGE: rcond, berr, err_bnds_norm, err_bnds_comp, info, x, params = NumRu::Lapack.dgerfsx( trans, equed, a, af, ipiv, r, c, b, x, params, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGERFSX( TRANS, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV, R, C, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, IWORK, INFO ) * Purpose * ======= * * DGERFSX improves the computed solution to a system of linear * equations and provides error bounds and backward error estimates * for the solution. In addition to normwise error bound, the code * provides maximum componentwise error bound if possible. See * comments for ERR_BNDS_NORM and ERR_BNDS_COMP for details of the * error bounds. * * The original system of linear equations may have been equilibrated * before calling this routine, as described by arguments EQUED, R * and C below. In this case, the solution and error bounds returned * are for the original unequilibrated system. * * Arguments * ========= * * Some optional parameters are bundled in the PARAMS array. These * settings determine how refinement is performed, but often the * defaults are acceptable. If the defaults are acceptable, users * can pass NPARAMS = 0 which prevents the source code from accessing * the PARAMS argument. * * TRANS (input) CHARACTER*1 * Specifies the form of the system of equations: * = 'N': A * X = B (No transpose) * = 'T': A**T * X = B (Transpose) * = 'C': A**H * X = B (Conjugate transpose = Transpose) * * EQUED (input) CHARACTER*1 * Specifies the form of equilibration that was done to A * before calling this routine. This is needed to compute * the solution and error bounds correctly. * = 'N': No equilibration * = 'R': Row equilibration, i.e., A has been premultiplied by * diag(R). * = 'C': Column equilibration, i.e., A has been postmultiplied * by diag(C). * = 'B': Both row and column equilibration, i.e., A has been * replaced by diag(R) * A * diag(C). * The right hand side B has been changed accordingly. * * N (input) INTEGER * The order of the matrix A. N >= 0. * * NRHS (input) INTEGER * The number of right hand sides, i.e., the number of columns * of the matrices B and X. NRHS >= 0. * * A (input) DOUBLE PRECISION array, dimension (LDA,N) * The original N-by-N matrix A. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * AF (input) DOUBLE PRECISION array, dimension (LDAF,N) * The factors L and U from the factorization A = P*L*U * as computed by DGETRF. * * LDAF (input) INTEGER * The leading dimension of the array AF. LDAF >= max(1,N). * * IPIV (input) INTEGER array, dimension (N) * The pivot indices from DGETRF; for 1<=i<=N, row i of the * matrix was interchanged with row IPIV(i). * * R (input) DOUBLE PRECISION array, dimension (N) * The row scale factors for A. If EQUED = 'R' or 'B', A is * multiplied on the left by diag(R); if EQUED = 'N' or 'C', R * is not accessed. * If R is accessed, each element of R should be a power of the radix * to ensure a reliable solution and error estimates. Scaling by * powers of the radix does not cause rounding errors unless the * result underflows or overflows. Rounding errors during scaling * lead to refining with a matrix that is not equivalent to the * input matrix, producing error estimates that may not be * reliable. * * C (input) DOUBLE PRECISION array, dimension (N) * The column scale factors for A. If EQUED = 'C' or 'B', A is * multiplied on the right by diag(C); if EQUED = 'N' or 'R', C * is not accessed. * If C is accessed, each element of C should be a power of the radix * to ensure a reliable solution and error estimates. Scaling by * powers of the radix does not cause rounding errors unless the * result underflows or overflows. Rounding errors during scaling * lead to refining with a matrix that is not equivalent to the * input matrix, producing error estimates that may not be * reliable. * * B (input) DOUBLE PRECISION array, dimension (LDB,NRHS) * The right hand side matrix B. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,N). * * X (input/output) DOUBLE PRECISION array, dimension (LDX,NRHS) * On entry, the solution matrix X, as computed by DGETRS. * On exit, the improved solution matrix X. * * LDX (input) INTEGER * The leading dimension of the array X. LDX >= max(1,N). * * RCOND (output) DOUBLE PRECISION * Reciprocal scaled condition number. This is an estimate of the * reciprocal Skeel condition number of the matrix A after * equilibration (if done). If this is less than the machine * precision (in particular, if it is zero), the matrix is singular * to working precision. Note that the error may still be small even * if this number is very small and the matrix appears ill- * conditioned. * * BERR (output) DOUBLE PRECISION array, dimension (NRHS) * Componentwise relative backward error. This is the * componentwise relative backward error of each solution vector X(j) * (i.e., the smallest relative change in any element of A or B that * makes X(j) an exact solution). * * N_ERR_BNDS (input) INTEGER * Number of error bounds to return for each right hand side * and each type (normwise or componentwise). See ERR_BNDS_NORM and * ERR_BNDS_COMP below. * * ERR_BNDS_NORM (output) DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS) * For each right-hand side, this array contains information about * various error bounds and condition numbers corresponding to the * normwise relative error, which is defined as follows: * * Normwise relative error in the ith solution vector: * max_j (abs(XTRUE(j,i) - X(j,i))) * ------------------------------ * max_j abs(X(j,i)) * * The array is indexed by the type of error information as described * below. There currently are up to three pieces of information * returned. * * The first index in ERR_BNDS_NORM(i,:) corresponds to the ith * right-hand side. * * The second index in ERR_BNDS_NORM(:,err) contains the following * three fields: * err = 1 "Trust/don't trust" boolean. Trust the answer if the * reciprocal condition number is less than the threshold * sqrt(n) * dlamch('Epsilon'). * * err = 2 "Guaranteed" error bound: The estimated forward error, * almost certainly within a factor of 10 of the true error * so long as the next entry is greater than the threshold * sqrt(n) * dlamch('Epsilon'). This error bound should only * be trusted if the previous boolean is true. * * err = 3 Reciprocal condition number: Estimated normwise * reciprocal condition number. Compared with the threshold * sqrt(n) * dlamch('Epsilon') to determine if the error * estimate is "guaranteed". These reciprocal condition * numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some * appropriately scaled matrix Z. * Let Z = S*A, where S scales each row by a power of the * radix so all absolute row sums of Z are approximately 1. * * See Lapack Working Note 165 for further details and extra * cautions. * * ERR_BNDS_COMP (output) DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS) * For each right-hand side, this array contains information about * various error bounds and condition numbers corresponding to the * componentwise relative error, which is defined as follows: * * Componentwise relative error in the ith solution vector: * abs(XTRUE(j,i) - X(j,i)) * max_j ---------------------- * abs(X(j,i)) * * The array is indexed by the right-hand side i (on which the * componentwise relative error depends), and the type of error * information as described below. There currently are up to three * pieces of information returned for each right-hand side. If * componentwise accuracy is not requested (PARAMS(3) = 0.0), then * ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most * the first (:,N_ERR_BNDS) entries are returned. * * The first index in ERR_BNDS_COMP(i,:) corresponds to the ith * right-hand side. * * The second index in ERR_BNDS_COMP(:,err) contains the following * three fields: * err = 1 "Trust/don't trust" boolean. Trust the answer if the * reciprocal condition number is less than the threshold * sqrt(n) * dlamch('Epsilon'). * * err = 2 "Guaranteed" error bound: The estimated forward error, * almost certainly within a factor of 10 of the true error * so long as the next entry is greater than the threshold * sqrt(n) * dlamch('Epsilon'). This error bound should only * be trusted if the previous boolean is true. * * err = 3 Reciprocal condition number: Estimated componentwise * reciprocal condition number. Compared with the threshold * sqrt(n) * dlamch('Epsilon') to determine if the error * estimate is "guaranteed". These reciprocal condition * numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some * appropriately scaled matrix Z. * Let Z = S*(A*diag(x)), where x is the solution for the * current right-hand side and S scales each row of * A*diag(x) by a power of the radix so all absolute row * sums of Z are approximately 1. * * See Lapack Working Note 165 for further details and extra * cautions. * * NPARAMS (input) INTEGER * Specifies the number of parameters set in PARAMS. If .LE. 0, the * PARAMS array is never referenced and default values are used. * * PARAMS (input / output) DOUBLE PRECISION array, dimension (NPARAMS) * Specifies algorithm parameters. If an entry is .LT. 0.0, then * that entry will be filled with default value used for that * parameter. Only positions up to NPARAMS are accessed; defaults * are used for higher-numbered parameters. * * PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative * refinement or not. * Default: 1.0D+0 * = 0.0 : No refinement is performed, and no error bounds are * computed. * = 1.0 : Use the double-precision refinement algorithm, * possibly with doubled-single computations if the * compilation environment does not support DOUBLE * PRECISION. * (other values are reserved for future use) * * PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual * computations allowed for refinement. * Default: 10 * Aggressive: Set to 100 to permit convergence using approximate * factorizations or factorizations other than LU. If * the factorization uses a technique other than * Gaussian elimination, the guarantees in * err_bnds_norm and err_bnds_comp may no longer be * trustworthy. * * PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code * will attempt to find a solution with small componentwise * relative error in the double-precision algorithm. Positive * is true, 0.0 is false. * Default: 1.0 (attempt componentwise convergence) * * WORK (workspace) DOUBLE PRECISION array, dimension (4*N) * * IWORK (workspace) INTEGER array, dimension (N) * * INFO (output) INTEGER * = 0: Successful exit. The solution to every right-hand side is * guaranteed. * < 0: If INFO = -i, the i-th argument had an illegal value * > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization * has been completed, but the factor U is exactly singular, so * the solution and error bounds could not be computed. RCOND = 0 * is returned. * = N+J: The solution corresponding to the Jth right-hand side is * not guaranteed. The solutions corresponding to other right- * hand sides K with K > J may not be guaranteed as well, but * only the first such right-hand side is reported. If a small * componentwise error is not requested (PARAMS(3) = 0.0) then * the Jth right-hand side is the first with a normwise error * bound that is not guaranteed (the smallest J such * that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0) * the Jth right-hand side is the first with either a normwise or * componentwise error bound that is not guaranteed (the smallest * J such that either ERR_BNDS_NORM(J,1) = 0.0 or * ERR_BNDS_COMP(J,1) = 0.0). See the definition of * ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information * about all of the right-hand sides check ERR_BNDS_NORM or * ERR_BNDS_COMP. * * ================================================================== *go to the page top

USAGE: tau, info, a = NumRu::Lapack.dgerq2( a, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGERQ2( M, N, A, LDA, TAU, WORK, INFO ) * Purpose * ======= * * DGERQ2 computes an RQ factorization of a real m by n matrix A: * A = R * Q. * * 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, if m <= n, the upper triangle of the subarray * A(1:m,n-m+1:n) contains the m by m upper triangular matrix R; * if m >= n, the elements on and above the (m-n)-th subdiagonal * contain the m by n upper trapezoidal matrix R; the remaining * elements, with the array TAU, represent the orthogonal matrix * Q as a product of elementary reflectors (see Further * Details). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) * The scalar factors of the elementary reflectors (see Further * Details). * * WORK (workspace) DOUBLE PRECISION array, dimension (M) * * 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 scalar, and v is a real vector with * v(n-k+i+1:n) = 0 and v(n-k+i) = 1; v(1:n-k+i-1) is stored on exit in * A(m-k+i,1:n-k+i-1), and tau in TAU(i). * * ===================================================================== *go to the page top

USAGE: tau, work, info, a = NumRu::Lapack.dgerqf( m, a, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGERQF( M, N, A, LDA, TAU, WORK, LWORK, INFO ) * Purpose * ======= * * DGERQF computes an RQ factorization of a real M-by-N matrix A: * A = R * Q. * * 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, * if m <= n, the upper triangle of the subarray * A(1:m,n-m+1:n) contains the M-by-M upper triangular matrix R; * if m >= n, the elements on and above the (m-n)-th subdiagonal * contain the M-by-N upper trapezoidal matrix R; * the remaining elements, with the array TAU, represent the * orthogonal matrix Q as a product of min(m,n) elementary * reflectors (see Further Details). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) * The scalar factors of the elementary reflectors (see Further * Details). * * 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 >= max(1,M). * For optimum performance LWORK >= M*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 scalar, and v is a real vector with * v(n-k+i+1:n) = 0 and v(n-k+i) = 1; v(1:n-k+i-1) is stored on exit in * A(m-k+i,1:n-k+i-1), and tau in TAU(i). * * ===================================================================== * * .. Local Scalars .. LOGICAL LQUERY INTEGER I, IB, IINFO, IWS, K, KI, KK, LDWORK, LWKOPT, $ MU, NB, NBMIN, NU, NX * .. * .. External Subroutines .. EXTERNAL DGERQ2, DLARFB, DLARFT, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. External Functions .. INTEGER ILAENV EXTERNAL ILAENV * ..go to the page top

USAGE: scale, rhs = NumRu::Lapack.dgesc2( a, rhs, ipiv, jpiv, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE ) * Purpose * ======= * * DGESC2 solves a system of linear equations * * A * X = scale* RHS * * with a general N-by-N matrix A using the LU factorization with * complete pivoting computed by DGETC2. * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. * * A (input) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the LU part of the factorization of the n-by-n * matrix A computed by DGETC2: A = P * L * U * Q * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1, N). * * RHS (input/output) DOUBLE PRECISION array, dimension (N). * On entry, the right hand side vector b. * On exit, the solution vector X. * * IPIV (input) INTEGER array, dimension (N). * The pivot indices; for 1 <= i <= N, row i of the * matrix has been interchanged with row IPIV(i). * * JPIV (input) INTEGER array, dimension (N). * The pivot indices; for 1 <= j <= N, column j of the * matrix has been interchanged with column JPIV(j). * * SCALE (output) DOUBLE PRECISION * On exit, SCALE contains the scale factor. SCALE is chosen * 0 <= SCALE <= 1 to prevent owerflow in the solution. * * Further Details * =============== * * Based on contributions by * Bo Kagstrom and Peter Poromaa, Department of Computing Science, * Umea University, S-901 87 Umea, Sweden. * * ===================================================================== *go to the page top

USAGE: s, u, vt, work, info, a = NumRu::Lapack.dgesdd( jobz, a, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, IWORK, INFO ) * Purpose * ======= * * DGESDD computes the singular value decomposition (SVD) of a real * M-by-N matrix A, optionally computing the left and right singular * vectors. If singular vectors are desired, it uses a * divide-and-conquer algorithm. * * The SVD is written * * A = U * SIGMA * transpose(V) * * where SIGMA is an M-by-N matrix which is zero except for its * min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and * V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA * are the singular values of A; they are real and non-negative, and * are returned in descending order. The first min(m,n) columns of * U and V are the left and right singular vectors of A. * * Note that the routine returns VT = V**T, not V. * * The divide and conquer algorithm makes very mild assumptions about * floating point arithmetic. It will work on machines with a guard * digit in add/subtract, or on those binary machines without guard * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or * Cray-2. It could conceivably fail on hexadecimal or decimal machines * without guard digits, but we know of none. * * Arguments * ========= * * JOBZ (input) CHARACTER*1 * Specifies options for computing all or part of the matrix U: * = 'A': all M columns of U and all N rows of V**T are * returned in the arrays U and VT; * = 'S': the first min(M,N) columns of U and the first * min(M,N) rows of V**T are returned in the arrays U * and VT; * = 'O': If M >= N, the first N columns of U are overwritten * on the array A and all rows of V**T are returned in * the array VT; * otherwise, all columns of U are returned in the * array U and the first M rows of V**T are overwritten * in the array A; * = 'N': no columns of U or rows of V**T are computed. * * M (input) INTEGER * The number of rows of the input matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the input matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the M-by-N matrix A. * On exit, * if JOBZ = 'O', A is overwritten with the first N columns * of U (the left singular vectors, stored * columnwise) if M >= N; * A is overwritten with the first M rows * of V**T (the right singular vectors, stored * rowwise) otherwise. * if JOBZ .ne. 'O', the contents of A are destroyed. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * S (output) DOUBLE PRECISION array, dimension (min(M,N)) * The singular values of A, sorted so that S(i) >= S(i+1). * * U (output) DOUBLE PRECISION array, dimension (LDU,UCOL) * UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N; * UCOL = min(M,N) if JOBZ = 'S'. * If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M * orthogonal matrix U; * if JOBZ = 'S', U contains the first min(M,N) columns of U * (the left singular vectors, stored columnwise); * if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced. * * LDU (input) INTEGER * The leading dimension of the array U. LDU >= 1; if * JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M. * * VT (output) DOUBLE PRECISION array, dimension (LDVT,N) * If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the * N-by-N orthogonal matrix V**T; * if JOBZ = 'S', VT contains the first min(M,N) rows of * V**T (the right singular vectors, stored rowwise); * if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced. * * LDVT (input) INTEGER * The leading dimension of the array VT. LDVT >= 1; if * JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N; * if JOBZ = 'S', LDVT >= min(M,N). * * 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 >= 1. * If JOBZ = 'N', * LWORK >= 3*min(M,N) + max(max(M,N),7*min(M,N)). * If JOBZ = 'O', * LWORK >= 3*min(M,N) + * max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)). * If JOBZ = 'S' or 'A' * LWORK >= 3*min(M,N) + * max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)). * For good performance, LWORK should generally be larger. * If LWORK = -1 but other input arguments are legal, WORK(1) * returns the optimal LWORK. * * IWORK (workspace) INTEGER array, dimension (8*min(M,N)) * * INFO (output) INTEGER * = 0: successful exit. * < 0: if INFO = -i, the i-th argument had an illegal value. * > 0: DBDSDC did not converge, updating process failed. * * Further Details * =============== * * Based on contributions by * Ming Gu and Huan Ren, Computer Science Division, University of * California at Berkeley, USA * * ===================================================================== *go to the page top

USAGE: ipiv, info, a, b = NumRu::Lapack.dgesv( a, b, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO ) * Purpose * ======= * * DGESV computes the solution to a real system of linear equations * A * X = B, * where A is an N-by-N matrix and X and B are N-by-NRHS matrices. * * The LU decomposition with partial pivoting and row interchanges is * used to factor A as * A = P * L * U, * where P is a permutation matrix, L is unit lower triangular, and U is * upper triangular. The factored form of A is then used to solve the * system of equations A * X = B. * * Arguments * ========= * * N (input) INTEGER * The number of linear equations, i.e., the order of the * matrix A. N >= 0. * * NRHS (input) INTEGER * The number of right hand sides, i.e., the number of columns * of the matrix B. NRHS >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the N-by-N coefficient matrix A. * On exit, the factors L and U from the factorization * A = P*L*U; the unit diagonal elements of L are not stored. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * IPIV (output) INTEGER array, dimension (N) * The pivot indices that define the permutation matrix P; * row i of the matrix was interchanged with row IPIV(i). * * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) * On entry, the N-by-NRHS matrix of right hand side matrix B. * On exit, if INFO = 0, the N-by-NRHS solution matrix X. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, U(i,i) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, so the solution could not be computed. * * ===================================================================== * * .. External Subroutines .. EXTERNAL DGETRF, DGETRS, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * ..go to the page top

USAGE: s, u, vt, work, info, a = NumRu::Lapack.dgesvd( jobu, jobvt, a, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO ) * Purpose * ======= * * DGESVD computes the singular value decomposition (SVD) of a real * M-by-N matrix A, optionally computing the left and/or right singular * vectors. The SVD is written * * A = U * SIGMA * transpose(V) * * where SIGMA is an M-by-N matrix which is zero except for its * min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and * V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA * are the singular values of A; they are real and non-negative, and * are returned in descending order. The first min(m,n) columns of * U and V are the left and right singular vectors of A. * * Note that the routine returns V**T, not V. * * Arguments * ========= * * JOBU (input) CHARACTER*1 * Specifies options for computing all or part of the matrix U: * = 'A': all M columns of U are returned in array U: * = 'S': the first min(m,n) columns of U (the left singular * vectors) are returned in the array U; * = 'O': the first min(m,n) columns of U (the left singular * vectors) are overwritten on the array A; * = 'N': no columns of U (no left singular vectors) are * computed. * * JOBVT (input) CHARACTER*1 * Specifies options for computing all or part of the matrix * V**T: * = 'A': all N rows of V**T are returned in the array VT; * = 'S': the first min(m,n) rows of V**T (the right singular * vectors) are returned in the array VT; * = 'O': the first min(m,n) rows of V**T (the right singular * vectors) are overwritten on the array A; * = 'N': no rows of V**T (no right singular vectors) are * computed. * * JOBVT and JOBU cannot both be 'O'. * * M (input) INTEGER * The number of rows of the input matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the input matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the M-by-N matrix A. * On exit, * if JOBU = 'O', A is overwritten with the first min(m,n) * columns of U (the left singular vectors, * stored columnwise); * if JOBVT = 'O', A is overwritten with the first min(m,n) * rows of V**T (the right singular vectors, * stored rowwise); * if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A * are destroyed. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * S (output) DOUBLE PRECISION array, dimension (min(M,N)) * The singular values of A, sorted so that S(i) >= S(i+1). * * U (output) DOUBLE PRECISION array, dimension (LDU,UCOL) * (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'. * If JOBU = 'A', U contains the M-by-M orthogonal matrix U; * if JOBU = 'S', U contains the first min(m,n) columns of U * (the left singular vectors, stored columnwise); * if JOBU = 'N' or 'O', U is not referenced. * * LDU (input) INTEGER * The leading dimension of the array U. LDU >= 1; if * JOBU = 'S' or 'A', LDU >= M. * * VT (output) DOUBLE PRECISION array, dimension (LDVT,N) * If JOBVT = 'A', VT contains the N-by-N orthogonal matrix * V**T; * if JOBVT = 'S', VT contains the first min(m,n) rows of * V**T (the right singular vectors, stored rowwise); * if JOBVT = 'N' or 'O', VT is not referenced. * * LDVT (input) INTEGER * The leading dimension of the array VT. LDVT >= 1; if * JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N). * * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) * On exit, if INFO = 0, WORK(1) returns the optimal LWORK; * if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged * superdiagonal elements of an upper bidiagonal matrix B * whose diagonal is in S (not necessarily sorted). B * satisfies A = U * B * VT, so it has the same singular values * as A, and singular vectors related by U and VT. * * LWORK (input) INTEGER * The dimension of the array WORK. * LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)). * For good performance, LWORK should 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. * * INFO (output) INTEGER * = 0: successful exit. * < 0: if INFO = -i, the i-th argument had an illegal value. * > 0: if DBDSQR did not converge, INFO specifies how many * superdiagonals of an intermediate bidiagonal form B * did not converge to zero. See the description of WORK * above for details. * * ===================================================================== *go to the page top

USAGE: sva, info, a, v, work = NumRu::Lapack.dgesvj( joba, jobu, jobv, m, a, mv, v, work, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK, LWORK, INFO ) * Purpose * ======= * * DGESVJ computes the singular value decomposition (SVD) of a real * M-by-N matrix A, where M >= N. The SVD of A is written as * [++] [xx] [x0] [xx] * A = U * SIGMA * V^t, [++] = [xx] * [ox] * [xx] * [++] [xx] * where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal * matrix, and V is an N-by-N orthogonal matrix. The diagonal elements * of SIGMA are the singular values of A. The columns of U and V are the * left and the right singular vectors of A, respectively. * * Further Details * ~~~~~~~~~~~~~~~ * The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane * rotations. The rotations are implemented as fast scaled rotations of * Anda and Park [1]. In the case of underflow of the Jacobi angle, a * modified Jacobi transformation of Drmac [4] is used. Pivot strategy uses * column interchanges of de Rijk [2]. The relative accuracy of the computed * singular values and the accuracy of the computed singular vectors (in * angle metric) is as guaranteed by the theory of Demmel and Veselic [3]. * The condition number that determines the accuracy in the full rank case * is essentially min_{D=diag} kappa(A*D), where kappa(.) is the * spectral condition number. The best performance of this Jacobi SVD * procedure is achieved if used in an accelerated version of Drmac and * Veselic [5,6], and it is the kernel routine in the SIGMA library [7]. * Some tunning parameters (marked with [TP]) are available for the * implementer. * The computational range for the nonzero singular values is the machine * number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even * denormalized singular values can be computed with the corresponding * gradual loss of accurate digits. * * Contributors * ~~~~~~~~~~~~ * Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany) * * References * ~~~~~~~~~~ * [1] A. A. Anda and H. Park: Fast plane rotations with dynamic scaling. * SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174. * [2] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the * singular value decomposition on a vector computer. * SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371. * [3] J. Demmel and K. Veselic: Jacobi method is more accurate than QR. * [4] Z. Drmac: Implementation of Jacobi rotations for accurate singular * value computation in floating point arithmetic. * SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222. * [5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I. * SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342. * LAPACK Working note 169. * [6] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II. * SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362. * LAPACK Working note 170. * [7] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, * QSVD, (H,K)-SVD computations. * Department of Mathematics, University of Zagreb, 2008. * * Bugs, Examples and Comments * ~~~~~~~~~~~~~~~~~~~~~~~~~~~ * Please report all bugs and send interesting test examples and comments to * drmac@math.hr. Thank you. * * Arguments * ========= * * JOBA (input) CHARACTER* 1 * Specifies the structure of A. * = 'L': The input matrix A is lower triangular; * = 'U': The input matrix A is upper triangular; * = 'G': The input matrix A is general M-by-N matrix, M >= N. * * JOBU (input) CHARACTER*1 * Specifies whether to compute the left singular vectors * (columns of U): * = 'U': The left singular vectors corresponding to the nonzero * singular values are computed and returned in the leading * columns of A. See more details in the description of A. * The default numerical orthogonality threshold is set to * approximately TOL=CTOL*EPS, CTOL=DSQRT(M), EPS=DLAMCH('E'). * = 'C': Analogous to JOBU='U', except that user can control the * level of numerical orthogonality of the computed left * singular vectors. TOL can be set to TOL = CTOL*EPS, where * CTOL is given on input in the array WORK. * No CTOL smaller than ONE is allowed. CTOL greater * than 1 / EPS is meaningless. The option 'C' * can be used if M*EPS is satisfactory orthogonality * of the computed left singular vectors, so CTOL=M could * save few sweeps of Jacobi rotations. * See the descriptions of A and WORK(1). * = 'N': The matrix U is not computed. However, see the * description of A. * * JOBV (input) CHARACTER*1 * Specifies whether to compute the right singular vectors, that * is, the matrix V: * = 'V' : the matrix V is computed and returned in the array V * = 'A' : the Jacobi rotations are applied to the MV-by-N * array V. In other words, the right singular vector * matrix V is not computed explicitly, instead it is * applied to an MV-by-N matrix initially stored in the * first MV rows of V. * = 'N' : the matrix V is not computed and the array V is not * referenced * * M (input) INTEGER * The number of rows of the input matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the input matrix A. * M >= N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the M-by-N matrix A. * On exit : * If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C' : * If INFO .EQ. 0 : * RANKA orthonormal columns of U are returned in the * leading RANKA columns of the array A. Here RANKA <= N * is the number of computed singular values of A that are * above the underflow threshold DLAMCH('S'). The singular * vectors corresponding to underflowed or zero singular * values are not computed. The value of RANKA is returned * in the array WORK as RANKA=NINT(WORK(2)). Also see the * descriptions of SVA and WORK. The computed columns of U * are mutually numerically orthogonal up to approximately * TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'), * see the description of JOBU. * If INFO .GT. 0 : * the procedure DGESVJ did not converge in the given number * of iterations (sweeps). In that case, the computed * columns of U may not be orthogonal up to TOL. The output * U (stored in A), SIGMA (given by the computed singular * values in SVA(1:N)) and V is still a decomposition of the * input matrix A in the sense that the residual * ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small. * * If JOBU .EQ. 'N' : * If INFO .EQ. 0 : * Note that the left singular vectors are 'for free' in the * one-sided Jacobi SVD algorithm. However, if only the * singular values are needed, the level of numerical * orthogonality of U is not an issue and iterations are * stopped when the columns of the iterated matrix are * numerically orthogonal up to approximately M*EPS. Thus, * on exit, A contains the columns of U scaled with the * corresponding singular values. * If INFO .GT. 0 : * the procedure DGESVJ did not converge in the given number * of iterations (sweeps). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * SVA (workspace/output) DOUBLE PRECISION array, dimension (N) * On exit : * If INFO .EQ. 0 : * depending on the value SCALE = WORK(1), we have: * If SCALE .EQ. ONE : * SVA(1:N) contains the computed singular values of A. * During the computation SVA contains the Euclidean column * norms of the iterated matrices in the array A. * If SCALE .NE. ONE : * The singular values of A are SCALE*SVA(1:N), and this * factored representation is due to the fact that some of the * singular values of A might underflow or overflow. * If INFO .GT. 0 : * the procedure DGESVJ did not converge in the given number of * iterations (sweeps) and SCALE*SVA(1:N) may not be accurate. * * MV (input) INTEGER * If JOBV .EQ. 'A', then the product of Jacobi rotations in DGESVJ * is applied to the first MV rows of V. See the description of JOBV. * * V (input/output) DOUBLE PRECISION array, dimension (LDV,N) * If JOBV = 'V', then V contains on exit the N-by-N matrix of * the right singular vectors; * If JOBV = 'A', then V contains the product of the computed right * singular vector matrix and the initial matrix in * the array V. * If JOBV = 'N', then V is not referenced. * * LDV (input) INTEGER * The leading dimension of the array V, LDV .GE. 1. * If JOBV .EQ. 'V', then LDV .GE. max(1,N). * If JOBV .EQ. 'A', then LDV .GE. max(1,MV) . * * WORK (input/workspace/output) DOUBLE PRECISION array, dimension max(4,M+N). * On entry : * If JOBU .EQ. 'C' : * WORK(1) = CTOL, where CTOL defines the threshold for convergence. * The process stops if all columns of A are mutually * orthogonal up to CTOL*EPS, EPS=DLAMCH('E'). * It is required that CTOL >= ONE, i.e. it is not * allowed to force the routine to obtain orthogonality * below EPS. * On exit : * WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N) * are the computed singular values of A. * (See description of SVA().) * WORK(2) = NINT(WORK(2)) is the number of the computed nonzero * singular values. * WORK(3) = NINT(WORK(3)) is the number of the computed singular * values that are larger than the underflow threshold. * WORK(4) = NINT(WORK(4)) is the number of sweeps of Jacobi * rotations needed for numerical convergence. * WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep. * This is useful information in cases when DGESVJ did * not converge, as it can be used to estimate whether * the output is stil useful and for post festum analysis. * WORK(6) = the largest absolute value over all sines of the * Jacobi rotation angles in the last sweep. It can be * useful for a post festum analysis. * * LWORK (input) INTEGER * length of WORK, WORK >= MAX(6,M+N) * * INFO (output) INTEGER * = 0 : successful exit. * < 0 : if INFO = -i, then the i-th argument had an illegal value * > 0 : DGESVJ did not converge in the maximal allowed number (30) * of sweeps. The output may still be useful. See the * description of WORK. * * ===================================================================== * * .. Local Parameters .. DOUBLE PRECISION ZERO, HALF, ONE, TWO PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, + TWO = 2.0D0 ) INTEGER NSWEEP PARAMETER ( NSWEEP = 30 ) * .. * .. Local Scalars .. DOUBLE PRECISION AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG, + BIGTHETA, CS, CTOL, EPSLN, LARGE, MXAAPQ, + MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL, + SKL, SFMIN, SMALL, SN, T, TEMP1, THETA, + THSIGN, TOL INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1, + ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34, + N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP, + SWBAND LOGICAL APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK, + RSVEC, UCTOL, UPPER * .. * .. Local Arrays .. DOUBLE PRECISION FASTR( 5 ) * .. * .. Intrinsic Functions .. INTRINSIC DABS, DMAX1, DMIN1, DBLE, MIN0, DSIGN, DSQRT * .. * .. External Functions .. * .. * from BLAS DOUBLE PRECISION DDOT, DNRM2 EXTERNAL DDOT, DNRM2 INTEGER IDAMAX EXTERNAL IDAMAX * from LAPACK DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. * .. * from BLAS EXTERNAL DAXPY, DCOPY, DROTM, DSCAL, DSWAP * from LAPACK EXTERNAL DLASCL, DLASET, DLASSQ, XERBLA * EXTERNAL DGSVJ0, DGSVJ1 * ..go to the page top

USAGE: x, rcond, ferr, berr, work, info, a, af, ipiv, equed, r, c, b = NumRu::Lapack.dgesvx( fact, trans, a, b, [:af => af, :ipiv => ipiv, :equed => equed, :r => r, :c => c, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGESVX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, IWORK, INFO ) * Purpose * ======= * * DGESVX uses the LU factorization to compute the solution to a real * system of linear equations * A * X = B, * where A is an N-by-N matrix and X and B are N-by-NRHS matrices. * * Error bounds on the solution and a condition estimate are also * provided. * * Description * =========== * * The following steps are performed: * * 1. If FACT = 'E', real scaling factors are computed to equilibrate * the system: * TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B * TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B * TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B * Whether or not the system will be equilibrated depends on the * scaling of the matrix A, but if equilibration is used, A is * overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N') * or diag(C)*B (if TRANS = 'T' or 'C'). * * 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the * matrix A (after equilibration if FACT = 'E') as * A = P * L * U, * where P is a permutation matrix, L is a unit lower triangular * matrix, and U is upper triangular. * * 3. If some U(i,i)=0, so that U is exactly singular, then the routine * returns with INFO = i. Otherwise, the factored form of A is used * to estimate the condition number of the matrix A. If the * reciprocal of the condition number is less than machine precision, * INFO = N+1 is returned as a warning, but the routine still goes on * to solve for X and compute error bounds as described below. * * 4. The system of equations is solved for X using the factored form * of A. * * 5. Iterative refinement is applied to improve the computed solution * matrix and calculate error bounds and backward error estimates * for it. * * 6. If equilibration was used, the matrix X is premultiplied by * diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so * that it solves the original system before equilibration. * * Arguments * ========= * * FACT (input) CHARACTER*1 * Specifies whether or not the factored form of the matrix A is * supplied on entry, and if not, whether the matrix A should be * equilibrated before it is factored. * = 'F': On entry, AF and IPIV contain the factored form of A. * If EQUED is not 'N', the matrix A has been * equilibrated with scaling factors given by R and C. * A, AF, and IPIV are not modified. * = 'N': The matrix A will be copied to AF and factored. * = 'E': The matrix A will be equilibrated if necessary, then * copied to AF and factored. * * TRANS (input) CHARACTER*1 * Specifies the form of the system of equations: * = 'N': A * X = B (No transpose) * = 'T': A**T * X = B (Transpose) * = 'C': A**H * X = B (Transpose) * * N (input) INTEGER * The number of linear equations, i.e., the order of the * matrix A. N >= 0. * * NRHS (input) INTEGER * The number of right hand sides, i.e., the number of columns * of the matrices B and X. NRHS >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the N-by-N matrix A. If FACT = 'F' and EQUED is * not 'N', then A must have been equilibrated by the scaling * factors in R and/or C. A is not modified if FACT = 'F' or * 'N', or if FACT = 'E' and EQUED = 'N' on exit. * * On exit, if EQUED .ne. 'N', A is scaled as follows: * EQUED = 'R': A := diag(R) * A * EQUED = 'C': A := A * diag(C) * EQUED = 'B': A := diag(R) * A * diag(C). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * AF (input or output) DOUBLE PRECISION array, dimension (LDAF,N) * If FACT = 'F', then AF is an input argument and on entry * contains the factors L and U from the factorization * A = P*L*U as computed by DGETRF. If EQUED .ne. 'N', then * AF is the factored form of the equilibrated matrix A. * * If FACT = 'N', then AF is an output argument and on exit * returns the factors L and U from the factorization A = P*L*U * of the original matrix A. * * If FACT = 'E', then AF is an output argument and on exit * returns the factors L and U from the factorization A = P*L*U * of the equilibrated matrix A (see the description of A for * the form of the equilibrated matrix). * * LDAF (input) INTEGER * The leading dimension of the array AF. LDAF >= max(1,N). * * IPIV (input or output) INTEGER array, dimension (N) * If FACT = 'F', then IPIV is an input argument and on entry * contains the pivot indices from the factorization A = P*L*U * as computed by DGETRF; row i of the matrix was interchanged * with row IPIV(i). * * If FACT = 'N', then IPIV is an output argument and on exit * contains the pivot indices from the factorization A = P*L*U * of the original matrix A. * * If FACT = 'E', then IPIV is an output argument and on exit * contains the pivot indices from the factorization A = P*L*U * of the equilibrated matrix A. * * EQUED (input or output) CHARACTER*1 * Specifies the form of equilibration that was done. * = 'N': No equilibration (always true if FACT = 'N'). * = 'R': Row equilibration, i.e., A has been premultiplied by * diag(R). * = 'C': Column equilibration, i.e., A has been postmultiplied * by diag(C). * = 'B': Both row and column equilibration, i.e., A has been * replaced by diag(R) * A * diag(C). * EQUED is an input argument if FACT = 'F'; otherwise, it is an * output argument. * * R (input or output) DOUBLE PRECISION array, dimension (N) * The row scale factors for A. If EQUED = 'R' or 'B', A is * multiplied on the left by diag(R); if EQUED = 'N' or 'C', R * is not accessed. R is an input argument if FACT = 'F'; * otherwise, R is an output argument. If FACT = 'F' and * EQUED = 'R' or 'B', each element of R must be positive. * * C (input or output) DOUBLE PRECISION array, dimension (N) * The column scale factors for A. If EQUED = 'C' or 'B', A is * multiplied on the right by diag(C); if EQUED = 'N' or 'R', C * is not accessed. C is an input argument if FACT = 'F'; * otherwise, C is an output argument. If FACT = 'F' and * EQUED = 'C' or 'B', each element of C must be positive. * * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) * On entry, the N-by-NRHS right hand side matrix B. * On exit, * if EQUED = 'N', B is not modified; * if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by * diag(R)*B; * if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is * overwritten by diag(C)*B. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,N). * * X (output) DOUBLE PRECISION array, dimension (LDX,NRHS) * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X * to the original system of equations. Note that A and B are * modified on exit if EQUED .ne. 'N', and the solution to the * equilibrated system is inv(diag(C))*X if TRANS = 'N' and * EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C' * and EQUED = 'R' or 'B'. * * LDX (input) INTEGER * The leading dimension of the array X. LDX >= max(1,N). * * RCOND (output) DOUBLE PRECISION * The estimate of the reciprocal condition number of the matrix * A after equilibration (if done). If RCOND is less than the * machine precision (in particular, if RCOND = 0), the matrix * is singular to working precision. This condition is * indicated by a return code of INFO > 0. * * FERR (output) DOUBLE PRECISION array, dimension (NRHS) * The estimated forward error bound for each solution vector * X(j) (the j-th column of the solution matrix X). * If XTRUE is the true solution corresponding to X(j), FERR(j) * is an estimated upper bound for the magnitude of the largest * element in (X(j) - XTRUE) divided by the magnitude of the * largest element in X(j). The estimate is as reliable as * the estimate for RCOND, and is almost always a slight * overestimate of the true error. * * BERR (output) DOUBLE PRECISION array, dimension (NRHS) * The componentwise relative backward error of each solution * vector X(j) (i.e., the smallest relative change in * any element of A or B that makes X(j) an exact solution). * * WORK (workspace/output) DOUBLE PRECISION array, dimension (4*N) * On exit, WORK(1) contains the reciprocal pivot growth * factor norm(A)/norm(U). The "max absolute element" norm is * used. If WORK(1) is much less than 1, then the stability * of the LU factorization of the (equilibrated) matrix A * could be poor. This also means that the solution X, condition * estimator RCOND, and forward error bound FERR could be * unreliable. If factorization fails with 0go to the page top0: if INFO = i, and i is * <= N: U(i,i) is exactly zero. The factorization has * been completed, but the factor U is exactly * singular, so the solution and error bounds * could not be computed. RCOND = 0 is returned. * = N+1: U is nonsingular, but RCOND is less than machine * precision, meaning that the matrix is singular * to working precision. Nevertheless, the * solution and error bounds are computed because * there are a number of situations where the * computed solution can be more accurate than the * value of RCOND would suggest. * * ===================================================================== *

USAGE: x, rcond, rpvgrw, berr, err_bnds_norm, err_bnds_comp, info, a, af, ipiv, equed, r, c, b, params = NumRu::Lapack.dgesvxx( fact, trans, a, af, ipiv, equed, r, c, b, params, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGESVXX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, EQUED, R, C, B, LDB, X, LDX, RCOND, RPVGRW, BERR, N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, IWORK, INFO ) * Purpose * ======= * * DGESVXX uses the LU factorization to compute the solution to a * double precision system of linear equations A * X = B, where A is an * N-by-N matrix and X and B are N-by-NRHS matrices. * * If requested, both normwise and maximum componentwise error bounds * are returned. DGESVXX will return a solution with a tiny * guaranteed error (O(eps) where eps is the working machine * precision) unless the matrix is very ill-conditioned, in which * case a warning is returned. Relevant condition numbers also are * calculated and returned. * * DGESVXX accepts user-provided factorizations and equilibration * factors; see the definitions of the FACT and EQUED options. * Solving with refinement and using a factorization from a previous * DGESVXX call will also produce a solution with either O(eps) * errors or warnings, but we cannot make that claim for general * user-provided factorizations and equilibration factors if they * differ from what DGESVXX would itself produce. * * Description * =========== * * The following steps are performed: * * 1. If FACT = 'E', double precision scaling factors are computed to equilibrate * the system: * * TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B * TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B * TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B * * Whether or not the system will be equilibrated depends on the * scaling of the matrix A, but if equilibration is used, A is * overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N') * or diag(C)*B (if TRANS = 'T' or 'C'). * * 2. If FACT = 'N' or 'E', the LU decomposition is used to factor * the matrix A (after equilibration if FACT = 'E') as * * A = P * L * U, * * where P is a permutation matrix, L is a unit lower triangular * matrix, and U is upper triangular. * * 3. If some U(i,i)=0, so that U is exactly singular, then the * routine returns with INFO = i. Otherwise, the factored form of A * is used to estimate the condition number of the matrix A (see * argument RCOND). If the reciprocal of the condition number is less * than machine precision, the routine still goes on to solve for X * and compute error bounds as described below. * * 4. The system of equations is solved for X using the factored form * of A. * * 5. By default (unless PARAMS(LA_LINRX_ITREF_I) is set to zero), * the routine will use iterative refinement to try to get a small * error and error bounds. Refinement calculates the residual to at * least twice the working precision. * * 6. If equilibration was used, the matrix X is premultiplied by * diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so * that it solves the original system before equilibration. * * Arguments * ========= * * Some optional parameters are bundled in the PARAMS array. These * settings determine how refinement is performed, but often the * defaults are acceptable. If the defaults are acceptable, users * can pass NPARAMS = 0 which prevents the source code from accessing * the PARAMS argument. * * FACT (input) CHARACTER*1 * Specifies whether or not the factored form of the matrix A is * supplied on entry, and if not, whether the matrix A should be * equilibrated before it is factored. * = 'F': On entry, AF and IPIV contain the factored form of A. * If EQUED is not 'N', the matrix A has been * equilibrated with scaling factors given by R and C. * A, AF, and IPIV are not modified. * = 'N': The matrix A will be copied to AF and factored. * = 'E': The matrix A will be equilibrated if necessary, then * copied to AF and factored. * * TRANS (input) CHARACTER*1 * Specifies the form of the system of equations: * = 'N': A * X = B (No transpose) * = 'T': A**T * X = B (Transpose) * = 'C': A**H * X = B (Conjugate Transpose = Transpose) * * N (input) INTEGER * The number of linear equations, i.e., the order of the * matrix A. N >= 0. * * NRHS (input) INTEGER * The number of right hand sides, i.e., the number of columns * of the matrices B and X. NRHS >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the N-by-N matrix A. If FACT = 'F' and EQUED is * not 'N', then A must have been equilibrated by the scaling * factors in R and/or C. A is not modified if FACT = 'F' or * 'N', or if FACT = 'E' and EQUED = 'N' on exit. * * On exit, if EQUED .ne. 'N', A is scaled as follows: * EQUED = 'R': A := diag(R) * A * EQUED = 'C': A := A * diag(C) * EQUED = 'B': A := diag(R) * A * diag(C). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * AF (input or output) DOUBLE PRECISION array, dimension (LDAF,N) * If FACT = 'F', then AF is an input argument and on entry * contains the factors L and U from the factorization * A = P*L*U as computed by DGETRF. If EQUED .ne. 'N', then * AF is the factored form of the equilibrated matrix A. * * If FACT = 'N', then AF is an output argument and on exit * returns the factors L and U from the factorization A = P*L*U * of the original matrix A. * * If FACT = 'E', then AF is an output argument and on exit * returns the factors L and U from the factorization A = P*L*U * of the equilibrated matrix A (see the description of A for * the form of the equilibrated matrix). * * LDAF (input) INTEGER * The leading dimension of the array AF. LDAF >= max(1,N). * * IPIV (input or output) INTEGER array, dimension (N) * If FACT = 'F', then IPIV is an input argument and on entry * contains the pivot indices from the factorization A = P*L*U * as computed by DGETRF; row i of the matrix was interchanged * with row IPIV(i). * * If FACT = 'N', then IPIV is an output argument and on exit * contains the pivot indices from the factorization A = P*L*U * of the original matrix A. * * If FACT = 'E', then IPIV is an output argument and on exit * contains the pivot indices from the factorization A = P*L*U * of the equilibrated matrix A. * * EQUED (input or output) CHARACTER*1 * Specifies the form of equilibration that was done. * = 'N': No equilibration (always true if FACT = 'N'). * = 'R': Row equilibration, i.e., A has been premultiplied by * diag(R). * = 'C': Column equilibration, i.e., A has been postmultiplied * by diag(C). * = 'B': Both row and column equilibration, i.e., A has been * replaced by diag(R) * A * diag(C). * EQUED is an input argument if FACT = 'F'; otherwise, it is an * output argument. * * R (input or output) DOUBLE PRECISION array, dimension (N) * The row scale factors for A. If EQUED = 'R' or 'B', A is * multiplied on the left by diag(R); if EQUED = 'N' or 'C', R * is not accessed. R is an input argument if FACT = 'F'; * otherwise, R is an output argument. If FACT = 'F' and * EQUED = 'R' or 'B', each element of R must be positive. * If R is output, each element of R is a power of the radix. * If R is input, each element of R should be a power of the radix * to ensure a reliable solution and error estimates. Scaling by * powers of the radix does not cause rounding errors unless the * result underflows or overflows. Rounding errors during scaling * lead to refining with a matrix that is not equivalent to the * input matrix, producing error estimates that may not be * reliable. * * C (input or output) DOUBLE PRECISION array, dimension (N) * The column scale factors for A. If EQUED = 'C' or 'B', A is * multiplied on the right by diag(C); if EQUED = 'N' or 'R', C * is not accessed. C is an input argument if FACT = 'F'; * otherwise, C is an output argument. If FACT = 'F' and * EQUED = 'C' or 'B', each element of C must be positive. * If C is output, each element of C is a power of the radix. * If C is input, each element of C should be a power of the radix * to ensure a reliable solution and error estimates. Scaling by * powers of the radix does not cause rounding errors unless the * result underflows or overflows. Rounding errors during scaling * lead to refining with a matrix that is not equivalent to the * input matrix, producing error estimates that may not be * reliable. * * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) * On entry, the N-by-NRHS right hand side matrix B. * On exit, * if EQUED = 'N', B is not modified; * if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by * diag(R)*B; * if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is * overwritten by diag(C)*B. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,N). * * X (output) DOUBLE PRECISION array, dimension (LDX,NRHS) * If INFO = 0, the N-by-NRHS solution matrix X to the original * system of equations. Note that A and B are modified on exit * if EQUED .ne. 'N', and the solution to the equilibrated system is * inv(diag(C))*X if TRANS = 'N' and EQUED = 'C' or 'B', or * inv(diag(R))*X if TRANS = 'T' or 'C' and EQUED = 'R' or 'B'. * * LDX (input) INTEGER * The leading dimension of the array X. LDX >= max(1,N). * * RCOND (output) DOUBLE PRECISION * Reciprocal scaled condition number. This is an estimate of the * reciprocal Skeel condition number of the matrix A after * equilibration (if done). If this is less than the machine * precision (in particular, if it is zero), the matrix is singular * to working precision. Note that the error may still be small even * if this number is very small and the matrix appears ill- * conditioned. * * RPVGRW (output) DOUBLE PRECISION * Reciprocal pivot growth. On exit, this contains the reciprocal * pivot growth factor norm(A)/norm(U). The "max absolute element" * norm is used. If this is much less than 1, then the stability of * the LU factorization of the (equilibrated) matrix A could be poor. * This also means that the solution X, estimated condition numbers, * and error bounds could be unreliable. If factorization fails with * 0go to the page top0 and <= N: U(INFO,INFO) is exactly zero. The factorization * has been completed, but the factor U is exactly singular, so * the solution and error bounds could not be computed. RCOND = 0 * is returned. * = N+J: The solution corresponding to the Jth right-hand side is * not guaranteed. The solutions corresponding to other right- * hand sides K with K > J may not be guaranteed as well, but * only the first such right-hand side is reported. If a small * componentwise error is not requested (PARAMS(3) = 0.0) then * the Jth right-hand side is the first with a normwise error * bound that is not guaranteed (the smallest J such * that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0) * the Jth right-hand side is the first with either a normwise or * componentwise error bound that is not guaranteed (the smallest * J such that either ERR_BNDS_NORM(J,1) = 0.0 or * ERR_BNDS_COMP(J,1) = 0.0). See the definition of * ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information * about all of the right-hand sides check ERR_BNDS_NORM or * ERR_BNDS_COMP. * * ================================================================== *

USAGE: ipiv, jpiv, info, a = NumRu::Lapack.dgetc2( a, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGETC2( N, A, LDA, IPIV, JPIV, INFO ) * Purpose * ======= * * DGETC2 computes an LU factorization with complete pivoting of the * n-by-n matrix A. The factorization has the form A = P * L * U * Q, * where P and Q are permutation matrices, L is lower triangular with * unit diagonal elements and U is upper triangular. * * This is the Level 2 BLAS algorithm. * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA, N) * On entry, the n-by-n matrix A to be factored. * On exit, the factors L and U from the factorization * A = P*L*U*Q; the unit diagonal elements of L are not stored. * If U(k, k) appears to be less than SMIN, U(k, k) is given the * value of SMIN, i.e., giving a nonsingular perturbed system. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * IPIV (output) INTEGER array, dimension(N). * The pivot indices; for 1 <= i <= N, row i of the * matrix has been interchanged with row IPIV(i). * * JPIV (output) INTEGER array, dimension(N). * The pivot indices; for 1 <= j <= N, column j of the * matrix has been interchanged with column JPIV(j). * * INFO (output) INTEGER * = 0: successful exit * > 0: if INFO = k, U(k, k) is likely to produce owerflow if * we try to solve for x in Ax = b. So U is perturbed to * avoid the overflow. * * Further Details * =============== * * Based on contributions by * Bo Kagstrom and Peter Poromaa, Department of Computing Science, * Umea University, S-901 87 Umea, Sweden. * * ===================================================================== *go to the page top

USAGE: ipiv, info, a = NumRu::Lapack.dgetf2( m, a, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGETF2( M, N, A, LDA, IPIV, INFO ) * Purpose * ======= * * DGETF2 computes an LU factorization of a general m-by-n matrix A * using partial pivoting with row interchanges. * * The factorization has the form * A = P * L * U * where P is a permutation matrix, L is lower triangular with unit * diagonal elements (lower trapezoidal if m > n), and U is upper * triangular (upper trapezoidal if m < n). * * This is the right-looking Level 2 BLAS version of the algorithm. * * 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 to be factored. * On exit, the factors L and U from the factorization * A = P*L*U; the unit diagonal elements of L are not stored. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * IPIV (output) INTEGER array, dimension (min(M,N)) * The pivot indices; for 1 <= i <= min(M,N), row i of the * matrix was interchanged with row IPIV(i). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * > 0: if INFO = k, U(k,k) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, and division by zero will occur if it is used * to solve a system of equations. * * ===================================================================== *go to the page top

USAGE: ipiv, info, a = NumRu::Lapack.dgetrf( m, a, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO ) * Purpose * ======= * * DGETRF computes an LU factorization of a general M-by-N matrix A * using partial pivoting with row interchanges. * * The factorization has the form * A = P * L * U * where P is a permutation matrix, L is lower triangular with unit * diagonal elements (lower trapezoidal if m > n), and U is upper * triangular (upper trapezoidal if m < n). * * This is the right-looking Level 3 BLAS version of the algorithm. * * 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 to be factored. * On exit, the factors L and U from the factorization * A = P*L*U; the unit diagonal elements of L are not stored. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * IPIV (output) INTEGER array, dimension (min(M,N)) * The pivot indices; for 1 <= i <= min(M,N), row i of the * matrix was interchanged with row IPIV(i). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, U(i,i) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, and division by zero will occur if it is used * to solve a system of equations. * * ===================================================================== *go to the page top

USAGE: work, info, a = NumRu::Lapack.dgetri( a, ipiv, [:lwork => lwork, :usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO ) * Purpose * ======= * * DGETRI computes the inverse of a matrix using the LU factorization * computed by DGETRF. * * This method inverts U and then computes inv(A) by solving the system * inv(A)*L = inv(U) for inv(A). * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the factors L and U from the factorization * A = P*L*U as computed by DGETRF. * On exit, if INFO = 0, the inverse of the original matrix A. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * IPIV (input) INTEGER array, dimension (N) * The pivot indices from DGETRF; for 1<=i<=N, row i of the * matrix was interchanged with row IPIV(i). * * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) * On exit, if INFO=0, then WORK(1) returns the optimal LWORK. * * LWORK (input) INTEGER * The dimension of the array WORK. LWORK >= max(1,N). * For optimal performance LWORK >= N*NB, where NB is * the optimal blocksize returned by ILAENV. * * 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 * > 0: if INFO = i, U(i,i) is exactly zero; the matrix is * singular and its inverse could not be computed. * * ===================================================================== *go to the page top

USAGE: info, b = NumRu::Lapack.dgetrs( trans, a, ipiv, b, [:usage => usage, :help => help]) FORTRAN MANUAL SUBROUTINE DGETRS( TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO ) * Purpose * ======= * * DGETRS solves a system of linear equations * A * X = B or A' * X = B * with a general N-by-N matrix A using the LU factorization computed * by DGETRF. * * Arguments * ========= * * TRANS (input) CHARACTER*1 * Specifies the form of the system of equations: * = 'N': A * X = B (No transpose) * = 'T': A'* X = B (Transpose) * = 'C': A'* X = B (Conjugate transpose = Transpose) * * N (input) INTEGER * The order of the matrix A. N >= 0. * * NRHS (input) INTEGER * The number of right hand sides, i.e., the number of columns * of the matrix B. NRHS >= 0. * * A (input) DOUBLE PRECISION array, dimension (LDA,N) * The factors L and U from the factorization A = P*L*U * as computed by DGETRF. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * IPIV (input) INTEGER array, dimension (N) * The pivot indices from DGETRF; for 1<=i<=N, row i of the * matrix was interchanged with row IPIV(i). * * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) * On entry, the right hand side matrix B. * On exit, the solution matrix X. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * * ===================================================================== *go to the page top

back to matrix types

back to data types