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
|
#ifndef LUDECOMP_H
#define LUDECOMP_H
/*
* Copyright 2000 Graeme W. Gill
* All rights reserved.
*
* This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
* see the License.txt file for licencing details.
*/
#ifdef __cplusplus
extern "C" {
#endif
/* NOTE:- lu decomp rearanges the rows of the matrix */
/* by swapping pointers rather than exchanging data, */
/* so the matrix must be addressed by the **pointer */
/* if it is re-used after an ludecomp!!! */
/* Solve the simultaneous linear equations A.X = B */
/* Return 1 if the matrix is singular, 0 if OK */
int
solve_se(
double **a, /* A[][] input matrix, returns LU decimposition of A */
double *b, /* B[] input array, returns solution X[] */
int n /* Dimensionality */
);
/* Solve the simultaneous linear equations A.X = B, with polishing */
/* Return 1 if the matrix is singular, 0 if OK */
int
polished_solve_se(
double **a, /* A[][] input matrix, returns LU decimposition of A */
double *b, /* B[] input array, returns solution X[] */
int n /* Dimensionality */
);
/* Decompose the square matrix A[][] into lower and upper triangles */
/* NOTE that rows get swaped by swapping matrix pointers! */
/* Return 1 if the matrix is singular. */
int
lu_decomp(
double **a, /* A input array, output upper and lower triangles. */
int n, /* Dimensionality */
int *pivx, /* Return pivoting row permutations record */
double *rip /* Row interchange parity, +/- 1.0, used for determinant */
);
/* Solve a set of simultaneous equations from the */
/* LU decomposition, by back substitution. */
void
lu_backsub(
double **a, /* A[][] LU decomposed matrix */
int n, /* Dimensionality */
int *pivx, /* Pivoting row permutations record */
double *b /* Input B[] vecor, return X[] */
);
/* Polish a solution for equations */
void
lu_polish(
double **a, /* Original A[][] matrix */
double **lua, /* LU decomposition of A[][] */
int n, /* Dimensionality */
double *b, /* B[] vector of equation */
double *x, /* X[] solution to be polished */
int *pivx /* Pivoting row permutations record */
);
/* Invert a matrix A using lu decomposition */
/* NOTE that it returns transposed inverse by normal convention. */
/* Use sym_matrix_trans() to fix this, or use matrix_trans_mult() */
/* Return 1 if the matrix is singular, 0 if OK */
int
lu_invert(
double **a, /* A[][] input matrix, returns inversion of A transposed */
int n /* Dimensionality */
);
/* Invert a matrix A using lu decomposition */
/* The normal convention inverse is returned */
/* Return 1 if the matrix is singular, 0 if OK */
int
lu_invert_normal(
double **a, /* A[][] input matrix, returns inversion of A */
int n /* Dimensionality */
);
/* Invert a matrix A using lu decomposition, and polish it. */
/* NOTE that it returns transposed inverse by normal convention. */
/* Use sym_matrix_trans() to fix this, or use matrix_trans_mult() */
/* Return 1 if the matrix is singular, 0 if OK */
int
lu_polished_invert(
double **a, /* A[][] input matrix, returns inversion of A */
int n /* Dimensionality */
);
/* Pseudo-Invert matrix A using lu decomposition */
/* Return 1 if the matrix is singular, 0 if OK */
int
lu_psinvert(
double **out, /* Output[0..N-1][0..M-1] */
double **in, /* Input[0..M-1][0..N-1] input matrix */
int m, /* In Rows */
int n /* In Columns */
);
/* - - - - - - - - - - - - - - - - - - - */
/* Use Cholesky decomposition on a symetric positive-definite matrix. */
/* Only the upper triangle of the matrix A is accessed. */
/* L returns the decomposition */
/* Return nz if A is not positive-definite */
int llt_decomp(double **L, double **A, int n);
/* Solve a set of simultaneous equations A.x = b from the */
/* LLt decomposition, by back substitution. */
void llt_backsub(
double **L, /* A[][] LLt decomposition in lower triangle */
int n, /* Dimensionality */
double *b, /* Input B[] */
double *x /* Return X[] (may be same as B[]) */
);
#ifdef __cplusplus
}
#endif
#endif /* LUDECOMP_H */
|
