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
|
/**
* @file msvmmaj_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 "msvmmaj_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 Cholesky factorization of a real symmetric positive
* definite matrix.
*
* @details
* This function uses the external LAPACK routine dpotrf.
*
* @param[in] UPLO which triangle of A is stored
* @param[in] N order of A
* @param[in,out] A double precision array of size (LDA, N). On
* exit contains the factor U or L of the Cholesky
* factorization
* @param[in] LDA leading dimension of A
* @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
* order i is not positive
* definite
*
* See the LAPACK documentation at:
* http://www.netlib.org/lapack/explore-html/d0/d8a/dpotrf_8f.html
*/
int dpotrf(char UPLO, int N, double *A, int LDA)
{
extern void dpotrf_(char *UPLO, int *N, double *A, int *LDA, int *INFOp);
int INFO;
dpotrf_(&UPLO, &N, A, &LDA, &INFO);
return INFO;
}
|
