Functions for QR decomposition and matrix inversion.
More...
Go to the source code of this file.
|
void | RSolve (double **a, int n, double *d, double *b) |
|
void | QRDecompose (double **a, int n, double *c, double *d, int *sing) |
|
void | QRSolve (double **a, int n, double *c, double *d, double *b) |
|
void | QRUpdate (double **r, double **qt, int n, double *u, double *v) |
|
void | rotate (double **r, double **qt, int n, int i, double a, double b) |
|
Functions for QR decomposition and matrix inversion.
- Author
- B. Vaidya (bvaid.nosp@m.ya@u.nosp@m.nito..nosp@m.it)
- Date
- September 11, 2014
Definition in file math_qr_decomp.c.
void QRDecompose |
( |
double ** |
a, |
|
|
int |
n, |
|
|
double * |
c, |
|
|
double * |
d, |
|
|
int * |
sing |
|
) |
| |
Definition at line 34 of file math_qr_decomp.c.
46 double scale, sigma, sum, tau, temp;
49 for (k = 0; k <
n-1; k++) {
51 for (i = k; i <
n; i++){
52 if ((temp = fabs (
a[i][k])) > scale)
60 for (i = k; i <
n; i++){
64 for (i = k; i <
n; i++) {
69 sigma = fabs(sqrt(sum));
71 sigma = -fabs(sqrt(sum));
79 for (i = k; i <
n; i++) {
80 sum += a[
i][
k] * a[
i][
j];
83 for (i = k; i <
n; i++) {
84 a[
i][
j] -= tau * a[
i][
k];
90 if (d[n-1] == 0.0) *sing = 1;
void QRSolve |
( |
double ** |
a, |
|
|
int |
n, |
|
|
double * |
c, |
|
|
double * |
d, |
|
|
double * |
b |
|
) |
| |
Definition at line 97 of file math_qr_decomp.c.
void RSolve(double **a, int n, double *d, double *b)
void QRUpdate |
( |
double ** |
r, |
|
|
double ** |
qt, |
|
|
int |
n, |
|
|
double * |
u, |
|
|
double * |
v |
|
) |
| |
Definition at line 121 of file math_qr_decomp.c.
130 for (k=
n-1;k>=0;k--){
134 for (i = k-1; i>=0; i--){
135 rotate(r,qt,
n, i, u[i], -u[i+1]);
136 if (u[i] == 0.0) u[
i] = fabs(u[i+1]);
137 else if (fabs(u[i]) > fabs(u[i+1]))
138 u[
i] = fabs(u[i])*sqrt(1.0 + (u[i+1]/u[i])*(u[i+1]/u[i]));
140 u[
i] = fabs(u[i+1])*sqrt(1.0 + (u[i]/u[i+1])*(u[i]/u[i+1]));
143 for (j=0;j<
n; j++) r[0][j] += u[0]*
v[j];
145 rotate(r,qt,n,i,r[i][i], -r[i+1][i]);
void rotate(double **r, double **qt, int n, int i, double a, double b)
void rotate |
( |
double ** |
r, |
|
|
double ** |
qt, |
|
|
int |
n, |
|
|
int |
i, |
|
|
double |
a, |
|
|
double |
b |
|
) |
| |
Definition at line 149 of file math_qr_decomp.c.
155 double c,
s, w, y, fact;
159 s = (b >= 0.0 ? 1.0 : -1.0);
160 }
else if (fabs(
a) > fabs(b)){
163 c = fabs(1.0/sqrt(1.0 + (fact*fact)));
165 c = -fabs(1.0/sqrt(1.0 + (fact*fact)));
171 s = fabs(1.0/sqrt(1.0 + (fact*fact)));
173 s = -fabs(1.0/sqrt(1.0 + (fact*fact)));
182 r[
i+1][
j] = s*y + c*w;
188 qt[
i][
j] = c*y - s*w;
189 qt[
i+1][
j] = s*y + c*w;
void RSolve |
( |
double ** |
a, |
|
|
int |
n, |
|
|
double * |
d, |
|
|
double * |
b |
|
) |
| |
Definition at line 12 of file math_qr_decomp.c.
26 for (j = i+1; j <
n; j++){
29 b[
i] = (b[
i] - sum)/d[i];