3 #include <grass/gmath.h>
5 static double at, bt, ct;
7 #define PYTHAG(a,b) ((at=fabs(a)) > (bt=fabs(b)) ? \
8 (ct=bt/at,at*sqrt(1.0+ct*ct)) : (bt ? (ct=at/bt,bt*sqrt(1.0+ct*ct)): 0.0))
10 static double maxarg1, maxarg2;
12 #define MAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
13 (maxarg1) : (maxarg2))
14 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
16 int G_svdcmp(
double **a,
int m,
int n,
double *
w,
double **v)
18 int flag, i, its, j, jj, k, ii = 0, nm = 0;
19 double c, f,
h,
s, x,
y, z;
20 double anorm = 0.0,
g = 0.0, scale = 0.0;
31 for (i = 0; i <=
n; i++) {
36 for (k = i; k <= m; k++)
37 scale += fabs(a[k][i]);
39 for (k = i; k <= m; k++) {
41 s += a[k][i] * a[k][i];
44 g = -
SIGN(sqrt(s), f);
48 for (j = ii; j <=
n; j++) {
49 for (s = 0.0, k = i; k <= m; k++)
50 s += a[k][i] * a[k][j];
52 for (k = i; k <= m; k++)
53 a[k][j] += f * a[k][i];
56 for (k = i; k <= m; k++)
62 if (i <= m && i != n) {
63 for (k = ii; k <=
n; k++)
64 scale += fabs(a[i][k]);
66 for (k = ii; k <=
n; k++) {
68 s += a[i][k] * a[i][k];
71 g = -
SIGN(sqrt(s), f);
74 for (k = ii; k <=
n; k++)
77 for (j = ii; j <= m; j++) {
78 for (s = 0.0, k = ii; k <=
n; k++)
79 s += a[j][k] * a[i][k];
80 for (k = ii; k <=
n; k++)
81 a[j][k] += s * rv1[k];
84 for (k = ii; k <=
n; k++)
88 anorm =
MAX(anorm, (fabs(w[i]) + fabs(rv1[i])));
90 for (i = n; i >= 0; i--) {
93 for (j = ii; j <=
n; j++)
94 v[j][i] = (a[i][j] / a[i][ii]) /
g;
95 for (j = ii; j <=
n; j++) {
96 for (s = 0.0, k = ii; k <=
n; k++)
97 s += a[i][k] * v[k][j];
98 for (k = ii; k <=
n; k++)
99 v[k][j] += s * v[k][i];
102 for (j = ii; j <=
n; j++)
103 v[i][j] = v[j][i] = 0.0;
109 for (i = n; i >= 0; i--) {
113 for (j = ii; j <=
n; j++)
118 for (j = ii; j <=
n; j++) {
119 for (s = 0.0, k = ii; k <= m; k++)
120 s += a[k][i] * a[k][j];
121 f = (s / a[i][i]) *
g;
122 for (k = i; k <= m; k++)
123 a[k][j] += f * a[k][i];
126 for (j = i; j <= m; j++)
130 for (j = i; j <= m; j++)
135 for (k = n; k >= 0; k--) {
136 for (its = 1; its <= 30; its++) {
138 for (ii = k; ii >= 0; ii--) {
140 if (fabs(rv1[ii]) + anorm == anorm) {
144 if (fabs(w[nm]) + anorm == anorm)
150 for (i = ii; i <= k; i++) {
152 if (fabs(f) + anorm != anorm) {
159 for (j = 0; j <= m; j++) {
162 a[j][nm] = y * c + z *
s;
163 a[j][i] = z * c - y *
s;
172 for (j = 0; j <=
n; j++)
173 v[j][k] = (-v[j][k]);
184 f = ((y - z) * (y + z) + (
g -
h) * (
g + h)) / (2.0 * h * y);
186 f = ((x - z) * (x + z) + h * ((y / (f +
SIGN(
g, f))) -
h)) / x;
188 for (j = ii; j <= nm; j++) {
202 for (jj = 0; jj <=
n; jj++) {
205 v[jj][j] = x * c + z *
s;
206 v[jj][i] = z * c - x *
s;
215 f = (c *
g) + (s * y);
216 x = (c *
y) - (s * g);
217 for (jj = 0; jj <= m; jj++) {
220 a[jj][j] = y * c + z *
s;
221 a[jj][i] = z * c - y *
s;
238 int m,
int n,
double b[],
double x[])
244 for (j = 0; j <
n; j++) {
247 for (i = 0; i < m; i++)
253 for (j = 0; j <
n; j++) {
255 for (i = 0; i <
n; i++)
256 s += v[j][i] * tmp[i];
272 for (i = 0; i <
n; i++)
276 for (i = 0; i <
n; i++)
int G_svbksb(double **u, double w[], double **v, int m, int n, double b[], double x[])
int G_svdcmp(double **a, int m, int n, double *w, double **v)
void G_free_vector(double *v)
Vector memory deallocation.
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
int G_svelim(double *w, int n)
tuple h
panel.defaultSize = wx.CheckBox(panel, id = wx.ID_ANY, label = _("Use default size")) panel...