Coverage for /usr/lib/python3/dist-packages/scipy/sparse/linalg/_eigen/lobpcg/lobpcg.py: 3%

1""" 

2Locally Optimal Block Preconditioned Conjugate Gradient Method (LOBPCG). 

3 

4References 

5---------- 

6.. [1] A. V. Knyazev (2001), 

7 Toward the Optimal Preconditioned Eigensolver: Locally Optimal 

8 Block Preconditioned Conjugate Gradient Method. 

9 SIAM Journal on Scientific Computing 23, no. 2, 

10 pp. 517-541. :doi:`10.1137/S1064827500366124` 

11 

12.. [2] A. V. Knyazev, I. Lashuk, M. E. Argentati, and E. Ovchinnikov (2007), 

13 Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) 

14 in hypre and PETSc. :arxiv:`0705.2626` 

15 

16.. [3] A. V. Knyazev's C and MATLAB implementations: 

17 https://github.com/lobpcg/blopex 

18""" 

19 

20import warnings 

21import numpy as np 

22from scipy.linalg import (inv, eigh, cho_factor, cho_solve, 

23 cholesky, LinAlgError) 

24from scipy.sparse.linalg import LinearOperator 

25from scipy.sparse import issparse 

26 

27__all__ = ["lobpcg"] 

28 

29 

30def _report_nonhermitian(M, name): 

31 """ 

32 Report if `M` is not a Hermitian matrix given its type. 

33 """ 

34 from scipy.linalg import norm 

35 

36 md = M - M.T.conj() 

37 nmd = norm(md, 1) 

38 tol = 10 * np.finfo(M.dtype).eps 

39 tol = max(tol, tol * norm(M, 1)) 

40 if nmd > tol: 

41 warnings.warn( 

42 f"Matrix {name} of the type {M.dtype} is not Hermitian: " 

43 f"condition: {nmd} < {tol} fails.", 

44 UserWarning, stacklevel=4 

45 ) 

46 

47def _as2d(ar): 

48 """ 

49 If the input array is 2D return it, if it is 1D, append a dimension, 

50 making it a column vector. 

51 """ 

52 if ar.ndim == 2: 

53 return ar 

54 else: # Assume 1! 

55 aux = np.array(ar, copy=False) 

56 aux.shape = (ar.shape[0], 1) 

57 return aux 

58 

59 

60def _makeMatMat(m): 

61 if m is None: 

62 return None 

63 elif callable(m): 

64 return lambda v: m(v) 

65 else: 

66 return lambda v: m @ v 

67 

68 

69def _matmul_inplace(x, y, verbosityLevel=0): 

70 """Perform 'np.matmul' in-place if possible. 

71 

72 If some sufficient conditions for inplace matmul are met, do so. 

73 Otherwise try inplace update and fall back to overwrite if that fails. 

74 """ 

75 if x.flags["CARRAY"] and x.shape[1] == y.shape[1] and x.dtype == y.dtype: 

76 # conditions where we can guarantee that inplace updates will work; 

77 # i.e. x is not a view/slice, x & y have compatible dtypes, and the 

78 # shape of the result of x @ y matches the shape of x. 

79 np.matmul(x, y, out=x) 

80 else: 

81 # ideally, we'd have an exhaustive list of conditions above when 

82 # inplace updates are possible; since we don't, we opportunistically 

83 # try if it works, and fall back to overwriting if necessary 

84 try: 

85 np.matmul(x, y, out=x) 

86 except Exception: 

87 if verbosityLevel: 

88 warnings.warn( 

89 "Inplace update of x = x @ y failed, " 

90 "x needs to be overwritten.", 

91 UserWarning, stacklevel=3 

92 ) 

93 x = x @ y 

94 return x 

95 

96 

97def _applyConstraints(blockVectorV, factYBY, blockVectorBY, blockVectorY): 

98 """Changes blockVectorV in-place.""" 

99 YBV = blockVectorBY.T.conj() @ blockVectorV 

100 tmp = cho_solve(factYBY, YBV) 

101 blockVectorV -= blockVectorY @ tmp 

102 

103 

104def _b_orthonormalize(B, blockVectorV, blockVectorBV=None, 

105 verbosityLevel=0): 

106 """in-place B-orthonormalize the given block vector using Cholesky.""" 

107 if blockVectorBV is None: 

108 if B is None: 

109 blockVectorBV = blockVectorV 

110 else: 

111 try: 

112 blockVectorBV = B(blockVectorV) 

113 except Exception as e: 

114 if verbosityLevel: 

115 warnings.warn( 

116 f"Secondary MatMul call failed with error\n" 

117 f"{e}\n", 

118 UserWarning, stacklevel=3 

119 ) 

120 return None, None, None 

121 if blockVectorBV.shape != blockVectorV.shape: 

122 raise ValueError( 

123 f"The shape {blockVectorV.shape} " 

124 f"of the orthogonalized matrix not preserved\n" 

125 f"and changed to {blockVectorBV.shape} " 

126 f"after multiplying by the secondary matrix.\n" 

127 ) 

128 

129 VBV = blockVectorV.T.conj() @ blockVectorBV 

130 try: 

131 # VBV is a Cholesky factor from now on... 

132 VBV = cholesky(VBV, overwrite_a=True) 

133 VBV = inv(VBV, overwrite_a=True) 

134 blockVectorV = _matmul_inplace( 

135 blockVectorV, VBV, 

136 verbosityLevel=verbosityLevel 

137 ) 

138 if B is not None: 

139 blockVectorBV = _matmul_inplace( 

140 blockVectorBV, VBV, 

141 verbosityLevel=verbosityLevel 

142 ) 

143 return blockVectorV, blockVectorBV, VBV 

144 except LinAlgError: 

145 if verbosityLevel: 

146 warnings.warn( 

147 "Cholesky has failed.", 

148 UserWarning, stacklevel=3 

149 ) 

150 return None, None, None 

151 

152 

153def _get_indx(_lambda, num, largest): 

154 """Get `num` indices into `_lambda` depending on `largest` option.""" 

155 ii = np.argsort(_lambda) 

156 if largest: 

157 ii = ii[:-num - 1:-1] 

158 else: 

159 ii = ii[:num] 

160 

161 return ii 

162 

163 

164def _handle_gramA_gramB_verbosity(gramA, gramB, verbosityLevel): 

165 if verbosityLevel: 

166 _report_nonhermitian(gramA, "gramA") 

167 _report_nonhermitian(gramB, "gramB") 

168 

169 

170def lobpcg( 

171 A, 

172 X, 

173 B=None, 

174 M=None, 

175 Y=None, 

176 tol=None, 

177 maxiter=None, 

178 largest=True, 

179 verbosityLevel=0, 

180 retLambdaHistory=False, 

181 retResidualNormsHistory=False, 

182 restartControl=20, 

183): 

184 """Locally Optimal Block Preconditioned Conjugate Gradient Method (LOBPCG). 

185 

186 LOBPCG is a preconditioned eigensolver for large real symmetric and complex 

187 Hermitian definite generalized eigenproblems. 

188 

189 Parameters 

190 ---------- 

191 A : {sparse matrix, ndarray, LinearOperator, callable object} 

192 The Hermitian linear operator of the problem, usually given by a 

193 sparse matrix. Often called the "stiffness matrix". 

194 X : ndarray, float32 or float64 

195 Initial approximation to the ``k`` eigenvectors (non-sparse). 

196 If `A` has ``shape=(n,n)`` then `X` must have ``shape=(n,k)``. 

197 B : {sparse matrix, ndarray, LinearOperator, callable object} 

198 Optional. By default ``B = None``, which is equivalent to identity. 

199 The right hand side operator in a generalized eigenproblem if present. 

200 Often called the "mass matrix". Must be Hermitian positive definite. 

201 M : {sparse matrix, ndarray, LinearOperator, callable object} 

202 Optional. By default ``M = None``, which is equivalent to identity. 

203 Preconditioner aiming to accelerate convergence. 

204 Y : ndarray, float32 or float64, default: None 

205 An ``n-by-sizeY`` ndarray of constraints with ``sizeY < n``. 

206 The iterations will be performed in the ``B``-orthogonal complement 

207 of the column-space of `Y`. `Y` must be full rank if present. 

208 tol : scalar, optional 

209 The default is ``tol=n*sqrt(eps)``. 

210 Solver tolerance for the stopping criterion. 

211 maxiter : int, default: 20 

212 Maximum number of iterations. 

213 largest : bool, default: True 

214 When True, solve for the largest eigenvalues, otherwise the smallest. 

215 verbosityLevel : int, optional 

216 By default ``verbosityLevel=0`` no output. 

217 Controls the solver standard/screen output. 

218 retLambdaHistory : bool, default: False 

219 Whether to return iterative eigenvalue history. 

220 retResidualNormsHistory : bool, default: False 

221 Whether to return iterative history of residual norms. 

222 restartControl : int, optional. 

223 Iterations restart if the residuals jump ``2**restartControl`` times 

224 compared to the smallest recorded in ``retResidualNormsHistory``. 

225 The default is ``restartControl=20``, making the restarts rare for 

226 backward compatibility. 

227 

228 Returns 

229 ------- 

230 lambda : ndarray of the shape ``(k, )``. 

231 Array of ``k`` approximate eigenvalues. 

232 v : ndarray of the same shape as ``X.shape``. 

233 An array of ``k`` approximate eigenvectors. 

234 lambdaHistory : ndarray, optional. 

235 The eigenvalue history, if `retLambdaHistory` is ``True``. 

236 ResidualNormsHistory : ndarray, optional. 

237 The history of residual norms, if `retResidualNormsHistory` 

238 is ``True``. 

239 

240 Notes 

241 ----- 

242 The iterative loop runs ``maxit=maxiter`` (20 if ``maxit=None``) 

243 iterations at most and finishes earler if the tolerance is met. 

244 Breaking backward compatibility with the previous version, LOBPCG 

245 now returns the block of iterative vectors with the best accuracy rather 

246 than the last one iterated, as a cure for possible divergence. 

247 

248 If ``X.dtype == np.float32`` and user-provided operations/multiplications 

249 by `A`, `B`, and `M` all preserve the ``np.float32`` data type, 

250 all the calculations and the output are in ``np.float32``. 

251 

252 The size of the iteration history output equals to the number of the best 

253 (limited by `maxit`) iterations plus 3: initial, final, and postprocessing. 

254 

255 If both `retLambdaHistory` and `retResidualNormsHistory` are ``True``, 

256 the return tuple has the following format 

257 ``(lambda, V, lambda history, residual norms history)``. 

258 

259 In the following ``n`` denotes the matrix size and ``k`` the number 

260 of required eigenvalues (smallest or largest). 

261 

262 The LOBPCG code internally solves eigenproblems of the size ``3k`` on every 

263 iteration by calling the dense eigensolver `eigh`, so if ``k`` is not 

264 small enough compared to ``n``, it makes no sense to call the LOBPCG code. 

265 Moreover, if one calls the LOBPCG algorithm for ``5k > n``, it would likely 

266 break internally, so the code calls the standard function `eigh` instead. 

267 It is not that ``n`` should be large for the LOBPCG to work, but rather the 

268 ratio ``n / k`` should be large. It you call LOBPCG with ``k=1`` 

269 and ``n=10``, it works though ``n`` is small. The method is intended 

270 for extremely large ``n / k``. 

271 

272 The convergence speed depends basically on three factors: 

273 

274 1. Quality of the initial approximations `X` to the seeking eigenvectors. 

275 Randomly distributed around the origin vectors work well if no better 

276 choice is known. 

277 

278 2. Relative separation of the desired eigenvalues from the rest 

279 of the eigenvalues. One can vary ``k`` to improve the separation. 

280 

281 3. Proper preconditioning to shrink the spectral spread. 

282 For example, a rod vibration test problem (under tests 

283 directory) is ill-conditioned for large ``n``, so convergence will be 

284 slow, unless efficient preconditioning is used. For this specific 

285 problem, a good simple preconditioner function would be a linear solve 

286 for `A`, which is easy to code since `A` is tridiagonal. 

287 

288 References 

289 ---------- 

290 .. [1] A. V. Knyazev (2001), 

291 Toward the Optimal Preconditioned Eigensolver: Locally Optimal 

292 Block Preconditioned Conjugate Gradient Method. 

293 SIAM Journal on Scientific Computing 23, no. 2, 

294 pp. 517-541. :doi:`10.1137/S1064827500366124` 

295 

296 .. [2] A. V. Knyazev, I. Lashuk, M. E. Argentati, and E. Ovchinnikov 

297 (2007), Block Locally Optimal Preconditioned Eigenvalue Xolvers 

298 (BLOPEX) in hypre and PETSc. :arxiv:`0705.2626` 

299 

300 .. [3] A. V. Knyazev's C and MATLAB implementations: 

301 https://github.com/lobpcg/blopex 

302 

303 Examples 

304 -------- 

305 Our first example is minimalistic - find the largest eigenvalue of 

306 a diagonal matrix by solving the non-generalized eigenvalue problem 

307 ``A x = lambda x`` without constraints or preconditioning. 

308 

309 >>> import numpy as np 

310 >>> from scipy.sparse import spdiags 

311 >>> from scipy.sparse.linalg import LinearOperator, aslinearoperator 

312 >>> from scipy.sparse.linalg import lobpcg 

313 

314 The square matrix size is 

315 

316 >>> n = 100 

317 

318 and its diagonal entries are 1, ..., 100 defined by 

319 

320 >>> vals = np.arange(1, n + 1).astype(np.int16) 

321 

322 The first mandatory input parameter in this test is 

323 the sparse diagonal matrix `A` 

324 of the eigenvalue problem ``A x = lambda x`` to solve. 

325 

326 >>> A = spdiags(vals, 0, n, n) 

327 >>> A = A.astype(np.int16) 

328 >>> A.toarray() 

329 array([[ 1, 0, 0, ..., 0, 0, 0], 

330 [ 0, 2, 0, ..., 0, 0, 0], 

331 [ 0, 0, 3, ..., 0, 0, 0], 

332 ..., 

333 [ 0, 0, 0, ..., 98, 0, 0], 

334 [ 0, 0, 0, ..., 0, 99, 0], 

335 [ 0, 0, 0, ..., 0, 0, 100]], dtype=int16) 

336 

337 The second mandatory input parameter `X` is a 2D array with the 

338 row dimension determining the number of requested eigenvalues. 

339 `X` is an initial guess for targeted eigenvectors. 

340 `X` must have linearly independent columns. 

341 If no initial approximations available, randomly oriented vectors 

342 commonly work best, e.g., with components normally distributed 

343 around zero or uniformly distributed on the interval [-1 1]. 

344 Setting the initial approximations to dtype ``np.float32`` 

345 forces all iterative values to dtype ``np.float32`` speeding up 

346 the run while still allowing accurate eigenvalue computations. 

347 

348 >>> k = 1 

349 >>> rng = np.random.default_rng() 

350 >>> X = rng.normal(size=(n, k)) 

351 >>> X = X.astype(np.float32) 

352 

353 >>> eigenvalues, _ = lobpcg(A, X, maxiter=60) 

354 >>> eigenvalues 

355 array([100.]) 

356 >>> eigenvalues.dtype 

357 dtype('float32') 

358 

359 LOBPCG needs only access the matrix product with `A` rather 

360 then the matrix itself. Since the matrix `A` is diagonal in 

361 this example, one can write a function of the product 

362 ``A @ X`` using the diagonal values ``vals`` only, e.g., by 

363 element-wise multiplication with broadcasting 

364 

365 >>> A_f = lambda X: vals[:, np.newaxis] * X 

366 

367 and use the handle ``A_f`` to this callable function as an input 

368 

369 >>> eigenvalues, _ = lobpcg(A_f, X, maxiter=60) 

370 >>> eigenvalues 

371 array([100.]) 

372 

373 The next example illustrates computing 3 smallest eigenvalues of 

374 the same matrix given by the function handle ``A_f`` with 

375 constraints and preconditioning. 

376 

377 >>> k = 3 

378 >>> X = rng.normal(size=(n, k)) 

379 

380 Constraints - an optional input parameter is a 2D array comprising 

381 of column vectors that the eigenvectors must be orthogonal to 

382 

383 >>> Y = np.eye(n, 3) 

384 

385 The preconditioner acts as the inverse of `A` in this example, but 

386 in the reduced precision ``np.float32`` even though the initial `X` 

387 and thus all iterates and the output are in full ``np.float64``. 

388 

389 >>> inv_vals = 1./vals 

390 >>> inv_vals = inv_vals.astype(np.float32) 

391 >>> M = lambda X: inv_vals[:, np.newaxis] * X 

392 

393 Let us now solve the eigenvalue problem for the matrix `A` first 

394 without preconditioning requesting 80 iterations 

395 

396 >>> eigenvalues, _ = lobpcg(A_f, X, Y=Y, largest=False, maxiter=80) 

397 >>> eigenvalues 

398 array([4., 5., 6.]) 

399 >>> eigenvalues.dtype 

400 dtype('float64') 

401 

402 With preconditioning we need only 20 iterations from the same `X` 

403 

404 >>> eigenvalues, _ = lobpcg(A_f, X, Y=Y, M=M, largest=False, maxiter=20) 

405 >>> eigenvalues 

406 array([4., 5., 6.]) 

407 

408 Note that the vectors passed in `Y` are the eigenvectors of the 3 

409 smallest eigenvalues. The results returned above are orthogonal to those. 

410 

411 Finally, the primary matrix `A` may be indefinite, e.g., after shifting 

412 ``vals`` by 50 from 1, ..., 100 to -49, ..., 50, we still can compute 

413 the 3 smallest or largest eigenvalues. 

414 

415 >>> vals = vals - 50 

416 >>> X = rng.normal(size=(n, k)) 

417 >>> eigenvalues, _ = lobpcg(A_f, X, largest=False, maxiter=99) 

418 >>> eigenvalues 

419 array([-49., -48., -47.]) 

420 >>> eigenvalues, _ = lobpcg(A_f, X, largest=True, maxiter=99) 

421 >>> eigenvalues 

422 array([50., 49., 48.]) 

423 

424 """ 

425 blockVectorX = X 

426 bestblockVectorX = blockVectorX 

427 blockVectorY = Y 

428 residualTolerance = tol 

429 if maxiter is None: 

430 maxiter = 20 

431 

432 bestIterationNumber = maxiter 

433 

434 sizeY = 0 

435 if blockVectorY is not None: 

436 if len(blockVectorY.shape) != 2: 

437 warnings.warn( 

438 f"Expected rank-2 array for argument Y, instead got " 

439 f"{len(blockVectorY.shape)}, " 

440 f"so ignore it and use no constraints.", 

441 UserWarning, stacklevel=2 

442 ) 

443 blockVectorY = None 

444 else: 

445 sizeY = blockVectorY.shape[1] 

446 

447 # Block size. 

448 if blockVectorX is None: 

449 raise ValueError("The mandatory initial matrix X cannot be None") 

450 if len(blockVectorX.shape) != 2: 

451 raise ValueError("expected rank-2 array for argument X") 

452 

453 n, sizeX = blockVectorX.shape 

454 

455 # Data type of iterates, determined by X, must be inexact 

456 if not np.issubdtype(blockVectorX.dtype, np.inexact): 

457 warnings.warn( 

458 f"Data type for argument X is {blockVectorX.dtype}, " 

459 f"which is not inexact, so casted to np.float32.", 

460 UserWarning, stacklevel=2 

461 ) 

462 blockVectorX = np.asarray(blockVectorX, dtype=np.float32) 

463 

464 if retLambdaHistory: 

465 lambdaHistory = np.zeros((maxiter + 3, sizeX), 

466 dtype=blockVectorX.dtype) 

467 if retResidualNormsHistory: 

468 residualNormsHistory = np.zeros((maxiter + 3, sizeX), 

469 dtype=blockVectorX.dtype) 

470 

471 if verbosityLevel: 

472 aux = "Solving " 

473 if B is None: 

474 aux += "standard" 

475 else: 

476 aux += "generalized" 

477 aux += " eigenvalue problem with" 

478 if M is None: 

479 aux += "out" 

480 aux += " preconditioning\n\n" 

481 aux += "matrix size %d\n" % n 

482 aux += "block size %d\n\n" % sizeX 

483 if blockVectorY is None: 

484 aux += "No constraints\n\n" 

485 else: 

486 if sizeY > 1: 

487 aux += "%d constraints\n\n" % sizeY 

488 else: 

489 aux += "%d constraint\n\n" % sizeY 

490 print(aux) 

491 

492 if (n - sizeY) < (5 * sizeX): 

493 warnings.warn( 

494 f"The problem size {n} minus the constraints size {sizeY} " 

495 f"is too small relative to the block size {sizeX}. " 

496 f"Using a dense eigensolver instead of LOBPCG iterations." 

497 f"No output of the history of the iterations.", 

498 UserWarning, stacklevel=2 

499 ) 

500 

501 sizeX = min(sizeX, n) 

502 

503 if blockVectorY is not None: 

504 raise NotImplementedError( 

505 "The dense eigensolver does not support constraints." 

506 ) 

507 

508 # Define the closed range of indices of eigenvalues to return. 

509 if largest: 

510 eigvals = (n - sizeX, n - 1) 

511 else: 

512 eigvals = (0, sizeX - 1) 

513 

514 try: 

515 if isinstance(A, LinearOperator): 

516 A = A(np.eye(n, dtype=int)) 

517 elif callable(A): 

518 A = A(np.eye(n, dtype=int)) 

519 if A.shape != (n, n): 

520 raise ValueError( 

521 f"The shape {A.shape} of the primary matrix\n" 

522 f"defined by a callable object is wrong.\n" 

523 ) 

524 elif issparse(A): 

525 A = A.toarray() 

526 else: 

527 A = np.asarray(A) 

528 except Exception as e: 

529 raise Exception( 

530 f"Primary MatMul call failed with error\n" 

531 f"{e}\n") 

532 

533 if B is not None: 

534 try: 

535 if isinstance(B, LinearOperator): 

536 B = B(np.eye(n, dtype=int)) 

537 elif callable(B): 

538 B = B(np.eye(n, dtype=int)) 

539 if B.shape != (n, n): 

540 raise ValueError( 

541 f"The shape {B.shape} of the secondary matrix\n" 

542 f"defined by a callable object is wrong.\n" 

543 ) 

544 elif issparse(B): 

545 B = B.toarray() 

546 else: 

547 B = np.asarray(B) 

548 except Exception as e: 

549 raise Exception( 

550 f"Secondary MatMul call failed with error\n" 

551 f"{e}\n") 

552 

553 try: 

554 vals, vecs = eigh(A, 

555 B, 

556 subset_by_index=eigvals, 

557 check_finite=False) 

558 if largest: 

559 # Reverse order to be compatible with eigs() in 'LM' mode. 

560 vals = vals[::-1] 

561 vecs = vecs[:, ::-1] 

562 

563 return vals, vecs 

564 except Exception as e: 

565 raise Exception( 

566 f"Dense eigensolver failed with error\n" 

567 f"{e}\n" 

568 ) 

569 

570 if (residualTolerance is None) or (residualTolerance <= 0.0): 

571 residualTolerance = np.sqrt(np.finfo(blockVectorX.dtype).eps) * n 

572 

573 A = _makeMatMat(A) 

574 B = _makeMatMat(B) 

575 M = _makeMatMat(M) 

576 

577 # Apply constraints to X. 

578 if blockVectorY is not None: 

579 

580 if B is not None: 

581 blockVectorBY = B(blockVectorY) 

582 if blockVectorBY.shape != blockVectorY.shape: 

583 raise ValueError( 

584 f"The shape {blockVectorY.shape} " 

585 f"of the constraint not preserved\n" 

586 f"and changed to {blockVectorBY.shape} " 

587 f"after multiplying by the secondary matrix.\n" 

588 ) 

589 else: 

590 blockVectorBY = blockVectorY 

591 

592 # gramYBY is a dense array. 

593 gramYBY = blockVectorY.T.conj() @ blockVectorBY 

594 try: 

595 # gramYBY is a Cholesky factor from now on... 

596 gramYBY = cho_factor(gramYBY, overwrite_a=True) 

597 except LinAlgError as e: 

598 raise ValueError("Linearly dependent constraints") from e 

599 

600 _applyConstraints(blockVectorX, gramYBY, blockVectorBY, blockVectorY) 

601 

602 ## 

603 # B-orthonormalize X. 

604 blockVectorX, blockVectorBX, _ = _b_orthonormalize( 

605 B, blockVectorX, verbosityLevel=verbosityLevel) 

606 if blockVectorX is None: 

607 raise ValueError("Linearly dependent initial approximations") 

608 

609 ## 

610 # Compute the initial Ritz vectors: solve the eigenproblem. 

611 blockVectorAX = A(blockVectorX) 

612 if blockVectorAX.shape != blockVectorX.shape: 

613 raise ValueError( 

614 f"The shape {blockVectorX.shape} " 

615 f"of the initial approximations not preserved\n" 

616 f"and changed to {blockVectorAX.shape} " 

617 f"after multiplying by the primary matrix.\n" 

618 ) 

619 

620 gramXAX = blockVectorX.T.conj() @ blockVectorAX 

621 

622 _lambda, eigBlockVector = eigh(gramXAX, check_finite=False) 

623 ii = _get_indx(_lambda, sizeX, largest) 

624 _lambda = _lambda[ii] 

625 if retLambdaHistory: 

626 lambdaHistory[0, :] = _lambda 

627 

628 eigBlockVector = np.asarray(eigBlockVector[:, ii]) 

629 blockVectorX = _matmul_inplace( 

630 blockVectorX, eigBlockVector, 

631 verbosityLevel=verbosityLevel 

632 ) 

633 blockVectorAX = _matmul_inplace( 

634 blockVectorAX, eigBlockVector, 

635 verbosityLevel=verbosityLevel 

636 ) 

637 if B is not None: 

638 blockVectorBX = _matmul_inplace( 

639 blockVectorBX, eigBlockVector, 

640 verbosityLevel=verbosityLevel 

641 ) 

642 

643 ## 

644 # Active index set. 

645 activeMask = np.ones((sizeX,), dtype=bool) 

646 

647 ## 

648 # Main iteration loop. 

649 

650 blockVectorP = None # set during iteration 

651 blockVectorAP = None 

652 blockVectorBP = None 

653 

654 smallestResidualNorm = np.abs(np.finfo(blockVectorX.dtype).max) 

655 

656 iterationNumber = -1 

657 restart = True 

658 forcedRestart = False 

659 explicitGramFlag = False 

660 while iterationNumber < maxiter: 

661 iterationNumber += 1 

662 

663 if B is not None: 

664 aux = blockVectorBX * _lambda[np.newaxis, :] 

665 else: 

666 aux = blockVectorX * _lambda[np.newaxis, :] 

667 

668 blockVectorR = blockVectorAX - aux 

669 

670 aux = np.sum(blockVectorR.conj() * blockVectorR, 0) 

671 residualNorms = np.sqrt(np.abs(aux)) 

672 if retResidualNormsHistory: 

673 residualNormsHistory[iterationNumber, :] = residualNorms 

674 residualNorm = np.sum(np.abs(residualNorms)) / sizeX 

675 

676 if residualNorm < smallestResidualNorm: 

677 smallestResidualNorm = residualNorm 

678 bestIterationNumber = iterationNumber 

679 bestblockVectorX = blockVectorX 

680 elif residualNorm > 2**restartControl * smallestResidualNorm: 

681 forcedRestart = True 

682 blockVectorAX = A(blockVectorX) 

683 if blockVectorAX.shape != blockVectorX.shape: 

684 raise ValueError( 

685 f"The shape {blockVectorX.shape} " 

686 f"of the restarted iterate not preserved\n" 

687 f"and changed to {blockVectorAX.shape} " 

688 f"after multiplying by the primary matrix.\n" 

689 ) 

690 if B is not None: 

691 blockVectorBX = B(blockVectorX) 

692 if blockVectorBX.shape != blockVectorX.shape: 

693 raise ValueError( 

694 f"The shape {blockVectorX.shape} " 

695 f"of the restarted iterate not preserved\n" 

696 f"and changed to {blockVectorBX.shape} " 

697 f"after multiplying by the secondary matrix.\n" 

698 ) 

699 

700 ii = np.where(residualNorms > residualTolerance, True, False) 

701 activeMask = activeMask & ii 

702 currentBlockSize = activeMask.sum() 

703 

704 if verbosityLevel: 

705 print(f"iteration {iterationNumber}") 

706 print(f"current block size: {currentBlockSize}") 

707 print(f"eigenvalue(s):\n{_lambda}") 

708 print(f"residual norm(s):\n{residualNorms}") 

709 

710 if currentBlockSize == 0: 

711 break 

712 

713 activeBlockVectorR = _as2d(blockVectorR[:, activeMask]) 

714 

715 if iterationNumber > 0: 

716 activeBlockVectorP = _as2d(blockVectorP[:, activeMask]) 

717 activeBlockVectorAP = _as2d(blockVectorAP[:, activeMask]) 

718 if B is not None: 

719 activeBlockVectorBP = _as2d(blockVectorBP[:, activeMask]) 

720 

721 if M is not None: 

722 # Apply preconditioner T to the active residuals. 

723 activeBlockVectorR = M(activeBlockVectorR) 

724 

725 ## 

726 # Apply constraints to the preconditioned residuals. 

727 if blockVectorY is not None: 

728 _applyConstraints(activeBlockVectorR, 

729 gramYBY, 

730 blockVectorBY, 

731 blockVectorY) 

732 

733 ## 

734 # B-orthogonalize the preconditioned residuals to X. 

735 if B is not None: 

736 activeBlockVectorR = activeBlockVectorR - ( 

737 blockVectorX @ 

738 (blockVectorBX.T.conj() @ activeBlockVectorR) 

739 ) 

740 else: 

741 activeBlockVectorR = activeBlockVectorR - ( 

742 blockVectorX @ 

743 (blockVectorX.T.conj() @ activeBlockVectorR) 

744 ) 

745 

746 ## 

747 # B-orthonormalize the preconditioned residuals. 

748 aux = _b_orthonormalize( 

749 B, activeBlockVectorR, verbosityLevel=verbosityLevel) 

750 activeBlockVectorR, activeBlockVectorBR, _ = aux 

751 

752 if activeBlockVectorR is None: 

753 warnings.warn( 

754 f"Failed at iteration {iterationNumber} with accuracies " 

755 f"{residualNorms}\n not reaching the requested " 

756 f"tolerance {residualTolerance}.", 

757 UserWarning, stacklevel=2 

758 ) 

759 break 

760 activeBlockVectorAR = A(activeBlockVectorR) 

761 

762 if iterationNumber > 0: 

763 if B is not None: 

764 aux = _b_orthonormalize( 

765 B, activeBlockVectorP, activeBlockVectorBP, 

766 verbosityLevel=verbosityLevel 

767 ) 

768 activeBlockVectorP, activeBlockVectorBP, invR = aux 

769 else: 

770 aux = _b_orthonormalize(B, activeBlockVectorP, 

771 verbosityLevel=verbosityLevel) 

772 activeBlockVectorP, _, invR = aux 

773 # Function _b_orthonormalize returns None if Cholesky fails 

774 if activeBlockVectorP is not None: 

775 activeBlockVectorAP = _matmul_inplace( 

776 activeBlockVectorAP, invR, 

777 verbosityLevel=verbosityLevel 

778 ) 

779 restart = forcedRestart 

780 else: 

781 restart = True 

782 

783 ## 

784 # Perform the Rayleigh Ritz Procedure: 

785 # Compute symmetric Gram matrices: 

786 

787 if activeBlockVectorAR.dtype == "float32": 

788 myeps = 1 

789 else: 

790 myeps = np.sqrt(np.finfo(activeBlockVectorR.dtype).eps) 

791 

792 if residualNorms.max() > myeps and not explicitGramFlag: 

793 explicitGramFlag = False 

794 else: 

795 # Once explicitGramFlag, forever explicitGramFlag. 

796 explicitGramFlag = True 

797 

798 # Shared memory assingments to simplify the code 

799 if B is None: 

800 blockVectorBX = blockVectorX 

801 activeBlockVectorBR = activeBlockVectorR 

802 if not restart: 

803 activeBlockVectorBP = activeBlockVectorP 

804 

805 # Common submatrices: 

806 gramXAR = np.dot(blockVectorX.T.conj(), activeBlockVectorAR) 

807 gramRAR = np.dot(activeBlockVectorR.T.conj(), activeBlockVectorAR) 

808 

809 gramDtype = activeBlockVectorAR.dtype 

810 if explicitGramFlag: 

811 gramRAR = (gramRAR + gramRAR.T.conj()) / 2 

812 gramXAX = np.dot(blockVectorX.T.conj(), blockVectorAX) 

813 gramXAX = (gramXAX + gramXAX.T.conj()) / 2 

814 gramXBX = np.dot(blockVectorX.T.conj(), blockVectorBX) 

815 gramRBR = np.dot(activeBlockVectorR.T.conj(), activeBlockVectorBR) 

816 gramXBR = np.dot(blockVectorX.T.conj(), activeBlockVectorBR) 

817 else: 

818 gramXAX = np.diag(_lambda).astype(gramDtype) 

819 gramXBX = np.eye(sizeX, dtype=gramDtype) 

820 gramRBR = np.eye(currentBlockSize, dtype=gramDtype) 

821 gramXBR = np.zeros((sizeX, currentBlockSize), dtype=gramDtype) 

822 

823 if not restart: 

824 gramXAP = np.dot(blockVectorX.T.conj(), activeBlockVectorAP) 

825 gramRAP = np.dot(activeBlockVectorR.T.conj(), activeBlockVectorAP) 

826 gramPAP = np.dot(activeBlockVectorP.T.conj(), activeBlockVectorAP) 

827 gramXBP = np.dot(blockVectorX.T.conj(), activeBlockVectorBP) 

828 gramRBP = np.dot(activeBlockVectorR.T.conj(), activeBlockVectorBP) 

829 if explicitGramFlag: 

830 gramPAP = (gramPAP + gramPAP.T.conj()) / 2 

831 gramPBP = np.dot(activeBlockVectorP.T.conj(), 

832 activeBlockVectorBP) 

833 else: 

834 gramPBP = np.eye(currentBlockSize, dtype=gramDtype) 

835 

836 gramA = np.block( 

837 [ 

838 [gramXAX, gramXAR, gramXAP], 

839 [gramXAR.T.conj(), gramRAR, gramRAP], 

840 [gramXAP.T.conj(), gramRAP.T.conj(), gramPAP], 

841 ] 

842 ) 

843 gramB = np.block( 

844 [ 

845 [gramXBX, gramXBR, gramXBP], 

846 [gramXBR.T.conj(), gramRBR, gramRBP], 

847 [gramXBP.T.conj(), gramRBP.T.conj(), gramPBP], 

848 ] 

849 ) 

850 

851 _handle_gramA_gramB_verbosity(gramA, gramB, verbosityLevel) 

852 

853 try: 

854 _lambda, eigBlockVector = eigh(gramA, 

855 gramB, 

856 check_finite=False) 

857 except LinAlgError as e: 

858 # raise ValueError("eigh failed in lobpcg iterations") from e 

859 if verbosityLevel: 

860 warnings.warn( 

861 f"eigh failed at iteration {iterationNumber} \n" 

862 f"with error {e} causing a restart.\n", 

863 UserWarning, stacklevel=2 

864 ) 

865 # try again after dropping the direction vectors P from RR 

866 restart = True 

867 

868 if restart: 

869 gramA = np.block([[gramXAX, gramXAR], [gramXAR.T.conj(), gramRAR]]) 

870 gramB = np.block([[gramXBX, gramXBR], [gramXBR.T.conj(), gramRBR]]) 

871 

872 _handle_gramA_gramB_verbosity(gramA, gramB, verbosityLevel) 

873 

874 try: 

875 _lambda, eigBlockVector = eigh(gramA, 

876 gramB, 

877 check_finite=False) 

878 except LinAlgError as e: 

879 # raise ValueError("eigh failed in lobpcg iterations") from e 

880 warnings.warn( 

881 f"eigh failed at iteration {iterationNumber} with error\n" 

882 f"{e}\n", 

883 UserWarning, stacklevel=2 

884 ) 

885 break 

886 

887 ii = _get_indx(_lambda, sizeX, largest) 

888 _lambda = _lambda[ii] 

889 eigBlockVector = eigBlockVector[:, ii] 

890 if retLambdaHistory: 

891 lambdaHistory[iterationNumber + 1, :] = _lambda 

892 

893 # Compute Ritz vectors. 

894 if B is not None: 

895 if not restart: 

896 eigBlockVectorX = eigBlockVector[:sizeX] 

897 eigBlockVectorR = eigBlockVector[sizeX: 

898 sizeX + currentBlockSize] 

899 eigBlockVectorP = eigBlockVector[sizeX + currentBlockSize:] 

900 

901 pp = np.dot(activeBlockVectorR, eigBlockVectorR) 

902 pp += np.dot(activeBlockVectorP, eigBlockVectorP) 

903 

904 app = np.dot(activeBlockVectorAR, eigBlockVectorR) 

905 app += np.dot(activeBlockVectorAP, eigBlockVectorP) 

906 

907 bpp = np.dot(activeBlockVectorBR, eigBlockVectorR) 

908 bpp += np.dot(activeBlockVectorBP, eigBlockVectorP) 

909 else: 

910 eigBlockVectorX = eigBlockVector[:sizeX] 

911 eigBlockVectorR = eigBlockVector[sizeX:] 

912 

913 pp = np.dot(activeBlockVectorR, eigBlockVectorR) 

914 app = np.dot(activeBlockVectorAR, eigBlockVectorR) 

915 bpp = np.dot(activeBlockVectorBR, eigBlockVectorR) 

916 

917 blockVectorX = np.dot(blockVectorX, eigBlockVectorX) + pp