jens/to_be_added
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
Files
master
to_be_added/lapack_cwrap.h
T
jens 05da8712f6 Initial import 
git-svn-id: http://moon:8086/svn/software/trunk/libsrc/_to_be_added@1 b431acfa-c32f-4a4a-93f1-934dc6c82436
2014-07-19 07:44:42 +00:00

1390 lines
55 KiB
C
Executable File
Raw Permalink Blame History

// ------------------------------------------------------------
// 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
Reference in New Issue View Git Blame Copy Permalink