New version of Fractalworks w full mulitprocessor support

[Duncan C]

Animation, you say?  Bring it on.  I'd love something new and fast for parameter-interpolated zoom animations on OSX.  Have you used Graham Anderson's Escape or Jesse Jones' Mandella?

I haven't tried either. Clearly I need to get my head out of my code long enough to survey the other apps that are out there. Can you point me to links to these two?


Duncan C

[Duncan C]

Properly done, these can easily accelerate your code by a factor of 8 or more.  Note that it's hard to get really good speed improvements on SSE as easily as you can with AltiVec, due to the limited number of vector registers (until you go to 64-bit x86 code, which you won't be able to reasonably do until Leopard ships), but depending on the image, it's not that hard to get 3x+ due to AltiVec (quadrium often gets 3.5 times the performance with AltiVec enabled, but only about 2x for SSE)

i got nearly a clean 4x speedup from using sse, and that's on the simple mandelbrot/julia - not something more compute-limited like quadrium; what gives?

Threading is also a must - basically split the image into parts, and have each part calculated by a different thread.  It's not all that hard to add in once you've got the basics of threading (and this is a fairly simple threading approach).  You pretty much get N times the performance then (where N is the number of cores your machine has).

that's of course assuming you have an equal workload for all pixels, otherwise you'll have threads twiddling thumbs after a while. to be fair you did mention that uniform subdivision is a simple approach, so this is obviously @ duncan and not at you.

The way I deal with that is to break the plot into a fairly large number of pieces, and queue the pieces up for computation. I break the plot into at least (num_processors x 2) pieces. I spawn as many compute threads as I have processors, and when a thread finishes, I spawn another one for the next piece of the plot. I'm still fine-tuning the optimum number of pieces. I do sometimes end up with one thread working alone at the end if it contains a large number of pixels at max iterations.

For brute-force computation, it would be pretty easy to break the plot into pices that have roughly the same number of iterations, by sub-sampling the region and estimating the total number of iterations required.

However that approach breaks down for a boundary following algorithm (see below), since contiguous regions with the same iteration value would only require that the perimeter pixels be computed. I haven't yet figured out how I'll break up plots using that approach. I might just go for equal area and call that "close enough."

My program already builds a histogram for a plot, tracking the number of pixels at each iteration value. I use the histogram in building a color table to display the plot.  I use a large color change between iteration values with lots of pixels, and small color change between iteration values with a small number of pixels. This makes color tables that adapt nicely to the complexity of the plot. Once I implement boundary following, I may add a second histogram that tracks "average number of contiuous pixels" for each iteration value. That should give me a better metric for building color tables than my current histogram does.

gcc offers two different complex number packages - one from the C99 (?) complex data type (which offers an intrinsic complex data type that allows to compiler to do things like basic operations), and another that involves the <complex> C++ template based library.  Neither are accelerated to support AltiVec/PPC (but due to that way that C++ works, you can extend that one via template specialization, but it's a fair amount of work).

i'm pretty sure core 2 duo has specific instructions for accelerating complex number computations (finally!), so it's reasonable to assume that a recent gcc will spit those instructions if compiled with the right march-flag (microarchitecture). of course that doesn't help our friend duncan, but i'd be interested to know if that's indeed the case.

finally, there are ways to speed up general fractal computations without resorting to special cases for various fractal types, but in the end your time will be more rewarded by just biting the bullet and going for the big speedup. put another way, if you ever want to catch up with quadrium, you'll have to get your hands dirty

I'm planning to tackle a boundary following algorithm for Mandelbrot and Julia sets. This produces HUGE speed improvements for all but the most convoluted of plots, since each portion of an iteration band usually has an area that's large compaired to it's perimeter. For plots where the iteration values change too rapidly to have regions with significant interiors, there's still a large payoff in using boundary following to identify Mandelbrot and Julia regions, since those are the most expensive of all to calculate.

There's a Java app I found that blows the doors off anything else I've seen, using double precision math in Java, with no machine-specific optimization. It is able to render several frames a second at 1000 iterations on my single processor G4. See this link: http://www.ibiblio.org/e-notes/MSet/Anim/ManJuOrb.htm. It's very impressive. The author even posted the source code, which should help me quite a bit.


Duncan C

[lycium]

I haven't tried either. Clearly I need to get my head out of my code long enough to survey the other apps that are out there. Can you point me to links to these two?

lazy! 10 seconds of googling reveals:

http://asia.cnet.com/downloads/mac/swinfo/0,39001908,20000561r-39035314s,00.htm
http://xahlee.org/PageTwo_dir/MathPrograms_dir/fractals.html

now to answer your other questions...

The way I deal with that is to break the plot into a fairly large number of pieces, and queue the pieces up for computation. I break the plot into at least (num_processors x 2) pieces. I spawn as many compute threads as I have processors, and when a thread finishes, I spawn another one for the next piece of the plot. I'm still fine-tuning the optimum number of pieces. I do sometimes end up with one thread working alone at the end if it contains a large number of pixels at max iterations.

more pieces == more overhead, it's not about that at all. it's about load balancing, which is best solved by a thread pool. see recent discussions between trifox and i about mutatorwhatitsthingy's multithreading on these forums.

For brute-force computation, it would be pretty easy to break the plot into pices that have roughly the same number of iterations, by sub-sampling the region and estimating the total number of iterations required.

that's not easy at all, and demands a priori knowledge of your fractal iteration. i say again: mandelbrots and julias belong in the 80s, you're fighting a losing battle trying to optimise that specific problem more than others have - rather make ANY fractal render more efficiently, and render something new!

I'm planning to tackle a boundary following algorithm for Mandelbrot and Julia sets. This produces HUGE speed improvements for all but the most convoluted of plots, since each portion of an iteration band usually has an area that's large compaired to it's perimeter. For plots where the iteration values change too rapidly to have regions with significant interiors, there's still a large payoff in using boundary following to identify Mandelbrot and Julia regions, since those are the most expensive of all to calculate.

let's look ahead a little bit...

say you get your boundary following algo working, massively outrageously optimised to the point where you can do 1920x1200 images at 120fps like my graphics card can. you're still rendering mandelbrots while everyone else is doing the cool stuff!!

i'd really like to know why you're focusing all this attention on such a ridiculously simple problem, with no intention of walking your own path

[Duncan C]

I haven't tried either. Clearly I need to get my head out of my code long enough to survey the other apps that are out there. Can you point me to links to these two?

lazy! 10 seconds of googling reveals:

http://asia.cnet.com/downloads/mac/swinfo/0,39001908,20000561r-39035314s,00.htm
http://xahlee.org/PageTwo_dir/MathPrograms_dir/fractals.html

Guilty as charged. Thanks for the links.

now to answer your other questions...

The way I deal with that is to break the plot into a fairly large number of pieces, and queue the pieces up for computation. I break the plot into at least (num_processors x 2) pieces. I spawn as many compute threads as I have processors, and when a thread finishes, I spawn another one for the next piece of the plot. I'm still fine-tuning the optimum number of pieces. I do sometimes end up with one thread working alone at the end if it contains a large number of pixels at max iterations.

more pieces == more overhead, it's not about that at all. it's about load balancing, which is best solved by a thread pool. see recent discussions between trifox and i about mutatorwhatitsthingy's multithreading on these forums.

Understood. That's what I was saying. There's an optimum number of pieces that gets the best balance between overhead and load balancing.

For brute-force computation, it would be pretty easy to break the plot into pices that have roughly the same number of iterations, by sub-sampling the region and estimating the total number of iterations required.

that's not easy at all, and demands a priori knowledge of your fractal iteration.

I'm surpised you say that. I haven't tackled it yet, but I would think sampling every n'th pixel in a plot and iterating it would give a good estimate of the "iteration depth" of the plot. If I wrote the iteration values into my plot data then the computation time would not be wasted.

i say again: mandelbrots and julias belong in the 80s, you're fighting a losing battle trying to optimise that specific problem more than others have - rather make ANY fractal render more efficiently, and render something new!

I'm planning to tackle a boundary following algorithm for Mandelbrot and Julia sets. This produces HUGE speed improvements for all but the most convoluted of plots, since each portion of an iteration band usually has an area that's large compaired to it's perimeter. For plots where the iteration values change too rapidly to have regions with significant interiors, there's still a large payoff in using boundary following to identify Mandelbrot and Julia regions, since those are the most expensive of all to calculate.

let's look ahead a little bit...

say you get your boundary following algo working, massively outrageously optimised to the point where you can do 1920x1200 images at 120fps like my graphics card can.

Are you saying your GRAPHICS CARD can do high iteration plots at 1920x1200 at 120 fps, or that it can display such plots at that rate?
I've heard about people writing code to use the multiprocessor ability of their graphics cards to do fractal rendering, but 1920x1200 sounds a bit over the top.

you're still rendering mandelbrots while everyone else is doing the cool stuff!!

i'd really like to know why you're focusing all this attention on such a ridiculously simple problem, with no intention of walking your own path

I'm just getting started. I'm tackling Mandelbrot and Julia sets first, to get warmed up. I also think there's a lot of potential for original work there, even still. It's a very rich subject. I have significant plans for animated "movies," for example. When I exhaust that, i'll move on to other formulas. Many of the tools and approaches I'm using will be quite applicable to other types of iterative fractals.

[lycium]

Are you saying your GRAPHICS CARD can do high iteration plots at 1920x1200 at 120 fps, or that it can display such plots at that rate?
I've heard about people writing code to use the multiprocessor ability of their graphics cards to do fractal rendering, but 1920x1200 sounds a bit over the top.

there's a screenshot in an older thread of mine: http://www.fractalforums.com/index.php?topic=386.0

yeah it's "only" single precision floating point (also used for colour-computations), but double precision can easily be used too. the take-away point is that no one cares because mandelbrots and julias have been done TO DEATH, and that application was written in a few seconds probably just to illustrate the computational power of modern gpus.

I'm just getting started. I'm tackling Mandelbrot and Julia sets first, to get warmed up. I also think there's a lot of potential for original work there, even still. It's a very rich subject. I have significant plans for animated "movies," for example. When I exhaust that, i'll move on to other formulas. Many of the tools and approaches I'm using will be quite applicable to other types of iterative fractals.

the "one step at a time" point is fair enough, but you didn't really answer my question: doing a single specific kind of fractal fast is sidestepping the greater goal of writing a good, generic fractal renderer - that code simply won't be used because as soon as you start exploring interesting other kinds of fractals 9/10th of your program is useless! edge tracing, symmetry exploitation, none of those techniques will help when you have fully flexible fractal iteration.

an excellent point of comparison at this stage is quadrium. he constructs his fractal iterations via a directed graph, entered via a gui, and has to iterate THAT fast! that's what i call a good learning experience, and certainly his results speak volumes: http://gandreas.deviantart.com/


edit: oh and you misunderstood what i said about load balancing; splitting the image into smaller pieces alone is a very poor way to distribute your computation. once you've split the image, you should dole out work to the worker threads as soon as they are done with their previous task, in a dynamic fashion. a simple and near-perfect alternative that only works in this instance, is to just spawn two threads to do the even/odd scanlines (with obvious generalisation for N cores).

[Duncan C]

snip...

 oh and you misunderstood what i said about load balancing; splitting the image into smaller pieces alone is a very poor way to distribute your computation. once you've split the image, you should dole out work to the worker threads as soon as they are done with their previous task, in a dynamic fashion. a simple and near-perfect alternative that only works in this instance, is to just spawn two threads to do the even/odd scanlines (with obvious generalisation for N cores).

I think I understood you.

I divide my plot into pieces, and create a single queue of "render tasks", each of which represents a piece of roughly the same size. (Much like a single line at a bank feeds customers to the tellers.)

I then start a rendering process where I take render tasks off my task queue and hand them to worker threads. My code keeps as many worker threads running as there are processors. When a thread completes, I take another render task off the queue and spawn another worker thread to render it. As a result, all the processors on the system stay busy until the queue of render tasks has run out. The worst case for my code is when the last task in the queue takes significantly longer than any other. In that case, the queue is empty and the other processors end up idle.

It would be better if each task required roughly the same time to complete, but I am going to re-factor my rendering soon in a way such that estimating the compute time will be pretty difficult. (boundary following.) I expect the variation in compute times between tasks to drop quite a bit, however, since large areas at max. iterations are the largest reason for slow rendering areas, and boundary following greatly reduces the number of points in a region that actually have to be computed.

I should probably do a little extra work and have my worker threads go into an idle state between tasks rather than shutting down. That would save the overhead of spawning a new thread for each task. My implementation was definitely "quick and dirty." Is that what you were referring to? Do you think the overhead of spawning a worker thread for each task has a major impact on efficiency?


Duncan C

[lycium]

clarified for sure, misunderstanding was on my part; yeah there is considerable overhead in spawning/terminating threads compared to just having a pool of them and re-using. it depends on your application though, like if you use lots of supersampling i doubt you can measure the difference (unless you're creating really a ton of threads); for realtime purposes (like generating fractals at 120 fps ) you'll want the pool.

still, in this simple instance where your workload is still relatively consistent, i think the best (measured by results/effort) approach would just be to spawn as many threads as you have processors*cores, and let them process the scanlines in an interleaved fashion.

[Duncan C]

clarified for sure, misunderstanding was on my part; yeah there is considerable overhead in spawning/terminating threads compared to just having a pool of them and re-using. it depends on your application though, like if you use lots of supersampling i doubt you can measure the difference (unless you're creating really a ton of threads); for realtime purposes (like generating fractals at 120 fps ) you'll want the pool.

still, in this simple instance where your workload is still relatively consistent, i think the best (measured by results/effort) approach would just be to spawn as many threads as you have processors*cores, and let them process the scanlines in an interleaved fashion.

I create less than 10 threads for a typical plot. I'll have to look at refactoring my code so that I keep the worker threads around until my task pool is empty. Given the small number of threads, however, I suspect the speed improvement will be very small.

Doing alternating scan-lines is not a good choice for me because my next step is the boundary following algorithm.

My goal is multiple frames per second on an average machine for a typical plot. I want my program to be able to do animated "walks" through the mandelbrot set, and through Julia sets. That's the real "meat" of my application, and where I feel there's room for innovation. Julia sets transform in fascinating ways as the seed point moves from one place to another, and the animations can be quite striking. If I can get the rendering speed fast enough, I'll be able to generate these walks "on the fly."

I've discovered through investigation that Julia sets have different structure on several different scales (at least 3) based on the location of their source point. The structure of the macro scale (width of 3 - 4) is dictated by the Julia's source-point on the macro scale (e.g. in the "cleft" of the Cardioid, in "seahorse valley", along the central spike at i=0, etc. )If you zoom in on the julia, at a width of .2 or so, the detail at that scale takes the structure of the second level detail from the Mandelbrot set. If you zoom in further on the Julia plot, it's third scale of detail (at a width of around .02 or so) takes the form of the next level of detail from the source point, etc.

I'm planning on creating animations that show a Mandelbrot plot at several scales, or levels of magnification, and corrisponding Julia sets, also at several levels of magnification. I'll create a Julia walk that moves through different structures and sub-structures of the Mandelbrot set and shows the change in structure of the Julia set, also at several different scales.

I need as much rendering speed as possible in order to get a reasonable frame rate for such an animation.

I also need at least double precision for my calculations, because the computations fall apart when the magnification gets too high. The different scales in Mandelbrot plots seem to require larges changes in magnification than different scales in Julias. After 4 or 5 scales on Mandelbrot plots, I run out of precision (my computations seem to fall apart at around 10^-14 or so.)


Duncan C

[lycium]

i look forward to seeing the results!

[gandreas]

Can you point me in the direction of some resources on Altavec/SSE optimization? I haven't tackled either before.

http://developer.apple.com/hardwaredrivers/ve/quickstart.html
http://developer.apple.com/documentation/Performance/Conceptual/Accelerate_sse_migration/migration_intro/chapter_1_section_1.html
http://developer.apple.com/hardwaredrivers/ve/sse.html

Also, isn't the optimization different for G4 vs G5 and Intel Core duo vs. Core 2 duo (SSE vs SSE2)?

Not all that much - especially with regards to the vector stuff (and the compiler can handle the instruction scheduling issues).  The G4 and G5 Altivec is essentially identical, and the baseline hardware requirements for OS X on Intel is SSE2 (and SSE3 may be available as well, but its not guaranteed).

It seems to me I'd get decent performance benifits by teaching my code to do simultaneous operations for complex add and multiply operations. (since Z^2 works out to (a+bi)^2, or a^2 +2bi - b^2. The a^, 2bi, and b^2 could all be performed at once. Adding two complex numbers, I should also be able to use Altevec/SSE to add the real and imaginary components simultaneously.)

You'd think so, but not really, because the amount of time it'll take to load up the two value and then write them back out (because the step before or afterwards needs them as singled values) will destroy any benefit.

Ideally, you want to do many similar operations on many parallel independent chunks of data - so you'll get better results if you calculate 4 pixels at once than if you try to split up parts of the calculation for only one 1 pixel.

Is it true that PPC G4 altavec is single precision floating point only? That's a deal-killer for me. My app is based on double precison. Single precision falls apart too fast for high "magnification" zooms. Double precision seems to fall apart at around 10^-14.

That is true - AltiVec only does single precision.  SSE can do double precision, but you only can process two elements at a time then (so you can only theoretically get a maximum of 2x improvement instead of 4x).

Also note that unless you're doing 64 bit mode (which requires Leopard to do anything other than a command line tool), you may not have enough SSE registers available to come close to that 4x (or 2x) improvement.

Note there are tricks that can be done with single precision to get more precision, but basically they involve using two of them to simulate double precision (basically, the number is broken down into two numbers such that their sum is the first number, for example, 123.456 would be broken down into 123e0 + 456e-3, allowing for twice the precision bits, though the same exponent range).

For other sorts of performance enhancements, you may want to look at XGrid, which would allow for distributed rendering (actually, the XGrid sample projects include one to draw a mandelbrot image, though not very well).  The XGrid sample may even have an AltiVec version as well...

[David Makin]

(my computations seem to fall apart at around 10^-14 or so.)
Duncan C

Hi Duncan - if you're using double and the calculations start falling apart at 10^-14 or so then I'm guessing you've implimented the same optimisation that I did in MMFrac i.e. you're using deltas for the x(real) and y(imag) loops.
To get greater precision using double you need to explicitly calculate start+pixelcount*delta for the real and imag loop values instead of using deltas.

[lycium]

should be roughly the same speed.

[Duncan C]

(my computations seem to fall apart at around 10^-14 or so.)
Duncan C

Hi Duncan - if you're using double and the calculations start falling apart at 10^-14 or so then I'm guessing you've implimented the same optimisation that I did in MMFrac i.e. you're using deltas for the x(real) and y(imag) loops.
To get greater precision using double you need to explicitly calculate start+pixelcount*delta for the real and imag loop values instead of using deltas.

Hmm. I do indeed use a precalcuated deta (or step) value to calculate my x and y value for a point. How much more precision do you think I'd gain if I refactored my code to calculate x and y position each time?

And I wonder if this would make a detectable difference in performance. Mulitplies usually take longer than adds. Of course, we're talking a single multiply vs. a single add per pixel, not per iteration.

Duncan C

P.S. I already ran into precision problems. My code does a "snap to grid" on my real and imaginary bounds so that the origin lands on a pixel. That lets me map symmetric regions of a plot without being off by a fraction of a step. Those calculations run into precision problems, and I had to resort to "long double" variables for that (once per plot)


Duncan C

[lycium]

How much more precision do you think I'd gain if I refactored my code to calculate x and y position each time?

lazy! test it out, that's < 1 line of code. or you could do a little thought experiment: how often have your x/y values been subjected to roundoff errors in the further reaches of your image (depending on how you split it)?

i guarantee you that the "snapping" code you have, besides precluding supersampling by the sound of things, is a zillion times slower than just 2 extra multiplies.

And I wonder if this would make a detectable difference in performance. Mulitplies usually take longer than adds.

lazy! read my post above your last with sse, afaik multiply and add have the same latency but you have less throughput with a multiply (half speed?).

[Duncan C]

How much more precision do you think I'd gain if I refactored my code to calculate x and y position each time?

lazy! test it out, that's < 1 line of code. or you could do a little thought experiment: how often have your x/y values been subjected to roundoff errors in the further reaches of your image (depending on how you split it)?

Actually, it's a design change to my app. It's built on passing step values to the plot code.

i guarantee you that the "snapping" code you have, besides precluding supersampling by the sound of things, is a zillion times slower than just 2 extra multiplies.

I think not. My "snap to grid" logic  takes place ONCE, at the beginning of a plot, and involves a half-dozen floating point calculations, once. 2 multiplies/pixel is a lot more calculations than a half-dozen calcuations for the whole plot. Also, my snap to grid code lets me recognize symmetric areas of my plot and skip ALL the calculations for symmetric areas, and instead just copy the pixels. That saves potentially millions of floating point calculations.

And I wonder if this would make a detectable difference in performance. Mulitplies usually take longer than adds.

lazy! read my post above your last with sse, afaik multiply and add have the same latency but you have less throughput with a multiply (half speed?).
[/quote]