(dot-eigenvectors ...) changes eigenvectors when mu_inv != null

[thchr]

Calling (dot-eigenvectors ev band-start) for some set of eigenvectors ev appears to effect a change to the computed eigenvectors when mu_inv is not null (i.e. when μ ≠ 1, equivalently when B ≠ H).
Accordingly, after calling (dot-eigenvectors ...) in a μ ≠ 1 case, subsequent calls to e.g. (get-eigenvector ...) yield unpredictable results. This does not happen when μ = 1.

The issue is illustrated by the (unfortunately, somewhat lengthy) example below:

(use-modules (ice-9 format)) ; allow control of print output formatting

; ----- basic, shared properties of both calculations (square lattice)
(set! geometry-lattice (make lattice (size 1 1 no-size)
                         (basis1 1 0)
                         (basis2 0 1)))

(set! k-points (list (vector3 0.125 0.125 0))) ; arbitrary k-point

(set! resolution 32)
(set! num-bands 6)

; ----- mu = 1 calculation & overlaps
(define matprops (make dielectric (epsilon 13)))
(set! geometry (list (make cylinder (center 0 0 0)
                                    (radius 0.13)
                                    (height infinity)
                                    (material matprops))))

(run-tm)

(define overlaps1 (dot-eigenvectors (get-eigenvectors 1 num-bands) 1))
(define overlaps2 (dot-eigenvectors (get-eigenvectors 1 num-bands) 1))

; ----- mu \noteq 1 calculation & overlaps
(define matprops (make medium-anisotropic (epsilon-diag 13 13 13)
                                          (epsilon-offdiag 0 0 0)
                                          (mu-diag 1 1 1)
                                          (mu-offdiag 0+0.40i 0 0)))
(set! geometry (list (make cylinder (center 0 0 0)
                                    (radius 0.13)
                                    (height infinity)
                                    (material matprops))))
(run-tm)

(define overlaps-mu1 (dot-eigenvectors (get-eigenvectors 1 num-bands) 1))
(define overlaps-mu2 (dot-eigenvectors (get-eigenvectors 1 num-bands) 1))


; ----- function for printing magnitude of sqmatrix object neatly
(define (print-sqmatrix-magnitude matrix size)
        (map (lambda (i1)
                (map (lambda (i2)
                        (print (format #f "~10,2,2g" (magnitude (sqmatrix-ref matrix i1 i2)) ", "))
                     ) (iota size))
                (print "\n")
             ) (iota size))
)
; ----- test orthonormalization, i.e. that <Bi|Hj>=delta_ij
(print "\n--- In the case mu = 1, dot-eigenvectors works perfectly ---\n")
(print "1st call:\n") (print-sqmatrix-magnitude overlaps1 num-bands)
(print "2nd call:\n") (print-sqmatrix-magnitude overlaps2 num-bands) (print "\n")

(print "\n--- In the case mu \\noteq 1, dot-eigenvectors CHANGES the  ---\n"
         "--- eigenvectors - subsequent calls to dot-eigenvectors do ---\n"
         "--- not yield the same results:                            ---\n")
(print "1st call:\n") (print-sqmatrix-magnitude overlaps-mu1 num-bands)
(print "2nd call:\n") (print-sqmatrix-magnitude overlaps-mu2 num-bands)

which yields the results (which are wrong only in the second call to (dot-eigenvectors ...) in the μ ≠ 1 case)

--- In the case mu = 1, dot-eigenvectors works perfectly ---
1st call:
  1.00      2.04E-14  2.29E-14  2.06E-14  7.40E-15  5.06E-15
  2.04E-14  1.00      3.01E-14  2.94E-14  1.82E-14  8.14E-14
  2.29E-14  3.01E-14  1.00      1.16E-14  1.64E-14  6.45E-14
  2.06E-14  2.94E-14  1.16E-14  1.00      2.48E-14  4.33E-14
  7.40E-15  1.82E-14  1.64E-14  2.48E-14  1.00      5.70E-14
  5.06E-15  8.14E-14  6.45E-14  4.33E-14  5.70E-14  1.00
2nd call:
  1.00      2.04E-14  2.29E-14  2.06E-14  7.40E-15  5.06E-15
  2.04E-14  1.00      3.01E-14  2.94E-14  1.82E-14  8.14E-14
  2.29E-14  3.01E-14  1.00      1.16E-14  1.64E-14  6.45E-14
  2.06E-14  2.94E-14  1.16E-14  1.00      2.48E-14  4.33E-14
  7.40E-15  1.82E-14  1.64E-14  2.48E-14  1.00      5.70E-14
  5.06E-15  8.14E-14  6.45E-14  4.33E-14  5.70E-14  1.00


--- In the case mu \noteq 1, dot-eigenvectors CHANGES the  ---
--- eigenvectors - subsequent calls to dot-eigenvectors do ---
--- not yield the same results:                            ---
1st call:
  1.00      9.66E-14  5.49E-14  5.42E-14  8.08E-14  7.39E-14
  9.67E-14  1.00      4.40E-13  1.33E-13  4.62E-13  2.18E-13
  5.49E-14  4.40E-13  1.00      5.60E-13  5.20E-13  4.44E-13
  5.42E-14  1.33E-13  5.60E-13   1.0      6.79E-13  5.78E-13
  8.08E-14  4.62E-13  5.20E-13  6.79E-13   1.0      2.95E-13
  7.39E-14  2.18E-13  4.44E-13  5.78E-13  2.95E-13   1.0
2nd call:                                                      <<< (issue here)
   1.0      9.36E-03  6.81E-02  9.62E-02  5.95E-02  8.86E-03   <<< (issue here)
  9.36E-03   1.1      6.49E-02  4.65E-02  3.46E-02  5.02E-02   <<< (issue here) 
  6.81E-02  6.49E-02  0.89      0.10      2.74E-02  3.31E-02   <<< (issue here)
  9.62E-02  4.65E-02  0.10       1.4      0.22      0.11       <<< (issue here)
  5.95E-02  3.46E-02  2.74E-02  0.22       1.1      4.07E-02   <<< (issue here)
  8.86E-03  5.02E-02  3.31E-02  0.11      4.07E-02   1.1       <<< (issue here)

Another consequence is that calling e.g. (get-bfield ...) or (get-hfield ...) after calling (dot-eigenvectors ...) will differ from calling them before (dot-eigenvectors ...).

I imagine the fix must lie somewhere in lines 293-319 of matrix-smob.c but, after much staring, I haven't been able to see what would cause this.

(In principle, I could make do with (integrate-fields ...), (get-hfield ...), & (get-bfield ...) combinations, but they are quite a bit slower)

Thanks!

[thchr]

I've had a chance to look at this some more.

It appears that Hblock points to H, or at least that assigning to Hblock.data effects the same change to H.data. This was somewhat surprising to me. Is this intended?
This behavior means that calls to e.g. maxwell_compute_H_from_B(mdata, H, Hblock, ...) will change H (this is the reason for the issue above). Thus, trying to evaluate anything which depends on H after such a call will give inconsistent results.

Working with the W[0] work matrix instead of Hblock in dot_eigenvectors(...) sidesteps this issue, in principle - but that then ignores the comments in the code about W[0] then not necessarily having sufficient space.

[stevengj]

If you have ≤ 11 bands (the eigensolver block size), then Hblock == H. See these lines.

W[0] should have the same size as Hblock, so it is the right thing. The comment you are referring to says that you can't dot all of the bands at once, you have to do them in blocks. W[0] has enough space for these blocks, though.

[stevengj]

I don't know why I used Hblock there and not W[0] … it would be great if you submitted a pull request to change this to W[0].

[thchr]

Thanks for explaining that. I checked for >11 bands, and the problem indeed goes away.

I'd be happy to submit a pull request. Before doing so, I was hoping to quickly clarify one thing:
... to avoid trouble in the case that I'm requesting a large number of bands simultaneously be converted from B to H --- would an "ugly" splitting condition, based on the number of bands requested, be acceptable? (i.e. use W[0] if num-bands ≤ 11 and Hblock otherwise?)

[thchr]

(or, I guess, more generally, if num_bands ≤ block_size)

[stevengj]

No, just use W[0] in all cases. It is the same size as Hblock.

[thchr]

Ah, I see. Okay, thanks.