|
39 | 39 | ## TODO it seems there is a lack of orthogonalization immediately after locking, maybe investigate this to save on some choleskys |
40 | 40 | ## TODO debug orthogonalizations when A=I |
41 | 41 |
|
42 | | -# TODO Use @debug for this |
43 | | -# vprintln(args...) = println(args...) # Uncomment for output |
44 | | -vprintln(args...) = nothing |
45 | | - |
46 | 42 | using LinearAlgebra |
47 | 43 | import LinearAlgebra: BlasFloat |
48 | 44 | import Base: * |
@@ -180,7 +176,7 @@ normest(M) = maximum(abs.(diag(M))) + norm(M - Diagonal(diag(M))) |
180 | 176 | success = true |
181 | 177 | catch err |
182 | 178 | @assert isa(err, PosDefException) |
183 | | - vprintln("fail") |
| 179 | + @debug "fail" |
184 | 180 | # see https://arxiv.org/pdf/1809.11085.pdf for a nice analysis |
185 | 181 | # We are not being very clever here; but this should very rarely happen |
186 | 182 | # so it should be OK |
@@ -220,7 +216,7 @@ normest(M) = maximum(abs.(diag(M))) + norm(M - Diagonal(diag(M))) |
220 | 216 | # condR = 1/LAPACK.trcon!('I', 'U', 'N', Array(R)) |
221 | 217 | condR = normest(R)*norminvR # in practice this seems to be an OK estimate |
222 | 218 |
|
223 | | - vprintln("Ortho(X) success? $success ", eps(real(T))*condR^2, " < $tol") |
| 219 | + @debug "Ortho(X) success? $success ", eps(real(T))*condR^2, " < $tol" |
224 | 220 |
|
225 | 221 | # a good a posteriori error is that X'X - I is eps()*κ(R)^2; |
226 | 222 | # in practice this seems to be sometimes very overconservative |
|
289 | 285 | niter > 10 && error("Ortho(X,Y) is failing badly, this should never happen") |
290 | 286 | niter += 1 |
291 | 287 | end |
292 | | - vprintln("ortho choleskys: ", ninners) # get how many Choleskys are performed |
| 288 | + @debug "ortho choleskys: ", ninners # get how many Choleskys are performed |
293 | 289 |
|
294 | 290 | # @assert (norm(BY'X)) < tol |
295 | 291 | # @assert (norm(X'X-I)) < tol |
|
400 | 396 | resid_history[i + nlocked, niter+1] = norm(new_R[:, i]) |
401 | 397 | end |
402 | 398 | end |
403 | | - vprintln(niter, " ", resid_history[:, niter+1]) |
| 399 | + @debug niter, " ", resid_history[:, niter+1] |
404 | 400 |
|
405 | 401 | # it is actually a good question of knowing when to |
406 | 402 | # precondition. Here seems sensible, but it could plausibly be |
|
418 | 414 | for i=nlocked+1:M |
419 | 415 | if resid_history[i, niter+1] < tol |
420 | 416 | nlocked += 1 |
421 | | - vprintln("locked $nlocked") |
| 417 | + @debug "locked $nlocked" |
422 | 418 | else |
423 | 419 | # We lock in order, assuming that the lowest |
424 | 420 | # eigenvectors converge first; might be tricky otherwise |
|
428 | 424 | end |
429 | 425 |
|
430 | 426 | if display_progress |
431 | | - println("Iter $niter, converged $(nlocked)/$(n_conv_check), resid ", |
432 | | - norm(resid_history[1:n_conv_check, niter+1])) |
| 427 | + @info "Iter $niter, converged $(nlocked)/$(n_conv_check), resid " * |
| 428 | + "$(norm(resid_history[1:n_conv_check, niter+1]))" |
433 | 429 | end |
434 | 430 |
|
435 | 431 | if nlocked >= n_conv_check # Converged! |
|
0 commit comments