38 VSinvUt_(A.m(), A.n()),
43 const label Um = U_.
m();
44 const label Un = U_.
n();
53 for (label i = 0; i < Um; i++)
61 for (label
k = i;
k < Un;
k++)
63 scale +=
mag(U_[
k][i]);
68 for (label
k = i;
k < Un;
k++)
71 s += U_[
k][i]*U_[
k][i];
79 for (label j = l-1; j < Um; j++)
82 for (label
k = i;
k < Un;
k++)
84 s += U_[
k][i]*U_[
k][j];
88 for (label
k = i;
k < A.
n();
k++)
90 U_[
k][j] += f*U_[
k][i];
94 for (label
k = i;
k < Un;
k++)
105 if (i+1 <= Un && i != Um)
107 for (label
k = l-1;
k < Um;
k++)
109 scale +=
mag(U_[i][
k]);
114 for (label
k=l-1;
k < Um;
k++)
117 s += U_[i][
k]*U_[i][
k];
120 scalar
f = U_[i][l-1];
125 for (label
k = l-1;
k < Um;
k++)
130 for (label j = l-1; j < Un; j++)
133 for (label
k = l-1;
k < Um;
k++)
135 s += U_[j][
k]*U_[i][
k];
138 for (label
k = l-1;
k < Um;
k++)
140 U_[j][
k] += s*rv1[
k];
143 for (label
k = l-1;
k < Um;
k++)
150 anorm =
max(anorm,
mag(S_[i]) +
mag(rv1[i]));
153 for (label i = Um-1; i >= 0; i--)
159 for (label j = l; j < Um; j++)
161 V_[j][i] = (U_[i][j]/U_[i][l])/g;
164 for (label j=l; j < Um; j++)
167 for (label
k = l;
k < Um;
k++)
169 s += U_[i][
k]*V_[
k][j];
172 for (label
k = l;
k < Um;
k++)
174 V_[
k][j] += s*V_[
k][i];
179 for (label j = l; j < Um;j++)
181 V_[i][j] = V_[j][i] = 0.0;
190 for (label i =
min(Um, Un) - 1; i >= 0; i--)
195 for (label j = l; j < Um; j++)
204 for (label j = l; j < Um; j++)
207 for (label
k = l;
k < Un;
k++)
209 s += U_[
k][i]*U_[
k][j];
212 scalar
f = (s/U_[i][i])*g;
214 for (label
k = i;
k < Un;
k++)
216 U_[
k][j] += f*U_[
k][i];
220 for (label j = i; j < Un; j++)
227 for (label j = i; j < Un; j++)
236 for (label
k = Um-1;
k >= 0;
k--)
238 for (label its = 0; its < 35; its++)
243 for (l =
k; l >= 0; l--)
246 if (
mag(rv1[l]) + anorm == anorm)
251 if (
mag(S_[nm]) + anorm == anorm)
break;
258 for (label i = l-1; i <
k+1; i++)
263 if (
mag(f) + anorm == anorm)
break;
272 for (label j = 0; j < Un; j++)
274 scalar
y = U_[j][nm];
276 U_[j][nm] = y*c + z*s;
277 U_[j][i] = z*c - y*s;
290 for (label j = 0; j < Um; j++) V_[j][
k] = -V_[j][
k];
299 "(scalarRectangularMatrix& A, const scalar minCondition)"
300 ) <<
"no convergence in 35 SVD iterations"
309 scalar
f = ((y - z)*(y + z) + (g -
h)*(g + h))/(2.0*h*y);
311 f = ((x - z)*(x + z) + h*((y/(f + sign(g, f))) -
h))/x;
315 for (label j = l; j <= nm; j++)
331 for (label jj = 0; jj < Um; jj++)
335 V_[jj][j] = x*c + z*s;
336 V_[jj][i] = z*c - x*s;
350 for (label jj=0; jj < Un; jj++)
354 U_[jj][j] = y*c + z*s;
355 U_[jj][i] = z*c - y*s;
365 const scalar minS = minCondition*S_[
findMax(S_)];
366 for (label i = 0; i < S_.size(); i++)