1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
|
/**
* @file gensvm_lapack.c
* @author Gertjan van den Burg
* @date August 9, 2013
* @brief Utility functions for interacting with LAPACK
*
* @details
* Functions in this file are auxiliary functions which make it easier
* to use LAPACK functions from liblapack.
*/
#include "gensvm_lapack.h"
/**
* @brief Solve AX = B where A is symmetric positive definite.
*
* @details
* Solve a linear system of equations AX = B where A is symmetric positive
* definite. This function uses the externel LAPACK routine dposv.
*
* @param[in] UPLO which triangle of A is stored
* @param[in] N order of A
* @param[in] NRHS number of columns of B
* @param[in,out] A double precision array of size (LDA, N). On
* exit contains the upper or lower factor of the
* Cholesky factorization of A.
* @param[in] LDA leading dimension of A
* @param[in,out] B double precision array of size (LDB, NRHS). On
* exit contains the N-by-NRHS solution matrix X.
* @param[in] LDB the leading dimension of B
* @returns info parameter which contains the status of the
* computation:
* - =0: success
* - <0: if -i, the i-th argument had
* an illegal value
* - >0: if i, the leading minor of A
* was not positive definite
*
* See the LAPACK documentation at:
* http://www.netlib.org/lapack/explore-html/dc/de9/group__double_p_osolve.html
*/
int dposv(char UPLO, int N, int NRHS, double *A, int LDA, double *B,
int LDB)
{
extern void dposv_(char *UPLO, int *Np, int *NRHSp, double *A,
int *LDAp, double *B, int *LDBp, int *INFOp);
int INFO;
dposv_(&UPLO, &N, &NRHS, A, &LDA, B, &LDB, &INFO);
return INFO;
}
/**
* @brief Solve a system of equations AX = B where A is symmetric.
*
* @details
* Solve a linear system of equations AX = B where A is symmetric. This
* function uses the external LAPACK routine dsysv.
*
* @param[in] UPLO which triangle of A is stored
* @param[in] N order of A
* @param[in] NRHS number of columns of B
* @param[in,out] A double precision array of size (LDA, N). On
* exit contains the block diagonal matrix D and
* the multipliers used to obtain the factor U or
* L from the factorization A = U*D*U**T or
* A = L*D*L**T.
* @param[in] LDA leading dimension of A
* @param[in] IPIV integer array containing the details of D
* @param[in,out] B double precision array of size (LDB, NRHS). On
* exit contains the N-by-NRHS matrix X
* @param[in] LDB leading dimension of B
* @param[out] WORK double precision array of size max(1,LWORK). On
* exit, WORK(1) contains the optimal LWORK
* @param[in] LWORK the length of WORK, can be used for determining
* the optimal blocksize for dsystrf.
* @returns info parameter which contains the status of the
* computation:
* - =0: success
* - <0: if -i, the i-th argument had an
* illegal value
* - >0: if i, D(i, i) is exactly zero,
* no solution can be computed.
*
* See the LAPACK documentation at:
* http://www.netlib.org/lapack/explore-html/d6/d0e/group__double_s_ysolve.html
*/
int dsysv(char UPLO, int N, int NRHS, double *A, int LDA, int *IPIV,
double *B, int LDB, double *WORK, int LWORK)
{
extern void dsysv_(char *UPLO, int *Np, int *NRHSp, double *A,
int *LDAp, int *IPIV, double *B, int *LDBp,
double *WORK, int *LWORK, int *INFOp);
int INFO;
dsysv_(&UPLO, &N, &NRHS, A, &LDA, IPIV, B, &LDB, WORK, &LWORK, &INFO);
return INFO;
}
/**
* @brief Compute the eigenvalues and optionally the eigenvectors of a
* symmetric matrix.
*
* @details
* See the LAPACK documentation at:
* http://www.netlib.org/lapack/explore-html/d2/d97/dsyevx_8f.html
*
*
*/
int dsyevx(char JOBZ, char RANGE, char UPLO, int N, double *A, int LDA,
double VL, double VU, int IL, int IU, double ABSTOL, int *M,
double *W, double *Z, int LDZ, double *WORK, int LWORK,
int *IWORK, int *IFAIL)
{
extern void dsyevx_(char *JOBZ, char *RANGE, char *UPLO, int *Np,
double *A, int *LDAp, double *VLp, double *VUp,
int *ILp, int *IUp, double *ABSTOLp, int *M,
double *W, double *Z, int *LDZp, double *WORK,
int *LWORKp, int *IWORK, int *IFAIL, int *INFOp);
int INFO;
dsyevx_(&JOBZ, &RANGE, &UPLO, &N, A, &LDA, &VL, &VU, &IL, &IU, &ABSTOL,
M, W, Z, &LDZ, WORK, &LWORK, IWORK, IFAIL, &INFO);
return INFO;
}
/**
* @brief Determine double precision machine parameters.
*
* @details
* See the LAPACK documentation at:
* http://www.netlib.org/lapack/explore-html/d5/dd4/dlamch_8f.html
*/
double dlamch(char CMACH)
{
extern double dlamch_(char *CMACH);
return dlamch_(&CMACH);
}
|
