[ODE] ODE announcement: new iterative solver

[Russ Smith]

hi all,

i have implemented a new, iterative solver for ODE that will make it as
fast as the other physics engines out there. the new solver (called
'QuickStep') operates in linear time, not O(n^3) like the old big-matrix
solver. it also uses substantially less memory. it's not quite ready to
check in yet - a day of two of cleanup and documenting are required. 
Alen Ladavac at croteam has generously shared a few tricks with me and
offered the benefit of their extensive testing with a similar algorithm:
thanks Alen!

the core of the new solver (i.e. where it spends all its time) is very
simple compared to the nightmare of complexity that is the old 'big
matrix' solver, so there are many optimization opportunities. to prove
it, the complete core of the new solver is at the bottom of this email.

after a period of testing, the new solver will become the default
method, so hopefully we can silence those complaints about "ODE goes
slow and crashes when i stack too many things".

i need some volunteers to test the new solver, and if you can provide
speed comparisons with other engines that would be great too.

russ.

--
the core of QuickStep right now is as follows (there's still a few
tweaks left to do though) :


for (int iteration=0; iteration < num_iterations; iteration++) {
    for (int i=0; i<m; i++) {
	J_ptr = J + i*12;
	iMJ_ptr = iMJ + i*12;
	int b1 = jb[i*2];
	int b2 = jb[i*2+1];
	dReal delta = b[i] - lambda[i]*Ad[i];
	dRealMutablePtr fc_ptr = fc + 6*b1;
	for (j=0; j<6; j++) delta -= fc_ptr[j] * J_ptr[j];
	J_ptr += 6;
	if (b2 >= 0) {
		fc_ptr = fc + 6*b2;
		for (j=0; j<6; j++) delta -= fc_ptr[j] * J_ptr[j];
	}
	dReal new_lambda = lambda[i] + delta;
	if (new_lambda < lo[i]) {
		delta = lo[i]-lambda[i];
		lambda[i] = lo[i];
	}
	else if (new_lambda > hi[i]) {
		delta = hi[i]-lambda[i];
		lambda[i] = hi[i];
	}
	else {
		lambda[i] = new_lambda;
	}
	fc_ptr = fc + 6*b1;
	for (j=0; j<6; j++) fc_ptr[j] += delta * iMJ_ptr[j];
	iMJ_ptr += 6;
	if (b2 >= 0) {
		fc_ptr = fc + 6*b2;
		for (j=0; j<6; j++) fc_ptr[j] += delta * iMJ_ptr[j];
	}
    }
}


-- 
Russell Smith
http://www.q12.org