# Determining the algorithmic complexity

A few of the iterative matrix algorithms (CG,GMRES etc.) I have authored are acting rather funny. They converge to the right answers but take abnormally long time to run. I am in the process of finding out why.

One of the first steps I thought is to find out the algorithmic complexity of the methods. For instance, I need to know if CG is indeed taking $O(N^2)$ as expected and similar.

How do I find out the exact algorithmic complexity of the methods? (I am looking for some way for to get the exact bound on the algorithm like $N^2 + 5N$).

And even before that, is this a valid first step for improving performance?

As of now, I am concentrating on single core performance (which by itself is pretty bad right now).

I am also looking at trying to determine the memory accesses. Is there any Free software (not necessarily Open Source) which would allow me to do so in a nice GUI? I use VTune for parallel but it is useless for serial.

1. I have tried googling but all algorithm problems end up in the search/sort portion of computing which is quite different from iterative matrix algorithms.

2. I have tried solving the algorithm for matrix sizes from 1 to 1000, plot the graph of average time taken for that size and curve fit it (quadratic). But I don't seem to be getting anywhere through this.

EDIT: Just to be clear, I want to verify that the complexity of the algorithm in practice is the same as the one predicted by theory. I want to verify that my algorithm of say, Matrix-Vector is indeed $O(N^2)$.

Further, I don't mind knowing per iteration complexity.

-
What is preventing you from computing the complexity of each of your function calls, e.g., mycall(n) := ... for i<- 1 to n do myroutine ..., has n*my_routine, so $Complexity_{mycall}(n)=nComplexity_{myroutine}(n)$? –  Deathbreath Mar 1 '12 at 19:30
To add to Deathbreath's comment also pay attention to the space performance...once the active set grows beyond the size of the cache there are concerns related to efficient use of that newly scarce resource, and if you have to start swapping main memory... –  dmckee Mar 2 '12 at 0:02
I agree with Deathbreath... krylov subspace methods aren't necessarily guaranteed except for specific classes of matrices. They often depend upon spectral properties, preconditioners, and the convergence criteria itself. I don't think we can specify much more than the complexity of each iteration. –  Paul Mar 2 '12 at 5:12
Deathbreath, @Paul, I don't mind looking at each iteration complexity. But I am interested in computing this and not calculating by hand. I want to know if the complexity actually being used in the algorithm is the same as the one that is predicted by the calculations you suggested. –  Inquest Mar 2 '12 at 14:24
@Nunoxic: What happens to your N vs. time plot if you divide by the iterations? How many different starting points do you use? –  Deathbreath Mar 2 '12 at 20:13

You have the code so implement a global counter for the number of floating point operations and increment the counter on every function call. There is no way around counting operations. You can use a similar global counter to measure the number of memory accesses.

To get the performance estimates, measure the time for a single iteration to get the FLOP/TIME. Then you can compare FLOP/TIME against your processor peak performance. Of course, you can do the same for the total number of data movements (in Bytes).

A piece of software that has counting flops functionality is PETSc, see their documentation for functions PetscLogFlops and the PetscGetFlops. The documentation provides good example codes and link to source code.

For an instructive discussion on what it involves to get close to peak FLOPS performance on current CPUs, see this thread.

For a book reference you can look into chapter 3 of Introduction to Algorithms.

-
Could you please provide a reference (even a textbook would do) for the things you suggested? –  Inquest Mar 4 '12 at 15:15
@Nunoxic, I added some references. –  fcruz Mar 4 '12 at 18:47

From a computer science point of view, there is no obvious way to automatically infer a program time/resource complexity. This is a hard problem in the realm of static analysis, and is mainly theoretically difficult because of the undecidability of the halting problem.

You can do the complexity analysis by hand though if your code is of manageable size. Inspect all your iterations and subroutines bottom-up, and calculate the complexity just like algorithm analysis books would do with psuedocode. But you have to make sure that the runtime semantics of your program are preserved by the compiler/optimizer.

From a practical point of view, you can use a profiler (e.g. gprof and valgrind) to diagnose your code performance while/after running. It can tell you about the pressure points in your code, and where your program spends most of its runtime. You can start from this info and diagnose your implementation.

-

One can compute the complexity of an iteration for iterative sovlers (e.g. number of operations), but this might fluctuate because of compiler magic. However, counting gets a good impression of what happends, and you can do it by yourself (or have the code couting it).

A complexity analysis for the total run will be hard, because the cost measure is somehow arbitrary, and, much more important, the total run time for, e.g., conjugate gradient will depend on the spectral properties of your matrix.

-
Please see edit. I don't mind the per iteration count. –  Inquest Mar 2 '12 at 14:27

Start by looking at the best out there, see: BLAS/LAPACK/ARPACK/UMFPACK/Trilinos/PETSc (in rough order of level of sophistication) These all rely on highly processor-specific, cache-optimized operations, and hence are nearly as fast as can be expected. (These are all c/c++/fortran monsters and there is a reason for that!)

I assume this is related to your previous questions about home-made linear algebra solvers, and I want to reiterate that making your own implementation, while a valuable learning experience, will almost certainly be slower. In my experience, you are doing really well if you get within a factor (i.e. 1/2) of the performance of these codes!

Mostly what causes big performance hits after you have validated the correctness of the algorithim is accidental allocation/copies/deallocation of memory. A quick run of Valgrind might help with this one. Other than that, make sure you are not storing intermediate products or doing extra work!

Lastly, it is worth note that iterative methods act finnicky given certain inputs. Start by making sure your implementation behaves appropriately given certain sorts of known inputs that have well-defined execution pathways. An example might be solving a diagonal, or Laplace matrix where the intermediate products have an easily graph-interpertable form.

-
The algorithm uses BLAS subroutines but I am concerned about how they are being used by me. Further, my iterative methods work for every matrix (that it is predicted to work on). The problem isn't convergence, it is slow convergence. –  Inquest Mar 2 '12 at 14:28
How well does it perform on a diagonal matrix? –  meawoppl Mar 2 '12 at 23:04
Fairly quick convergence. A 3000x3000 matrix (dense and random) takes about 0.3 seconds to converge (On my laptop's 2.13GHz Intel i3 330M processor) –  Inquest Mar 3 '12 at 4:59