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
|
#ifndef SVD_H
#define SVD_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
/*
U[] decomposes A[]'s columns into orthogonal, singular vectors.
U[]'s columns are vectors that sum to 1.0, i.e. they leave a vectors normal unchanged.
The inverse of U[] is its transpose.
U[]'s columns corresponding to non-zero W[] are the orthonormal vectors that span
the (input) RANGE space. Columns corresponding to zero W[] are zero.
W[] will return singlular values, i.e. the weighting of the singular vectors.
It's inverse is is the reciprical of its elements.
V[] composes the singular vectors back into A[]'s rows.
V[]'s columns and rows are orthonormal.
V[]'s columns corresponding to non-zero W[] are the orthonormal vectors that span
the (output) RANGE space.
v[]'s columns corresponding to zero W[] are the (output) othonormal vectors that span
the NULL space.
The inverse of V[] is its transpose.
To re-form, A = U.W.Vt, i.e. multiply by transpose of V.
i.e. mult. input vector[m] by U[] converts to [n] compact, orthogonal
basis vectors. W then scales them appropiately, setting null space
vectors to 0. V[] then transforms from the orthogonal basis vectors
to A[]'s output space.
To reveal NULL space, make sure n >= m, since U[] vectors corrsponding
to zero's are set to zero.
*/
/* Compute Singular Value Decomposition of A = U.W.Vt */
/* Return status value: */
/* 0 - no error */
/* 1 - m < n error */
int svdecomp(
double **a, /* A[0..m-1][0..n-1], return U[0..m-1][0..n-1] */
double *w, /* return W[0..n-1] = singular values */
double **v, /* return V[0..n-1][0..n-1] (not transpose!) */
int m, /* Number of equations */
int n /* Number of unknowns */
);
/* Threshold the singular values W[] */
void svdthresh(
double w[], /* Singular values */
int n /* Number of unknowns */
);
/* Threshold the singular values W[] to give a specific dof */
void svdsetthresh(
double w[], /* Singular values */
int n, /* Number of unknowns */
int dof /* Expected degree of freedom */
);
/* Use output of svdcmp() to solve overspecified and/or */
/* singular equation A.x = b */
int svdbacksub(
double **u, /* U[0..m-1][0..n-1] U, W, V SVD decomposition of A[][] */
double *w, /* W[0..n-1] */
double **v, /* V[0..n-1][0..n-1] (not transpose!) */
double b[], /* B[0..m-1] Right hand side of equation */
double x[], /* X[0..n-1] Return solution. (May be the same as b[]) */
int m, /* Number of equations */
int n /* Number of unknowns */
);
/* Solve the equation A.x = b using SVD */
/* (The w[] values are thresholded for best accuracy) */
/* Return non-zero if no solution found */
/* !!! Note that A[][] will be changed !!! */
int svdsolve(
double **a, /* A[0..m-1][0..n-1] input A[][], will return U[][] */
double b[], /* B[0..m-1] Right hand side of equation, return solution */
int m, /* Number of equations */
int n /* Number of unknowns */
);
/* Solve the equation A.x = b using SVD */
/* The top s out of n singular values will be used */
/* Return non-zero if no solution found */
/* !!! Note that A[][] will be changed !!! */
int svdsolve_s(
double **a, /* A[0..m-1][0..n-1] input A[][], will return U[][] */
double b[], /* B[0..m-1] Right hand side of equation, return solution */
int m, /* Number of equations */
int n, /* Number of unknowns */
int s /* Number of unknowns */
);
/* Solve the equation A.x = b using Direct calculation, LU or SVD as appropriate */
/* Return non-zero if no solution found */
/* !!! Note that A[][] will be changed !!! */
int gen_solve_se(
double **a, /* A[0..m-1][0..n-1] input A[][], will return U[][] */
double b[], /* B[0..m-1] Right hand side of equation, return solution */
int m, /* Number of equations */
int n /* Number of unknowns */
);
#ifdef __cplusplus
}
#endif
#endif /* SVD_H */
|
