Loading...
Searching...
No Matches
qrbdv.c
Go to the documentation of this file.
/* qrbdv.c CCMATH mathematics library source code.
*
* Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
* This code may be redistributed under the terms of the GNU library
* public license (LGPL). ( See the lgpl.license file for details.)
* ------------------------------------------------------------------------
*/
#include "ccmath.h"
int qrbdv(double *dm, double *em, double *um, int mm, double *vm, int m)
{
int i, j, k, n, jj, nm;
double u, x, y, a, b, c, s, t, w, *p, *q;
for (j = 1, t = fabs(dm[0]); j < m; ++j)
if ((s = fabs(dm[j]) + fabs(em[j - 1])) > t)
t = s;
t *= 1.e-15;
n = 100 * m;
nm = m;
for (j = 0; m > 1 && j < n; ++j) {
for (k = m - 1; k > 0; --k) {
if (fabs(em[k - 1]) < t)
break;
if (fabs(dm[k - 1]) < t) {
for (i = k, s = 1., c = 0.; i < m; ++i) {
a = s * em[i - 1];
b = dm[i];
em[i - 1] *= c;
dm[i] = u = sqrt(a * a + b * b);
s = -a / u;
c = b / u;
for (jj = 0, p = um + k - 1; jj < mm; ++jj, p += mm) {
q = p + i - k + 1;
w = c * *p + s * *q;
*q = c * *q - s * *p;
*p = w;
}
}
break;
}
}
y = dm[k];
x = dm[m - 1];
u = em[m - 2];
a = (y + x) * (y - x) - u * u;
s = y * em[k];
b = s + s;
u = sqrt(a * a + b * b);
if (u != 0.) {
c = sqrt((u + a) / (u + u));
if (c != 0.)
s /= (c * u);
else
s = 1.;
for (i = k; i < m - 1; ++i) {
b = em[i];
if (i > k) {
a = s * em[i];
b *= c;
em[i - 1] = u = sqrt(x * x + a * a);
c = x / u;
s = a / u;
}
a = c * y + s * b;
b = c * b - s * y;
for (jj = 0, p = vm + i; jj < nm; ++jj, p += nm) {
w = c * *p + s * *(p + 1);
*(p + 1) = c * *(p + 1) - s * *p;
*p = w;
}
s *= dm[i + 1];
dm[i] = u = sqrt(a * a + s * s);
y = c * dm[i + 1];
c = a / u;
s /= u;
x = c * b + s * y;
y = c * y - s * b;
for (jj = 0, p = um + i; jj < mm; ++jj, p += mm) {
w = c * *p + s * *(p + 1);
*(p + 1) = c * *(p + 1) - s * *p;
*p = w;
}
}
}
em[m - 2] = x;
dm[m - 1] = y;
if (fabs(x) < t)
--m;
if (m == k + 1)
--m;
}
return j;
}
ccmath.h
b
double b
Definition driver/set_window.c:5
t
double t
Definition driver/set_window.c:5
qrbdv
int qrbdv(double *dm, double *em, double *um, int mm, double *vm, int m)
Definition qrbdv.c:9
x
#define x