3

Given a complex square matrix G and a square matrix M, which I can compute quickly from G, I need to calculate efficiently the matrix G' defined by the matrix equation G' ⋅ M⊤ = G in my C++ program. Afterwards the original G can be discarded.

Since my experience with numerical linear algebra is limited I have so far relied on the Armadillo library, which has very nice syntax and can provide good performance. My approach here would then be to take the transpose of both sides of the equation and solve the problem M ⋅ (G')⊤ = G⊤ by calling 

using namespace arma;
G = trans(solve(M, trans(G));

But if I have understood the heavily template-based code correctly, this would involve actually performing the transpositions and copying data around the LAPACK routine cgesv. Of course these are O(N2) operations compared to the O(N3) of the actual linear solver, but I'd still rather avoid them in this case.

I have never used LAPACK directly up to now, but would there be the possibility to gain performance by directly calling a LAPACK routine without having to actually carry out the transpositions? Or would you suggest a different approach?

rerx
  • 1,133
  • 8
  • 19
  • How big are those matrices? If you're talking about something larger than 10x10, the transposition will be negligible compared to linear solver. – Piotr Kołaczkowski Sep 25 '13 at 18:32
  • There is no good reason to resort to direct LAPACK calls for this. You said it yourself, the transposition is extremely cheap compared to the solver, so you're getting a performance bonus that may be too small to measure in return for making your code a lot less readable. Using LAPACK directly from C/C++ is rather messy. – us2012 Sep 25 '13 at 18:47
  • 1
    Also take a look at the Eigen library: http://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html I *think* they use some template magic that means the transposition is never actually computed, it just provides different accessor routines to the next operation (some kind of decomposition for the solver). I don't know much about armadillo, but maybe they do that as well? – us2012 Sep 25 '13 at 18:53
  • The matrices would be between 64x64 and 1024x1024 for now, potentially bigger later, if I manage to further optimize some other parts of the algorithms. You are right for sure that performance gains would be marginal at this point. This is more a matter of principle (and maybe wasted energy...): It just feels so _unnecessary_ to do the transpositions. – rerx Sep 25 '13 at 19:11
  • @us2012 Eigen is definitely a promising library as well. For now I have relied on Armadillo mostly because the major parts of the computations can be done by the Intel MKL, which provides some degree of parallelization and _should_ be the best-optimized code for the cluster of Xeon CPUs my programs run on. – rerx Sep 25 '13 at 19:16
  • 2
    A matter of what principle exactly? Premature optimization? I'm only partly joking ;) . There's something else that catches my eye though, you say that `M` can be computed quickly from `G`. Could this not lead to an easier solution of your system? – us2012 Sep 25 '13 at 20:11
  • With "quickly" I meant O(N^2). I have spent some thought on whether there is an easier solution, but could not find anything better than O(N^3). This might be another question, though, maybe on the mathematics sister site.. – rerx Sep 25 '13 at 20:59

1 Answers1

0

I have once provided a C++ code for lapack. Take a look at my answer which provides an example on how to use lapack from the original Fortran library.

Understanding LAPACK calls in C++ with a simple example

Community
  • 1
  • 1
The Quantum Physicist
  • 24,987
  • 19
  • 103
  • 189
  • While this is related, I was rather looking for lapack advise for this specific problem – rerx Sep 26 '13 at 07:27
  • In that link you have a manual and you have an example on how to solve general matrices. With all due respect, but what do you need exactly? Should I solve it for you? What do you need more than an example with different dimensionality? I'm very puzzled by your answer. – The Quantum Physicist Sep 27 '13 at 08:49