Sum along diagonals

[mole6e23(at)hotmail.com (Todd Clements)]

Hi there..

I thought I'd ask this question, more out of curiosity as to a possible
result than for really needing to optimize it (this code is used
approximately three times for every blue moon). I'm always impressed with
some of the results from this group, most recently:

a=a[*,where(histogram(b,MIN=0,MAX=(size(a,/DIMENSIONS))[1]-1,BINSIZE=1) eq 0)]

I wouldn't have thought of compacting all that into one line!

Every once in a while (not often enough to make me worry about optimizing
too much), I want to take a not necessarily square matrix and get the sum
along the diagonals, such as the following, with the theoretical function
sum_diag:

IDL> blah = indgen( 4, 4 )
IDL> print, blah
       0       1       2       3
       4       5       6       7
       8       9      10      11
      12      13      14      15
IDL> print, sum_diag( blah )
       0       5      15      30      30      25      15

which is the series [0, 4+1, 8+5+2, 12+9+6+3, ... ]

Of course, to be difficult, I'd like it to work for non-square matrices as well:

IDL> blah = indgen( 5, 3 )
IDL> print,blah
       0       1       2       3       4
       5       6       7       8       9
      10      11      12      13      14

and the result would be the series [0, 5+1, 10+6+2, 11+7+3, ... ]

The best I've found is to count the diagonals, loop through those, and
then loop again through the possible row and column indices in that array.
It's pretty ugly, and very time consuming for a 512x512 array (most of
what I want). It just occured to me now that for the squre matrix, you
could take "blah + transpose( blah )" and now you only have to worry about
the lower half of the matrix (and the doubled diagonal elements), so it
might be a bit better, but still messy.

So if you're bored, I'd love to see a better solution than a bunch of
long, nested loops, even if just only for the square case (since that is
90% of the cases in my data)

Enjoy!
Todd