9 int qrbdu1(
double *dm,
double *em,
double *um,
int mm,
double *vm,
int m)
11 int i, j, k, n, jj, nm;
13 double u,
x, y, a,
b, c, s,
t, w, *p, *q;
15 for (j = 1,
t = fabs(dm[0]); j < m; ++j)
16 if ((s = fabs(dm[j]) + fabs(em[j - 1])) >
t)
21 for (j = 0; m > 1 && j < n; ++j) {
22 for (k = m - 1; k > 0; --k) {
23 if (fabs(em[k - 1]) <
t)
25 if (fabs(dm[k - 1]) <
t) {
26 for (i = k, s = 1., c = 0.; i < m; ++i) {
30 dm[i] = u = sqrt(a * a +
b *
b);
33 for (jj = 0, p = um + k - 1; jj < mm; ++jj, p += nm) {
46 a = (y +
x) * (y -
x) - u * u;
49 u = sqrt(a * a +
b *
b);
51 c = sqrt((u + a) / (u + u));
56 for (i = k; i < m - 1; ++i) {
61 em[i - 1] = u = sqrt(
x *
x + a * a);
67 for (jj = 0, p = vm + i; jj < nm; ++jj, p += nm) {
68 w = c * *p + s * *(p + 1);
69 *(p + 1) = c * *(p + 1) - s * *p;
73 dm[i] = u = sqrt(a * a + s * s);
79 for (jj = 0, p = um + i; jj < mm; ++jj, p += nm) {
80 w = c * *p + s * *(p + 1);
81 *(p + 1) = c * *(p + 1) - s * *p;
int qrbdu1(double *dm, double *em, double *um, int mm, double *vm, int m)