52 return std::fabs(a) * std::sqrt(1 + c*c);
106 for (k = 0; k < std::max(nct,nrt); k++) {
113 for (i = k; i <
m_; i++) {
114 s_[k] = hypot(
s_[k],A[i][k]);
120 for (i = k; i <
m_; i++) {
127 for (j = k+1; j <
n_; j++) {
128 if ((k < nct) && (
s_[k] != 0.0)) {
133 for (i = k; i <
m_; i++) {
134 t += A[i][k]*A[i][j];
137 for (i = k; i <
m_; i++) {
138 A[i][j] +=
t*A[i][k];
152 for (i = k; i <
m_; i++) {
162 for (i = k+1; i <
n_; i++) {
163 e[k] = hypot(e[k],e[i]);
169 for (i = k+1; i <
n_; i++) {
175 if ((k+1 <
m_) && (e[k] != 0.0)) {
179 for (i = k+1; i <
m_; i++) {
182 for (j = k+1; j <
n_; j++) {
183 for (i = k+1; i <
m_; i++) {
184 work[i] += e[j]*A[i][j];
187 for (j = k+1; j <
n_; j++) {
188 Real t = -e[j]/e[k+1];
189 for (i = k+1; i <
m_; i++) {
190 A[i][j] +=
t*work[i];
198 for (i = k+1; i <
n_; i++) {
207 s_[nct] = A[nct][nct];
210 e[nrt] = A[nrt][
n_-1];
216 for (j = nct; j <
n_; j++) {
217 for (i = 0; i <
m_; i++) {
222 for (k = nct-1; k >= 0; --k) {
224 for (j = k+1; j <
n_; ++j) {
226 for (i = k; i <
m_; i++) {
227 t +=
U_[i][k]*
U_[i][j];
230 for (i = k; i <
m_; i++) {
231 U_[i][j] +=
t*
U_[i][k];
234 for (i = k; i <
m_; i++ ) {
235 U_[i][k] = -
U_[i][k];
237 U_[k][k] = 1.0 +
U_[k][k];
238 for (i = 0; i < k-1; i++) {
242 for (i = 0; i <
m_; i++) {
251 for (k =
n_-1; k >= 0; --k) {
252 if ((k < nrt) && (e[k] != 0.0)) {
253 for (j = k+1; j <
n_; ++j) {
255 for (i = k+1; i <
n_; i++) {
256 t +=
V_[i][k]*
V_[i][j];
259 for (i = k+1; i <
n_; i++) {
260 V_[i][j] +=
t*
V_[i][k];
264 for (i = 0; i <
n_; i++) {
274 Real eps = std::pow(2.0,-52.0);
291 for (k = p-2; k >= -1; --k) {
295 if (std::fabs(e[k]) <= eps*(std::fabs(
s_[k]) +
296 std::fabs(
s_[k+1]))) {
305 for (ks = p-1; ks >= k; --ks) {
309 Real t = (ks != p ?
Real(std::fabs(e[ks])) : 0.) +
310 (ks != k+1 ?
Real(std::fabs(e[ks-1])) : 0.);
311 if (std::fabs(
s_[ks]) <= eps*
t) {
318 }
else if (ks == p-1) {
336 for (j = p-2; j >= k; --j) {
345 for (i = 0; i <
n_; i++) {
346 t = cs*
V_[i][j] + sn*
V_[i][p-1];
347 V_[i][p-1] = -sn*
V_[i][j] + cs*
V_[i][p-1];
359 for (j = k; j < p; j++) {
366 for (i = 0; i <
m_; i++) {
367 t = cs*
U_[i][j] + sn*
U_[i][k-1];
368 U_[i][k-1] = -sn*
U_[i][j] + cs*
U_[i][k-1];
380 Real scale = std::max(
383 std::max(std::fabs(
s_[p-1]),
389 Real spm1 =
s_[p-2]/scale;
390 Real epm1 = e[p-2]/scale;
392 Real ek = e[k]/scale;
393 Real b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0;
394 Real c = (sp*epm1)*(sp*epm1);
396 if ((
b != 0.0) || (c != 0.0)) {
397 shift = std::sqrt(
b*
b + c);
401 shift = c/(
b + shift);
403 Real f = (sk + sp)*(sk - sp) + shift;
408 for (j = k; j < p-1; j++) {
415 f = cs*
s_[j] + sn*e[j];
416 e[j] = cs*e[j] - sn*
s_[j];
418 s_[j+1] = cs*
s_[j+1];
419 for (i = 0; i <
n_; i++) {
420 t = cs*
V_[i][j] + sn*
V_[i][j+1];
421 V_[i][j+1] = -sn*
V_[i][j] + cs*
V_[i][j+1];
428 f = cs*e[j] + sn*
s_[j+1];
429 s_[j+1] = -sn*e[j] + cs*
s_[j+1];
433 for (i = 0; i <
m_; i++) {
434 t = cs*
U_[i][j] + sn*
U_[i][j+1];
435 U_[i][j+1] = -sn*
U_[i][j] + cs*
U_[i][j+1];
453 for (i = 0; i <= pp; i++) {
454 V_[i][k] = -
V_[i][k];
461 if (
s_[k] >=
s_[k+1]) {
466 for (i = 0; i <
n_; i++) {
471 for (i = 0; i <
m_; i++) {
530 const Size numericalRank = this->
rank();
531 for (
Size i=0; i<numericalRank; i++)
1-D array used in linear algebra.
Matrix used in linear algebra.
const Array & singularValues() const
Array solveFor(const Array &) const
ext::function< Real(Real)> b
QL_INTEGER Integer
integer number
std::size_t Size
size of a container
void swap(Array &v, Array &w) noexcept
Matrix transpose(const Matrix &m)
Matrix inverse(const Matrix &m)
ext::shared_ptr< YieldTermStructure > r
singular value decomposition