Coverage for /usr/lib/python3/dist-packages/mpmath/matrices/linalg.py: 6%

1"""

2Linear algebra

3--------------

5Linear equations

6................

8Basic linear algebra is implemented; you can for example solve the linear

9equation system::

11 x + 2*y = -10

12 3*x + 4*y = 10

14using ``lu_solve``::

16 >>> from mpmath import *

17 >>> mp.pretty = False

18 >>> A = matrix([[1, 2], [3, 4]])

19 >>> b = matrix([-10, 10])

20 >>> x = lu_solve(A, b)

21 >>> x

22 matrix(

23 [['30.0'],

24 ['-20.0']])

26If you don't trust the result, use ``residual`` to calculate the residual ||A*x-b||::

28 >>> residual(A, x, b)

29 matrix(

30 [['3.46944695195361e-18'],

31 ['3.46944695195361e-18']])

32 >>> str(eps)

33 '2.22044604925031e-16'

35As you can see, the solution is quite accurate. The error is caused by the

36inaccuracy of the internal floating point arithmetic. Though, it's even smaller

37than the current machine epsilon, which basically means you can trust the

38result.

40If you need more speed, use NumPy, or ``fp.lu_solve`` for a floating-point computation.

42 >>> fp.lu_solve(A, b)

43 matrix(

44 [['30.0'],

45 ['-20.0']])

47``lu_solve`` accepts overdetermined systems. It is usually not possible to solve

48such systems, so the residual is minimized instead. Internally this is done

49using Cholesky decomposition to compute a least squares approximation. This means

50that that ``lu_solve`` will square the errors. If you can't afford this, use

51``qr_solve`` instead. It is twice as slow but more accurate, and it calculates

52the residual automatically.

55Matrix factorization

56....................

58The function ``lu`` computes an explicit LU factorization of a matrix::

60 >>> P, L, U = lu(matrix([[0,2,3],[4,5,6],[7,8,9]]))

61 >>> print(P)

62 [0.0 0.0 1.0]

63 [1.0 0.0 0.0]

64 [0.0 1.0 0.0]

65 >>> print(L)

66 [ 1.0 0.0 0.0]

67 [ 0.0 1.0 0.0]

68 [0.571428571428571 0.214285714285714 1.0]

69 >>> print(U)

70 [7.0 8.0 9.0]

71 [0.0 2.0 3.0]

72 [0.0 0.0 0.214285714285714]

73 >>> print(P.T*L*U)

74 [0.0 2.0 3.0]

75 [4.0 5.0 6.0]

76 [7.0 8.0 9.0]

78Interval matrices

79-----------------

81Matrices may contain interval elements. This allows one to perform

82basic linear algebra operations such as matrix multiplication

83and equation solving with rigorous error bounds::

85 >>> a = iv.matrix([['0.1','0.3','1.0'],

86 ... ['7.1','5.5','4.8'],

87 ... ['3.2','4.4','5.6']])

88 >>>

89 >>> b = iv.matrix(['4','0.6','0.5'])

90 >>> c = iv.lu_solve(a, b)

91 >>> print(c)

92 [ [5.2582327113062568605927528666, 5.25823271130625686059275702219]]

93 [[-13.1550493962678375411635581388, -13.1550493962678375411635540152]]

94 [ [7.42069154774972557628979076189, 7.42069154774972557628979190734]]

95 >>> print(a*c)

96 [ [3.99999999999999999999999844904, 4.00000000000000000000000155096]]

97 [[0.599999999999999999999968898009, 0.600000000000000000000031763736]]

98 [[0.499999999999999999999979320485, 0.500000000000000000000020679515]]

99"""

100

101# TODO:

102# *implement high-level qr()

103# *test unitvector

104# *iterative solving

105

106from copy import copy

107

108from ..libmp.backend import xrange

109

110class LinearAlgebraMethods(object):

111

112 def LU_decomp(ctx, A, overwrite=False, use_cache=True):

113 """

114 LU-factorization of a n*n matrix using the Gauss algorithm.

115 Returns L and U in one matrix and the pivot indices.

116

117 Use overwrite to specify whether A will be overwritten with L and U.

118 """

119 if not A.rows == A.cols:

120 raise ValueError('need n*n matrix')

121 # get from cache if possible

122 if use_cache and isinstance(A, ctx.matrix) and A._LU:

123 return A._LU

124 if not overwrite:

125 orig = A

126 A = A.copy()

127 tol = ctx.absmin(ctx.mnorm(A,1) * ctx.eps) # each pivot element has to be bigger

128 n = A.rows

129 p = [None]*(n - 1)

130 for j in xrange(n - 1):

131 # pivoting, choose max(abs(reciprocal row sum)*abs(pivot element))

132 biggest = 0

133 for k in xrange(j, n):

134 s = ctx.fsum([ctx.absmin(A[k,l]) for l in xrange(j, n)])

135 if ctx.absmin(s) <= tol:

136 raise ZeroDivisionError('matrix is numerically singular')

137 current = 1/s * ctx.absmin(A[k,j])

138 if current > biggest: # TODO: what if equal?

139 biggest = current

140 p[j] = k

141 # swap rows according to p

142 ctx.swap_row(A, j, p[j])

143 if ctx.absmin(A[j,j]) <= tol:

144 raise ZeroDivisionError('matrix is numerically singular')

145 # calculate elimination factors and add rows

146 for i in xrange(j + 1, n):

147 A[i,j] /= A[j,j]

148 for k in xrange(j + 1, n):

149 A[i,k] -= A[i,j]*A[j,k]

150 if ctx.absmin(A[n - 1,n - 1]) <= tol:

151 raise ZeroDivisionError('matrix is numerically singular')

152 # cache decomposition

153 if not overwrite and isinstance(orig, ctx.matrix):

154 orig._LU = (A, p)

155 return A, p

156

157 def L_solve(ctx, L, b, p=None):

158 """

159 Solve the lower part of a LU factorized matrix for y.

160 """

161 if L.rows != L.cols:

162 raise RuntimeError("need n*n matrix")

163 n = L.rows

164 if len(b) != n:

165 raise ValueError("Value should be equal to n")

166 b = copy(b)

167 if p: # swap b according to p

168 for k in xrange(0, len(p)):

169 ctx.swap_row(b, k, p[k])

170 # solve

171 for i in xrange(1, n):

172 for j in xrange(i):

173 b[i] -= L[i,j] * b[j]

174 return b

175

176 def U_solve(ctx, U, y):

177 """

178 Solve the upper part of a LU factorized matrix for x.

179 """

180 if U.rows != U.cols:

181 raise RuntimeError("need n*n matrix")

182 n = U.rows

183 if len(y) != n:

184 raise ValueError("Value should be equal to n")

185 x = copy(y)

186 for i in xrange(n - 1, -1, -1):

187 for j in xrange(i + 1, n):

188 x[i] -= U[i,j] * x[j]

189 x[i] /= U[i,i]

190 return x

191

192 def lu_solve(ctx, A, b, **kwargs):

193 """

194 Ax = b => x

195

196 Solve a determined or overdetermined linear equations system.

197 Fast LU decomposition is used, which is less accurate than QR decomposition

198 (especially for overdetermined systems), but it's twice as efficient.

199 Use qr_solve if you want more precision or have to solve a very ill-

200 conditioned system.

201

202 If you specify real=True, it does not check for overdeterminded complex

203 systems.

204 """

205 prec = ctx.prec

206 try:

207 ctx.prec += 10

208 # do not overwrite A nor b

209 A, b = ctx.matrix(A, **kwargs).copy(), ctx.matrix(b, **kwargs).copy()

210 if A.rows < A.cols:

211 raise ValueError('cannot solve underdetermined system')

212 if A.rows > A.cols:

213 # use least-squares method if overdetermined

214 # (this increases errors)

215 AH = A.H

216 A = AH * A

217 b = AH * b

218 if (kwargs.get('real', False) or

219 not sum(type(i) is ctx.mpc for i in A)):

220 # TODO: necessary to check also b?

221 x = ctx.cholesky_solve(A, b)

222 else:

223 x = ctx.lu_solve(A, b)

224 else:

225 # LU factorization

226 A, p = ctx.LU_decomp(A)

227 b = ctx.L_solve(A, b, p)

228 x = ctx.U_solve(A, b)

229 finally:

230 ctx.prec = prec

231 return x

232

233 def improve_solution(ctx, A, x, b, maxsteps=1):

234 """

235 Improve a solution to a linear equation system iteratively.

236

237 This re-uses the LU decomposition and is thus cheap.

238 Usually 3 up to 4 iterations are giving the maximal improvement.

239 """

240 if A.rows != A.cols:

241 raise RuntimeError("need n*n matrix") # TODO: really?

242 for _ in xrange(maxsteps):

243 r = ctx.residual(A, x, b)

244 if ctx.norm(r, 2) < 10*ctx.eps:

245 break

246 # this uses cached LU decomposition and is thus cheap

247 dx = ctx.lu_solve(A, -r)

248 x += dx

249 return x

250

251 def lu(ctx, A):

252 """

253 A -> P, L, U

254

255 LU factorisation of a square matrix A. L is the lower, U the upper part.

256 P is the permutation matrix indicating the row swaps.

257

258 P*A = L*U

259

260 If you need efficiency, use the low-level method LU_decomp instead, it's

261 much more memory efficient.

262 """

263 # get factorization

264 A, p = ctx.LU_decomp(A)

265 n = A.rows

266 L = ctx.matrix(n)

267 U = ctx.matrix(n)

268 for i in xrange(n):

269 for j in xrange(n):

270 if i > j:

271 L[i,j] = A[i,j]

272 elif i == j:

273 L[i,j] = 1

274 U[i,j] = A[i,j]

275 else:

276 U[i,j] = A[i,j]

277 # calculate permutation matrix

278 P = ctx.eye(n)

279 for k in xrange(len(p)):

280 ctx.swap_row(P, k, p[k])

281 return P, L, U

282

283 def unitvector(ctx, n, i):

284 """

285 Return the i-th n-dimensional unit vector.

286 """

287 assert 0 < i <= n, 'this unit vector does not exist'

288 return [ctx.zero]*(i-1) + [ctx.one] + [ctx.zero]*(n-i)

289

290 def inverse(ctx, A, **kwargs):

291 """

292 Calculate the inverse of a matrix.

293

294 If you want to solve an equation system Ax = b, it's recommended to use

295 solve(A, b) instead, it's about 3 times more efficient.

296 """

297 prec = ctx.prec

298 try:

299 ctx.prec += 10

300 # do not overwrite A

301 A = ctx.matrix(A, **kwargs).copy()

302 n = A.rows

303 # get LU factorisation

304 A, p = ctx.LU_decomp(A)

305 cols = []

306 # calculate unit vectors and solve corresponding system to get columns

307 for i in xrange(1, n + 1):

308 e = ctx.unitvector(n, i)

309 y = ctx.L_solve(A, e, p)

310 cols.append(ctx.U_solve(A, y))

311 # convert columns to matrix

312 inv = []

313 for i in xrange(n):

314 row = []

315 for j in xrange(n):

316 row.append(cols[j][i])

317 inv.append(row)

318 result = ctx.matrix(inv, **kwargs)

319 finally:

320 ctx.prec = prec

321 return result

322

323 def householder(ctx, A):

324 """

325 (A|b) -> H, p, x, res

326

327 (A|b) is the coefficient matrix with left hand side of an optionally

328 overdetermined linear equation system.

329 H and p contain all information about the transformation matrices.

330 x is the solution, res the residual.

331 """

332 if not isinstance(A, ctx.matrix):

333 raise TypeError("A should be a type of ctx.matrix")

334 m = A.rows

335 n = A.cols

336 if m < n - 1:

337 raise RuntimeError("Columns should not be less than rows")

338 # calculate Householder matrix

339 p = []

340 for j in xrange(0, n - 1):

341 s = ctx.fsum(abs(A[i,j])**2 for i in xrange(j, m))

342 if not abs(s) > ctx.eps:

343 raise ValueError('matrix is numerically singular')

344 p.append(-ctx.sign(ctx.re(A[j,j])) * ctx.sqrt(s))

345 kappa = ctx.one / (s - p[j] * A[j,j])

346 A[j,j] -= p[j]

347 for k in xrange(j+1, n):

348 y = ctx.fsum(ctx.conj(A[i,j]) * A[i,k] for i in xrange(j, m)) * kappa

349 for i in xrange(j, m):

350 A[i,k] -= A[i,j] * y

351 # solve Rx = c1

352 x = [A[i,n - 1] for i in xrange(n - 1)]

353 for i in xrange(n - 2, -1, -1):

354 x[i] -= ctx.fsum(A[i,j] * x[j] for j in xrange(i + 1, n - 1))

355 x[i] /= p[i]

356 # calculate residual

357 if not m == n - 1:

358 r = [A[m-1-i, n-1] for i in xrange(m - n + 1)]

359 else:

360 # determined system, residual should be 0

361 r = [0]*m # maybe a bad idea, changing r[i] will change all elements

362 return A, p, x, r

363

364 #def qr(ctx, A):

365 # """

366 # A -> Q, R

367 #

368 # QR factorisation of a square matrix A using Householder decomposition.

369 # Q is orthogonal, this leads to very few numerical errors.

370 #

371 # A = Q*R

372 # """

373 # H, p, x, res = householder(A)

374 # TODO: implement this

375

376 def residual(ctx, A, x, b, **kwargs):

377 """

378 Calculate the residual of a solution to a linear equation system.

379

380 r = A*x - b for A*x = b

381 """

382 oldprec = ctx.prec

383 try:

384 ctx.prec *= 2

385 A, x, b = ctx.matrix(A, **kwargs), ctx.matrix(x, **kwargs), ctx.matrix(b, **kwargs)

386 return A*x - b

387 finally:

388 ctx.prec = oldprec

389

390 def qr_solve(ctx, A, b, norm=None, **kwargs):

391 """

392 Ax = b => x, ||Ax - b||

393

394 Solve a determined or overdetermined linear equations system and

395 calculate the norm of the residual (error).

396 QR decomposition using Householder factorization is applied, which gives very

397 accurate results even for ill-conditioned matrices. qr_solve is twice as

398 efficient.

399 """

400 if norm is None:

401 norm = ctx.norm

402 prec = ctx.prec

403 try:

404 ctx.prec += 10

405 # do not overwrite A nor b

406 A, b = ctx.matrix(A, **kwargs).copy(), ctx.matrix(b, **kwargs).copy()

407 if A.rows < A.cols:

408 raise ValueError('cannot solve underdetermined system')

409 H, p, x, r = ctx.householder(ctx.extend(A, b))

410 res = ctx.norm(r)

411 # calculate residual "manually" for determined systems

412 if res == 0:

413 res = ctx.norm(ctx.residual(A, x, b))

414 return ctx.matrix(x, **kwargs), res

415 finally:

416 ctx.prec = prec

417

418 def cholesky(ctx, A, tol=None):

419 r"""

420 Cholesky decomposition of a symmetric positive-definite matrix `A`.

421 Returns a lower triangular matrix `L` such that `A = L \times L^T`.

422 More generally, for a complex Hermitian positive-definite matrix,

423 a Cholesky decomposition satisfying `A = L \times L^H` is returned.

424

425 The Cholesky decomposition can be used to solve linear equation

426 systems twice as efficiently as LU decomposition, or to

427 test whether `A` is positive-definite.

428

429 The optional parameter ``tol`` determines the tolerance for

430 verifying positive-definiteness.

431

432 **Examples**

433

434 Cholesky decomposition of a positive-definite symmetric matrix::

435

436 >>> from mpmath import *

437 >>> mp.dps = 25; mp.pretty = True

438 >>> A = eye(3) + hilbert(3)

439 >>> nprint(A)

440 [ 2.0 0.5 0.333333]

441 [ 0.5 1.33333 0.25]

442 [0.333333 0.25 1.2]

443 >>> L = cholesky(A)

444 >>> nprint(L)

445 [ 1.41421 0.0 0.0]

446 [0.353553 1.09924 0.0]

447 [0.235702 0.15162 1.05899]

448 >>> chop(A - L*L.T)

449 [0.0 0.0 0.0]

450 [0.0 0.0 0.0]

451 [0.0 0.0 0.0]

452

453 Cholesky decomposition of a Hermitian matrix::

454

455 >>> A = eye(3) + matrix([[0,0.25j,-0.5j],[-0.25j,0,0],[0.5j,0,0]])

456 >>> L = cholesky(A)

457 >>> nprint(L)

458 [ 1.0 0.0 0.0]

459 [(0.0 - 0.25j) (0.968246 + 0.0j) 0.0]

460 [ (0.0 + 0.5j) (0.129099 + 0.0j) (0.856349 + 0.0j)]

461 >>> chop(A - L*L.H)

462 [0.0 0.0 0.0]

463 [0.0 0.0 0.0]

464 [0.0 0.0 0.0]

465

466 Attempted Cholesky decomposition of a matrix that is not positive

467 definite::

468

469 >>> A = -eye(3) + hilbert(3)

470 >>> L = cholesky(A)

471 Traceback (most recent call last):

472 ...

473 ValueError: matrix is not positive-definite

474

475 **References**

476

477 1. [Wikipedia]_ http://en.wikipedia.org/wiki/Cholesky_decomposition

478

479 """

480 if not isinstance(A, ctx.matrix):

481 raise RuntimeError("A should be a type of ctx.matrix")

482 if not A.rows == A.cols:

483 raise ValueError('need n*n matrix')

484 if tol is None:

485 tol = +ctx.eps

486 n = A.rows

487 L = ctx.matrix(n)

488 for j in xrange(n):

489 c = ctx.re(A[j,j])

490 if abs(c-A[j,j]) > tol:

491 raise ValueError('matrix is not Hermitian')

492 s = c - ctx.fsum((L[j,k] for k in xrange(j)),

493 absolute=True, squared=True)

494 if s < tol:

495 raise ValueError('matrix is not positive-definite')

496 L[j,j] = ctx.sqrt(s)

497 for i in xrange(j, n):

498 it1 = (L[i,k] for k in xrange(j))

499 it2 = (L[j,k] for k in xrange(j))

500 t = ctx.fdot(it1, it2, conjugate=True)

501 L[i,j] = (A[i,j] - t) / L[j,j]

502 return L

503

504 def cholesky_solve(ctx, A, b, **kwargs):

505 """

506 Ax = b => x

507

508 Solve a symmetric positive-definite linear equation system.

509 This is twice as efficient as lu_solve.

510

511 Typical use cases:

512 * A.T*A

513 * Hessian matrix

514 * differential equations

515 """

516 prec = ctx.prec

517 try:

518 ctx.prec += 10

519 # do not overwrite A nor b

520 A, b = ctx.matrix(A, **kwargs).copy(), ctx.matrix(b, **kwargs).copy()

521 if A.rows != A.cols:

522 raise ValueError('can only solve determined system')

523 # Cholesky factorization

524 L = ctx.cholesky(A)

525 # solve

526 n = L.rows

527 if len(b) != n:

528 raise ValueError("Value should be equal to n")

529 for i in xrange(n):

530 b[i] -= ctx.fsum(L[i,j] * b[j] for j in xrange(i))

531 b[i] /= L[i,i]

532 x = ctx.U_solve(L.T, b)

533 return x

534 finally:

535 ctx.prec = prec

536

537 def det(ctx, A):

538 """

539 Calculate the determinant of a matrix.

540 """

541 prec = ctx.prec

542 try:

543 # do not overwrite A

544 A = ctx.matrix(A).copy()

545 # use LU factorization to calculate determinant

546 try:

547 R, p = ctx.LU_decomp(A)

548 except ZeroDivisionError:

549 return 0

550 z = 1

551 for i, e in enumerate(p):

552 if i != e:

553 z *= -1

554 for i in xrange(A.rows):

555 z *= R[i,i]

556 return z

557 finally:

558 ctx.prec = prec

559

560 def cond(ctx, A, norm=None):

561 """

562 Calculate the condition number of a matrix using a specified matrix norm.

563

564 The condition number estimates the sensitivity of a matrix to errors.

565 Example: small input errors for ill-conditioned coefficient matrices

566 alter the solution of the system dramatically.

567

568 For ill-conditioned matrices it's recommended to use qr_solve() instead

569 of lu_solve(). This does not help with input errors however, it just avoids

570 to add additional errors.

571

572 Definition: cond(A) = ||A|| * ||A**-1||

573 """

574 if norm is None:

575 norm = lambda x: ctx.mnorm(x,1)

576 return norm(A) * norm(ctx.inverse(A))

577

578 def lu_solve_mat(ctx, a, b):

579 """Solve a * x = b where a and b are matrices."""

580 r = ctx.matrix(a.rows, b.cols)

581 for i in range(b.cols):

582 c = ctx.lu_solve(a, b.column(i))

583 for j in range(len(c)):

584 r[j, i] = c[j]

585 return r

586

587 def qr(ctx, A, mode = 'full', edps = 10):

588 """

589 Compute a QR factorization $A = QR$ where

590 A is an m x n matrix of real or complex numbers where m >= n

591

592 mode has following meanings:

593 (1) mode = 'raw' returns two matrixes (A, tau) in the

594 internal format used by LAPACK

595 (2) mode = 'skinny' returns the leading n columns of Q

596 and n rows of R

597 (3) Any other value returns the leading m columns of Q

598 and m rows of R

599

600 edps is the increase in mp precision used for calculations

601

602 **Examples**

603

604 >>> from mpmath import *

605 >>> mp.dps = 15

606 >>> mp.pretty = True

607 >>> A = matrix([[1, 2], [3, 4], [1, 1]])

608 >>> Q, R = qr(A)

609 >>> Q

610 [-0.301511344577764 0.861640436855329 0.408248290463863]

611 [-0.904534033733291 -0.123091490979333 -0.408248290463863]

612 [-0.301511344577764 -0.492365963917331 0.816496580927726]

613 >>> R

614 [-3.3166247903554 -4.52267016866645]

615 [ 0.0 0.738548945875996]

616 [ 0.0 0.0]

617 >>> Q * R

618 [1.0 2.0]

619 [3.0 4.0]

620 [1.0 1.0]

621 >>> chop(Q.T * Q)

622 [1.0 0.0 0.0]

623 [0.0 1.0 0.0]

624 [0.0 0.0 1.0]

625 >>> B = matrix([[1+0j, 2-3j], [3+j, 4+5j]])

626 >>> Q, R = qr(B)

627 >>> nprint(Q)

628 [ (-0.301511 + 0.0j) (0.0695795 - 0.95092j)]

629 [(-0.904534 - 0.301511j) (-0.115966 + 0.278318j)]

630 >>> nprint(R)

631 [(-3.31662 + 0.0j) (-5.72872 - 2.41209j)]

632 [ 0.0 (3.91965 + 0.0j)]

633 >>> Q * R

634 [(1.0 + 0.0j) (2.0 - 3.0j)]

635 [(3.0 + 1.0j) (4.0 + 5.0j)]

636 >>> chop(Q.T * Q.conjugate())

637 [1.0 0.0]

638 [0.0 1.0]

639

640 """

641

642 # check values before continuing

643 assert isinstance(A, ctx.matrix)

644 m = A.rows

645 n = A.cols

646 assert n > 1

647 assert m >= n

648 assert edps >= 0

649

650 # check for complex data type

651 cmplx = any(type(x) is ctx.mpc for x in A)

652

653 # temporarily increase the precision and initialize

654 with ctx.extradps(edps):

655 tau = ctx.matrix(n,1)

656 A = A.copy()

657

658 # ---------------

659 # FACTOR MATRIX A

660 # ---------------

661 if cmplx:

662 one = ctx.mpc('1.0', '0.0')

663 zero = ctx.mpc('0.0', '0.0')

664 rzero = ctx.mpf('0.0')

665

666 # main loop to factor A (complex)

667 for j in xrange(0, n):

668 alpha = A[j,j]

669 alphr = ctx.re(alpha)

670 alphi = ctx.im(alpha)

671

672 if (m-j) >= 2:

673 xnorm = ctx.fsum( A[i,j]*ctx.conj(A[i,j]) for i in xrange(j+1, m) )

674 xnorm = ctx.re( ctx.sqrt(xnorm) )

675 else:

676 xnorm = rzero

677

678 if (xnorm == rzero) and (alphi == rzero):

679 tau[j] = zero

680 continue

681

682 if alphr < rzero:

683 beta = ctx.sqrt(alphr**2 + alphi**2 + xnorm**2)

684 else:

685 beta = -ctx.sqrt(alphr**2 + alphi**2 + xnorm**2)

686

687 tau[j] = ctx.mpc( (beta - alphr) / beta, -alphi / beta )

688 t = -ctx.conj(tau[j])

689 za = one / (alpha - beta)

690

691 for i in xrange(j+1, m):

692 A[i,j] *= za

693

694 A[j,j] = one

695 for k in xrange(j+1, n):

696 y = ctx.fsum(A[i,j] * ctx.conj(A[i,k]) for i in xrange(j, m))

697 temp = t * ctx.conj(y)

698 for i in xrange(j, m):

699 A[i,k] += A[i,j] * temp

700

701 A[j,j] = ctx.mpc(beta, '0.0')

702 else:

703 one = ctx.mpf('1.0')

704 zero = ctx.mpf('0.0')

705

706 # main loop to factor A (real)

707 for j in xrange(0, n):

708 alpha = A[j,j]

709

710 if (m-j) > 2:

711 xnorm = ctx.fsum( (A[i,j])**2 for i in xrange(j+1, m) )

712 xnorm = ctx.sqrt(xnorm)

713 elif (m-j) == 2:

714 xnorm = abs( A[m-1,j] )

715 else:

716 xnorm = zero

717

718 if xnorm == zero:

719 tau[j] = zero

720 continue

721

722 if alpha < zero:

723 beta = ctx.sqrt(alpha**2 + xnorm**2)

724 else:

725 beta = -ctx.sqrt(alpha**2 + xnorm**2)

726

727 tau[j] = (beta - alpha) / beta

728 t = -tau[j]

729 da = one / (alpha - beta)

730

731 for i in xrange(j+1, m):

732 A[i,j] *= da

733

734 A[j,j] = one

735 for k in xrange(j+1, n):

736 y = ctx.fsum( A[i,j] * A[i,k] for i in xrange(j, m) )

737 temp = t * y

738 for i in xrange(j,m):

739 A[i,k] += A[i,j] * temp

740

741 A[j,j] = beta

742

743 # return factorization in same internal format as LAPACK

744 if (mode == 'raw') or (mode == 'RAW'):

745 return A, tau

746

747 # ----------------------------------

748 # FORM Q USING BACKWARD ACCUMULATION

749 # ----------------------------------

750

751 # form R before the values are overwritten

752 R = A.copy()

753 for j in xrange(0, n):

754 for i in xrange(j+1, m):

755 R[i,j] = zero

756

757 # set the value of p (number of columns of Q to return)

758 p = m

759 if (mode == 'skinny') or (mode == 'SKINNY'):

760 p = n

761

762 # add columns to A if needed and initialize

763 A.cols += (p-n)

764 for j in xrange(0, p):

765 A[j,j] = one

766 for i in xrange(0, j):

767 A[i,j] = zero

768

769 # main loop to form Q

770 for j in xrange(n-1, -1, -1):

771 t = -tau[j]

772 A[j,j] += t

773

774 for k in xrange(j+1, p):

775 if cmplx:

776 y = ctx.fsum(A[i,j] * ctx.conj(A[i,k]) for i in xrange(j+1, m))

777 temp = t * ctx.conj(y)

778 else:

779 y = ctx.fsum(A[i,j] * A[i,k] for i in xrange(j+1, m))

780 temp = t * y

781 A[j,k] = temp

782 for i in xrange(j+1, m):

783 A[i,k] += A[i,j] * temp

784

785 for i in xrange(j+1, m):

786 A[i, j] *= t

787

788 return A, R[0:p,0:n]

789

790 # ------------------

791 # END OF FUNCTION QR

792 # ------------------