// ------------------------------------------------------------ // lapack_cwrap.h // Wrapping LAPACK FORTRAN-function to C-Functions // // 26.02.2005, J.Ahrensfeld // ------------------------------------------------------------ #ifndef LAPACK_CWRAP_H #define LAPACK_CWRAP_H // ------------------------------------------------------------ // LAPACK functions // ------------------------------------------------------------ int dgtsv(int N, int NRHS, double *DL, double *D, double *DU, double *B, int ldb); // ------------------------------------------------------------ int dgesv(int N, int NHRS, double *A, int LDA, int *IPIV, double *B, int LDB); // ------------------------------------------------------------ // SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO ) // // -- LAPACK driver routine (version 3.0) -- // Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., // Courant Institute, Argonne National Lab, and Rice University // March 31, 1993 // // .. Scalar Arguments .. // INTEGER INFO, LDA, LDB, N, NRHS // .. // .. Array Arguments .. // INTEGER IPIV( * ) // DOUBLE PRECISION A( LDA, * ), B( LDB, * ) // .. // // Purpose // ======= // // DGESV computes the solution to a real system of linear equations // A * X = B, // where A is an N-by-N matrix and X and B are N-by-NRHS matrices. // // The LU decomposition with partial pivoting and row interchanges is // used to factor A as // A = P * L * U, // where P is a permutation matrix, L is unit lower triangular, and U is // upper triangular. The factored form of A is then used to solve the // system of equations A * X = B. // // Arguments // ========= // // N (input) INTEGER // The number of linear equations, i.e., the order of the // matrix A. N >= 0. // // NRHS (input) INTEGER // The number of right hand sides, i.e., the number of columns // of the matrix B. NRHS >= 0. // // A (input/output) DOUBLE PRECISION array, dimension (LDA,N) // On entry, the N-by-N coefficient matrix A. // On exit, the factors L and U from the factorization // A = P*L*U; the unit diagonal elements of L are not stored. // // LDA (input) INTEGER // The leading dimension of the array A. LDA >= max(1,N). // // IPIV (output) INTEGER array, dimension (N) // The pivot indices that define the permutation matrix P; // row i of the matrix was interchanged with row IPIV(i). // // B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) // On entry, the N-by-NRHS matrix of right hand side matrix B. // On exit, if INFO = 0, the N-by-NRHS solution matrix X. // // LDB (input) INTEGER // The leading dimension of the array B. LDB >= max(1,N). // // INFO (output) INTEGER // = 0: successful exit // < 0: if INFO = -i, the i-th argument had an illegal value // > 0: if INFO = i, U(i,i) is exactly zero. The factorization // has been completed, but the factor U is exactly // singular, so the solution could not be computed. // ------------------------------------------------------------ int dsyev(char JOBZ, char UPLO, int N, double *A, int LDA, double *W, double *WORK, int LWORK); // ------------------------------------------------------------ // SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO ) // // DSYEV computes all eigenvalues and, optionally, eigenvectors of a // real symmetric matrix A. // // Arguments // ========= // // JOBZ (input) CHARACTER*1 // = 'N': Compute eigenvalues only; // = 'V': Compute eigenvalues and eigenvectors. // // UPLO (input) CHARACTER*1 // = 'U': Upper triangle of A is stored; // = 'L': Lower triangle of A is stored. // // N (input) INTEGER // The order of the matrix A. N >= 0. // // A (input/output) DOUBLE PRECISION 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) INTEGER // The leading dimension of the array A. LDA >= max(1,N). // // W (output) DOUBLE PRECISION array, dimension (N) // If INFO = 0, the eigenvalues in ascending order. // // WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) // On exit, if INFO = 0, WORK(1) returns the optimal LWORK. // // LWORK (input) INTEGER // The length of the array WORK. LWORK >= max(1,3*N-1). // For optimal efficiency, LWORK >= (NB+2)*N, // where NB is the blocksize for DSYTRD returned by ILAENV. // // If LWORK = -1, then a workspace query is assumed; the routine // only calculates the optimal size of the WORK array, returns // this value as the first entry of the WORK array, and no error // message related to LWORK is issued by XERBLA. // // INFO (output) INTEGER // = 0: successful exit // < 0: if INFO = -i, the i-th argument had an illegal value // > 0: if INFO = i, the algorithm failed to converge; i // off-diagonal elements of an intermediate tridiagonal // form did not converge to zero. // ------------------------------------------------------------ int dgebrd(int M, int N, double *A, int LDA, double *D, double *E, double *TAUQ, double *TAUP, double *WORK, int LWORK); // ------------------------------------------------------------ // SUBROUTINE DGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, // $ 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 // // .. Scalar Arguments .. // INTEGER INFO, LDA, LWORK, M, N // .. // .. Array Arguments .. // DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), TAUP( * ), // $ TAUQ( * ), WORK( * ) // .. // // Purpose // ======= // // DGEBRD reduces a general real M-by-N matrix A to upper or lower // bidiagonal form B by an orthogonal transformation: Q**T * A * P = B. // // If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal. // // Arguments // ========= // // M (input) INTEGER // The number of rows in the matrix A. M >= 0. // // N (input) INTEGER // The number of columns in the matrix A. N >= 0. // // A (input/output) DOUBLE PRECISION array, dimension (LDA,N) // On entry, the M-by-N general matrix to be reduced. // On exit, // if m >= n, the diagonal and the first superdiagonal are // overwritten with the upper bidiagonal matrix B; the // elements below the diagonal, with the array TAUQ, represent // the orthogonal matrix Q as a product of elementary // reflectors, and the elements above the first superdiagonal, // with the array TAUP, represent the orthogonal matrix P as // a product of elementary reflectors; // if m < n, the diagonal and the first subdiagonal are // overwritten with the lower bidiagonal matrix B; the // elements below the first subdiagonal, with the array TAUQ, // represent the orthogonal matrix Q as a product of // elementary reflectors, and the elements above the diagonal, // with the array TAUP, represent the orthogonal matrix P as // a product of elementary reflectors. // See Further Details. // // LDA (input) INTEGER // The leading dimension of the array A. LDA >= max(1,M). // // D (output) DOUBLE PRECISION array, dimension (min(M,N)) // The diagonal elements of the bidiagonal matrix B: // D(i) = A(i,i). // // E (output) DOUBLE PRECISION array, dimension (min(M,N)-1) // The off-diagonal elements of the bidiagonal matrix B: // if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1; // if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1. // // TAUQ (output) DOUBLE PRECISION array dimension (min(M,N)) // The scalar factors of the elementary reflectors which // represent the orthogonal matrix Q. See Further Details. // // TAUP (output) DOUBLE PRECISION array, dimension (min(M,N)) // The scalar factors of the elementary reflectors which // represent the orthogonal matrix P. See Further Details. // // WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) // On exit, if INFO = 0, WORK(1) returns the optimal LWORK. // // LWORK (input) INTEGER // The length of the array WORK. LWORK >= max(1,M,N). // For optimum performance LWORK >= (M+N)*NB, where NB // is the optimal blocksize. // // If LWORK = -1, then a workspace query is assumed; the routine // only calculates the optimal size of the WORK array, returns // this value as the first entry of the WORK array, and no error // message related to LWORK is issued by XERBLA. // // INFO (output) INTEGER // = 0: successful exit // < 0: if INFO = -i, the i-th argument had an illegal value. // // ------------------------------------------------------------ int dbdsqr(char UPLO, int N, int NCVT, int NRU, int NCC, double *D, double *E, double *VT, int LDVT, double *U, int LDU, double *C, int LDC, double *WORK); // ------------------------------------------------------------ // SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, // $ LDU, C, LDC, WORK, 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, 1999 // // .. Scalar Arguments .. // CHARACTER UPLO // INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU // .. // .. Array Arguments .. // DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ), // $ VT( LDVT, * ), WORK( * ) // .. // // Purpose // ======= // // 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) CHARACTER*1 // = 'U': B is upper bidiagonal; // = 'L': B is lower bidiagonal. // // N (input) INTEGER // The order of the matrix B. N >= 0. // // NCVT (input) INTEGER // The number of columns of the matrix VT. NCVT >= 0. // // NRU (input) INTEGER // The number of rows of the matrix U. NRU >= 0. // // NCC (input) INTEGER // The number of columns of the matrix C. NCC >= 0. // // D (input/output) DOUBLE PRECISION 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 PRECISION 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 PRECISION 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) INTEGER // The leading dimension of the array VT. // LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0. // // U (input/output) DOUBLE PRECISION 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) INTEGER // The leading dimension of the array U. LDU >= max(1,NRU). // // C (input/output) DOUBLE PRECISION 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) INTEGER // The leading dimension of the array C. // LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0. // // WORK (workspace) DOUBLE PRECISION array, dimension (4*N) // // INFO (output) INTEGER // = 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. // // ------------------------------------------------------------ int dgesvd(char JOBU, char JOBVT, int M, int N, double *A, int LDA, double *S, double *U, int LDU, double *VT, int LDVT, double *WORK, int LWORK); // ------------------------------------------------------------ // SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, // $ WORK, LWORK, INFO ) // // -- LAPACK driver routine (version 3.0) -- // Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., // Courant Institute, Argonne National Lab, and Rice University // October 31, 1999 // // .. Scalar Arguments .. // CHARACTER JOBU, JOBVT // INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N // .. // .. Array Arguments .. // DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ), // $ VT( LDVT, * ), WORK( * ) // .. // // Purpose // ======= // // DGESVD computes the singular value decomposition (SVD) of a real // M-by-N matrix A, optionally computing the left and/or right singular // vectors. The SVD is written // // A = U * SIGMA * transpose(V) // // where SIGMA is an M-by-N matrix which is zero except for its // min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and // V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA // are the singular values of A; they are real and non-negative, and // are returned in descending order. The first min(m,n) columns of // U and V are the left and right singular vectors of A. // // Note that the routine returns V**T, not V. // // Arguments // ========= // // JOBU (input) CHARACTER*1 // Specifies options for computing all or part of the matrix U: // = 'A': all M columns of U are returned in array U: // = 'S': the first min(m,n) columns of U (the left singular // vectors) are returned in the array U; // = 'O': the first min(m,n) columns of U (the left singular // vectors) are overwritten on the array A; // = 'N': no columns of U (no left singular vectors) are // computed. // // JOBVT (input) CHARACTER*1 // Specifies options for computing all or part of the matrix // V**T: // = 'A': all N rows of V**T are returned in the array VT; // = 'S': the first min(m,n) rows of V**T (the right singular // vectors) are returned in the array VT; // = 'O': the first min(m,n) rows of V**T (the right singular // vectors) are overwritten on the array A; // = 'N': no rows of V**T (no right singular vectors) are // computed. // // JOBVT and JOBU cannot both be 'O'. // // M (input) INTEGER // The number of rows of the input matrix A. M >= 0. // // N (input) INTEGER // The number of columns of the input matrix A. N >= 0. // // A (input/output) DOUBLE PRECISION array, dimension (LDA,N) // On entry, the M-by-N matrix A. // On exit, // if JOBU = 'O', A is overwritten with the first min(m,n) // columns of U (the left singular vectors, // stored columnwise); // if JOBVT = 'O', A is overwritten with the first min(m,n) // rows of V**T (the right singular vectors, // stored rowwise); // if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A // are destroyed. // // LDA (input) INTEGER // The leading dimension of the array A. LDA >= max(1,M). // // S (output) DOUBLE PRECISION array, dimension (min(M,N)) // The singular values of A, sorted so that S(i) >= S(i+1). // // U (output) DOUBLE PRECISION array, dimension (LDU,UCOL) // (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'. // If JOBU = 'A', U contains the M-by-M orthogonal matrix U; // if JOBU = 'S', U contains the first min(m,n) columns of U // (the left singular vectors, stored columnwise); // if JOBU = 'N' or 'O', U is not referenced. // // LDU (input) INTEGER // The leading dimension of the array U. LDU >= 1; if // JOBU = 'S' or 'A', LDU >= M. // // VT (output) DOUBLE PRECISION array, dimension (LDVT,N) // If JOBVT = 'A', VT contains the N-by-N orthogonal matrix // V**T; // if JOBVT = 'S', VT contains the first min(m,n) rows of // V**T (the right singular vectors, stored rowwise); // if JOBVT = 'N' or 'O', VT is not referenced. // // LDVT (input) INTEGER // The leading dimension of the array VT. LDVT >= 1; if // JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N). // // WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) // On exit, if INFO = 0, WORK(1) returns the optimal LWORK; // if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged // superdiagonal elements of an upper bidiagonal matrix B // whose diagonal is in S (not necessarily sorted). B // satisfies A = U * B * VT, so it has the same singular values // as A, and singular vectors related by U and VT. // // LWORK (input) INTEGER // The dimension of the array WORK. LWORK >= 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) INTEGER // = 0: successful exit. // < 0: if INFO = -i, the i-th argument had an illegal value. // > 0: if DBDSQR did not converge, INFO specifies how many // superdiagonals of an intermediate bidiagonal form B // did not converge to zero. See the description of WORK // above for details. // int dgecon(char NORM, int N, double *A, int LDA, double ANORM, double *RCOND, double *DWORK, int *IWORK); // SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, // $ 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 // // .. Scalar Arguments .. // CHARACTER NORM // INTEGER INFO, LDA, N // DOUBLE PRECISION ANORM, RCOND // .. // .. Array Arguments .. // INTEGER IWORK( * ) // DOUBLE PRECISION A( LDA, * ), WORK( * ) // .. // // Purpose // ======= // // DGECON estimates the reciprocal of the condition number of a general // real matrix A, in either the 1-norm or the infinity-norm, using // the LU factorization computed by DGETRF. // // An estimate is obtained for norm(inv(A)), and the reciprocal of the // condition number is computed as // RCOND = 1 / ( norm(A) * norm(inv(A)) ). // // Arguments // ========= // // NORM (input) CHARACTER*1 // Specifies whether the 1-norm condition number or the // infinity-norm condition number is required: // = '1' or 'O': 1-norm; // = 'I': Infinity-norm. // // N (input) INTEGER // The order of the matrix A. N >= 0. // // A (input) DOUBLE PRECISION array, dimension (LDA,N) // The factors L and U from the factorization A = P*L*U // as computed by DGETRF. // // LDA (input) INTEGER // The leading dimension of the array A. LDA >= max(1,N). // // ANORM (input) DOUBLE PRECISION // If NORM = '1' or 'O', the 1-norm of the original matrix A. // If NORM = 'I', the infinity-norm of the original matrix A. // // RCOND (output) DOUBLE PRECISION // The reciprocal of the condition number of the matrix A, // computed as RCOND = 1/(norm(A) * norm(inv(A))). // // WORK (workspace) DOUBLE PRECISION array, dimension (4*N) // // IWORK (workspace) INTEGER array, dimension (N) // // INFO (output) INTEGER // = 0: successful exit // ------------------------------------------------------------ // BLAS functions LEVEL 1 // ------------------------------------------------------------ void dswap (int n, double *dx, int incx, double *dy, int incy); void dcopy (int n, double *dx, int incx, double *dy, int incy); void dscal (int N, double da, double *dx, int incx); void daxpy (int n, double da, double *dx, int incx, double *dy, int incy); double ddot (int n, double *dx, int incx, double *dy, int incy); double dnrm2 (int n, double *dx, int incx); double dasum (int n, double *dx, int incx); int idamax (int n, double *dx, int incx); // ------------------------------------------------------------ // BLAS functions LEVEL 2 // ------------------------------------------------------------ // ------------------------------------------------------------ void dgemv(char TRANS, int M, int N, double ALPHA, double *A, int LDA, double *X, int INCX, double BETA, double *Y, int INCY); // ------------------------------------------------------------ // SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY ) // // DGEMV performs one of the matrix-vector operations // // y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, // // where alpha and beta are scalars, x and y are vectors and A is an // m by n matrix. // // Parameters // ========== // // TRANS - CHARACTER*1. // On entry, TRANS specifies the operation to be performed as // follows: // // TRANS = 'N' or 'n' y := alpha*A*x + beta*y. // // TRANS = 'T' or 't' y := alpha*A'*x + beta*y. // // TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. // // Unchanged on exit. // // M - INTEGER. // On entry, M specifies the number of rows of the matrix A. // M must be at least zero. // Unchanged on exit. // // N - INTEGER. // On entry, N specifies the number of columns of the matrix A. // N must be at least zero. // Unchanged on exit. // // ALPHA - DOUBLE PRECISION. // On entry, ALPHA specifies the scalar alpha. // Unchanged on exit. // // A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). // Before entry, the leading m by n part of the array A must // contain the matrix of coefficients. // Unchanged on exit. // // LDA - INTEGER. // On entry, LDA specifies the first dimension of A as declared // in the calling (sub) program. LDA must be at least // max( 1, m ). // Unchanged on exit. // // X - DOUBLE PRECISION array of DIMENSION at least // ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' // and at least // ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. // Before entry, the incremented array X must contain the // vector x. // Unchanged on exit. // // INCX - INTEGER. // On entry, INCX specifies the increment for the elements of // X. INCX must not be zero. // Unchanged on exit. // // BETA - DOUBLE PRECISION. // On entry, BETA specifies the scalar beta. When BETA is // supplied as zero then Y need not be set on input. // Unchanged on exit. // // Y - DOUBLE PRECISION array of DIMENSION at least // ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' // and at least // ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. // Before entry with BETA non-zero, the incremented array Y // must contain the vector y. On exit, Y is overwritten by the // updated vector y. // // INCY - INTEGER. // On entry, INCY specifies the increment for the elements of // Y. INCY must not be zero. // Unchanged on exit. // ------------------------------------------------------------ // BLAS functions LEVEL 2 // ------------------------------------------------------------ void dger(int m,int n, double alpha, double *x, int incx, double *y, int incy, double *A, int lda); // ------------------------------------------------------------ // SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) // DGER performs the rank 1 operation // A := alpha*x*y' + A, // where alpha is a scalar, x is an m element vector, y is an n element // vector and A is an m by n matrix. // // Parameters // ========== // // M - INTEGER. // On entry, M specifies the number of rows of the matrix A. // M must be at least zero. // Unchanged on exit. // // N - INTEGER. // On entry, N specifies the number of columns of the matrix A. // N must be at least zero. // Unchanged on exit. // // ALPHA - DOUBLE PRECISION. // On entry, ALPHA specifies the scalar alpha. // Unchanged on exit. // // X - DOUBLE PRECISION array of dimension at least // ( 1 + ( m - 1 )//abs( INCX ) ). // Before entry, the incremented array X must contain the m // element vector x. // Unchanged on exit. // // INCX - INTEGER. // On entry, INCX specifies the increment for the elements of // X. INCX must not be zero. // Unchanged on exit. // // Y - DOUBLE PRECISION array of dimension at least // ( 1 + ( n - 1 )//abs( INCY ) ). // Before entry, the incremented array Y must contain the n // element vector y. // Unchanged on exit. // // INCY - INTEGER. // On entry, INCY specifies the increment for the elements of // Y. INCY must not be zero. // Unchanged on exit. // // A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). // Before entry, the leading m by n part of the array A must // contain the matrix of coefficients. On exit, A is // overwritten by the updated matrix. // // LDA - INTEGER. // On entry, LDA specifies the first dimension of A as declared // in the calling (sub) program. LDA must be at least // max( 1, m ). // Unchanged on exit. // ------------------------------------------------------------ void dsyr(char UPLO, int n, double alpha, double *x, int incx, double *A, int lda); // ------------------------------------------------------------ // SUBROUTINE DSYR ( UPLO, N, ALPHA, X, INCX, A, LDA ) // .. Scalar Arguments .. // DOUBLE PRECISION ALPHA // INTEGER INCX, LDA, N // CHARACTER*1 UPLO // .. Array Arguments .. // DOUBLE PRECISION A( LDA, * ), X( * ) // .. // // Purpose // ======= // // DSYR performs the symmetric rank 1 operation // // A := alpha*x*x' + A, // // where alpha is a real scalar, x is an n element vector and A is an // n by n symmetric matrix. // // Parameters // ========== // // UPLO - CHARACTER*1. // On entry, UPLO specifies whether the upper or lower // triangular part of the array A is to be referenced as // follows: // // UPLO = 'U' or 'u' Only the upper triangular part of A // is to be referenced. // // UPLO = 'L' or 'l' Only the lower triangular part of A // is to be referenced. // // Unchanged on exit. // // N - INTEGER. // On entry, N specifies the order of the matrix A. // N must be at least zero. // Unchanged on exit. // // ALPHA - DOUBLE PRECISION. // On entry, ALPHA specifies the scalar alpha. // Unchanged on exit. // // X - DOUBLE PRECISION array of dimension at least // ( 1 + ( n - 1 )*abs( INCX ) ). // Before entry, the incremented array X must contain the n // element vector x. // Unchanged on exit. // // INCX - INTEGER. // On entry, INCX specifies the increment for the elements of // X. INCX must not be zero. // Unchanged on exit. // // A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). // Before entry with UPLO = 'U' or 'u', the leading n by n // upper triangular part of the array A must contain the upper // triangular part of the symmetric matrix and the strictly // lower triangular part of A is not referenced. On exit, the // upper triangular part of the array A is overwritten by the // upper triangular part of the updated matrix. // Before entry with UPLO = 'L' or 'l', the leading n by n // lower triangular part of the array A must contain the lower // triangular part of the symmetric matrix and the strictly // upper triangular part of A is not referenced. On exit, the // lower triangular part of the array A is overwritten by the // lower triangular part of the updated matrix. // // LDA - INTEGER. // On entry, LDA specifies the first dimension of A as declared // in the calling (sub) program. LDA must be at least // max( 1, n ). // Unchanged on exit. // ------------------------------------------------------------ // BLAS functions LEVEL 3 // ------------------------------------------------------------ void dsyrk(char UPLO, char TRANS, int N, int K, double ALPHA, double *A, int LDA, double BETA, double *C, int LDC); // ------------------------------------------------------------ // SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA, // $ BETA, C, LDC ) // .. Scalar Arguments .. // CHARACTER*1 UPLO, TRANS // INTEGER N, K, LDA, LDC // DOUBLE PRECISION ALPHA, BETA // .. Array Arguments .. // DOUBLE PRECISION A( LDA, * ), C( LDC, * ) // .. // // Purpose // ======= // // DSYRK performs one of the symmetric rank k operations // // C := alpha*A*A' + beta*C, // // or // // C := alpha*A'*A + beta*C, // // where alpha and beta are scalars, C is an n by n symmetric matrix // and A is an n by k matrix in the first case and a k by n matrix // in the second case. // // Parameters // ========== // // UPLO - CHARACTER*1. // On entry, UPLO specifies whether the upper or lower // triangular part of the array C is to be referenced as // follows: // // UPLO = 'U' or 'u' Only the upper triangular part of C // is to be referenced. // // UPLO = 'L' or 'l' Only the lower triangular part of C // is to be referenced. // // Unchanged on exit. // // TRANS - CHARACTER*1. // On entry, TRANS specifies the operation to be performed as // follows: // // TRANS = 'N' or 'n' C := alpha*A*A' + beta*C. // // TRANS = 'T' or 't' C := alpha*A'*A + beta*C. // // TRANS = 'C' or 'c' C := alpha*A'*A + beta*C. // // Unchanged on exit. // // N - INTEGER. // On entry, N specifies the order of the matrix C. N must be // at least zero. // Unchanged on exit. // // K - INTEGER. // On entry with TRANS = 'N' or 'n', K specifies the number // of columns of the matrix A, and on entry with // TRANS = 'T' or 't' or 'C' or 'c', K specifies the number // of rows of the matrix A. K must be at least zero. // Unchanged on exit. // // ALPHA - DOUBLE PRECISION. // On entry, ALPHA specifies the scalar alpha. // Unchanged on exit. // // A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is // k when TRANS = 'N' or 'n', and is n otherwise. // Before entry with TRANS = 'N' or 'n', the leading n by k // part of the array A must contain the matrix A, otherwise // the leading k by n part of the array A must contain the // matrix A. // Unchanged on exit. // // LDA - INTEGER. // On entry, LDA specifies the first dimension of A as declared // in the calling (sub) program. When TRANS = 'N' or 'n' // then LDA must be at least max( 1, n ), otherwise LDA must // be at least max( 1, k ). // Unchanged on exit. // // BETA - DOUBLE PRECISION. // On entry, BETA specifies the scalar beta. // Unchanged on exit. // // C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). // Before entry with UPLO = 'U' or 'u', the leading n by n // upper triangular part of the array C must contain the upper // triangular part of the symmetric matrix and the strictly // lower triangular part of C is not referenced. On exit, the // upper triangular part of the array C is overwritten by the // upper triangular part of the updated matrix. // Before entry with UPLO = 'L' or 'l', the leading n by n // lower triangular part of the array C must contain the lower // triangular part of the symmetric matrix and the strictly // upper triangular part of C is not referenced. On exit, the // lower triangular part of the array C is overwritten by the // lower triangular part of the updated matrix. // // LDC - INTEGER. // On entry, LDC specifies the first dimension of C as declared // in the calling (sub) program. LDC must be at least // max( 1, n ). // Unchanged on exit. // ------------------------------------------------------------ void dgemm(char TRANSA, char TRANSB, int M, int N, int K, double ALPHA, double *A, int LDA, double *B, int LDB, double BETA, double *C, int LDC); // ------------------------------------------------------------ // SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC ) // // DGEMM performs one of the matrix-matrix operations // // C := alpha*op( A )*op( B ) + beta*C, // // where op( X ) is one of // // op( X ) = X or op( X ) = X', // // alpha and beta are scalars, and A, B and C are matrices, with op( A ) // an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. // // Parameters // ========== // // TRANSA - CHARACTER*1. // On entry, TRANSA specifies the form of op( A ) to be used in // the matrix multiplication as follows: // // TRANSA = 'N' or 'n', op( A ) = A. // // TRANSA = 'T' or 't', op( A ) = A'. // // TRANSA = 'C' or 'c', op( A ) = A'. // // Unchanged on exit. // // TRANSB - CHARACTER*1. // On entry, TRANSB specifies the form of op( B ) to be used in // the matrix multiplication as follows: // // TRANSB = 'N' or 'n', op( B ) = B. // // TRANSB = 'T' or 't', op( B ) = B'. // // TRANSB = 'C' or 'c', op( B ) = B'. // // Unchanged on exit. // // M - INTEGER. // On entry, M specifies the number of rows of the matrix // op( A ) and of the matrix C. M must be at least zero. // Unchanged on exit. // // N - INTEGER. // On entry, N specifies the number of columns of the matrix // op( B ) and the number of columns of the matrix C. N must be // at least zero. // Unchanged on exit. // // K - INTEGER. // On entry, K specifies the number of columns of the matrix // op( A ) and the number of rows of the matrix op( B ). K must // be at least zero. // Unchanged on exit. // // ALPHA - DOUBLE PRECISION. // On entry, ALPHA specifies the scalar alpha. // Unchanged on exit. // // A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is // k when TRANSA = 'N' or 'n', and is m otherwise. // Before entry with TRANSA = 'N' or 'n', the leading m by k // part of the array A must contain the matrix A, otherwise // the leading k by m part of the array A must contain the // matrix A. // Unchanged on exit. // // LDA - INTEGER. // On entry, LDA specifies the first dimension of A as declared // in the calling (sub) program. When TRANSA = 'N' or 'n' then // LDA must be at least max( 1, m ), otherwise LDA must be at // least max( 1, k ). // Unchanged on exit. // // B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is // n when TRANSB = 'N' or 'n', and is k otherwise. // Before entry with TRANSB = 'N' or 'n', the leading k by n // part of the array B must contain the matrix B, otherwise // the leading n by k part of the array B must contain the // matrix B. // Unchanged on exit. // // LDB - INTEGER. // On entry, LDB specifies the first dimension of B as declared // in the calling (sub) program. When TRANSB = 'N' or 'n' then // LDB must be at least max( 1, k ), otherwise LDB must be at // least max( 1, n ). // Unchanged on exit. // // BETA - DOUBLE PRECISION. // On entry, BETA specifies the scalar beta. When BETA is // supplied as zero then C need not be set on input. // Unchanged on exit. // // C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). // Before entry, the leading m by n part of the array C must // contain the matrix C, except when beta is zero, in which // case C need not be set on entry. // On exit, the array C is overwritten by the m by n matrix // ( alpha*op( A )*op( B ) + beta*C ). // // LDC - INTEGER. // On entry, LDC specifies the first dimension of C as declared // in the calling (sub) program. LDC must be at least // max( 1, m ). // Unchanged on exit. // ------------------------------------------------------------ void dtrmm(char SIDE, char UPLO, char TRANSA, char DIAG, int M, int N, double ALPHA, double *A, int LDA, double *B, int LDB); // ------------------------------------------------------------ // SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, // B, LDB ) // .. Scalar Arguments .. // CHARACTER*1 SIDE, UPLO, TRANSA, DIAG // INTEGER M, N, LDA, LDB // DOUBLE PRECISION ALPHA // .. Array Arguments .. // DOUBLE PRECISION A( LDA, * ), B( LDB, * ) // .. // // Purpose // ======= // // DTRMM performs one of the matrix-matrix operations // // B := alpha*op( A )*B, or B := alpha*B*op( A ), // // where alpha is a scalar, B is an m by n matrix, A is a unit, or // non-unit, upper or lower triangular matrix and op( A ) is one of // // op( A ) = A or op( A ) = A'. // // Parameters // ========== // // SIDE - CHARACTER*1. // On entry, SIDE specifies whether op( A ) multiplies B from // the left or right as follows: // // SIDE = 'L' or 'l' B := alpha*op( A )*B. // // SIDE = 'R' or 'r' B := alpha*B*op( A ). // // Unchanged on exit. // // UPLO - CHARACTER*1. // On entry, UPLO specifies whether the matrix A is an upper or // lower triangular matrix as follows: // // UPLO = 'U' or 'u' A is an upper triangular matrix. // // UPLO = 'L' or 'l' A is a lower triangular matrix. // // Unchanged on exit. // // TRANSA - CHARACTER*1. // On entry, TRANSA specifies the form of op( A ) to be used in // the matrix multiplication as follows: // // TRANSA = 'N' or 'n' op( A ) = A. // // TRANSA = 'T' or 't' op( A ) = A'. // // TRANSA = 'C' or 'c' op( A ) = A'. // // Unchanged on exit. // // DIAG - CHARACTER*1. // On entry, DIAG specifies whether or not A is unit triangular // as follows: // // DIAG = 'U' or 'u' A is assumed to be unit triangular. // // DIAG = 'N' or 'n' A is not assumed to be unit // triangular. // // Unchanged on exit. // // M - INTEGER. // On entry, M specifies the number of rows of B. M must be at // least zero. // Unchanged on exit. // // N - INTEGER. // On entry, N specifies the number of columns of B. N must be // at least zero. // Unchanged on exit. // // ALPHA - DOUBLE PRECISION. // On entry, ALPHA specifies the scalar alpha. When alpha is // zero then A is not referenced and B need not be set before // entry. // Unchanged on exit. // // A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m // when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. // Before entry with UPLO = 'U' or 'u', the leading k by k // upper triangular part of the array A must contain the upper // triangular matrix and the strictly lower triangular part of // A is not referenced. // Before entry with UPLO = 'L' or 'l', the leading k by k // lower triangular part of the array A must contain the lower // triangular matrix and the strictly upper triangular part of // A is not referenced. // Note that when DIAG = 'U' or 'u', the diagonal elements of // A are not referenced either, but are assumed to be unity. // Unchanged on exit. // // LDA - INTEGER. // On entry, LDA specifies the first dimension of A as declared // in the calling (sub) program. When SIDE = 'L' or 'l' then // LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' // then LDA must be at least max( 1, n ). // Unchanged on exit. // // B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). // Before entry, the leading m by n part of the array B must // contain the matrix B, and on exit is overwritten by the // transformed matrix. // // LDB - INTEGER. // On entry, LDB specifies the first dimension of B as declared // in the calling (sub) program. LDB must be at least // max( 1, m ). // Unchanged on exit. // ------------------------------------------------------------ void dtrsm(char SIDE, char UPLO, char TRANSA, char DIAG, int M, int N, double ALPHA, double *A, int LDA, double *B, int LDB); // ------------------------------------------------------------ // SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, // $ B, LDB ) // .. Scalar Arguments .. // CHARACTER*1 SIDE, UPLO, TRANSA, DIAG // INTEGER M, N, LDA, LDB // DOUBLE PRECISION ALPHA // .. Array Arguments .. // DOUBLE PRECISION A( LDA, * ), B( LDB, * ) // .. // // Purpose // ======= // // DTRSM solves one of the matrix equations // // op( A )*X = alpha*B, or X*op( A ) = alpha*B, // // where alpha is a scalar, X and B are m by n matrices, A is a unit, or // non-unit, upper or lower triangular matrix and op( A ) is one of // // op( A ) = A or op( A ) = A'. // // The matrix X is overwritten on B. // // Parameters // ========== // // SIDE - CHARACTER*1. // On entry, SIDE specifies whether op( A ) appears on the left // or right of X as follows: // // SIDE = 'L' or 'l' op( A )*X = alpha*B. // // SIDE = 'R' or 'r' X*op( A ) = alpha*B. // // Unchanged on exit. // // UPLO - CHARACTER*1. // On entry, UPLO specifies whether the matrix A is an upper or // lower triangular matrix as follows: // // UPLO = 'U' or 'u' A is an upper triangular matrix. // // UPLO = 'L' or 'l' A is a lower triangular matrix. // // Unchanged on exit. // // TRANSA - CHARACTER*1. // On entry, TRANSA specifies the form of op( A ) to be used in // the matrix multiplication as follows: // // TRANSA = 'N' or 'n' op( A ) = A. // // TRANSA = 'T' or 't' op( A ) = A'. // // TRANSA = 'C' or 'c' op( A ) = A'. // // Unchanged on exit. // // DIAG - CHARACTER*1. // On entry, DIAG specifies whether or not A is unit triangular // as follows: // // DIAG = 'U' or 'u' A is assumed to be unit triangular. // // DIAG = 'N' or 'n' A is not assumed to be unit // triangular. // // Unchanged on exit. // // M - INTEGER. // On entry, M specifies the number of rows of B. M must be at // least zero. // Unchanged on exit. // // N - INTEGER. // On entry, N specifies the number of columns of B. N must be // at least zero. // Unchanged on exit. // // ALPHA - DOUBLE PRECISION. // On entry, ALPHA specifies the scalar alpha. When alpha is // zero then A is not referenced and B need not be set before // entry. // Unchanged on exit. // // A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m // when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. // Before entry with UPLO = 'U' or 'u', the leading k by k // upper triangular part of the array A must contain the upper // triangular matrix and the strictly lower triangular part of // A is not referenced. // Before entry with UPLO = 'L' or 'l', the leading k by k // lower triangular part of the array A must contain the lower // triangular matrix and the strictly upper triangular part of // A is not referenced. // Note that when DIAG = 'U' or 'u', the diagonal elements of // A are not referenced either, but are assumed to be unity. // Unchanged on exit. // // LDA - INTEGER. // On entry, LDA specifies the first dimension of A as declared // in the calling (sub) program. When SIDE = 'L' or 'l' then // LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' // then LDA must be at least max( 1, n ). // Unchanged on exit. // // B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). // Before entry, the leading m by n part of the array B must // contain the right-hand side matrix B, and on exit is // overwritten by the solution matrix X. // // LDB - INTEGER. // On entry, LDB specifies the first dimension of B as declared // in the calling (sub) program. LDB must be at least // max( 1, m ). // Unchanged on exit. // ------------------------------------------------------------ int dgetri(int N, double *A, int LDA, int *IPIV, double *WORK, int LWORK); // ------------------------------------------------------------ // SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, 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 // // .. Scalar Arguments .. // INTEGER INFO, LDA, LWORK, N // .. // .. Array Arguments .. // INTEGER IPIV( * ) // DOUBLE PRECISION A( LDA, * ), WORK( * ) // .. // // Purpose // ======= // // DGETRI computes the inverse of a matrix using the LU factorization // computed by DGETRF. // // This method inverts U and then computes inv(A) by solving the system // inv(A)*L = inv(U) for inv(A). // // Arguments // ========= // // N (input) INTEGER // The order of the matrix A. N >= 0. // // A (input/output) DOUBLE PRECISION array, dimension (LDA,N) // On entry, the factors L and U from the factorization // A = P*L*U as computed by DGETRF. // On exit, if INFO = 0, the inverse of the original matrix A. // // LDA (input) INTEGER // The leading dimension of the array A. LDA >= max(1,N). // // IPIV (input) INTEGER array, dimension (N) // The pivot indices from DGETRF; for 1<=i<=N, row i of the // matrix was interchanged with row IPIV(i). // // WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) // On exit, if INFO=0, then WORK(1) returns the optimal LWORK. // // LWORK (input) INTEGER // The dimension of the array WORK. LWORK >= max(1,N). // For optimal performance LWORK >= N*NB, where NB is // the optimal blocksize returned by ILAENV. // // If LWORK = -1, then a workspace query is assumed; the routine // only calculates the optimal size of the WORK array, returns // this value as the first entry of the WORK array, and no error // message related to LWORK is issued by XERBLA. // // INFO (output) INTEGER // = 0: successful exit // < 0: if INFO = -i, the i-th argument had an illegal value // > 0: if INFO = i, U(i,i) is exactly zero; the matrix is // singular and its inverse could not be computed. // ------------------------------------------------------------ int dpotri(char UPLO, int N, double *A, int LDA); // ------------------------------------------------------------ // SUBROUTINE DPOTRI( UPLO, N, A, LDA, INFO ) // // -- LAPACK routine (version 3.0) -- // Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., // Courant Institute, Argonne National Lab, and Rice University // March 31, 1993 // // .. Scalar Arguments .. // CHARACTER UPLO // INTEGER INFO, LDA, N // .. // .. Array Arguments .. // DOUBLE PRECISION A( LDA, * ) // .. // // Purpose // ======= // // DPOTRI computes the inverse of a real symmetric positive definite // matrix A using the Cholesky factorization A = U**T*U or A = L*L**T // computed by DPOTRF. // // Arguments // ========= // // UPLO (input) CHARACTER*1 // = 'U': Upper triangle of A is stored; // = 'L': Lower triangle of A is stored. // // N (input) INTEGER // The order of the matrix A. N >= 0. // // A (input/output) DOUBLE PRECISION array, dimension (LDA,N) // On entry, the triangular factor U or L from the Cholesky // factorization A = U**T*U or A = L*L**T, as computed by // DPOTRF. // On exit, the upper or lower triangle of the (symmetric) // inverse of A, overwriting the input factor U or L. // // LDA (input) INTEGER // The leading dimension of the array A. LDA >= max(1,N). // // INFO (output) INTEGER // = 0: successful exit // < 0: if INFO = -i, the i-th argument had an illegal value // > 0: if INFO = i, the (i,i) element of the factor U or L is // zero, and the inverse could not be computed. #endif // LAPACK_CWRAP_H