278
A. Yu. Morozov and D. L. Reviznikov
r = 0, jr = m / 2
ii = {1, 2, ..., n}
jj = {1, 2, ..., m}
jj = jj / jr
while r < min(n, m):
ir = -1
for i = 1, ..., n:
U[i][r] = F[i][jr];
for k = 1, ..., r:
U[i][r] -= U[i][k] * V[k][jr]
if (i in ii) and (ir < 0 or |U[i][r]| > |U[ir][r]|):
ir = i
if ir < 0 or |U[ir][r]| < eps:
break
ii = ii / ir
jr = -1
for j = 1, ..., m:
V[r][j] = F[ir][j]
for k = 1, ..., r:
V[r][j] -= U[ir][k] * V[k][j]
V[r][j] /= U[ir][r]
if (j in jj) and (jr < 0 or |V[r][j]| > |V[r][jr]|):
jr = j
if jr < 0 or |V[r][jr]| < eps:
r += 1
break
jj = jj / jr
r += 1
return U, V, r
Fig. 19.2 Pseudo-code of the algorithm
Fig. 19.3 Rows and
columns, from which the
entire matrix is restored