Some doubts about types

[stevengj]

@baptiste

My biggest doubt about this approach though, is that ultimately A appears in a 3Nx3N linear system A x = b (coupled-dipole equations), and if A is a NxN matrix of SMatrices I don’t really see how it’s going to work (ultimately some Lapack routine will need the 3Nx3N matrix, won’t it?).

You can do Ajk = ...some SMatrix.. and then do A[3j-2:3j, 3k-3:3k] = Ajk to write the Smatrix into the “flat” array A in-place.

@baptiste

α_jj = transpose(Rm) * Diagonal(Alpha[jj]) * Rm

You can do this by dot(Rm, Alpha[jj] .* Rm) or Rm' * (Alpha[jj] .* Rm)without constructing the Diagonal matrix.

Note also that you can just do foo' instead of transpose(foo) — the two are equivalent for real vectors like your coordinate vectors here.

@baptiste

Ajk = SMatrix(...)

If you’ve done things correctly, you shouldn’t need the SMatrix(...) constructor here… the matrix expression should already give an SMatrix (because it’s built out of SMatrix expressions like rjkhat * rjkhat').

(You can add a type assertion ::SMatrix at the end of the expression if you want to check that you did this correctly. A type assertion will throw an error if the expression is not the expected type, but is eliminated at compile-time if the expression has the correct type and type-inference succeeds.)