|
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 |
|
411 | 407 | resid_history[i + nlocked, niter+1] = norm(new_R[:, i]) |
412 | 408 | end |
413 | 409 | end |
414 | | - vprintln(niter, " ", resid_history[:, niter+1]) |
| 410 | + @debug niter, " ", resid_history[:, niter+1] |
415 | 411 |
|
416 | 412 | # it is actually a good question of knowing when to |
417 | 413 | # precondition. Here seems sensible, but it could plausibly be |
|
429 | 425 | for i=nlocked+1:M |
430 | 426 | if resid_history[i, niter+1] < tol |
431 | 427 | nlocked += 1 |
432 | | - vprintln("locked $nlocked") |
| 428 | + @debug "locked $nlocked" |
433 | 429 | else |
434 | 430 | # We lock in order, assuming that the lowest |
435 | 431 | # eigenvectors converge first; might be tricky otherwise |
|
439 | 435 | end |
440 | 436 |
|
441 | 437 | if display_progress |
442 | | - println("Iter $niter, converged $(nlocked)/$(n_conv_check), resid ", |
443 | | - norm(resid_history[1:n_conv_check, niter+1])) |
| 438 | + @info "Iter $niter, converged $(nlocked)/$(n_conv_check), resid " * |
| 439 | + "$(norm(resid_history[1:n_conv_check, niter+1]))" |
444 | 440 | end |
445 | 441 |
|
446 | 442 | if nlocked >= n_conv_check # Converged! |
|
0 commit comments