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
122
123
124
125
126
127
128
129
130
131
|
template <typename Scalar>
void ei_qrfac(int m, int n, Scalar *a, int
lda, int pivot, int *ipvt, Scalar *rdiag)
{
/* System generated locals */
int a_dim1, a_offset;
/* Local variables */
int i, j, k, jp1;
Scalar sum;
int kmax;
Scalar temp;
int minmn;
Scalar ajnorm;
Matrix< Scalar, Dynamic, 1 > wa(n+1);
/* Parameter adjustments */
--rdiag;
a_dim1 = lda;
a_offset = 1 + a_dim1 * 1;
a -= a_offset;
--ipvt;
/* Function Body */
const Scalar epsmch = epsilon<Scalar>();
/* compute the initial column norms and initialize several arrays. */
for (j = 1; j <= n; ++j) {
rdiag[j] = Map< Matrix< Scalar, Dynamic, 1 > >(&a[j * a_dim1 + 1],m).blueNorm();
wa[j] = rdiag[j];
if (pivot)
ipvt[j] = j;
}
/* reduce a to r with householder transformations. */
minmn = std::min(m,n);
for (j = 1; j <= minmn; ++j) {
if (! (pivot)) {
goto L40;
}
/* bring the column of largest norm into the pivot position. */
kmax = j;
for (k = j; k <= n; ++k) {
if (rdiag[k] > rdiag[kmax]) {
kmax = k;
}
/* L20: */
}
if (kmax == j) {
goto L40;
}
for (i = 1; i <= m; ++i) {
temp = a[i + j * a_dim1];
a[i + j * a_dim1] = a[i + kmax * a_dim1];
a[i + kmax * a_dim1] = temp;
/* L30: */
}
rdiag[kmax] = rdiag[j];
wa[kmax] = wa[j];
k = ipvt[j];
ipvt[j] = ipvt[kmax];
ipvt[kmax] = k;
L40:
/* compute the householder transformation to reduce the */
/* j-th column of a to a multiple of the j-th unit vector. */
ajnorm = Map< Matrix< Scalar, Dynamic, 1 > >(&a[j + j * a_dim1],m-j+1).blueNorm();
if (ajnorm == 0.) {
goto L100;
}
if (a[j + j * a_dim1] < 0.) {
ajnorm = -ajnorm;
}
for (i = j; i <= m; ++i) {
a[i + j * a_dim1] /= ajnorm;
/* L50: */
}
a[j + j * a_dim1] += 1.;
/* apply the transformation to the remaining columns */
/* and update the norms. */
jp1 = j + 1;
if (n < jp1) {
goto L100;
}
for (k = jp1; k <= n; ++k) {
sum = 0.;
for (i = j; i <= m; ++i) {
sum += a[i + j * a_dim1] * a[i + k * a_dim1];
/* L60: */
}
temp = sum / a[j + j * a_dim1];
for (i = j; i <= m; ++i) {
a[i + k * a_dim1] -= temp * a[i + j * a_dim1];
/* L70: */
}
if (! (pivot) || rdiag[k] == 0.) {
goto L80;
}
temp = a[j + k * a_dim1] / rdiag[k];
/* Computing MAX */
/* Computing 2nd power */
rdiag[k] *= ei_sqrt((std::max(Scalar(0.), Scalar(1.)-ei_abs2(temp))));
/* Computing 2nd power */
if (Scalar(.05) * ei_abs2(rdiag[k] / wa[k]) > epsmch) {
goto L80;
}
rdiag[k] = Map< Matrix< Scalar, Dynamic, 1 > >(&a[jp1 + k * a_dim1],m-j).blueNorm();
wa[k] = rdiag[k];
L80:
/* L90: */
;
}
L100:
rdiag[j] = -ajnorm;
/* L110: */
}
return;
/* last card of subroutine qrfac. */
} /* qrfac_ */
|
