#ifndef _NUMclapack_h_ #define _NUMclapack_h_ /* NUMclapack.h * * Copyright (C) 1994-2003 David Weenink * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or (at * your option) any later version. * * This program is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ /* djmw 20020923 GPL header djmw 20030310 Latest modification */ /* The following routines all use FORTRAN column-major storage of matrices. A consequence is that all matrices must have been allocated as a single block of mxn elements. All matrices are passed as a vector of mxn elements. matrix[i][j] => vector[(j-1)*m + i] "Fortran" matrix[i][j] => vector[(i-1)*n + j] "C" The consequence is that you have to transpose C matrices before you pass them to a CLAPACK routine. Sometimes you can avoid transposition by considering the solution of the transposed problem (e.g. See code in SVD_compute). */ int NUMlapack_dbdsqr(char *uplo, long *n, long *ncvt, long *nru, long *ncc, double *d, double *e, double *vt, long *ldvt, double *u, long *ldu, double *c, long *ldc, double *work, long *info); /* Purpose ======= NUMlapack_dbdsqr computes the singular value decomposition (SVD) of a real N-by-N (upper or lower) bidiagonal matrix B: B = Q * S * P' (P' denotes the transpose of P), where S is a diagonal matrix with non-negative diagonal elements (the singular values of B), and Q and P are orthogonal matrices. The routine computes S, and optionally computes U * Q, P' * VT, or Q' * C, for given real input matrices U, VT, and C. See "Computing Small Singular Values of Bidiagonal Matrices With Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11, no. 5, pp. 873-912, Sept 1990) and "Accurate singular values and differential qd algorithms," by B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics Department, University of California at Berkeley, July 1992 for a detailed description of the algorithm. Arguments ========= UPLO (input) char* = 'U': B is upper bidiagonal; = 'L': B is lower bidiagonal. N (input) long The order of the matrix B. N >= 0. NCVT (input) long The number of columns of the matrix VT. NCVT >= 0. NRU (input) long The number of rows of the matrix U. NRU >= 0. NCC (input) long The number of columns of the matrix C. NCC >= 0. D (input/output) double array, dimension (N) On entry, the n diagonal elements of the bidiagonal matrix B. On exit, if INFO=0, the singular values of B in decreasing order. E (input/output) double array, dimension (N) On entry, the elements of E contain the offdiagonal elements of the bidiagonal matrix whose SVD is desired. On normal exit (INFO = 0), E is destroyed. If the algorithm does not converge (INFO > 0), D and E will contain the diagonal and superdiagonal elements of a bidiagonal matrix orthogonally equivalent to the one given as input. E(N) is used for workspace. VT (input/output) double array, dimension (LDVT, NCVT) On entry, an N-by-NCVT matrix VT. On exit, VT is overwritten by P' * VT. VT is not referenced if NCVT = 0. LDVT (input) long The leading dimension of the array VT. LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0. U (input/output) double array, dimension (LDU, N) On entry, an NRU-by-N matrix U. On exit, U is overwritten by U * Q. U is not referenced if NRU = 0. LDU (input) long The leading dimension of the array U. LDU >= max(1,NRU). C (input/output) double array, dimension (LDC, NCC) On entry, an N-by-NCC matrix C. On exit, C is overwritten by Q' * C. C is not referenced if NCC = 0. LDC (input) long The leading dimension of the array C. LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0. WORK (workspace) double array, dimension (4*N) INFO (output) long = 0: successful exit < 0: If INFO = -i, the i-th argument had an illegal value > 0: the algorithm did not converge; D and E contain the elements of a bidiagonal matrix which is orthogonally similar to the input matrix B; if INFO = i, i elements of E have not converged to zero. Internal Parameters =================== TOLMUL double, default = max(10,min(100,EPS**(-1/8))) TOLMUL controls the convergence criterion of the QR loop. If it is positive, TOLMUL*EPS is the desired relative precision in the computed singular values. If it is negative, abs(TOLMUL*EPS*sigma_max) is the desired absolute accuracy in the computed singular values (corresponds to relative accuracy abs(TOLMUL*EPS) in the largest singular value. abs(TOLMUL) should be between 1 and 1/EPS, and preferably between 10 (for fast convergence) and .1/EPS (for there to be some accuracy in the results). Default is to lose at either one eighth or 2 of the available decimal digits in each computed singular value (whichever is smaller). MAXITR long, default = 6 MAXITR controls the maximum number of passes of the algorithm through its inner loop. The algorithms stops (and so fails to converge) if the number of passes through the inner loop exceeds MAXITR*N**2. ===================================================================== */ int NUMlapack_dgebd2(long *m, long *n, double *a, long *lda, double *d, double *e, double *tauq, double *taup, double *work, long *info); /* Purpose ======= NUMlapack_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) long The number of rows in the matrix A. M >= 0. N (input) long The number of columns in the matrix A. N >= 0. A (input/output) double 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) long The leading dimension of the array A. LDA >= max(1,M). D (output) double array, dimension (min(M,N)) The diagonal elements of the bidiagonal matrix B: D(i) = A(i,i). E (output) double 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 array dimension (min(M,N)) The scalar factors of the elementary reflectors which represent the orthogonal matrix Q. See Further Details. TAUP (output) double array, dimension (min(M,N)) The scalar factors of the elementary reflectors which represent the orthogonal matrix P. See Further Details. WORK (workspace) double array, dimension (max(M,N)) INFO (output) long = 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). ===================================================================== */ int NUMlapack_dgebrd(long *m, long *n, double *a, long *lda, double *d, double *e, double *tauq, double *taup, double *work, long *lwork, long *info); /* Purpose ======= NUMlapack_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) long The number of rows in the matrix A. M >= 0. N (input) long The number of columns in the matrix A. N >= 0. A (input/output) double 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) long The leading dimension of the array A. LDA >= max(1,M). D (output) double array, dimension (min(M,N)) The diagonal elements of the bidiagonal matrix B: D(i) = A(i,i). E (output) double 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 array dimension (min(M,N)) The scalar factors of the elementary reflectors which represent the orthogonal matrix Q. See Further Details. TAUP (output) double 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 array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) long 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) long = 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). ===================================================================== */ int NUMlapack_dgebak (char *job, char *side, long *n, long *ilo, long *ihi, double *scale, long *m, double *v, long *ldv, long *info); /* -- LAPACK routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University September 30, 1994 Purpose ======= NUMlapack_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 NUMlapack_dgebal. Arguments ========= JOB (input) char* 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 NUMlapack_dgebal. SIDE (input) char* = 'R': V contains right eigenvectors; = 'L': V contains left eigenvectors. N (input) long The number of rows of the matrix V. N >= 0. ILO (input) long IHI (input) long The integers ILO and IHI determined by NUMlapack_dgebal. 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. SCALE (input) double array, dimension (N) Details of the permutation and scaling factors, as returned by NUMlapack_dgebal. M (input) long The number of columns of the matrix V. M >= 0. V (input/output) double array, dimension (LDV,M) On entry, the matrix of right or left eigenvectors to be transformed, as returned by DHSEIN or NUMlapack_dtrevc. On exit, V is overwritten by the transformed eigenvectors. LDV (input) long The leading dimension of the array V. LDV >= max(1,N). INFO (output) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value. ===================================================================== */ int NUMlapack_dgebal (char *job, long *n, double *a, long *lda, long *ilo, long *ihi, double *scale, long *info); /* -- LAPACK routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University June 30, 1999 Purpose ======= NUMlapack_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) char* 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) long The order of the matrix A. N >= 0. A (input/output) double 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) long The leading dimension of the array A. LDA >= max(1,N). ILO (output) long IHI (output) long 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 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) long = 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 ===================================================================== */ int NUMlapack_dgeev (char *jobvl, char *jobvr, long *n, double *a, long *lda, double *wr, double *wi, double *vl, long *ldvl, double *vr, long *ldvr, double *work, long *lwork, long *info); /* -- LAPACK driver routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University December 8, 1999 Purpose ======= NUMlapack_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) char* = 'N': left eigenvectors of A are not computed; = 'V': left eigenvectors of A are computed. JOBVR (input) char* = 'N': right eigenvectors of A are not computed; = 'V': right eigenvectors of A are computed. N (input) long The order of the matrix A. N >= 0. A (input/output) double array, dimension (LDA,N) On entry, the N-by-N matrix A. On exit, A has been overwritten. LDA (input) long The leading dimension of the array A. LDA >= max(1,N). WR (output) double array, dimension (N) WI (output) double 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 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) long The leading dimension of the array VL. LDVL >= 1; if JOBVL = 'V', LDVL >= N. VR (output) double 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) long The leading dimension of the array VR. LDVR >= 1; if JOBVR = 'V', LDVR >= N. WORK (workspace/output) double array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) long 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) long = 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. ===================================================================== */ int NUMlapack_dgehd2 (long *n, long *ilo, long *ihi, double *a, long *lda, double *tau, double *work, long *info); /* -- LAPACK routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University October 31, 1992 Purpose ======= NUMlapack_dgehd2 reduces a real general matrix A to upper Hessenberg form H by an orthogonal similarity transformation: Q' * A * Q = H . Arguments ========= N (input) long The order of the matrix A. N >= 0. ILO (input) long IHI (input) long 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 NUMlapack_dgebal; otherwise they should be set to 1 and N respectively. See Further Details. 1 <= ILO <= IHI <= max(1,N). A (input/output) double 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) long The leading dimension of the array A. LDA >= max(1,N). TAU (output) double array, dimension (N-1) The scalar factors of the elementary reflectors (see Further Details). WORK (workspace) double array, dimension (N) INFO (output) long = 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). ===================================================================== */ int NUMlapack_dgehrd (long *n, long *ilo, long *ihi, double *a, long *lda, double *tau, double *work, long *lwork, long *info); /* -- LAPACK routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University June 30, 1999 Purpose ======= NUMlapack_dgehrd reduces a real general matrix A to upper Hessenberg form H by an orthogonal similarity transformation: Q' * A * Q = H . Arguments ========= N (input) long The order of the matrix A. N >= 0. ILO (input) long IHI (input) long 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 NUMlapack_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 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) long The leading dimension of the array A. LDA >= max(1,N). TAU (output) double 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 array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) long 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) long = 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). ===================================================================== */ int NUMlapack_dgelq2 (long *m, long *n, double *a, long *lda, double *tau, double *work, long *info); /* Purpose ======= NUMlapack_dgelq2 computes an LQ factorization of a real m by n matrix A: A = L * Q. Arguments ========= M (input) long The number of rows of the matrix A. M >= 0. N (input) long The number of columns of the matrix A. N >= 0. A (input/output) double 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) long The leading dimension of the array A. LDA >= max(1,M). TAU (output) double array, dimension (min(M,N)) The scalar factors of the elementary reflectors (see Further Details). WORK (workspace) double array, dimension (M) INFO (output) long = 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). ===================================================================== */ int NUMlapack_dgelqf (long *m, long *n, double *a, long *lda, double *tau, double *work, long *lwork, long *info); /* Purpose ======= NUMlapack_dgelqf computes an LQ factorization of a real M-by-N matrix A: A = L * Q. Arguments ========= M (input) long The number of rows of the matrix A. M >= 0. N (input) long The number of columns of the matrix A. N >= 0. A (input/output) double 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) long The leading dimension of the array A. LDA >= max(1,M). TAU (output) double array, dimension (min(M,N)) The scalar factors of the elementary reflectors (see Further Details). WORK (workspace/output) double array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) long 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) long = 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). ===================================================================== */ int NUMlapack_dgelss (long *m, long *n, long *nrhs, double *a, long *lda, double *b, long *ldb, double *s, double *rcond, long *rank, double *work, long *lwork, long *info); /* Purpose ======= NUMlapack_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) long The number of rows of the matrix A. M >= 0. N (input) long The number of columns of the matrix A. N >= 0. NRHS (input) long The number of right hand sides, i.e., the number of columns of the matrices B and X. NRHS >= 0. A (input/output) double 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) long The leading dimension of the array A. LDA >= max(1,M). B (input/output) double 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) long The leading dimension of the array B. LDB >= max(1,max(M,N)). S (output) double 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 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) long The effective rank of A, i.e., the number of singular values which are greater than RCOND*S(1). WORK (workspace/output) double array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) long 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) long = 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. ===================================================================== */ int NUMlapack_dgeqpf (long *m, long *n, double *a, long *lda, long *jpvt, double *tau, double *work, long *info); /* Purpose ======= This routine is deprecated and has been replaced by routine DGEQP3. NUMlapack_dgeqpf computes a QR factorization with column pivoting of a real M-by-N matrix A: A*P = Q*R. Arguments ========= M (input) long The number of rows of the matrix A. M >= 0. N (input) long The number of columns of the matrix A. N >= 0 A (input/output) double 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) long The leading dimension of the array A. LDA >= max(1,M). JPVT (input/output) long 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 array, dimension (min(M,N)) The scalar factors of the elementary reflectors. WORK (workspace) double array, dimension (3*N) INFO (output) long = 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. ===================================================================== */ int NUMlapack_dgeqr2 (long *m, long *n, double *a, long *lda, double *tau, double *work, long *info); /* Purpose ======= NUMlapack_dgeqr2 computes a QR factorization of a real m by n matrix A: A = Q * R. Arguments ========= M (input) long The number of rows of the matrix A. M >= 0. N (input) long The number of columns of the matrix A. N >= 0. A (input/output) double 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) long The leading dimension of the array A. LDA >= max(1,M). TAU (output) double array, dimension (min(M,N)) The scalar factors of the elementary reflectors (see Further Details). WORK (workspace) double array, dimension (N) INFO (output) long = 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). ===================================================================== */ int NUMlapack_dgeqrf(long *m, long *n, double *a, long *lda, double *tau, double *work, long *lwork, long *info); /* Purpose ======= NUMlapack_dgeqrf computes a QR factorization of a real M-by-N matrix A: A = Q * R. Arguments ========= M (input) long The number of rows of the matrix A. M >= 0. N (input) long The number of columns of the matrix A. N >= 0. A (input/output) double 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) long The leading dimension of the array A. LDA >= max(1,M). TAU (output) double array, dimension (min(M,N)) The scalar factors of the elementary reflectors (see Further Details). WORK (workspace/output) double array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) long 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) long = 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). ===================================================================== */ int NUMlapack_dgerq2(long *m, long *n, double *a, long *lda, double *tau, double *work, long *info); /* Purpose ======= NUMlapack_dgerq2 computes an RQ factorization of a real m by n matrix A: A = R * Q. Arguments ========= M (input) long The number of rows of the matrix A. M >= 0. N (input) long The number of columns of the matrix A. N >= 0. A (input/output) double 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) long The leading dimension of the array A. LDA >= max(1,M). TAU (output) double array, dimension (min(M,N)) The scalar factors of the elementary reflectors (see Further Details). WORK (workspace) double array, dimension (M) INFO (output) long = 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). ===================================================================== */ int NUMlapack_dgesv (long *n, long *nrhs, double *a, long *lda, long *ipiv, double *b, long *ldb, long *info); /* Purpose ======= NUMlapack_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) long The number of linear equations, i.e., the order of the matrix A. N >= 0. NRHS (input) long The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0. A (input/output) double 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) long The leading dimension of the array A. LDA >= max(1,N). IPIV (output) long 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 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) long The leading dimension of the array B. LDB >= max(1,N). INFO (output) long = 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. ===================================================================== */ int NUMlapack_dgesvd (char *jobu, char *jobvt, long *m, long *n, double *a, long *lda, double *s, double *u, long *ldu, double *vt, long *ldvt, double *work, long *lwork, long *info); /* Purpose ======= NUMlapack_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) char* 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) char* 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) long The number of rows of the input matrix A. M >= 0. N (input) long The number of columns of the input matrix A. N >= 0. A (input/output) double 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) long The leading dimension of the array A. LDA >= max(1,M). S (output) double array, dimension (min(M,N)) The singular values of A, sorted so that S(i) >= S(i+1). U (output) double 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) long The leading dimension of the array U. LDU >= 1; if JOBU = 'S' or 'A', LDU >= M. VT (output) double 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) long 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 array, dimension (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) long The dimension of the array WORK. LWORK >= 1. LWORK >= MAX(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) long = 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. ===================================================================== */ int NUMlapack_dgetf2 (long *m, long *n, double *a, long *lda, long *ipiv, long *info); /* Purpose ======= NUMlapack_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) long The number of rows of the matrix A. M >= 0. N (input) long The number of columns of the matrix A. N >= 0. A (input/output) double 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) long The leading dimension of the array A. LDA >= max(1,M). IPIV (output) long 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) long = 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. ===================================================================== */ int NUMlapack_dgetri (long *n, double *a, long *lda, long *ipiv, double *work, long *lwork, long *info); /* Purpose ======= NUMlapack_dgetri computes the inverse of a matrix using the LU factorization computed by NUMlapack_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) long The order of the matrix A. N >= 0. A (input/output) double 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) long The leading dimension of the array A. LDA >= max(1,N). IPIV (input) long 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 array, dimension (LWORK) On exit, if INFO=0, then WORK(1) returns the optimal LWORK. LWORK (input) long The dimension of the array WORK. LWORK >= max(1,N). For optimal performance LWORK >= N*NB, where NB is the optimal blocksize returned by NUMlapack_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) long = 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. ===================================================================== */ int NUMlapack_dgetrf (long *m, long *n, double *a, long *lda, long *ipiv, long *info); /* Purpose ======= NUMlapack_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) long The number of rows of the matrix A. M >= 0. N (input) long The number of columns of the matrix A. N >= 0. A (input/output) double 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) long The leading dimension of the array A. LDA >= max(1,M). IPIV (output) long 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) long = 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. ===================================================================== */ int NUMlapack_dgetrs (char *trans, long *n, long *nrhs, double *a, long *lda, long *ipiv, double *b, long *ldb, long *info); /* Purpose ======= NUMlapack_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) char* 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) long The order of the matrix A. N >= 0. NRHS (input) long The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0. A (input) double array, dimension (LDA,N) The factors L and U from the factorization A = P*L*U as computed by DGETRF. LDA (input) long The leading dimension of the array A. LDA >= max(1,N). IPIV (input) long 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 array, dimension (LDB,NRHS) On entry, the right hand side matrix B. On exit, the solution matrix X. LDB (input) long The leading dimension of the array B. LDB >= max(1,N). INFO (output) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value ===================================================================== */ int NUMlapack_dggsvd (char *jobu, char *jobv, char *jobq, long *m, long *n, long *p, long *k, long *l, double *a, long *lda, double *b, long *ldb, double *alpha, double *beta, double *u, long *ldu, double *v, long *ldv, double *q, long *ldq, double *work, long *iwork, long *info); /* Purpose ======= NUMlapack_dggsvd computes the generalized singular value decomposition (GSVD) of an M-by-N real matrix A and P-by-N real matrix B: U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R ) where U, V and Q are orthogonal matrices, and Z' is the transpose of Z. Let K+L = the effective numerical rank of the matrix (A',B')', then R is a K+L-by-K+L nonsingular upper triangular matrix, D1 and D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and of the following structures, respectively: If M-K-L >= 0, K L D1 = K ( I 0 ) L ( 0 C ) M-K-L ( 0 0 ) K L D2 = L ( 0 S ) P-L ( 0 0 ) N-K-L K L ( 0 R ) = K ( 0 R11 R12 ) L ( 0 0 R22 ) where C = diag( ALPHA(K+1), ... , ALPHA(K+L) ), S = diag( BETA(K+1), ... , BETA(K+L) ), C**2 + S**2 = I. R is stored in A(1:K+L,N-K-L+1:N) on exit. If M-K-L < 0, K M-K K+L-M D1 = K ( I 0 0 ) M-K ( 0 C 0 ) K M-K K+L-M D2 = M-K ( 0 S 0 ) K+L-M ( 0 0 I ) P-L ( 0 0 0 ) N-K-L K M-K K+L-M ( 0 R ) = K ( 0 R11 R12 R13 ) M-K ( 0 0 R22 R23 ) K+L-M ( 0 0 0 R33 ) where C = diag( ALPHA(K+1), ... , ALPHA(M) ), S = diag( BETA(K+1), ... , BETA(M) ), C**2 + S**2 = I. (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored ( 0 R22 R23 ) in B(M-K+1:L,N+M-K-L+1:N) on exit. The routine computes C, S, R, and optionally the orthogonal transformation matrices U, V and Q. In particular, if B is an N-by-N nonsingular matrix, then the GSVD of A and B implicitly gives the SVD of A*inv(B): A*inv(B) = U*(D1*inv(D2))*V'. If ( A',B')' has orthonormal columns, then the GSVD of A and B is also equal to the CS decomposition of A and B. Furthermore, the GSVD can be used to derive the solution of the eigenvalue problem: A'*A x = lambda* B'*B x. In some literature, the GSVD of A and B is presented in the form U'*A*X = ( 0 D1 ), V'*B*X = ( 0 D2 ) where U and V are orthogonal and X is nonsingular, D1 and D2 are ``diagonal''. The former GSVD form can be converted to the latter form by taking the nonsingular matrix X as X = Q*( I 0 ) ( 0 inv(R) ). Arguments ========= JOBU (input) char* = 'U': Orthogonal matrix U is computed; = 'N': U is not computed. JOBV (input) char* = 'V': Orthogonal matrix V is computed; = 'N': V is not computed. JOBQ (input) char* = 'Q': Orthogonal matrix Q is computed; = 'N': Q is not computed. M (input) long The number of rows of the matrix A. M >= 0. N (input) long The number of columns of the matrices A and B. N >= 0. P (input) long The number of rows of the matrix B. P >= 0. K (output) long L (output) long On exit, K and L specify the dimension of the subblocks described in the Purpose section. K + L = effective numerical rank of (A',B')'. A (input/output) double array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, A contains the triangular matrix R, or part of R. See Purpose for details. LDA (input) long The leading dimension of the array A. LDA >= max(1,M). B (input/output) double array, dimension (LDB,N) On entry, the P-by-N matrix B. On exit, B contains the triangular matrix R if M-K-L < 0. See Purpose for details. LDB (input) long The leading dimension of the array B. LDA >= max(1,P). ALPHA (output) double array, dimension (N) BETA (output) double array, dimension (N) On exit, ALPHA and BETA contain the generalized singular value pairs of A and B; ALPHA(1:K) = 1, BETA(1:K) = 0, and if M-K-L >= 0, ALPHA(K+1:K+L) = C, BETA(K+1:K+L) = S, or if M-K-L < 0, ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0 BETA(K+1:M) =S, BETA(M+1:K+L) =1 and ALPHA(K+L+1:N) = 0 BETA(K+L+1:N) = 0 U (output) double array, dimension (LDU,M) If JOBU = 'U', U contains the M-by-M orthogonal matrix U. If JOBU = 'N', U is not referenced. LDU (input) long The leading dimension of the array U. LDU >= max(1,M) if JOBU = 'U'; LDU >= 1 otherwise. V (output) double array, dimension (LDV,P) If JOBV = 'V', V contains the P-by-P orthogonal matrix V. If JOBV = 'N', V is not referenced. LDV (input) long The leading dimension of the array V. LDV >= max(1,P) if JOBV = 'V'; LDV >= 1 otherwise. Q (output) double array, dimension (LDQ,N) If JOBQ = 'Q', Q contains the N-by-N orthogonal matrix Q. If JOBQ = 'N', Q is not referenced. LDQ (input) long The leading dimension of the array Q. LDQ >= max(1,N) if JOBQ = 'Q'; LDQ >= 1 otherwise. WORK (workspace) double array, dimension (max(3*N,M,P)+N) IWORK (workspace/output) long array, dimension (N) On exit, IWORK stores the sorting information. More precisely, the following loop will sort ALPHA for I = K+1, min(M,K+L) swap ALPHA(I) and ALPHA(IWORK(I)) endfor such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N). INFO (output)long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value. > 0: if INFO = 1, the Jacobi-type procedure failed to converge. For further details, see subroutine NUMlapack_dtgsja. Internal Parameters =================== TOLA double TOLB double TOLA and TOLB are the thresholds to determine the effective rank of (A',B')'. Generally, they are set to TOLA = MAX(M,N)*norm(A)*MAZHEPS, TOLB = MAX(P,N)*norm(B)*MAZHEPS. The size of TOLA and TOLB may affect the size of backward errors of the decomposition. ===================================================================== */ int NUMlapack_dggsvp (char *jobu, char *jobv, char *jobq, long *m, long *p, long *n, double *a, long *lda, double *b, long *ldb, double *tola, double *tolb, long *k, long *l, double *u, long *ldu, double *v, long *ldv, double *q, long *ldq, long *iwork, double *tau, double *work, long *info); /* Purpose ======= NUMlapack_dggsvp computes orthogonal matrices U, V and Q such that N-K-L K L U'*A*Q = K ( 0 A12 A13 ) if M-K-L >= 0; L ( 0 0 A23 ) M-K-L ( 0 0 0 ) N-K-L K L = K ( 0 A12 A13 ) if M-K-L < 0; M-K ( 0 0 A23 ) N-K-L K L V'*B*Q = L ( 0 0 B13 ) P-L ( 0 0 0 ) where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0, otherwise A23 is (M-K)-by-L upper trapezoidal. K+L = the effective numerical rank of the (M+P)-by-N matrix (A',B')'. Z' denotes the transpose of Z. This decomposition is the preprocessing step for computing the Generalized Singular Value Decomposition (GSVD), see subroutine NUMlapack_dggsvd. Arguments ========= JOBU (input) char* = 'U': Orthogonal matrix U is computed; = 'N': U is not computed. JOBV (input) char* = 'V': Orthogonal matrix V is computed; = 'N': V is not computed. JOBQ (input) char* = 'Q': Orthogonal matrix Q is computed; = 'N': Q is not computed. M (input) long The number of rows of the matrix A. M >= 0. P (input) long The number of rows of the matrix B. P >= 0. N (input) long The number of columns of the matrices A and B. N >= 0. A (input/output) double array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, A contains the triangular (or trapezoidal) matrix described in the Purpose section. LDA (input) long The leading dimension of the array A. LDA >= max(1,M). B (input/output) double array, dimension (LDB,N) On entry, the P-by-N matrix B. On exit, B contains the triangular matrix described in the Purpose section. LDB (input) long The leading dimension of the array B. LDB >= max(1,P). TOLA (input) double TOLB (input) double TOLA and TOLB are the thresholds to determine the effective numerical rank of matrix B and a subblock of A. Generally, they are set to TOLA = MAX(M,N)*norm(A)*MAZHEPS, TOLB = MAX(P,N)*norm(B)*MAZHEPS. The size of TOLA and TOLB may affect the size of backward errors of the decomposition. K (output) long L (output) long On exit, K and L specify the dimension of the subblocks described in Purpose. K + L = effective numerical rank of (A',B')'. U (output) double array, dimension (LDU,M) If JOBU = 'U', U contains the orthogonal matrix U. If JOBU = 'N', U is not referenced. LDU (input) long The leading dimension of the array U. LDU >= max(1,M) if JOBU = 'U'; LDU >= 1 otherwise. V (output) double array, dimension (LDV,M) If JOBV = 'V', V contains the orthogonal matrix V. If JOBV = 'N', V is not referenced. LDV (input) long The leading dimension of the array V. LDV >= max(1,P) if JOBV = 'V'; LDV >= 1 otherwise. Q (output) double array, dimension (LDQ,N) If JOBQ = 'Q', Q contains the orthogonal matrix Q. If JOBQ = 'N', Q is not referenced. LDQ (input) long The leading dimension of the array Q. LDQ >= max(1,N) if JOBQ = 'Q'; LDQ >= 1 otherwise. IWORK (workspace) long array, dimension (N) TAU (workspace) double array, dimension (N) WORK (workspace) double array, dimension (max(3*N,M,P)) INFO (output) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value. Further Details =============== The subroutine uses LAPACK subroutine DGEQPF for the QR factorization with column pivoting to detect the effective numerical rank of the a matrix. It may be replaced by a better rank determination strategy. ===================================================================== */ int NUMlapack_dhseqr (char *job, char *compz, long *n, long *ilo, long *ihi, double *h, long *ldh, double *wr, double *wi, double *z, long *ldz, double *work, long *lwork, long *info); /* -- LAPACK routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University June 30, 1999 Purpose ======= NUMlapack_dhseqr computes the eigenvalues of a real upper Hessenberg matrix H and, optionally, the matrices T and Z from the Schur decomposition H = Z T Z**T, where T is an upper quasi-triangular matrix (the Schur form), and Z is the orthogonal matrix of Schur vectors. Optionally Z may be postmultiplied into an input orthogonal matrix Q, so that this routine can give the Schur factorization of a matrix A which has been reduced to the Hessenberg form H by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T. Arguments ========= JOB (input) char* = 'E': compute eigenvalues only; = 'S': compute eigenvalues and the Schur form T. COMPZ (input) char* = 'N': no Schur vectors are computed; = 'I': Z is initialized to the unit matrix and the matrix Z of Schur vectors of H is returned; = 'V': Z must contain an orthogonal matrix Q on entry, and the product Q*Z is returned. N (input) long The order of the matrix H. N >= 0. ILO (input) long IHI (input) long It is assumed that H 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 NUMlapack_dgebal, and then passed to SGEHRD when the matrix output by NUMlapack_dgebal is reduced to Hessenberg form. Otherwise ILO and IHI should be set to 1 and N respectively. 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. H (input/output) double array, dimension (LDH,N) On entry, the upper Hessenberg matrix H. On exit, if JOB = 'S', H contains the upper quasi-triangular matrix T from the Schur decomposition (the Schur form); 2-by-2 diagonal blocks (corresponding to complex conjugate pairs of eigenvalues) are returned in standard form, with H(i,i) = H(i+1,i+1) and H(i+1,i)*H(i,i+1) < 0. If JOB = 'E', the contents of H are unspecified on exit. LDH (input) long The leading dimension of the array H. LDH >= max(1,N). WR (output) double array, dimension (N) WI (output) double array, dimension (N) The real and imaginary parts, respectively, of the computed eigenvalues. If two eigenvalues are computed as a complex conjugate pair, they are stored in consecutive elements of WR and WI, say the i-th and (i+1)th, with WI(i) > 0 and WI(i+1) < 0. If JOB = 'S', the eigenvalues are stored in the same order as on the diagonal of the Schur form returned in H, with WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i). Z (input/output) double array, dimension (LDZ,N) If COMPZ = 'N': Z is not referenced. If COMPZ = 'I': on entry, Z need not be set, and on exit, Z contains the orthogonal matrix Z of the Schur vectors of H. If COMPZ = 'V': on entry Z must contain an N-by-N matrix Q, which is assumed to be equal to the unit matrix except for the submatrix Z(ILO:IHI,ILO:IHI); on exit Z contains Q*Z. Normally Q is the orthogonal matrix generated by NUMlapack_dorghr after the call to NUMlapack_dgehrd which formed the Hessenberg matrix H. LDZ (input) long The leading dimension of the array Z. LDZ >= max(1,N) if COMPZ = 'I' or 'V'; LDZ >= 1 otherwise. WORK (workspace/output) double array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) long The dimension of the array WORK. LWORK >= max(1,N). 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) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, NUMlapack_dhseqr failed to compute all of the eigenvalues in a total of 30*(IHI-ILO+1) iterations; elements 1:ilo-1 and i+1:n of WR and WI contain those eigenvalues which have been successfully computed. ===================================================================== */ int NUMlapack_dlabad (double *smal, double *large); /* Purpose ======= NUMlapack_dlabad takes as input the values computed by DLAMCH for underflow and overflow, and returns the square root of each of these values if the log of LARGE is sufficiently large. This subroutine is intended to identify machines with a large exponent range, such as the Crays, and redefine the underflow and overflow limits to be the square roots of the values computed by DLAMCH. This subroutine is needed because DLAMCH does not compensate for poor arithmetic in the upper half of the exponent range, as is found on a Cray. Arguments ========= smal (input/output) double On entry, the underflow threshold as computed by DLAMCH. On exit, if LOG10(LARGE) is sufficiently large, the square root of smal, otherwise unchanged. LARGE (input/output) double On entry, the overflow threshold as computed by DLAMCH. On exit, if LOG10(LARGE) is sufficiently large, the square root of LARGE, otherwise unchanged. ===================================================================== If it looks like we're on a Cray, take the square root of smal and LARGE to avoid overflow and underflow problems. */ int NUMlapack_dlabrd (long *m, long *n, long *nb, double *a, long *lda, double *d, double *e, double *tauq, double *taup, double *x, long *ldx, double *y, long *ldy); /* Purpose ======= NUMlapack_dlabrd reduces the first NB rows and columns of a real general m by n matrix A to upper or lower bidiagonal form by an orthogonal transformation Q' * A * P, and returns the matrices X and Y which are needed to apply the transformation to the unreduced part of A. If m >= n, A is reduced to upper bidiagonal form; if m < n, to lower bidiagonal form. This is an auxiliary routine called by DGEBRD Arguments ========= M (input) long The number of rows in the matrix A. N (input) long The number of columns in the matrix A. NB (input) long The number of leading rows and columns of A to be reduced. A (input/output) double array, dimension (LDA,N) On entry, the m by n general matrix to be reduced. On exit, the first NB rows and columns of the matrix are overwritten; the rest of the array is unchanged. If m >= n, elements on and below the diagonal in the first NB columns, with the array TAUQ, represent the orthogonal matrix Q as a product of elementary reflectors; and elements above the diagonal in the first NB rows, with the array TAUP, represent the orthogonal matrix P as a product of elementary reflectors. If m < n, elements below the diagonal in the first NB columns, with the array TAUQ, represent the orthogonal matrix Q as a product of elementary reflectors, and elements on and above the diagonal in the first NB rows, with the array TAUP, represent the orthogonal matrix P as a product of elementary reflectors. See Further Details. LDA (input) long The leading dimension of the array A. LDA >= max(1,M). D (output) double array, dimension (NB) The diagonal elements of the first NB rows and columns of the reduced matrix. D(i) = A(i,i). E (output) double array, dimension (NB) The off-diagonal elements of the first NB rows and columns of the reduced matrix. TAUQ (output) double array dimension (NB) The scalar factors of the elementary reflectors which represent the orthogonal matrix Q. See Further Details. TAUP (output) double array, dimension (NB) The scalar factors of the elementary reflectors which represent the orthogonal matrix P. See Further Details. X (output) double array, dimension (LDX,NB) The m-by-nb matrix X required to update the unreduced part of A. LDX (input) long The leading dimension of the array X. LDX >= M. Y (output) double array, dimension (LDY,NB) The n-by-nb matrix Y required to update the unreduced part of A. LDY (output) long The leading dimension of the array Y. LDY >= N. Further Details =============== The matrices Q and P are represented as products of elementary reflectors: Q = H(1) H(2) . . . H(nb) and P = G(1) G(2) . . . G(nb) 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. If m >= n, v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in A(i:m,i); u(1:i) = 0, u(i+1) = 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). If m < n, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i:n) is stored on exit in A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i). The elements of the vectors v and u together form the m-by-nb matrix V and the nb-by-n matrix U' which are needed, with X and Y, to apply the transformation to the unreduced part of the matrix, using a block update of the form: A := A - V*Y' - X*U'. The contents of A on exit are illustrated by the following examples with nb = 2: m = 6 and n = 5 (m > n): m = 5 and n = 6 (m < n): ( 1 1 u1 u1 u1 ) ( 1 u1 u1 u1 u1 u1 ) ( v1 1 1 u2 u2 ) ( 1 1 u2 u2 u2 u2 ) ( v1 v2 a a a ) ( v1 1 a a a a ) ( v1 v2 a a a ) ( v1 v2 a a a a ) ( v1 v2 a a a ) ( v1 v2 a a a a ) ( v1 v2 a a a ) where a denotes an element of the original matrix which is unchanged, vi denotes an element of the vector defining H(i), and ui an element of the vector defining G(i). ===================================================================== */ int NUMlapack_dlacpy (char *uplo, long *m, long *n, double *a, long *lda, double *b, long *ldb); /* Purpose ======= NUMlapack_dlacpy copies all or part of a two-dimensional matrix A to another matrix B. Arguments ========= UPLO (input) char* Specifies the part of the matrix A to be copied to B. = 'U': Upper triangular part = 'L': Lower triangular part Otherwise: All of the matrix A M (input) long The number of rows of the matrix A. M >= 0. N (input) long The number of columns of the matrix A. N >= 0. A (input) double array, dimension (LDA,N) The m by n matrix A. If UPLO = 'U', only the upper triangle or trapezoid is accessed; if UPLO = 'L', only the lower triangle or trapezoid is accessed. LDA (input) long The leading dimension of the array A. LDA >= max(1,M). B (output) double array, dimension (LDB,N) On exit, B = A in the locations specified by UPLO. LDB (input) long The leading dimension of the array B. LDB >= max(1,M). ===================================================================== */ int NUMlapack_dladiv (double *a, double *b, double *c, double *d, double *p, double *q); /* -- LAPACK auxiliary routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University October 31, 1992 Purpose ======= NUMlapack_dladiv performs complex division in real arithmetic a + i*b p + i*q = --------- c + i*d The algorithm is due to Robert L. Smith and can be found in D. Knuth, The art of Computer Programming, Vol.2, p.195 Arguments ========= A (input) double B (input) double C (input) double D (input) double The scalars a, b, c, and d in the above expression. P (output) double Q (output) double The scalars p and q in the above expression. ===================================================================== */ int NUMlapack_dlae2 (double *a, double *b, double *c, double *rt1, double *rt2); /* Purpose ======= NUMlapack_dlae2 computes the eigenvalues of a 2-by-2 symmetric matrix [ A B ] [ B C ]. On return, RT1 is the eigenvalue of larger absolute value, and RT2 is the eigenvalue of smaller absolute value. Arguments ========= A (input) double The (1,1) element of the 2-by-2 matrix. B (input) double The (1,2) and (2,1) elements of the 2-by-2 matrix. C (input) double The (2,2) element of the 2-by-2 matrix. RT1 (output) double The eigenvalue of larger absolute value. RT2 (output) double The eigenvalue of smaller absolute value. Further Details =============== RT1 is accurate to a few ulps barring over/underflow. RT2 may be inaccurate if there is massive cancellation in the determinant A*C-B*B; higher precision or correctly rounded or correctly truncated arithmetic would be needed to compute RT2 accurately in all cases. Overflow is possible only if RT1 is within a factor of 5 of overflow. Underflow is harmless if the input data is 0 or exceeds underflow_threshold / macheps. ===================================================================== */ int NUMlapack_dlaev2 (double *a, double *b, double *c, double *rt1, double *rt2, double *cs1, double *sn1); /* Purpose ======= NUMlapack_dlaev2 computes the eigendecomposition of a 2-by-2 symmetric matrix [ A B ] [ B C ]. On return, RT1 is the eigenvalue of larger absolute value, RT2 is the eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right eigenvector for RT1, giving the decomposition [ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ] [-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ]. Arguments ========= A (input) double The (1,1) element of the 2-by-2 matrix. B (input) double The (1,2) element and the conjugate of the (2,1) element of the 2-by-2 matrix. C (input) double The (2,2) element of the 2-by-2 matrix. RT1 (output) double The eigenvalue of larger absolute value. RT2 (output) double The eigenvalue of smaller absolute value. CS1 (output) double SN1 (output) double The vector (CS1, SN1) is a unit right eigenvector for RT1. Further Details =============== RT1 is accurate to a few ulps barring over/underflow. RT2 may be inaccurate if there is massive cancellation in the determinant A*C-B*B; higher precision or correctly rounded or correctly truncated arithmetic would be needed to compute RT2 accurately in all cases. CS1 and SN1 are accurate to a few ulps barring over/underflow. Overflow is possible only if RT1 is within a factor of 5 of overflow. Underflow is harmless if the input data is 0 or exceeds underflow_threshold / macheps. ===================================================================== */ int NUMlapack_dlags2 (long *upper, double *a1, double *a2, double *a3, double *b1, double *b2, double *b3, double *csu, double *snu, double *csv, double *snv, double *csq, double *snq); /* Purpose ======= NUMlapack_dlags2 computes 2-by-2 orthogonal matrices U, V and Q, such that if ( UPPER ) then U'*A*Q = U'*( A1 A2 )*Q = ( x 0 ) ( 0 A3 ) ( x x ) and V'*B*Q = V'*( B1 B2 )*Q = ( x 0 ) ( 0 B3 ) ( x x ) or if ( .NOT.UPPER ) then U'*A*Q = U'*( A1 0 )*Q = ( x x ) ( A2 A3 ) ( 0 x ) and V'*B*Q = V'*( B1 0 )*Q = ( x x ) ( B2 B3 ) ( 0 x ) The rows of the transformed A and B are parallel, where U = ( CSU SNU ), V = ( CSV SNV ), Q = ( CSQ SNQ ) ( -SNU CSU ) ( -SNV CSV ) ( -SNQ CSQ ) Z' denotes the transpose of Z. Arguments ========= UPPER (input) long* (boolean) = TRUE: the input matrices A and B are upper triangular. = FALSE: the input matrices A and B are lower triangular. A1 (input) double A2 (input) double A3 (input) double On entry, A1, A2 and A3 are elements of the input 2-by-2 upper (lower) triangular matrix A. B1 (input) double B2 (input) double B3 (input) double On entry, B1, B2 and B3 are elements of the input 2-by-2 upper (lower) triangular matrix B. CSU (output) double SNU (output) double The desired orthogonal matrix U. CSV (output) double SNV (output) double The desired orthogonal matrix V. CSQ (output) double SNQ (output) double The desired orthogonal matrix Q. ===================================================================== */ int NUMlapack_dlahqr (int * wantt, int * wantz, long *n, long *ilo, long *ihi, double *h, long *ldh, double *wr, double *wi, long *iloz, long *ihiz, double *z, long *ldz, long *info); /* -- LAPACK auxiliary routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University June 30, 1999 Purpose ======= NUMlapack_dlahqr is an auxiliary routine called by NUMlapack_dhseqr to update the eigenvalues and Schur decomposition already computed by NUMlapack_dhseqr, by dealing with the Hessenberg submatrix in rows and columns ILO to IHI. Arguments ========= WANTT (input) int = .TRUE. : the full Schur form T is required; = .FALSE.: only eigenvalues are required. WANTZ (input) int = .TRUE. : the matrix of Schur vectors Z is required; = .FALSE.: Schur vectors are not required. N (input) long The order of the matrix H. N >= 0. ILO (input) long IHI (input) long It is assumed that H is already upper quasi-triangular in rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1). NUMlapack_dlahqr works primarily with the Hessenberg submatrix in rows and columns ILO to IHI, but applies transformations to all of H if WANTT is .TRUE.. 1 <= ILO <= max(1,IHI); IHI <= N. H (input/output) double array, dimension (LDH,N) On entry, the upper Hessenberg matrix H. On exit, if WANTT is .TRUE., H is upper quasi-triangular in rows and columns ILO:IHI, with any 2-by-2 diagonal blocks in standard form. If WANTT is .FALSE., the contents of H are unspecified on exit. LDH (input) long The leading dimension of the array H. LDH >= max(1,N). WR (output) double array, dimension (N) WI (output) double array, dimension (N) The real and imaginary parts, respectively, of the computed eigenvalues ILO to IHI are stored in the corresponding elements of WR and WI. If two eigenvalues are computed as a complex conjugate pair, they are stored in consecutive elements of WR and WI, say the i-th and (i+1)th, with WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the eigenvalues are stored in the same order as on the diagonal of the Schur form returned in H, with WR(i) = H(i,i), and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i). ILOZ (input) long IHIZ (input) long Specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 <= ILOZ <= ILO; IHI <= IHIZ <= N. Z (input/output) double array, dimension (LDZ,N) If WANTZ is .TRUE., on entry Z must contain the current matrix Z of transformations accumulated by NUMlapack_dhseqr, and on exit Z has been updated; transformations are applied only to the submatrix Z(ILOZ:IHIZ,ILO:IHI). If WANTZ is .FALSE., Z is not referenced. LDZ (input) long The leading dimension of the array Z. LDZ >= max(1,N). INFO (output) long = 0: successful exit > 0: NUMlapack_dlahqr failed to compute all the eigenvalues ILO to IHI in a total of 30*(IHI-ILO+1) iterations; if INFO = i, elements i+1:ihi of WR and WI contain those eigenvalues which have been successfully computed. Further Details =============== 2-96 Based on modifications by David Day, Sandia National Laboratory, USA ===================================================================== */ int NUMlapack_dlahrd (long *n, long *k, long *nb, double *a, long *lda, double *tau, double *t, long *ldt, double *y, long *ldy); /* -- LAPACK auxiliary routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University June 30, 1999 Purpose ======= NUMlapack_dlahrd reduces the first NB columns of a real general n-by-(n-k+1) matrix A so that elements below the k-th subdiagonal are zero. The reduction is performed by an orthogonal similarity transformation Q' * A * Q. The routine returns the matrices V and T which determine Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T. This is an auxiliary routine called by NUMlapack_dgehrd. Arguments ========= N (input) long The order of the matrix A. K (input) long The offset for the reduction. Elements below the k-th subdiagonal in the first NB columns are reduced to zero. NB (input) long The number of columns to be reduced. A (input/output) double array, dimension (LDA,N-K+1) On entry, the n-by-(n-k+1) general matrix A. On exit, the elements on and above the k-th subdiagonal in the first NB columns are overwritten with the corresponding elements of the reduced matrix; the elements below the k-th subdiagonal, with the array TAU, represent the matrix Q as a product of elementary reflectors. The other columns of A are unchanged. See Further Details. LDA (input) long The leading dimension of the array A. LDA >= max(1,N). TAU (output) double array, dimension (NB) The scalar factors of the elementary reflectors. See Further Details. T (output) double array, dimension (LDT,NB) The upper triangular matrix T. LDT (input) long The leading dimension of the array T. LDT >= NB. Y (output) double array, dimension (LDY,NB) The n-by-nb matrix Y. LDY (input) long The leading dimension of the array Y. LDY >= N. Further Details =============== The matrix Q is represented as a product of nb elementary reflectors Q = H(1) H(2) . . . H(nb). 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+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in A(i+k+1:n,i), and tau in TAU(i). The elements of the vectors v together form the (n-k+1)-by-nb matrix V which is needed, with T and Y, to apply the transformation to the unreduced part of the matrix, using an update of the form: A := (I - V*T*V') * (A - Y*V'). The contents of A on exit are illustrated by the following example with n = 7, k = 3 and nb = 2: ( a h a a a ) ( a h a a a ) ( a h a a a ) ( h h a a a ) ( v1 h a a a ) ( v1 v2 a a a ) ( v1 v2 a 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). ===================================================================== */ int NUMlapack_dlaln2 (int * ltrans, long *na, long *nw, double *smin, double *ca, double *a, long *lda, double *d1, double *d2, double *b, long *ldb, double *wr, double *wi, double *x, long *ldx, double *scale, double *xnorm, long *info); /* -- LAPACK auxiliary routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University October 31, 1992 Purpose ======= NUMlapack_dlaln2 solves a system of the form (ca A - w D ) X = s B or (ca A' - w D) X = s B with possible scaling ("s") and perturbation of A. (A' means A-transpose.) A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA real diagonal matrix, w is a real or complex value, and X and B are NA x 1 matrices -- real if w is real, complex if w is complex. NA may be 1 or 2. If w is complex, X and B are represented as NA x 2 matrices, the first column of each being the real part and the second being the imaginary part. "s" is a scaling factor (.LE. 1), computed by NUMlapack_dlaln2, which is so chosen that X can be computed without overflow. X is further scaled if necessary to assure that norm(ca A - w D)*norm(X) is less than overflow. If both singular values of (ca A - w D) are less than SMIN, SMIN*identity will be used instead of (ca A - w D). If only one singular value is less than SMIN, one element of (ca A - w D) will be perturbed enough to make the smallest singular value roughly SMIN. If both singular values are at least SMIN, (ca A - w D) will not be perturbed. In any case, the perturbation will be at most some small multiple of max( SMIN, ulp*norm(ca A - w D) ). The singular values are computed by infinity-norm approximations, and thus will only be correct to a factor of 2 or so. Note: all input quantities are assumed to be smaller than overflow by a reasonable factor. (See BIGNUM.) Arguments ========== LTRANS (input) int =.TRUE.: A-transpose will be used. =.FALSE.: A will be used (not transposed.) NA (input) long The size of the matrix A. It may (only) be 1 or 2. NW (input) long 1 if "w" is real, 2 if "w" is complex. It may only be 1 or 2. SMIN (input) double The desired lower bound on the singular values of A. This should be a safe distance away from underflow or overflow, say, between (underflow/machine precision) and (machine precision * overflow ). (See BIGNUM and ULP.) CA (input) double The coefficient c, which A is multiplied by. A (input) double array, dimension (LDA,NA) The NA x NA matrix A. LDA (input) long The leading dimension of A. It must be at least NA. D1 (input) double The 1,1 element in the diagonal matrix D. D2 (input) double The 2,2 element in the diagonal matrix D. Not used if NW=1. B (input) double array, dimension (LDB,NW) The NA x NW matrix B (right-hand side). If NW=2 ("w" is complex), column 1 contains the real part of B and column 2 contains the imaginary part. LDB (input) long The leading dimension of B. It must be at least NA. WR (input) double The real part of the scalar "w". WI (input) double The imaginary part of the scalar "w". Not used if NW=1. X (output) double array, dimension (LDX,NW) The NA x NW matrix X (unknowns), as computed by NUMlapack_dlaln2. If NW=2 ("w" is complex), on exit, column 1 will contain the real part of X and column 2 will contain the imaginary part. LDX (input) long The leading dimension of X. It must be at least NA. SCALE (output) double The scale factor that B must be multiplied by to insure that overflow does not occur when computing X. Thus, (ca A - w D) X will be SCALE*B, not B (ignoring perturbations of A.) It will be at most 1. XNORM (output) double The infinity-norm of X, when X is regarded as an NA x NW real matrix. INFO (output) long An error flag. It will be set to zero if no error occurs, a negative number if an argument is in error, or a positive number if ca A - w D had to be perturbed. The possible values are: = 0: No error occurred, and (ca A - w D) did not have to be perturbed. = 1: (ca A - w D) had to be perturbed to make its smallest (or only) singular value greater than SMIN. NOTE: In the interests of speed, this routine does not check the inputs for errors. ===================================================================== */ double NUMlapack_dlange (char *norm, long *m, long *n, double *a, long *lda, double *work); /* Purpose ======= NUMlapack_dlange returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real matrix A. Description =========== DLANGE returns the value DLANGE = ( max(abs(A(i,j))), NORM = 'M' or 'm' ( ( norm1(A), NORM = '1', 'O' or 'o' ( ( normI(A), NORM = 'I' or 'i' ( ( normF(A), NORM = 'F', 'f', 'E' or 'e' where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Note that max(abs(A(i,j))) is not a matrix norm. Arguments ========= NORM (input) char* Specifies the value to be returned in DLANGE as described above. M (input) long The number of rows of the matrix A. M >= 0. When M = 0, DLANGE is set to zero. N (input) long The number of columns of the matrix A. N >= 0. When N = 0, DLANGE is set to zero. A (input) double array, dimension (LDA,N) The m by n matrix A. LDA (input) long The leading dimension of the array A. LDA >= max(M,1). WORK (workspace) double array, dimension (LWORK), where LWORK >= M when NORM = 'I'; otherwise, WORK is not referenced. ===================================================================== */ double NUMlapack_dlanhs (char *norm, long *n, double *a, long *lda, double *work); /* -- LAPACK auxiliary routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University October 31, 1992 Purpose ======= NUMlapack_dlanhs returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a Hessenberg matrix A. Description =========== NUMlapack_dlanhs returns the value NUMlapack_dlanhs = ( max(abs(A(i,j))), NORM = 'M' or 'm' ( ( norm1(A), NORM = '1', 'O' or 'o' ( ( normI(A), NORM = 'I' or 'i' ( ( normF(A), NORM = 'F', 'f', 'E' or 'e' where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Note that max(abs(A(i,j))) is not a matrix norm. Arguments ========= NORM (input) char* Specifies the value to be returned in NUMlapack_dlanhs as described above. N (input) long The order of the matrix A. N >= 0. When N = 0, NUMlapack_dlanhs is set to zero. A (input) double array, dimension (LDA,N) The n by n upper Hessenberg matrix A; the part of A below the first sub-diagonal is not referenced. LDA (input) long The leading dimension of the array A. LDA >= max(N,1). WORK (workspace) double array, dimension (LWORK), where LWORK >= N when NORM = 'I'; otherwise, WORK is not referenced. ===================================================================== */ double NUMlapack_dlanst (char *norm, long *n, double *d, double *e); /* Purpose ======= NUMlapack_dlanst returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric tridiagonal matrix A. Description =========== NUMlapack_dlanst returns the value NUMlapack_dlanst = ( max(abs(A(i,j))), NORM = 'M' or 'm' ( ( norm1(A), NORM = '1', 'O' or 'o' ( ( normI(A), NORM = 'I' or 'i' ( ( normF(A), NORM = 'F', 'f', 'E' or 'e' where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Note that max(abs(A(i,j))) is not a matrix norm. Arguments ========= NORM (input) char* Specifies the value to be returned in DLANST as described above. N (input) long The order of the matrix A. N >= 0. When N = 0, DLANST is set to zero. D (input) double array, dimension (N) The diagonal elements of A. E (input) double array, dimension (N-1) The (n-1) sub-diagonal or super-diagonal elements of A. ===================================================================== */ double NUMlapack_dlansy (char *norm, char *uplo, long *n, double *a, long *lda, double *work); /* Purpose ======= NUMlapack_dlansy returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric matrix A. Description =========== NUMlapack_dlansy returns the value ( max(abs(A(i,j))), NORM = 'M' or 'm' ( ( norm1(A), NORM = '1', 'O' or 'o' ( ( normI(A), NORM = 'I' or 'i' ( ( normF(A), NORM = 'F', 'f', 'E' or 'e' where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Note that max(abs(A(i,j))) is not a matrix norm. Arguments ========= NORM (input) char* Specifies the value to be returned in DLANSY as described above. UPLO (input) char* Specifies whether the upper or lower triangular part of the symmetric matrix A is to be referenced. = 'U': Upper triangular part of A is referenced = 'L': Lower triangular part of A is referenced N (input) long The order of the matrix A. N >= 0. When N = 0, DLANSY is set to zero. A (input) double array, dimension (LDA,N) The symmetric matrix A. If UPLO = 'U', the leading n by n upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading n by n lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. LDA (input) long The leading dimension of the array A. LDA >= max(N,1). WORK (workspace) double array, dimension (LWORK), where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, WORK is not referenced. ===================================================================== */ int NUMlapack_dlanv2 (double *a, double *b, double *c, double *d, double *rt1r, double *rt1i, double *rt2r, double *rt2i, double *cs, double *sn); /* -- LAPACK driver routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University June 30, 1999 Purpose ======= NUMlapack_dlanv2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form: [ A B ] = [ CS -SN ] [ AA BB ] [ CS SN ] [ C D ] [ SN CS ] [ CC DD ] [-SN CS ] where either 1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or 2) AA = DD and BB*CC < 0, so that AA + or - sqrt(BB*CC) are complex conjugate eigenvalues. Arguments ========= A (input/output) double B (input/output) double C (input/output) double D (input/output) double On entry, the elements of the input matrix. On exit, they are overwritten by the elements of the standardised Schur form. RT1R (output) double RT1I (output) double RT2R (output) double RT2I (output) double The real and imaginary parts of the eigenvalues. If the eigenvalues are a complex conjugate pair, RT1I > 0. CS (output) double SN (output) double Parameters of the rotation matrix. Further Details =============== Modified by V. Sima, Research Institute for Informatics, Bucharest, Romania, to reduce the risk of cancellation errors, when computing real eigenvalues, and to ensure, if possible, that abs(RT1R) >= abs(RT2R). ===================================================================== */ int NUMlapack_dlapll(long *n, double *x, long *incx, double *y, long *incy, double *ssmin); /* Purpose ======= Given two column vectors X and Y, let A = ( X Y ). The subroutine first computes the QR factorization of A = Q*R, and then computes the SVD of the 2-by-2 upper triangular matrix R. The smaller singular value of R is returned in SSMIN, which is used as the measurement of the linear dependency of the vectors X and Y. Arguments ========= N (input) long The length of the vectors X and Y. X (input/output) double array, dimension (1+(N-1)*INCX) On entry, X contains the N-vector X. On exit, X is overwritten. INCX (input) long The increment between successive elements of X. INCX > 0. Y (input/output) double array, dimension (1+(N-1)*INCY) On entry, Y contains the N-vector Y. On exit, Y is overwritten. INCY (input) long The increment between successive elements of Y. INCY > 0. SSMIN (output) double The smallest singular value of the N-by-2 matrix A = ( X Y ). ===================================================================== */ double NUMlapack_dlapy2 (double *x, double *y); /* Purpose ======= NUMlapack_dlapy2 returns sqrt(x**2+y**2), taking care not to cause unnecessary overflow. Arguments ========= X (input) double Y (input) double X and Y specify the values x and y. ===================================================================== */ int NUMlapack_dlapmt (long *forwrd, long *m, long *n, double *x, long *ldx, long *k); /* Purpose ======= NUMlapack_dlapmt rearranges the columns of the M by N matrix X as specified by the permutation K(1),K(2),...,K(N) of the integers 1,...,N. If FORWRD = TRUE, forward permutation: X(*,K(J)) is moved X(*,J) for J = 1,2,...,N. If FORWRD = FALSE, backward permutation: X(*,J) is moved to X(*,K(J)) for J = 1,2,...,N. Arguments ========= FORWRD (input) long* (boolean) = TRUE, forward permutation = FALSE, backward permutation M (input) long The number of rows of the matrix X. M >= 0. N (input) long The number of columns of the matrix X. N >= 0. X (input/output) double array, dimension (LDX,N) On entry, the M by N matrix X. On exit, X contains the permuted matrix X. LDX (input) long The leading dimension of the array X, LDX >= MAX(1,M). K (input) long array, dimension (N) On entry, K contains the permutation vector. ===================================================================== */ int NUMlapack_dlarf (char *side, long *m, long *n, double *v, long *incv, double *tau, double *c, long *ldc, double *work); /* Purpose ======= NUMlapack_dlarf applies a real elementary reflector H to a real m by n matrix C, from either the left or the right. H is represented in the form H = I - tau * v * v' where tau is a real scalar and v is a real vector. If tau = 0, then H is taken to be the unit matrix. Arguments ========= SIDE (input) char* = 'L': form H * C = 'R': form C * H M (input) long The number of rows of the matrix C. N (input) long The number of columns of the matrix C. V (input) double array, dimension (1 + (M-1)*abs(INCV)) if SIDE = 'L' or (1 + (N-1)*abs(INCV)) if SIDE = 'R' The vector v in the representation of H. V is not used if TAU = 0. INCV (input) long The increment between elements of v. INCV <> 0. TAU (input) double The value tau in the representation of H. C (input/output) double array, dimension (LDC,N) On entry, the m by n matrix C. On exit, C is overwritten by the matrix H * C if SIDE = 'L', or C * H if SIDE = 'R'. LDC (input) long The leading dimension of the array C. LDC >= max(1,M). WORK (workspace) double array, dimension (N) if SIDE = 'L' or (M) if SIDE = 'R' ===================================================================== */ int NUMlapack_dlarfb (char *side, char *trans, char *direct, char *storev, long *m, long *n, long *k, double *v, long *ldv, double *t, long *ldt, double *c, long *ldc, double *work, long *ldwork); /* Purpose ======= NUMlapack_dlarfb applies a real block reflector H or its transpose H' to a real m by n matrix C, from either the left or the right. Arguments ========= SIDE (input) char* = 'L': apply H or H' from the Left = 'R': apply H or H' from the Right TRANS (input) char* = 'N': apply H (No transpose) = 'T': apply H' (Transpose) DIRECT (input) char* Indicates how H is formed from a product of elementary reflectors = 'F': H = H(1) H(2) . . . H(k) (Forward) = 'B': H = H(k) . . . H(2) H(1) (Backward) STOREV (input) char* Indicates how the vectors which define the elementary reflectors are stored: = 'C': Columnwise = 'R': Rowwise M (input) long The number of rows of the matrix C. N (input) long The number of columns of the matrix C. K (input) long The order of the matrix T (= the number of elementary reflectors whose product defines the block reflector). V (input) double array, dimension (LDV,K) if STOREV = 'C' (LDV,M) if STOREV = 'R' and SIDE = 'L' (LDV,N) if STOREV = 'R' and SIDE = 'R' The matrix V. See further details. LDV (input) long The leading dimension of the array V. If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); if STOREV = 'R', LDV >= K. T (input) double array, dimension (LDT,K) The triangular k by k matrix T in the representation of the block reflector. LDT (input) long The leading dimension of the array T. LDT >= K. C (input/output) double array, dimension (LDC,N) On entry, the m by n matrix C. On exit, C is overwritten by H*C or H'*C or C*H or C*H'. LDC (input) long The leading dimension of the array C. LDA >= max(1,M). WORK (workspace) double array, dimension (LDWORK,K) LDWORK (input) long The leading dimension of the array WORK. If SIDE = 'L', LDWORK >= max(1,N); if SIDE = 'R', LDWORK >= max(1,M). ===================================================================== */ int NUMlapack_dlarfg (long *n, double *alpha, double *x, long *incx, double *tau); /* Purpose ======= NUMlapack_dlarfg generates a real elementary reflector H of order n, such that H * ( alpha ) = ( beta ), H' * H = I. ( x ) ( 0 ) where alpha and beta are scalars, and x is an (n-1)-element real vector. H is represented in the form H = I - tau * ( 1 ) * ( 1 v' ) , ( v ) where tau is a real scalar and v is a real (n-1)-element vector. If the elements of x are all zero, then tau = 0 and H is taken to be the unit matrix. Otherwise 1 <= tau <= 2. Arguments ========= N (input) long The order of the elementary reflector. ALPHA (input/output) double On entry, the value alpha. On exit, it is overwritten with the value beta. X (input/output) double array, dimension (1+(N-2)*abs(INCX)) On entry, the vector x. On exit, it is overwritten with the vector v. INCX (input) long The increment between elements of X. INCX > 0. TAU (output) double The value tau. ===================================================================== */ int NUMlapack_dlarft (char *direct, char *storev, long *n, long *k, double *v, long *ldv, double *tau, double *t, long *ldt); /* Purpose ======= NUMlapack_dlarft forms the triangular factor T of a real block reflector H of order n, which is defined as a product of k elementary reflectors. If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular; If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular. If STOREV = 'C', the vector which defines the elementary reflector H(i) is stored in the i-th column of the array V, and H = I - V * T * V' If STOREV = 'R', the vector which defines the elementary reflector H(i) is stored in the i-th row of the array V, and H = I - V' * T * V Arguments ========= DIRECT (input) char* Specifies the order in which the elementary reflectors are multiplied to form the block reflector: = 'F': H = H(1) H(2) . . . H(k) (Forward) = 'B': H = H(k) . . . H(2) H(1) (Backward) STOREV (input) char* Specifies how the vectors which define the elementary reflectors are stored (see also Further Details): = 'C': columnwise = 'R': rowwise N (input) long The order of the block reflector H. N >= 0. K (input) long The order of the triangular factor T (= the number of elementary reflectors). K >= 1. V (input/output) double array, dimension (LDV,K) if STOREV = 'C' (LDV,N) if STOREV = 'R' The matrix V. See further details. LDV (input) long The leading dimension of the array V. If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. TAU (input) double array, dimension (K) TAU(i) must contain the scalar factor of the elementary reflector H(i). T (output) double array, dimension (LDT,K) The k by k triangular factor T of the block reflector. If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is lower triangular. The rest of the array is not used. LDT (input) long The leading dimension of the array T. LDT >= K. Further Details =============== The shape of the matrix V and the storage of the vectors which define the H(i) is best illustrated by the following example with n = 5 and k = 3. The elements equal to 1 are not stored; the corresponding array elements are modified but restored on exit. The rest of the array is not used. DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) ( v1 1 ) ( 1 v2 v2 v2 ) ( v1 v2 1 ) ( 1 v3 v3 ) ( v1 v2 v3 ) ( v1 v2 v3 ) DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': V = ( v1 v2 v3 ) V = ( v1 v1 1 ) ( v1 v2 v3 ) ( v2 v2 v2 1 ) ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) ( 1 v3 ) ( 1 ) ===================================================================== */ int NUMlapack_dlarfx (char *side, long *m, long *n, double *v, double *tau, double *c, long *ldc, double *work); /* -- LAPACK auxiliary routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University February 29, 1992 Purpose ======= NUMlapack_dlarfx applies a real elementary reflector H to a real m by n matrix C, from either the left or the right. H is represented in the form H = I - tau * v * v' where tau is a real scalar and v is a real vector. If tau = 0, then H is taken to be the unit matrix This version uses inline code if H has order < 11. Arguments ========= SIDE (input) char* = 'L': form H * C = 'R': form C * H M (input) long The number of rows of the matrix C. N (input) long The number of columns of the matrix C. V (input) double array, dimension (M) if SIDE = 'L' or (N) if SIDE = 'R' The vector v in the representation of H. TAU (input) double The value tau in the representation of H. C (input/output) double array, dimension (LDC,N) On entry, the m by n matrix C. On exit, C is overwritten by the matrix H * C if SIDE = 'L', or C * H if SIDE = 'R'. LDC (input) long The leading dimension of the array C. LDA >= (1,M). WORK (workspace) double array, dimension (N) if SIDE = 'L' or (M) if SIDE = 'R' WORK is not referenced if H has order < 11. ===================================================================== */ int NUMlapack_dlartg (double *f, double *g, double *cs, double *sn, double *r); /* Purpose ======= NUMlapack_dlartg generate a plane rotation so that [ CS SN ] . [ F ] = [ R ] where CS**2 + SN**2 = 1. [ -SN CS ] [ G ] [ 0 ] This is a slower, more accurate version of the BLAS1 routine DROTG, with the following other differences: F and G are unchanged on return. If G=0, then CS=1 and SN=0. If F=0 and (G .ne. 0), then CS=0 and SN=1 without doing any floating point operations (saves work in DBDSQR when there are zeros on the diagonal). If F exceeds G in magnitude, CS will be positive. Arguments ========= F (input) double The first component of vector to be rotated. G (input) double The second component of vector to be rotated. CS (output) double The cosine of the rotation. SN (output) double The sine of the rotation. R (output) double The nonzero component of the rotated vector. ===================================================================== */ int NUMlapack_dlas2 (double *f, double *g, double *h, double *ssmin, double *ssmax); /* Purpose ======= NUMlapack_dlas2 computes the singular values of the 2-by-2 matrix [ F G ] [ 0 H ]. On return, SSMIN is the smaller singular value and SSMAX is the larger singular value. Arguments ========= F (input) double The (1,1) element of the 2-by-2 matrix. G (input) double The (1,2) element of the 2-by-2 matrix. H (input) double The (2,2) element of the 2-by-2 matrix. SSMIN (output) double The smaller singular value. SSMAX (output) double The larger singular value. Further Details =============== Barring over/underflow, all output quantities are correct to within a few units in the last place (ulps), even in the absence of a guard digit in addition/subtraction. In IEEE arithmetic, the code works correctly if one matrix element is infinite. Overflow will not occur unless the largest singular value itself overflows, or is within a few ulps of overflow. (On machines with partial overflow, like the Cray, overflow may occur if the largest singular value is within a factor of 2 of overflow.) Underflow is harmless if underflow is gradual. Otherwise, results may correspond to a matrix modified by perturbations of size near the underflow threshold. ==================================================================== */ int NUMlapack_dlascl (char *type, long *kl, long *ku, double *cfrom, double *cto, long *m, long *n, double *a, long *lda, long *info); /* Purpose ======= NUMlapack_dlascl multiplies the M by N real matrix A by the real scalar CTO/CFROM. This is done without over/underflow as long as the final result CTO*A(I,J)/CFROM does not over/underflow. TYPE specifies that A may be full, upper triangular, lower triangular, upper Hessenberg, or banded. Arguments ========= TYPE (input) char* TYPE indices the storage type of the input matrix. = 'G': A is a full matrix. = 'L': A is a lower triangular matrix. = 'U': A is an upper triangular matrix. = 'H': A is an upper Hessenberg matrix. = 'B': A is a symmetric band matrix with lower bandwidth KL and upper bandwidth KU and with the only the lower half stored. = 'Q': A is a symmetric band matrix with lower bandwidth KL and upper bandwidth KU and with the only the upper half stored. = 'Z': A is a band matrix with lower bandwidth KL and upper bandwidth KU. KL (input) long The lower bandwidth of A. Referenced only if TYPE = 'B', 'Q' or 'Z'. KU (input) long The upper bandwidth of A. Referenced only if TYPE = 'B', 'Q' or 'Z'. CFROM (input) double CTO (input) double The matrix A is multiplied by CTO/CFROM. A(I,J) is computed without over/underflow if the final result CTO*A(I,J)/CFROM can be represented without over/underflow. CFROM must be nonzero. M (input) long The number of rows of the matrix A. M >= 0. N (input) long The number of columns of the matrix A. N >= 0. A (input/output) double array, dimension (LDA,M) The matrix to be multiplied by CTO/CFROM. See TYPE for the storage type. LDA (input) long The leading dimension of the array A. LDA >= max(1,M). INFO (output) long 0 - successful exit <0 - if INFO = -i, the i-th argument had an illegal value. ===================================================================== */ int NUMlapack_dlaset (char *uplo, long *m, long *n, double *alpha, double *beta, double *a, long *lda); /* Purpose ======= NUMlapack_dlaset initializes an m-by-n matrix A to BETA on the diagonal and ALPHA on the offdiagonals. Arguments ========= UPLO (input) char* Specifies the part of the matrix A to be set. = 'U': Upper triangular part is set; the strictly lower triangular part of A is not changed. = 'L': Lower triangular part is set; the strictly upper triangular part of A is not changed. Otherwise: All of the matrix A is set. M (input) long The number of rows of the matrix A. M >= 0. N (input) long The number of columns of the matrix A. N >= 0. ALPHA (input) double The constant to which the offdiagonal elements are to be set. BETA (input) double The constant to which the diagonal elements are to be set. A (input/output) double array, dimension (LDA,N) On exit, the leading m-by-n submatrix of A is set as follows: if UPLO = 'U', A(i,j) = ALPHA, 1<=i<=j-1, 1<=j<=n, if UPLO = 'L', A(i,j) = ALPHA, j+1<=i<=m, 1<=j<=n, otherwise, A(i,j) = ALPHA, 1<=i<=m, 1<=j<=n, i.ne.j, and, for all UPLO, A(i,i) = BETA, 1<=i<=min(m,n). LDA (input) long The leading dimension of the array A. LDA >= max(1,M). ===================================================================== */ int NUMlapack_dlasq1 (long *n, double *d, double *e, double *work, long *info); /* Purpose ======= NUMlapack_dlasq1 computes the singular values of a real N-by-N bidiagonal matrix with diagonal D and off-diagonal E. The singular values are computed to high relative accuracy, in the absence of denormalization, underflow and overflow. The algorithm was first presented in "Accurate singular values and differential qd algorithms" by K. V. Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230, 1994, and the present implementation is described in "An implementation of the dqds Algorithm (Positive Case)", LAPACK Working Note. Arguments ========= N (input) long The number of rows and columns in the matrix. N >= 0. D (input/output) double array, dimension (N) On entry, D contains the diagonal elements of the bidiagonal matrix whose SVD is desired. On normal exit, D contains the singular values in decreasing order. E (input/output) double array, dimension (N) On entry, elements E(1:N-1) contain the off-diagonal elements of the bidiagonal matrix whose SVD is desired. On exit, E is overwritten. WORK (workspace) double array, dimension (4*N) INFO (output) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: the algorithm failed = 1, a split was marked by a positive value in E = 2, current block of Z not diagonalized after 30*N iterations (in inner while loop) = 3, termination criterion of outer while loop not met (program created more than N unreduced blocks) ===================================================================== */ int NUMlapack_dlasq2 (long *n, double *z, long *info); /* Purpose ======= NUMlapack_dlasq2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd array Z to high relative accuracy are computed to high relative accuracy, in the absence of denormalization, underflow and overflow. To see the relation of Z to the tridiagonal matrix, let L be a unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and let U be an upper bidiagonal matrix with 1's above and diagonal Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the symmetric tridiagonal to which it is similar. Note : NUMlapack_dlasq2 defines a long* (boolean) variable, IEEE, which is true on machines which follow ieee-754 floating-point standard in their handling of infinities and NaNs, and false otherwise. This variable is passed to NUMlapack_dlasq3. Arguments ========= N (input) long The number of rows and columns in the matrix. N >= 0. Z (workspace) double array, dimension ( 4*N ) On entry Z holds the qd array. On exit, entries 1 to N hold the eigenvalues in decreasing order, Z( 2*N+1 ) holds the trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 ) holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of shifts that failed. INFO (output) long = 0: successful exit < 0: if the i-th argument is a scalar and had an illegal value, then INFO = -i, if the i-th argument is an array and the j-entry had an illegal value, then INFO = -(i*100+j) > 0: the algorithm failed = 1, a split was marked by a positive value in E = 2, current block of Z not diagonalized after 30*N iterations (in inner while loop) = 3, termination criterion of outer while loop not met (program created more than N unreduced blocks) Further Details =============== Local Variables: I0:N0 defines a current unreduced segment of Z. The shifts are accumulated in SIGMA. Iteration count is in ITER. Ping-pong is controlled by PP (alternates between 0 and 1). ===================================================================== */ int NUMlapack_dlasq3 (long *i0, long *n0, double *z, long *pp, double *dmin, double *sigma, double *desig, double *qmax, long *nfail, long *iter, long *ndiv, long *ieee); /* Purpose ======= NUMlapack_dlasq3 checks for deflation, computes a shift (TAU) and calls dqds. In case of failure it changes shifts, and tries again until output is positive. Arguments ========= I0 (input) long First index. N0 (input) long Last index. Z (input) double array, dimension ( 4*N ) Z holds the qd array. PP (input) long PP=0 for ping, PP=1 for pong. DMIN (output) double Minimum value of d. SIGMA (output) double Sum of shifts used in current segment. DESIG (input/output) double Lower order part of SIGMA QMAX (input) double Maximum value of q. NFAIL (output) long Number of times shift was too big. ITER (output) long Number of iterations. NDIV (output) long Number of divisions. TTYPE (output) long Shift type. IEEE (input) long* (boolean) Flag for IEEE or non IEEE arithmetic (passed to DLASQ5). ===================================================================== */ int NUMlapack_dlasq4 (long *i0, long *n0, double *z, long *pp, long *n0in, double *dmin, double *dmin1, double *dmin2, double *dn, double *dn1, double *dn2, double *tau, long *ttype); /* Purpose ======= NUMlapack_dlasq4 computes an approximation TAU to the smallest eigenvalue using values of d from the previous transform. I0 (input) long First index. N0 (input) long Last index. Z (input) double array, dimension ( 4*N ) Z holds the qd array. PP (input) long PP=0 for ping, PP=1 for pong. NOIN (input) long The value of N0 at start of EIGTEST. DMIN (input) double Minimum value of d. DMIN1 (input) double Minimum value of d, excluding D( N0 ). DMIN2 (input) double Minimum value of d, excluding D( N0 ) and D( N0-1 ). DN (input) double d(N) DN1 (input) double d(N-1) DN2 (input) double d(N-2) TAU (output) double This is the shift. TTYPE (output) long Shift type. Further Details =============== CNST1 = 9/16 ===================================================================== */ int NUMlapack_dlasq5 (long *i0, long *n0, double *z, long *pp, double *tau, double *dmin, double *dmin1, double *dmin2, double *dn, double *dnm1, double *dnm2, long *ieee); /* Purpose ======= NUMlapack_dlasq5 computes one dqds transform in ping-pong form, one version for IEEE machines another for non IEEE machines. Arguments ========= I0 (input) long First index. N0 (input) long Last index. Z (input) double array, dimension ( 4*N ) Z holds the qd array. EMIN is stored in Z(4*N0) to avoid an extra argument. PP (input) long PP=0 for ping, PP=1 for pong. TAU (input) double This is the shift. DMIN (output) double Minimum value of d. DMIN1 (output) double Minimum value of d, excluding D( N0 ). DMIN2 (output) double Minimum value of d, excluding D( N0 ) and D( N0-1 ). DN (output) double d(N0), the last value of d. DNM1 (output) double d(N0-1). DNM2 (output) double d(N0-2). IEEE (input) long* (boolean) Flag for IEEE or non IEEE arithmetic. ===================================================================== */ int NUMlapack_dlasq6 (long *i0, long *n0, double *z, long *pp, double *dmin, double *dmin1, double *dmin2, double *dn, double *dnm1, double *dnm2); /* Purpose ======= NUMlapack_dlasq6 computes one dqd (shift equal to zero) transform in ping-pong form, with protection against underflow and overflow. Arguments ========= I0 (input) long First index. N0 (input) long Last index. Z (input) double array, dimension ( 4*N ) Z holds the qd array. EMIN is stored in Z(4*N0) to avoid an extra argument. PP (input) long PP=0 for ping, PP=1 for pong. DMIN (output) double Minimum value of d. DMIN1 (output) double Minimum value of d, excluding D( N0 ). DMIN2 (output) double Minimum value of d, excluding D( N0 ) and D( N0-1 ). DN (output) double d(N0), the last value of d. DNM1 (output) double d(N0-1). DNM2 (output) double d(N0-2). ===================================================================== */ int NUMlapack_dlasr (char *side, char *pivot, char *direct, long *m, long *n, double *c, double *s, double *a, long *lda); /* Purpose ======= NUMlapack_dlasr performs the transformation A := P*A, when SIDE = 'L' or 'l' ( Left-hand side ) A := A*P', when SIDE = 'R' or 'r' ( Right-hand side ) where A is an m by n real matrix and P is an orthogonal matrix, consisting of a sequence of plane rotations determined by the parameters PIVOT and DIRECT as follows ( z = m when SIDE = 'L' or 'l' and z = n when SIDE = 'R' or 'r' ): When DIRECT = 'F' or 'f' ( Forward sequence ) then P = P( z - 1 )*...*P( 2 )*P( 1 ), and when DIRECT = 'B' or 'b' ( Backward sequence ) then P = P( 1 )*P( 2 )*...*P( z - 1 ), where P( k ) is a plane rotation matrix for the following planes: when PIVOT = 'V' or 'v' ( Variable pivot ), the plane ( k, k + 1 ) when PIVOT = 'T' or 't' ( Top pivot ), the plane ( 1, k + 1 ) when PIVOT = 'B' or 'b' ( Bottom pivot ), the plane ( k, z ) c( k ) and s( k ) must contain the cosine and sine that define the matrix P( k ). The two by two plane rotation part of the matrix P( k ), R( k ), is assumed to be of the form R( k ) = ( c( k ) s( k ) ). ( -s( k ) c( k ) ) This version vectorises across rows of the array A when SIDE = 'L'. Arguments ========= SIDE (input) char* Specifies whether the plane rotation matrix P is applied to A on the left or the right. = 'L': Left, compute A := P*A = 'R': Right, compute A:= A*P' DIRECT (input) char* Specifies whether P is a forward or backward sequence of plane rotations. = 'F': Forward, P = P( z - 1 )*...*P( 2 )*P( 1 ) = 'B': Backward, P = P( 1 )*P( 2 )*...*P( z - 1 ) PIVOT (input) char* Specifies the plane for which P(k) is a plane rotation matrix. = 'V': Variable pivot, the plane (k,k+1) = 'T': Top pivot, the plane (1,k+1) = 'B': Bottom pivot, the plane (k,z) M (input) long The number of rows of the matrix A. If m <= 1, an immediate return is effected. N (input) long The number of columns of the matrix A. If n <= 1, an immediate return is effected. C, S (input) double arrays, dimension (M-1) if SIDE = 'L' (N-1) if SIDE = 'R' c(k) and s(k) contain the cosine and sine that define the matrix P(k). The two by two plane rotation part of the matrix P(k), R(k), is assumed to be of the form R( k ) = ( c( k ) s( k ) ). ( -s( k ) c( k ) ) A (input/output) double array, dimension (LDA,N) The m by n matrix A. On exit, A is overwritten by P*A if SIDE = 'R' or by A*P' if SIDE = 'L'. LDA (input) long The leading dimension of the array A. LDA >= max(1,M). ===================================================================== */ int NUMlapack_dlasrt (char *id, long *n, double *d, long *info); /* Purpose ======= Sort the numbers in D in increasing order (if ID = 'I') or in decreasing order (if ID = 'D' ). Use Quick Sort, reverting to Insertion sort on arrays of size <= 20. Dimension of STACK limits N to about 2**32. Arguments ========= ID (input) char* = 'I': sort D in increasing order; = 'D': sort D in decreasing order. N (input) long The length of the array D. D (input/output) double array, dimension (N) On entry, the array to be sorted. On exit, D has been sorted into increasing order (D(1) <= ... <= D(N) ) or into decreasing order (D(1) >= ... >= D(N) ), depending on ID. INFO (output) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value ===================================================================== */ int NUMlapack_dlassq (long *n, double *x, long *incx, double *scale, double *sumsq); /* Purpose ======= NUMlapack_dlassq returns the values scl and smsq such that ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is assumed to be non-negative and scl returns the value scl = max( scale, abs( x( i ) ) ). scale and sumsq must be supplied in SCALE and SUMSQ and scl and smsq are overwritten on SCALE and SUMSQ respectively. The routine makes only one pass through the vector x. Arguments ========= N (input) long The number of elements to be used from the vector X. X (input) double array, dimension (N) The vector for which a scaled sum of squares is computed. x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n. INCX (input) long The increment between successive values of the vector X. INCX > 0. SCALE (input/output) double On entry, the value scale in the equation above. On exit, SCALE is overwritten with scl , the scaling factor for the sum of squares. SUMSQ (input/output) double On entry, the value sumsq in the equation above. On exit, SUMSQ is overwritten with smsq , the basic sum of squares from which scl has been factored out. ===================================================================== */ int NUMlapack_dlasv2 (double *f, double *g, double *h, double *ssmin, double *ssmax, double *snr, double *csr, double *snl, double *csl); /* Purpose ======= NUMlapack_dlasv2 computes the singular value decomposition of a 2-by-2 triangular matrix [ F G ] [ 0 H ]. On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and right singular vectors for abs(SSMAX), giving the decomposition [ CSL SNL ] [ F G ] [ CSR -SNR ] = [ SSMAX 0 ] [-SNL CSL ] [ 0 H ] [ SNR CSR ] [ 0 SSMIN ]. Arguments ========= F (input) double The (1,1) element of the 2-by-2 matrix. G (input) double The (1,2) element of the 2-by-2 matrix. H (input) double The (2,2) element of the 2-by-2 matrix. SSMIN (output) double abs(SSMIN) is the smaller singular value. SSMAX (output) double abs(SSMAX) is the larger singular value. SNL (output) double CSL (output) double The vector (CSL, SNL) is a unit left singular vector for the singular value abs(SSMAX). SNR (output) double CSR (output) double The vector (CSR, SNR) is a unit right singular vector for the singular value abs(SSMAX). Further Details =============== Any input parameter may be aliased with any output parameter. Barring over/underflow and assuming a guard digit in subtraction, all output quantities are correct to within a few units in the last place (ulps). In IEEE arithmetic, the code works correctly if one matrix element is infinite. Overflow will not occur unless the largest singular value itself overflows or is within a few ulps of overflow. (On machines with partial overflow, like the Cray, overflow may occur if the largest singular value is within a factor of 2 of overflow.) Underflow is harmless if underflow is gradual. Otherwise, results may correspond to a matrix modified by perturbations of size near the underflow threshold. ===================================================================== */ int NUMlapack_dlaswp (long *n, double *a, long *lda, long *k1, long *k2, long *ipiv, long *incx); /* Purpose ======= NUMlapack_dlaswp performs a series of row interchanges on the matrix A. One row interchange is initiated for each of rows K1 through K2 of A. Arguments ========= N (input) long The number of columns of the matrix A. A (input/output) double array, dimension (LDA,N) On entry, the matrix of column dimension N to which the row interchanges will be applied. On exit, the permuted matrix. LDA (input) long The leading dimension of the array A. K1 (input) long The first element of IPIV for which a row interchange will be done. K2 (input) long The last element of IPIV for which a row interchange will be done. IPIV (input) long array, dimension (M*abs(INCX)) The vector of pivot indices. Only the elements in positions K1 through K2 of IPIV are accessed. IPIV(K) = L implies rows K and L are to be interchanged. INCX (input) long The increment between successive values of IPIV. If IPIV is negative, the pivots are applied in reverse order. ===================================================================== */ int NUMlapack_dlatrd (char *uplo, long *n, long *nb, double *a, long *lda, double *e, double *tau, double *w, long *ldw); /* Purpose ======= NUMlapack_dlatrd reduces NB rows and columns of a real symmetric matrix A to symmetric tridiagonal form by an orthogonal similarity transformation Q' * A * Q, and returns the matrices V and W which are needed to apply the transformation to the unreduced part of A. If UPLO = 'U', DLATRD reduces the last NB rows and columns of a matrix, of which the upper triangle is supplied; if UPLO = 'L', DLATRD reduces the first NB rows and columns of a matrix, of which the lower triangle is supplied. This is an auxiliary routine called by NUMlapack_dsytrd. Arguments ========= UPLO (input) CHARACTER Specifies whether the upper or lower triangular part of the symmetric matrix A is stored: = 'U': Upper triangular = 'L': Lower triangular N (input) long The order of the matrix A. NB (input) long The number of rows and columns to be reduced. A (input/output) double array, dimension (LDA,N) On entry, the symmetric matrix A. If UPLO = 'U', the leading n-by-n upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading n-by-n lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. On exit: if UPLO = 'U', the last NB columns have been reduced to tridiagonal form, with the diagonal elements overwriting the diagonal elements of A; the elements above the diagonal with the array TAU, represent the orthogonal matrix Q as a product of elementary reflectors; if UPLO = 'L', the first NB columns have been reduced to tridiagonal form, with the diagonal elements overwriting the diagonal elements of A; 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) long The leading dimension of the array A. LDA >= (1,N). E (output) double array, dimension (N-1) If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal elements of the last NB columns of the reduced matrix; if UPLO = 'L', E(1:nb) contains the subdiagonal elements of the first NB columns of the reduced matrix. TAU (output) double array, dimension (N-1) The scalar factors of the elementary reflectors, stored in TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'. See Further Details. W (output) double array, dimension (LDW,NB) The n-by-nb matrix W required to update the unreduced part of A. LDW (input) long The leading dimension of the array W. LDW >= max(1,N). Further Details =============== If UPLO = 'U', the matrix Q is represented as a product of elementary reflectors Q = H(n) H(n-1) . . . H(n-nb+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(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i), and tau in TAU(i-1). If UPLO = 'L', the matrix Q is represented as a product of elementary reflectors Q = H(1) H(2) . . . H(nb). 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 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i), and tau in TAU(i). The elements of the vectors v together form the n-by-nb matrix V which is needed, with W, to apply the transformation to the unreduced part of the matrix, using a symmetric rank-2k update of the form: A := A - V*W' - W*V'. The contents of A on exit are illustrated by the following examples with n = 5 and nb = 2: if UPLO = 'U': if UPLO = 'L': ( a a a v4 v5 ) ( d ) ( a a v4 v5 ) ( 1 d ) ( a 1 v5 ) ( v1 1 a ) ( d 1 ) ( v1 v2 a a ) ( d ) ( v1 v2 a a a ) where d denotes a diagonal element of the reduced matrix, a denotes an element of the original matrix that is unchanged, and vi denotes an element of the vector defining H(i). ===================================================================== */ int NUMlapack_dorg2l (long *m, long *n, long *k, double * a, long *lda, double *tau, double *work, long *info); /* Purpose ======= NUMlapack_dorg2l generates an m by n real matrix Q with orthonormal columns, which is defined as the last n columns of a product of k elementary reflectors of order m Q = H(k) . . . H(2) H(1) as returned by NUMlapack_dgeqlf. Arguments ========= M (input) long The number of rows of the matrix Q. M >= 0. N (input) long The number of columns of the matrix Q. M >= N >= 0. K (input) long The number of elementary reflectors whose product defines the matrix Q. N >= K >= 0. A (input/output) double array, dimension (LDA,N) On entry, the (n-k+i)-th column must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by DGEQLF in the last k columns of its array argument A. On exit, the m by n matrix Q. LDA (input) long The first dimension of the array A. LDA >= max(1,M). TAU (input) double array, dimension (K) TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by DGEQLF. WORK (workspace) double array, dimension (N) INFO (output) long = 0: successful exit < 0: if INFO = -i, the i-th argument has an illegal value ===================================================================== */ int NUMlapack_dorg2r (long *m, long *n, long *k, double *a, long *lda, double *tau, double *work, long *info); /* Purpose ======= NUMlapack_dorg2r generates an m by n real matrix Q with orthonormal columns, which is defined as the first n columns of a product of k elementary reflectors of order m Q = H(1) H(2) . . . H(k) as returned by NUMlapack_dgeqrf. Arguments ========= M (input) long The number of rows of the matrix Q. M >= 0. N (input) long The number of columns of the matrix Q. M >= N >= 0. K (input) long The number of elementary reflectors whose product defines the matrix Q. N >= K >= 0. A (input/output) double array, dimension (LDA,N) On entry, the i-th column must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by DGEQRF in the first k columns of its array argument A. On exit, the m-by-n matrix Q. LDA (input) long The first dimension of the array A. LDA >= max(1,M). TAU (input) double array, dimension (K) TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by DGEQRF. WORK (workspace) double array, dimension (N) INFO (output) long = 0: successful exit < 0: if INFO = -i, the i-th argument has an illegal value ===================================================================== */ int NUMlapack_dorgbr (char *vect, long *m, long *n, long *k, double *a, long *lda, double *tau, double *work, long *lwork, long *info); /* Purpose ======= NUMlapack_dorgbr generates one of the real orthogonal matrices Q or P**T determined by DGEBRD when reducing a real matrix A to bidiagonal form: A = Q * B * P**T. Q and P**T are defined as products of elementary reflectors H(i) or G(i) respectively. If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q is of order M: if m >= k, Q = H(1) H(2) . . . H(k) and NUMlapack_dorgbr returns the first n columns of Q, where m >= n >= k; if m < k, Q = H(1) H(2) . . . H(m-1) and NUMlapack_dorgbr returns Q as an M-by-M matrix. If VECT = 'P', A is assumed to have been a K-by-N matrix, and P**T is of order N: if k < n, P**T = G(k) . . . G(2) G(1) and NUMlapack_dorgbr returns the first m rows of P**T, where n >= m >= k; if k >= n, P**T = G(n-1) . . . G(2) G(1) and NUMlapack_dorgbr returns P**T as an N-by-N matrix. Arguments ========= VECT (input) char* Specifies whether the matrix Q or the matrix P**T is required, as defined in the transformation applied by DGEBRD: = 'Q': generate Q; = 'P': generate P**T. M (input) long The number of rows of the matrix Q or P**T to be returned. M >= 0. N (input) long The number of columns of the matrix Q or P**T to be returned. N >= 0. If VECT = 'Q', M >= N >= min(M,K); if VECT = 'P', N >= M >= min(N,K). K (input) long If VECT = 'Q', the number of columns in the original M-by-K matrix reduced by DGEBRD. If VECT = 'P', the number of rows in the original K-by-N matrix reduced by DGEBRD. K >= 0. A (input/output) double array, dimension (LDA,N) On entry, the vectors which define the elementary reflectors, as returned by DGEBRD. On exit, the M-by-N matrix Q or P**T. LDA (input) long The leading dimension of the array A. LDA >= max(1,M). TAU (input) double array, dimension (min(M,K)) if VECT = 'Q' (min(N,K)) if VECT = 'P' TAU(i) must contain the scalar factor of the elementary reflector H(i) or G(i), which determines Q or P**T, as returned by DGEBRD in its array argument TAUQ or TAUP. WORK (workspace/output) double array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) long The dimension of the array WORK. LWORK >= max(1,min(M,N)). For optimum performance LWORK >= min(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) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value ===================================================================== */ int NUMlapack_dorghr (long *n, long *ilo, long *ihi, double *a, long *lda, double *tau, double *work, long *lwork, long *info); /* -- LAPACK routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University June 30, 1999 Purpose ======= NUMlapack_dorghr generates a real orthogonal matrix Q which is defined as the product of IHI-ILO elementary reflectors of order N, as returned by NUMlapack_dgehrd: Q = H(ilo) H(ilo+1) . . . H(ihi-1). Arguments ========= N (input) long The order of the matrix Q. N >= 0. ILO (input) long IHI (input) long ILO and IHI must have the same values as in the previous call of NUMlapack_dgehrd. Q is equal to the unit matrix except in the submatrix Q(ilo+1:ihi,ilo+1:ihi). 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. A (input/output) double array, dimension (LDA,N) On entry, the vectors which define the elementary reflectors, as returned by NUMlapack_dgehrd. On exit, the N-by-N orthogonal matrix Q. LDA (input) long The leading dimension of the array A. LDA >= max(1,N). TAU (input) double array, dimension (N-1) TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by NUMlapack_dgehrd. WORK (workspace/output) double array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) long The dimension of the array WORK. LWORK >= IHI-ILO. For optimum performance LWORK >= (IHI-ILO)*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) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value ===================================================================== */ int NUMlapack_dorgl2 (long *m, long *n, long *k, double *a, long *lda, double *tau, double *work, long *info); /* Purpose ======= NUMlapack_dorgl2 generates an m by n real matrix Q with orthonormal rows, which is defined as the first m rows of a product of k elementary reflectors of order n Q = H(k) . . . H(2) H(1) as returned by NUMlapack_dgelqf. Arguments ========= M (input) long The number of rows of the matrix Q. M >= 0. N (input) long The number of columns of the matrix Q. N >= M. K (input) long The number of elementary reflectors whose product defines the matrix Q. M >= K >= 0. A (input/output) double array, dimension (LDA,N) On entry, the i-th row must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by NUMlapack_dgelqf in the first k rows of its array argument A. On exit, the m-by-n matrix Q. LDA (input) long The first dimension of the array A. LDA >= max(1,M). TAU (input) double array, dimension (K) TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by NUMlapack_dgelqf. WORK (workspace) double array, dimension (M) INFO (output) long = 0: successful exit < 0: if INFO = -i, the i-th argument has an illegal value ===================================================================== */ int NUMlapack_dorglq (long *m, long *n, long *k, double *a, long *lda, double *tau, double *work, long *lwork, long *info); /* Purpose ======= NUMlapack_dorglq generates an M-by-N real matrix Q with orthonormal rows, which is defined as the first M rows of a product of K elementary reflectors of order N Q = H(k) . . . H(2) H(1) as returned by DGELQf. Arguments ========= M (input) long The number of rows of the matrix Q. M >= 0. N (input) long The number of columns of the matrix Q. N >= M. K (input) long The number of elementary reflectors whose product defines the matrix Q. M >= K >= 0. A (input/output) double array, dimension (LDA,N) On entry, the i-th row must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by NUMlapack_dgelqf in the first k rows of its array argument A. On exit, the M-by-N matrix Q. LDA (input) long The first dimension of the array A. LDA >= max(1,M). TAU (input) double array, dimension (K) TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by NUMlapack_dgelqf. WORK (workspace/output) double array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) long 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) long = 0: successful exit < 0: if INFO = -i, the i-th argument has an illegal value ===================================================================== */ int NUMlapack_dorgql (long *m, long *n, long *k, double *a, long *lda, double *tau, double *work, long *lwork, long *info); /* Purpose ======= NUMlapack_dorgql generates an M-by-N real matrix Q with orthonormal columns, which is defined as the last N columns of a product of K elementary reflectors of order M Q = H(k) . . . H(2) H(1) as returned by DGEQLF. Arguments ========= M (input) long The number of rows of the matrix Q. M >= 0. N (input) long The number of columns of the matrix Q. M >= N >= 0. K (input) long The number of elementary reflectors whose product defines the matrix Q. N >= K >= 0. A (input/output) double array, dimension (LDA,N) On entry, the (n-k+i)-th column must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by DGEQLF in the last k columns of its array argument A. On exit, the M-by-N matrix Q. LDA (input) long The first dimension of the array A. LDA >= max(1,M). TAU (input) double array, dimension (K) TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by DGEQLF. WORK (workspace/output) double array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) long 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) long = 0: successful exit < 0: if INFO = -i, the i-th argument has an illegal value ===================================================================== */ int NUMlapack_dorgqr (long *m, long *n, long *k, double *a, long *lda, double *tau, double *work, long *lwork, long *info); /* Purpose ======= NUMlapack_dorgqr generates an M-by-N real matrix Q with orthonormal columns, which is defined as the first N columns of a product of K elementary reflectors of order M Q = H(1) H(2) . . . H(k) as returned by DGEQRF. Arguments ========= M (input) long The number of rows of the matrix Q. M >= 0. N (input) long The number of columns of the matrix Q. M >= N >= 0. K (input) long The number of elementary reflectors whose product defines the matrix Q. N >= K >= 0. A (input/output) double array, dimension (LDA,N) On entry, the i-th column must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by DGEQRF in the first k columns of its array argument A. On exit, the M-by-N matrix Q. LDA (input) long The first dimension of the array A. LDA >= max(1,M). TAU (input) double array, dimension (K) TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by DGEQRF. WORK (workspace/output) double array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) long 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) long = 0: successful exit < 0: if INFO = -i, the i-th argument has an illegal value ===================================================================== */ int NUMlapack_dorgtr (char *uplo, long *n, double *a, long *lda, double *tau, double *work, long *lwork, long *info); /* Purpose ======= NUMlapack_dorgtr generates a real orthogonal matrix Q which is defined as the product of n-1 elementary reflectors of order N, as returned by NUMlapack_dsytrd: if UPLO = 'U', Q = H(n-1) . . . H(2) H(1), if UPLO = 'L', Q = H(1) H(2) . . . H(n-1). Arguments ========= UPLO (input) char* = 'U': Upper triangle of A contains elementary reflectors from NUMlapack_dsytrd; = 'L': Lower triangle of A contains elementary reflectors from NUMlapack_dsytrd. N (input) long The order of the matrix Q. N >= 0. A (input/output) double array, dimension (LDA,N) On entry, the vectors which define the elementary reflectors, as returned by NUMlapack_dsytrd. On exit, the N-by-N orthogonal matrix Q. LDA (input) long The leading dimension of the array A. LDA >= max(1,N). TAU (input) double array, dimension (N-1) TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by NUMlapack_dsytrd. WORK (workspace/output) double array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) long The dimension of the array WORK. LWORK >= max(1,N-1). For optimum performance LWORK >= (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) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value ===================================================================== */ int NUMlapack_dorm2r (char *side, char *trans, long *m, long *n, long *k, double *a, long *lda, double *tau, double *c, long *ldc, double *work, long *info); /* Purpose ======= NUMlapack_dorm2r overwrites the general real m by n matrix C with Q * C if SIDE = 'L' and TRANS = 'N', or Q'* C if SIDE = 'L' and TRANS = 'T', or C * Q if SIDE = 'R' and TRANS = 'N', or C * Q' if SIDE = 'R' and TRANS = 'T', where Q is a real orthogonal matrix defined as the product of k elementary reflectors Q = H(1) H(2) . . . H(k) as returned by DGEQRF. Q is of order m if SIDE = 'L' and of order n if SIDE = 'R'. Arguments ========= SIDE (input) char* = 'L': apply Q or Q' from the Left = 'R': apply Q or Q' from the Right TRANS (input) char* = 'N': apply Q (No transpose) = 'T': apply Q' (Transpose) M (input) long The number of rows of the matrix C. M >= 0. N (input) long The number of columns of the matrix C. N >= 0. K (input) long The number of elementary reflectors whose product defines the matrix Q. If SIDE = 'L', M >= K >= 0; if SIDE = 'R', N >= K >= 0. A (input) double array, dimension (LDA,K) The i-th column must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by DGEQRF in the first k columns of its array argument A. A is modified by the routine but restored on exit. LDA (input) long The leading dimension of the array A. If SIDE = 'L', LDA >= max(1,M); if SIDE = 'R', LDA >= max(1,N). TAU (input) double array, dimension (K) TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by DGEQRF. C (input/output) double array, dimension (LDC,N) On entry, the m by n matrix C. On exit, C is overwritten by Q*C or Q'*C or C*Q' or C*Q. LDC (input) long The leading dimension of the array C. LDC >= max(1,M). WORK (workspace) double array, dimension (N) if SIDE = 'L', (M) if SIDE = 'R' INFO (output) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value ===================================================================== */ int NUMlapack_dormbr (char *vect, char *side, char *trans, long *m, long *n, long *k, double *a, long *lda, double *tau, double *c, long *ldc, double *work, long *lwork, long *info); /* Purpose ======= If VECT = 'Q', NUMlapack_dormbr overwrites the general real M-by-N matrix C with SIDE = 'L' SIDE = 'R' TRANS = 'N': Q * C C * Q TRANS = 'T': Q**T * C C * Q**T If VECT = 'P', NUMlapack_dormbr overwrites the general real M-by-N matrix C with SIDE = 'L' SIDE = 'R' TRANS = 'N': P * C C * P TRANS = 'T': P**T * C C * P**T Here Q and P**T are the orthogonal matrices determined by DGEBRD when reducing a real matrix A to bidiagonal form: A = Q * B * P**T. Q and P**T are defined as products of elementary reflectors H(i) and G(i) respectively. Let nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Thus nq is the order of the orthogonal matrix Q or P**T that is applied. If VECT = 'Q', A is assumed to have been an NQ-by-K matrix: if nq >= k, Q = H(1) H(2) . . . H(k); if nq < k, Q = H(1) H(2) . . . H(nq-1). If VECT = 'P', A is assumed to have been a K-by-NQ matrix: if k < nq, P = G(1) G(2) . . . G(k); if k >= nq, P = G(1) G(2) . . . G(nq-1). Arguments ========= VECT (input) char* = 'Q': apply Q or Q**T; = 'P': apply P or P**T. SIDE (input) char* = 'L': apply Q, Q**T, P or P**T from the Left; = 'R': apply Q, Q**T, P or P**T from the Right. TRANS (input) char* = 'N': No transpose, apply Q or P; = 'T': Transpose, apply Q**T or P**T. M (input) long The number of rows of the matrix C. M >= 0. N (input) long The number of columns of the matrix C. N >= 0. K (input) long If VECT = 'Q', the number of columns in the original matrix reduced by DGEBRD. If VECT = 'P', the number of rows in the original matrix reduced by DGEBRD. K >= 0. A (input) double array, dimension (LDA,min(nq,K)) if VECT = 'Q' (LDA,nq) if VECT = 'P' The vectors which define the elementary reflectors H(i) and G(i), whose products determine the matrices Q and P, as returned by DGEBRD. LDA (input) long The leading dimension of the array A. If VECT = 'Q', LDA >= max(1,nq); if VECT = 'P', LDA >= max(1,min(nq,K)). TAU (input) double array, dimension (min(nq,K)) TAU(i) must contain the scalar factor of the elementary reflector H(i) or G(i) which determines Q or P, as returned by DGEBRD in the array argument TAUQ or TAUP. C (input/output) double array, dimension (LDC,N) On entry, the M-by-N matrix C. On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q or P*C or P**T*C or C*P or C*P**T. LDC (input) long The leading dimension of the array C. LDC >= max(1,M). WORK (workspace/output) double array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) long The dimension of the array WORK. If SIDE = 'L', LWORK >= max(1,N); if SIDE = 'R', LWORK >= max(1,M). For optimum performance LWORK >= N*NB if SIDE = 'L', and LWORK >= M*NB if SIDE = 'R', 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) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value ===================================================================== */ int NUMlapack_dorml2 (char *side, char *trans, long *m, long *n, long *k, double *a, long *lda, double *tau, double *c, long *ldc, double *work, long *info); /* Purpose ======= NUMlapack_dorml2 overwrites the general real m by n matrix C with Q * C if SIDE = 'L' and TRANS = 'N', or Q'* C if SIDE = 'L' and TRANS = 'T', or C * Q if SIDE = 'R' and TRANS = 'N', or C * Q' if SIDE = 'R' and TRANS = 'T', where Q is a real orthogonal matrix defined as the product of k elementary reflectors Q = H(k) . . . H(2) H(1) as returned by NUMlapack_dgelqf. Q is of order m if SIDE = 'L' and of order n if SIDE = 'R'. Arguments ========= SIDE (input) char* = 'L': apply Q or Q' from the Left = 'R': apply Q or Q' from the Right TRANS (input) char* = 'N': apply Q (No transpose) = 'T': apply Q' (Transpose) M (input) long The number of rows of the matrix C. M >= 0. N (input) long The number of columns of the matrix C. N >= 0. K (input) long The number of elementary reflectors whose product defines the matrix Q. If SIDE = 'L', M >= K >= 0; if SIDE = 'R', N >= K >= 0. A (input) double array, dimension (LDA,M) if SIDE = 'L', (LDA,N) if SIDE = 'R' The i-th row must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by NUMlapack_dgelqf in the first k rows of its array argument A. A is modified by the routine but restored on exit. LDA (input) long The leading dimension of the array A. LDA >= max(1,K). TAU (input) double array, dimension (K) TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by NUMlapack_dgelqf. C (input/output) double array, dimension (LDC,N) On entry, the m by n matrix C. On exit, C is overwritten by Q*C or Q'*C or C*Q' or C*Q. LDC (input) long The leading dimension of the array C. LDC >= max(1,M). WORK (workspace) double array, dimension (N) if SIDE = 'L', (M) if SIDE = 'R' INFO (output) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value ===================================================================== */ int NUMlapack_dormlq (char *side, char *trans, long *m, long *n, long *k, double *a, long *lda, double *tau, double *c, long *ldc, double *work, long *lwork, long *info); /* Purpose ======= NUMlapack_dormlq overwrites the general real M-by-N matrix C with SIDE = 'L' SIDE = 'R' TRANS = 'N': Q * C C * Q TRANS = 'T': Q**T * C C * Q**T where Q is a real orthogonal matrix defined as the product of k elementary reflectors Q = H(k) . . . H(2) H(1) as returned by NUMlapack_dgelqf. Q is of order M if SIDE = 'L' and of order N if SIDE = 'R'. Arguments ========= SIDE (input) char* = 'L': apply Q or Q**T from the Left; = 'R': apply Q or Q**T from the Right. TRANS (input) char* = 'N': No transpose, apply Q; = 'T': Transpose, apply Q**T. M (input) long The number of rows of the matrix C. M >= 0. N (input) long The number of columns of the matrix C. N >= 0. K (input) long The number of elementary reflectors whose product defines the matrix Q. If SIDE = 'L', M >= K >= 0; if SIDE = 'R', N >= K >= 0. A (input) double array, dimension (LDA,M) if SIDE = 'L', (LDA,N) if SIDE = 'R' The i-th row must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by NUMlapack_dgelqf in the first k rows of its array argument A. A is modified by the routine but restored on exit. LDA (input) long The leading dimension of the array A. LDA >= max(1,K). TAU (input) double array, dimension (K) TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by NUMlapack_dgelqf. C (input/output) double array, dimension (LDC,N) On entry, the M-by-N matrix C. On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q. LDC (input) long The leading dimension of the array C. LDC >= max(1,M). WORK (workspace/output) double array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) long The dimension of the array WORK. If SIDE = 'L', LWORK >= max(1,N); if SIDE = 'R', LWORK >= max(1,M). For optimum performance LWORK >= N*NB if SIDE = 'L', and LWORK >= M*NB if SIDE = 'R', 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) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value ===================================================================== */ int NUMlapack_dormqr (char *side, char *trans, long *m, long *n, long *k, double *a, long *lda, double *tau, double *c, long *ldc, double *work, long *lwork, long *info); /* Purpose ======= NUMlapack_dormqr overwrites the general real M-by-N matrix C with SIDE = 'L' SIDE = 'R' TRANS = 'N': Q * C C * Q TRANS = 'T': Q**T * C C * Q**T where Q is a real orthogonal matrix defined as the product of k elementary reflectors Q = H(1) H(2) . . . H(k) as returned by DGEQRF. Q is of order M if SIDE = 'L' and of order N if SIDE = 'R'. Arguments ========= SIDE (input) char* = 'L': apply Q or Q**T from the Left; = 'R': apply Q or Q**T from the Right. TRANS (input) char* = 'N': No transpose, apply Q; = 'T': Transpose, apply Q**T. M (input) long The number of rows of the matrix C. M >= 0. N (input) long The number of columns of the matrix C. N >= 0. K (input) long The number of elementary reflectors whose product defines the matrix Q. If SIDE = 'L', M >= K >= 0; if SIDE = 'R', N >= K >= 0. A (input) double array, dimension (LDA,K) The i-th column must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by DGEQRF in the first k columns of its array argument A. A is modified by the routine but restored on exit. LDA (input) long The leading dimension of the array A. If SIDE = 'L', LDA >= max(1,M); if SIDE = 'R', LDA >= max(1,N). TAU (input) double array, dimension (K) TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by DGEQRF. C (input/output) double array, dimension (LDC,N) On entry, the M-by-N matrix C. On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q. LDC (input) long The leading dimension of the array C. LDC >= max(1,M). WORK (workspace/output) double array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) long The dimension of the array WORK. If SIDE = 'L', LWORK >= max(1,N); if SIDE = 'R', LWORK >= max(1,M). For optimum performance LWORK >= N*NB if SIDE = 'L', and LWORK >= M*NB if SIDE = 'R', 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) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value ===================================================================== */ int NUMlapack_dormr2 (char *side, char *trans, long *m, long *n, long *k, double *a, long *lda, double *tau, double *c, long *ldc, double *work, long *info); /* Purpose ======= NUMlapack_dormr2 overwrites the general real m by n matrix C with Q * C if SIDE = 'L' and TRANS = 'N', or Q'* C if SIDE = 'L' and TRANS = 'T', or C * Q if SIDE = 'R' and TRANS = 'N', or C * Q' if SIDE = 'R' and TRANS = 'T', where Q is a real orthogonal matrix defined as the product of k elementary reflectors Q = H(1) H(2) . . . H(k) as returned by DGERQF. Q is of order m if SIDE = 'L' and of order n if SIDE = 'R'. Arguments ========= SIDE (input) char* = 'L': apply Q or Q' from the Left = 'R': apply Q or Q' from the Right TRANS (input) char* = 'N': apply Q (No transpose) = 'T': apply Q' (Transpose) M (input) long The number of rows of the matrix C. M >= 0. N (input) long The number of columns of the matrix C. N >= 0. K (input) long The number of elementary reflectors whose product defines the matrix Q. If SIDE = 'L', M >= K >= 0; if SIDE = 'R', N >= K >= 0. A (input) double array, dimension (LDA,M) if SIDE = 'L', (LDA,N) if SIDE = 'R' The i-th row must contain the vector which defines the elementary reflector H(i), for i = 1,2,...,k, as returned by DGERQF in the last k rows of its array argument A. A is modified by the routine but restored on exit. LDA (input) long The leading dimension of the array A. LDA >= max(1,K). TAU (input) double array, dimension (K) TAU(i) must contain the scalar factor of the elementary reflector H(i), as returned by DGERQF. C (input/output) double array, dimension (LDC,N) On entry, the m by n matrix C. On exit, C is overwritten by Q*C or Q'*C or C*Q' or C*Q. LDC (input) long The leading dimension of the array C. LDC >= max(1,M). WORK (workspace) double array, dimension (N) if SIDE = 'L', (M) if SIDE = 'R' INFO (output) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value ===================================================================== */ int NUMlapack_dpotf2 (char *uplo, long *n, double *a, long *lda, long *info); /* -- LAPACK routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University February 29, 1992 Purpose ======= NUMlapack_dpotf2 computes the Cholesky factorization of a real symmetric positive definite matrix A. The factorization has the form A = U' * U , if UPLO = 'U', or A = L * L', if UPLO = 'L', where U is an upper triangular matrix and L is lower triangular. This is the unblocked version of the algorithm, calling Level 2 BLAS. Arguments ========= UPLO (input) char * Specifies whether the upper or lower triangular part of the symmetric matrix A is stored. = 'U': Upper triangular = 'L': Lower triangular N (input) long The order of the matrix A. N >= 0. A (input/output) double array, dimension (LDA,N) On entry, the symmetric matrix A. If UPLO = 'U', the leading n by n upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading n by n lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. On exit, if INFO = 0, the factor U or L from the Cholesky factorization A = U'*U or A = L*L'. LDA (input) long The leading dimension of the array A. LDA >= max(1,N). INFO (output) long = 0: successful exit < 0: if INFO = -k, the k-th argument had an illegal value > 0: if INFO = k, the leading minor of order k is not positive definite, and the factorization could not be completed. */ int NUMlapack_drscl (long *n, double *sa, double *sx, long *incx); /* Purpose ======= NUMlapack_drscl multiplies an n-element real vector x by the real scalar 1/a. This is done without overflow or underflow as long as the final result x/a does not overflow or underflow. Arguments ========= N (input) long The number of components of the vector x. SA (input) double The scalar a which is used to divide each component of x. SA must be >= 0, or the subroutine will divide by zero. SX (input/output) double array, dimension (1+(N-1)*abs(INCX)) The n-element vector x. INCX (input) long The increment between successive values of the vector SX. > 0: SX(1) = X(1) and SX(1+(i-1)*INCX) = x(i), 1< i<= n ===================================================================== */ int NUMlapack_dsteqr (char *compz, long *n, double *d, double *e, double *z, long *ldz, double *work, long *info); /* Purpose ======= NUMlapack_dsteqr computes all eigenvalues and, optionally, eigenvectors of a symmetric tridiagonal matrix using the implicit QL or QR method. The eigenvectors of a full or band symmetric matrix can also be found if NUMlapack_dsytrd or DSPTRD or DSBTRD has been used to reduce this matrix to tridiagonal form. Arguments ========= COMPZ (input) char* = 'N': Compute eigenvalues only. = 'V': Compute eigenvalues and eigenvectors of the original symmetric matrix. On entry, Z must contain the orthogonal matrix used to reduce the original matrix to tridiagonal form. = 'I': Compute eigenvalues and eigenvectors of the tridiagonal matrix. Z is initialized to the identity matrix. N (input) long The order of the matrix. N >= 0. D (input/output) double array, dimension (N) On entry, the diagonal elements of the tridiagonal matrix. On exit, if INFO = 0, the eigenvalues in ascending order. E (input/output) double array, dimension (N-1) On entry, the (n-1) subdiagonal elements of the tridiagonal matrix. On exit, E has been destroyed. Z (input/output) double array, dimension (LDZ, N) On entry, if COMPZ = 'V', then Z contains the orthogonal matrix used in the reduction to tridiagonal form. On exit, if INFO = 0, then if COMPZ = 'V', Z contains the orthonormal eigenvectors of the original symmetric matrix, and if COMPZ = 'I', Z contains the orthonormal eigenvectors of the symmetric tridiagonal matrix. If COMPZ = 'N', then Z is not referenced. LDZ (input) long The leading dimension of the array Z. LDZ >= 1, and if eigenvectors are desired, then LDZ >= max(1,N). WORK (workspace) double array, dimension (max(1,2*N-2)) If COMPZ = 'N', then WORK is not referenced. INFO (output) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: the algorithm has failed to find all the eigenvalues in a total of 30*N iterations; if INFO = i, then i elements of E have not converged to zero; on exit, D and E contain the elements of a symmetric tridiagonal matrix which is orthogonally similar to the original matrix. ===================================================================== */ int NUMlapack_dsterf (long *n, double *d, double *e, long *info); /* Purpose ======= NUMlapack_dsterf computes all eigenvalues of a symmetric tridiagonal matrix using the Pal-Walker-Kahan variant of the QL or QR algorithm. Arguments ========= N (input) long The order of the matrix. N >= 0. D (input/output) double array, dimension (N) On entry, the n diagonal elements of the tridiagonal matrix. On exit, if INFO = 0, the eigenvalues in ascending order. E (input/output) double array, dimension (N-1) On entry, the (n-1) subdiagonal elements of the tridiagonal matrix. On exit, E has been destroyed. INFO (output) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: the algorithm failed to find all of the eigenvalues in a total of 30*N iterations; if INFO = i, then i elements of E have not converged to zero. ===================================================================== */ int NUMlapack_dsyev (char *jobz, char *uplo, long *n, double *a, long *lda, double *w, double *work, long *lwork, long *info); /* Purpose ======= NUMlapack_dsyev computes all eigenvalues and, optionally, eigenvectors of a real symmetric matrix A. Arguments ========= JOBZ (input) char* = 'N': Compute eigenvalues only; = 'V': Compute eigenvalues and eigenvectors. UPLO (input) char* = 'U': Upper triangle of A is stored; = 'L': Lower triangle of A is stored. N (input) long The order of the matrix A. N >= 0. A (input/output) double array, dimension (LDA, N) On entry, the symmetric matrix A. If UPLO = 'U', the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A. If UPLO = 'L', the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A. On exit, if JOBZ = 'V', then if INFO = 0, A contains the orthonormal eigenvectors of the matrix A. If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') or the upper triangle (if UPLO='U') of A, including the diagonal, is destroyed. LDA (input) long The leading dimension of the array A. LDA >= max(1,N). W (output) double array, dimension (N) If INFO = 0, the eigenvalues in ascending order. WORK (workspace/output) double array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) long The length of the array WORK. LWORK >= max(1,3*N-1). For optimal efficiency, LWORK >= (NB+2)*N, where NB is the blocksize for NUMlapack_dsytrd returned by NUMlapack_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) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, the algorithm failed to converge; i off-diagonal elements of an intermediate tridiagonal form did not converge to zero. ===================================================================== */ int NUMlapack_dsytd2 (char *uplo, long *n, double *a, long *lda, double *d, double *e, double *tau, long *info); /* Purpose ======= NUMlapack_dsytd2 reduces a real symmetric matrix A to symmetric tridiagonal form T by an orthogonal similarity transformation: Q' * A * Q = T. Arguments ========= UPLO (input) char* Specifies whether the upper or lower triangular part of the symmetric matrix A is stored: = 'U': Upper triangular = 'L': Lower triangular N (input) long The order of the matrix A. N >= 0. A (input/output) double array, dimension (LDA,N) On entry, the symmetric matrix A. If UPLO = 'U', the leading n-by-n upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading n-by-n lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. On exit, if UPLO = 'U', the diagonal and first superdiagonal of A are overwritten by the corresponding elements of the tridiagonal matrix T, and the elements above the first superdiagonal, with the array TAU, represent the orthogonal matrix Q as a product of elementary reflectors; if UPLO = 'L', the diagonal and first subdiagonal of A are over- written by the corresponding elements of the tridiagonal matrix T, 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) long The leading dimension of the array A. LDA >= max(1,N). D (output) double array, dimension (N) The diagonal elements of the tridiagonal matrix T: D(i) = A(i,i). E (output) double array, dimension (N-1) The off-diagonal elements of the tridiagonal matrix T: E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. TAU (output) double array, dimension (N-1) The scalar factors of the elementary reflectors (see Further Details). INFO (output) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value. Further Details =============== If UPLO = 'U', the matrix Q is represented as a product of elementary reflectors Q = H(n-1) . . . H(2) H(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(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in A(1:i-1,i+1), and tau in TAU(i). If UPLO = 'L', the matrix Q is represented as a product of elementary reflectors Q = H(1) H(2) . . . H(n-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 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i), and tau in TAU(i). The contents of A on exit are illustrated by the following examples with n = 5: if UPLO = 'U': if UPLO = 'L': ( d e v2 v3 v4 ) ( d ) ( d e v3 v4 ) ( e d ) ( d e v4 ) ( v1 e d ) ( d e ) ( v1 v2 e d ) ( d ) ( v1 v2 v3 e d ) where d and e denote diagonal and off-diagonal elements of T, and vi denotes an element of the vector defining H(i). ===================================================================== */ int NUMlapack_dsytrd (char *uplo, long *n, double *a, long *lda, double *d, double *e, double *tau, double *work, long *lwork, long *info); /* Purpose ======= NUMlapack_dsytrd reduces a real symmetric matrix A to real symmetric tridiagonal form T by an orthogonal similarity transformation: Q**T * A * Q = T. Arguments ========= UPLO (input) char* = 'U': Upper triangle of A is stored; = 'L': Lower triangle of A is stored. N (input) long The order of the matrix A. N >= 0. A (input/output) double array, dimension (LDA,N) On entry, the symmetric matrix A. If UPLO = 'U', the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. On exit, if UPLO = 'U', the diagonal and first superdiagonal of A are overwritten by the corresponding elements of the tridiagonal matrix T, and the elements above the first superdiagonal, with the array TAU, represent the orthogonal matrix Q as a product of elementary reflectors; if UPLO = 'L', the diagonal and first subdiagonal of A are over- written by the corresponding elements of the tridiagonal matrix T, 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) long The leading dimension of the array A. LDA >= max(1,N). D (output) double array, dimension (N) The diagonal elements of the tridiagonal matrix T: D(i) = A(i,i). E (output) double array, dimension (N-1) The off-diagonal elements of the tridiagonal matrix T: E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. TAU (output) double array, dimension (N-1) The scalar factors of the elementary reflectors (see Further Details). WORK (workspace/output) double array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK (input) long The dimension of the array WORK. LWORK >= 1. 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) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value Further Details =============== If UPLO = 'U', the matrix Q is represented as a product of elementary reflectors Q = H(n-1) . . . H(2) H(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(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in A(1:i-1,i+1), and tau in TAU(i). If UPLO = 'L', the matrix Q is represented as a product of elementary reflectors Q = H(1) H(2) . . . H(n-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 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i), and tau in TAU(i). The contents of A on exit are illustrated by the following examples with n = 5: if UPLO = 'U': if UPLO = 'L': ( d e v2 v3 v4 ) ( d ) ( d e v3 v4 ) ( e d ) ( d e v4 ) ( v1 e d ) ( d e ) ( v1 v2 e d ) ( d ) ( v1 v2 v3 e d ) where d and e denote diagonal and off-diagonal elements of T, and vi denotes an element of the vector defining H(i). ===================================================================== */ int NUMlapack_dtgsja(char *jobu, char *jobv, char *jobq, long *m, long *p, long *n, long *k, long *l, double *a, long *lda, double *b, long *ldb, double *tola, double *tolb, double *alpha, double *beta, double *u, long *ldu, double *v, long *ldv, double *q, long *ldq, double *work, long *ncycle, long *info); /* Purpose ======= NUMlapack_dtgsja computes the generalized singular value decomposition (GSVD) of two real upper triangular (or trapezoidal) matrices A and B. On entry, it is assumed that matrices A and B have the following forms, which may be obtained by the preprocessing subroutine NUMlapack_dggsvp from a general M-by-N matrix A and P-by-N matrix B: N-K-L K L A = K ( 0 A12 A13 ) if M-K-L >= 0; L ( 0 0 A23 ) M-K-L ( 0 0 0 ) N-K-L K L A = K ( 0 A12 A13 ) if M-K-L < 0; M-K ( 0 0 A23 ) N-K-L K L B = L ( 0 0 B13 ) P-L ( 0 0 0 ) where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0, otherwise A23 is (M-K)-by-L upper trapezoidal. On exit, U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R ), where U, V and Q are orthogonal matrices, Z' denotes the transpose of Z, R is a nonsingular upper triangular matrix, and D1 and D2 are ``diagonal'' matrices, which are of the following structures: If M-K-L >= 0, K L D1 = K ( I 0 ) L ( 0 C ) M-K-L ( 0 0 ) K L D2 = L ( 0 S ) P-L ( 0 0 ) N-K-L K L ( 0 R ) = K ( 0 R11 R12 ) K L ( 0 0 R22 ) L where C = diag( ALPHA(K+1), ... , ALPHA(K+L) ), S = diag( BETA(K+1), ... , BETA(K+L) ), C**2 + S**2 = I. R is stored in A(1:K+L,N-K-L+1:N) on exit. If M-K-L < 0, K M-K K+L-M D1 = K ( I 0 0 ) M-K ( 0 C 0 ) K M-K K+L-M D2 = M-K ( 0 S 0 ) K+L-M ( 0 0 I ) P-L ( 0 0 0 ) N-K-L K M-K K+L-M ( 0 R ) = K ( 0 R11 R12 R13 ) M-K ( 0 0 R22 R23 ) K+L-M ( 0 0 0 R33 ) where C = diag( ALPHA(K+1), ... , ALPHA(M) ), S = diag( BETA(K+1), ... , BETA(M) ), C**2 + S**2 = I. R = ( R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N) and R33 is stored ( 0 R22 R23 ) in B(M-K+1:L,N+M-K-L+1:N) on exit. The computation of the orthogonal transformation matrices U, V or Q is optional. These matrices may either be formed explicitly, or they may be postmultiplied into input matrices U1, V1, or Q1. Arguments ========= JOBU (input) char* = 'U': U must contain an orthogonal matrix U1 on entry, and the product U1*U is returned; = 'I': U is initialized to the unit matrix, and the orthogonal matrix U is returned; = 'N': U is not computed. JOBV (input) char* = 'V': V must contain an orthogonal matrix V1 on entry, and the product V1*V is returned; = 'I': V is initialized to the unit matrix, and the orthogonal matrix V is returned; = 'N': V is not computed. JOBQ (input) char* = 'Q': Q must contain an orthogonal matrix Q1 on entry, and the product Q1*Q is returned; = 'I': Q is initialized to the unit matrix, and the orthogonal matrix Q is returned; = 'N': Q is not computed. M (input) long The number of rows of the matrix A. M >= 0. P (input) long The number of rows of the matrix B. P >= 0. N (input) long The number of columns of the matrices A and B. N >= 0. K (input) long L (input) long K and L specify the subblocks in the input matrices A and B: A23 = A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,N-L+1:N) of A and B, whose GSVD is going to be computed by NUMlapack_dtgsja. See Further details. A (input/output) double array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, A(N-K+1:N,1:MIN(K+L,M) ) contains the triangular matrix R or part of R. See Purpose for details. LDA (input) long The leading dimension of the array A. LDA >= max(1,M). B (input/output) double array, dimension (LDB,N) On entry, the P-by-N matrix B. On exit, if necessary, B(M-K+1:L,N+M-K-L+1:N) contains a part of R. See Purpose for details. LDB (input) long The leading dimension of the array B. LDB >= max(1,P). TOLA (input) double TOLB (input) double TOLA and TOLB are the convergence criteria for the Jacobi- Kogbetliantz iteration procedure. Generally, they are the same as used in the preprocessing step, say TOLA = max(M,N)*norm(A)*MAZHEPS, TOLB = max(P,N)*norm(B)*MAZHEPS. ALPHA (output) double array, dimension (N) BETA (output) double array, dimension (N) On exit, ALPHA and BETA contain the generalized singular value pairs of A and B; ALPHA(1:K) = 1, BETA(1:K) = 0, and if M-K-L >= 0, ALPHA(K+1:K+L) = diag(C), BETA(K+1:K+L) = diag(S), or if M-K-L < 0, ALPHA(K+1:M)= C, ALPHA(M+1:K+L)= 0 BETA(K+1:M) = S, BETA(M+1:K+L) = 1. Furthermore, if K+L < N, ALPHA(K+L+1:N) = 0 and BETA(K+L+1:N) = 0. U (input/output) double array, dimension (LDU,M) On entry, if JOBU = 'U', U must contain a matrix U1 (usually the orthogonal matrix returned by NUMlapack_dggsvp). On exit, if JOBU = 'I', U contains the orthogonal matrix U; if JOBU = 'U', U contains the product U1*U. If JOBU = 'N', U is not referenced. LDU (input) long The leading dimension of the array U. LDU >= max(1,M) if JOBU = 'U'; LDU >= 1 otherwise. V (input/output) double array, dimension (LDV,P) On entry, if JOBV = 'V', V must contain a matrix V1 (usually the orthogonal matrix returned by NUMlapack_dggsvp). On exit, if JOBV = 'I', V contains the orthogonal matrix V; if JOBV = 'V', V contains the product V1*V. If JOBV = 'N', V is not referenced. LDV (input) long The leading dimension of the array V. LDV >= max(1,P) if JOBV = 'V'; LDV >= 1 otherwise. Q (input/output) double array, dimension (LDQ,N) On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually the orthogonal matrix returned by NUMlapack_dggsvp). On exit, if JOBQ = 'I', Q contains the orthogonal matrix Q; if JOBQ = 'Q', Q contains the product Q1*Q. If JOBQ = 'N', Q is not referenced. LDQ (input) long The leading dimension of the array Q. LDQ >= max(1,N) if JOBQ = 'Q'; LDQ >= 1 otherwise. WORK (workspace) double array, dimension (2*N) NCYCLE (output) long The number of cycles required for convergence. INFO (output) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value. = 1: the procedure does not converge after MAXIT cycles. Internal Parameters =================== MAXIT long MAXIT specifies the total loops that the iterative procedure may take. If after MAXIT cycles, the routine fails to converge, we return INFO = 1. Further Details =============== NUMlapack_dtgsja essentially uses a variant of Kogbetliantz algorithm to reduce min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and L-by-L matrix B13 to the form: U1'*A13*Q1 = C1*R1; V1'*B13*Q1 = S1*R1, where U1, V1 and Q1 are orthogonal matrix, and Z' is the transpose of Z. C1 and S1 are diagonal matrices satisfying C1**2 + S1**2 = I, and R1 is an L-by-L nonsingular upper triangular matrix. ===================================================================== */ int NUMlapack_dtrevc (char *side, char *howmny, int * select, long *n, double *t, long *ldt, double *vl, long *ldvl, double *vr, long *ldvr, long *mm, long *m, double *work, long *info); /* -- LAPACK routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University June 30, 1999 Purpose ======= NUMlapack_dtrevc computes some or all of the right and/or left eigenvectors of a real upper quasi-triangular matrix T. The right eigenvector x and the left eigenvector y of T corresponding to an eigenvalue w are defined by: T*x = w*x, y'*T = w*y' where y' denotes the conjugate transpose of the vector y. If all eigenvectors are requested, the routine may either return the matrices X and/or Y of right or left eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an input orthogonal matrix. If T was obtained from the real-Schur factorization of an original matrix A = Q*T*Q', then Q*X and Q*Y are the matrices of right or left eigenvectors of A. T must be in Schur canonical form (as returned by NUMlapack_dhseqr), that is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block has its diagonal elements equal and its off-diagonal elements of opposite sign. Corresponding to each 2-by-2 diagonal block is a complex conjugate pair of eigenvalues and eigenvectors; only one eigenvector of the pair is computed, namely the one corresponding to the eigenvalue with positive imaginary part. Arguments ========= SIDE (input) char* = 'R': compute right eigenvectors only; = 'L': compute left eigenvectors only; = 'B': compute both right and left eigenvectors. HOWMNY (input) char* = 'A': compute all right and/or left eigenvectors; = 'B': compute all right and/or left eigenvectors, and backtransform them using the input matrices supplied in VR and/or VL; = 'S': compute selected right and/or left eigenvectors, specified by the int array SELECT. SELECT (input/output) int array, dimension (N) If HOWMNY = 'S', SELECT specifies the eigenvectors to be computed. If HOWMNY = 'A' or 'B', SELECT is not referenced. To select the real eigenvector corresponding to a real eigenvalue w(j), SELECT(j) must be set to .TRUE.. To select the complex eigenvector corresponding to a complex conjugate pair w(j) and w(j+1), either SELECT(j) or SELECT(j+1) must be set to .TRUE.; then on exit SELECT(j) is .TRUE. and SELECT(j+1) is .FALSE.. N (input) long The order of the matrix T. N >= 0. T (input) double array, dimension (LDT,N) The upper quasi-triangular matrix T in Schur canonical form. LDT (input) long The leading dimension of the array T. LDT >= max(1,N). VL (input/output) double array, dimension (LDVL,MM) On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must contain an N-by-N matrix Q (usually the orthogonal matrix Q of Schur vectors returned by NUMlapack_dhseqr). On exit, if SIDE = 'L' or 'B', VL contains: if HOWMNY = 'A', the matrix Y of left eigenvectors of T; VL has the same quasi-lower triangular form as T'. If T(i,i) is a real eigenvalue, then the i-th column VL(i) of VL is its corresponding eigenvector. If T(i:i+1,i:i+1) is a 2-by-2 block whose eigenvalues are complex-conjugate eigenvalues of T, then VL(i)+sqrt(-1)*VL(i+1) is the complex eigenvector corresponding to the eigenvalue with positive real part. if HOWMNY = 'B', the matrix Q*Y; if HOWMNY = 'S', the left eigenvectors of T specified by SELECT, stored consecutively in the columns of VL, in the same order as their eigenvalues. A complex eigenvector corresponding to a complex eigenvalue is stored in two consecutive columns, the first holding the real part, and the second the imaginary part. If SIDE = 'R', VL is not referenced. LDVL (input) long The leading dimension of the array VL. LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise. VR (input/output) double array, dimension (LDVR,MM) On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must contain an N-by-N matrix Q (usually the orthogonal matrix Q of Schur vectors returned by NUMlapack_dhseqr). On exit, if SIDE = 'R' or 'B', VR contains: if HOWMNY = 'A', the matrix X of right eigenvectors of T; VR has the same quasi-upper triangular form as T. If T(i,i) is a real eigenvalue, then the i-th column VR(i) of VR is its corresponding eigenvector. If T(i:i+1,i:i+1) is a 2-by-2 block whose eigenvalues are complex-conjugate eigenvalues of T, then VR(i)+sqrt(-1)*VR(i+1) is the complex eigenvector corresponding to the eigenvalue with positive real part. if HOWMNY = 'B', the matrix Q*X; if HOWMNY = 'S', the right eigenvectors of T specified by SELECT, stored consecutively in the columns of VR, in the same order as their eigenvalues. A complex eigenvector corresponding to a complex eigenvalue is stored in two consecutive columns, the first holding the real part and the second the imaginary part. If SIDE = 'L', VR is not referenced. LDVR (input) long The leading dimension of the array VR. LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise. MM (input) long The number of columns in the arrays VL and/or VR. MM >= M. M (output) long The number of columns in the arrays VL and/or VR actually used to store the eigenvectors. If HOWMNY = 'A' or 'B', M is set to N. Each selected real eigenvector occupies one column and each selected complex eigenvector occupies two columns. WORK (workspace) double 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 algorithm used in this program is basically backward (forward) substitution, with scaling to make the the code robust against possible overflow. Each eigenvector is normalized so that the element of largest magnitude has magnitude 1; here the magnitude of a complex number (x,y) is taken to be |x| + |y|. ===================================================================== */ int NUMlapack_dtrti2 (char *uplo, char *diag, long *n, double *a, long *lda, long *info); /* Purpose ======= NUMlapack_dtrti2 computes the inverse of a real upper or lower triangular matrix. This is the Level 2 BLAS version of the algorithm. Arguments ========= UPLO (input) char* Specifies whether the matrix A is upper or lower triangular. = 'U': Upper triangular = 'L': Lower triangular DIAG (input) char* Specifies whether or not the matrix A is unit triangular. = 'N': Non-unit triangular = 'U': Unit triangular N (input) long The order of the matrix A. N >= 0. A (input/output) double array, dimension (LDA,N) On entry, the triangular matrix A. If UPLO = 'U', the leading n by n upper triangular part of the array A contains the upper triangular matrix, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading n by n lower triangular part of the array A contains the lower triangular matrix, and the strictly upper triangular part of A is not referenced. If DIAG = 'U', the diagonal elements of A are also not referenced and are assumed to be 1. On exit, the (triangular) inverse of the original matrix, in the same storage format. LDA (input) long The leading dimension of the array A. LDA >= max(1,N). INFO (output) long = 0: successful exit < 0: if INFO = -k, the k-th argument had an illegal value ===================================================================== */ int NUMlapack_dtrtri (char *uplo, char *diag, long *n, double * a, long *lda, long *info); /* Purpose ======= NUMlapack_dtrtri computes the inverse of a real upper or lower triangular matrix A. This is the Level 3 BLAS version of the algorithm. Arguments ========= UPLO (input) char* = 'U': A is upper triangular; = 'L': A is lower triangular. DIAG (input) char* = 'N': A is non-unit triangular; = 'U': A is unit triangular. N (input) long The order of the matrix A. N >= 0. A (input/output) double array, dimension (LDA,N) On entry, the triangular matrix A. If UPLO = 'U', the leading N-by-N upper triangular part of the array A contains the upper triangular matrix, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading N-by-N lower triangular part of the array A contains the lower triangular matrix, and the strictly upper triangular part of A is not referenced. If DIAG = 'U', the diagonal elements of A are also not referenced and are assumed to be 1. On exit, the (triangular) inverse of the original matrix, in the same storage format. LDA (input) long The leading dimension of the array A. LDA >= max(1,N). INFO (output) long = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, A(i,i) is exactly zero. The triangular matrix is singular and its inverse can not be computed. ===================================================================== */ long NUMlapack_ieeeck (long *ispec, float *zero, float *one); /* Purpose ======= NUMlapack_ieeeck is called from the NUMlapack_ilaenv to verify that Infinity and possibly NaN arithmetic is safe (i.e. will not trap). Arguments ========= ISPEC (input) long Specifies whether to test just for inifinity arithmetic or whether to test for infinity and NaN arithmetic. = 0: Verify infinity arithmetic only. = 1: Verify infinity and NaN arithmetic. ZERO (input) REAL Must contain the value 0.0 This is passed to prevent the compiler from optimizing away this code. ONE (input) REAL Must contain the value 1.0 This is passed to prevent the compiler from optimizing away this code. RETURN VALUE: long = 0: Arithmetic failed to produce the correct answers = 1: Arithmetic produced the correct answers */ long NUMlapack_ilaenv (long *ispec, char *name, char *opts, long *n1, long *n2, long *n3, long *n4, long name_len, long opts_len); /* Purpose ======= NUMlapack_ilaenv is called from the LAPACK routines to choose problem-dependent parameters for the local environment. See ISPEC for a description of the parameters. This version provides a set of parameters which should give good, but not optimal, performance on many of the currently available computers. Users are encouraged to modify this subroutine to set the tuning parameters for their particular machine using the option and problem size information in the arguments. This routine will not function correctly if it is converted to all lower case. Converting it to all upper case is allowed. Arguments ========= ISPEC (input) long Specifies the parameter to be returned as the value of NUMlapack_ilaenv. = 1: the optimal blocksize; if this value is 1, an unblocked algorithm will give the best performance. = 2: the minimum block size for which the block routine should be used; if the usable block size is less than this value, an unblocked routine should be used. = 3: the crossover point (in a block routine, for N less than this value, an unblocked routine should be used) = 4: the number of shifts, used in the nonsymmetric eigenvalue routines = 5: the minimum column dimension for blocking to be used; rectangular blocks must have dimension at least k by m, where k is given by NUMlapack_ilaenv(2,...) and m by NUMlapack_ilaenv(5,...) = 6: the crossover point for the SVD (when reducing an m by n matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds this value, a QR factorization is used first to reduce the matrix to a triangular form.) = 7: the number of processors = 8: the crossover point for the multishift QR and QZ methods for nonsymmetric eigenvalue problems. = 9: maximum size of the subproblems at the bottom of the computation tree in the divide-and-conquer algorithm (used by xGELSD and xGESDD) =10: ieee NaN arithmetic can be trusted not to trap =11: infinity arithmetic can be trusted not to trap NAME (input) CHARACTER*(*) The name of the calling subroutine, in either upper case or lower case. OPTS (input) CHARACTER*(*) The character options to the subroutine NAME, concatenated into a single character string. For example, UPLO = 'U', TRANS = 'T', and DIAG = 'N' for a triangular routine would be specified as OPTS = 'UTN'. N1 (input) long N2 (input) long N3 (input) long N4 (input) long Problem dimensions for the subroutine NAME; these may not all be required. (NUMlapack_ilaenv) (output) long >= 0: the value of the parameter specified by ISPEC < 0: if NUMlapack_ilaenv = -k, the k-th argument had an illegal value. Further Details =============== The following conventions have been used when calling NUMlapack_ilaenv from the LAPACK routines: 1) OPTS is a concatenation of all of the character options to subroutine NAME, in the same order that they appear in the argument list for NAME, even if they are not used in determining the value of the parameter specified by ISPEC. 2) The problem dimensions N1, N2, N3, N4 are specified in the order that they appear in the argument list for NAME. N1 is used first, N2 second, and so on, and unused problem dimensions are passed a value of -1. 3) The parameter value returned by NUMlapack_ilaenv is checked for validity in the calling subroutine. For example, NUMlapack_ilaenv is used to retrieve the optimal blocksize for STRTRI as follows: NB = NUMlapack_ilaenv( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 ) IF( NB.LE.1 ) NB = MAX( 1, N ) ===================================================================== */ #endif /* _NUMclapack_h_ */