Commit b78fc0b5 authored by David Saxton's avatar David Saxton

Speed up differential equation calculating by several hundred percent. The...

Speed up differential equation calculating by several hundred percent. The main lesson learnt today is to be wary of Qt's implicit sharing: Although it speeds up some usages, repeatedly freeing memory (on operator=) and then allocating memory (on detach()) can be overall a lot slower.

svn path=/trunk/KDE/kdeedu/kmplot/; revision=653656
parent bba8a84c
......@@ -281,6 +281,7 @@ Equation::Equation( Type type, Function * parent )
: m_type( type ),
m_parent( parent )
{
m_usesParameter = false;
mptr = 0;
if ( type == Differential || type == Cartesian )
......@@ -433,11 +434,9 @@ void Equation::updateVariables()
n += '\'';
}
}
}
bool Equation::usesParameter( ) const
{
//BEGIN Update whether we accept a parameter or not
int expectedNumVariables = 0;
switch ( m_type )
......@@ -457,11 +456,13 @@ bool Equation::usesParameter( ) const
expectedNumVariables = order()+1;
break;
case Constant:
return false;
case Constant:
expectedNumVariables = 0;
break;
}
return variables().size() > expectedNumVariables;
m_usesParameter = (variables().size() > expectedNumVariables);
//END Update whether we accept a parameter or not
}
......
......@@ -222,6 +222,7 @@ class DifferentialStates
const DifferentialState & operator[] ( int i ) const { return m_data[i]; }
void remove ( int i ) { m_data.remove(i); }
void remove ( int i, int count ) { m_data.remove( i, count ); }
void removeAll() { m_data.clear(); }
protected:
QVector<DifferentialState> m_data;
......@@ -290,7 +291,7 @@ class Equation
* \return whether the function accepts a parameter in addition to the x
* (and possibly y) variables.
*/
bool usesParameter() const;
bool usesParameter() const { return m_usesParameter; }
/**
* \return the name of the parameter variable (or a blank string if a
* parameter is not used).
......@@ -350,6 +351,7 @@ class Equation
*/
void updateVariables();
bool m_usesParameter;
const Type m_type;
QString m_fstr;
Function * m_parent;
......
......@@ -26,6 +26,7 @@
#include "vector.h"
#include <assert.h>
#include <string.h>
//BEGIN class Vector
Vector Vector::operator + ( const Vector & other ) const
......@@ -80,6 +81,22 @@ Vector & Vector::operator *= ( double x )
}
Vector & Vector::operator = ( const Vector & other )
{
// I don't much like Qt's copy-on-write behaviour.
// At least in KmPlot, it makes a small number of cases marginally faster,
// and a lot of cases a lot slower. This is because allocating and freeing
// memory is a lot slower than copying memory.
//
// In fact, it can make drawing differential equations several times slower.
resize( other.size() );
memcpy( m_data.data(), other.m_data.data(), size() * sizeof(double) );
return *this;
}
Vector & Vector::operator = ( const QVector<Value> & other )
{
int size = other.size();
......@@ -88,4 +105,35 @@ Vector & Vector::operator = ( const QVector<Value> & other )
(*this)[i] = other[i].value();
return *this;
}
void Vector::combine( const Vector & a, double k, const Vector & b )
{
assert( a.size() == b.size() );
int n = a.size();
resize( n );
double *d1 = m_data.data();
const double *d2 = a.m_data.data();
const double *d3 = b.m_data.data();
for ( int i = 0; i < n; ++i )
d1[i] = d2[i] + k * d3[i];
}
void Vector::addRK4( double dx, const Vector & k1, const Vector & k2, const Vector & k3, const Vector & k4 )
{
double *d = m_data.data();
const double *d1 = k1.m_data.data();
const double *d2 = k2.m_data.data();
const double *d3 = k3.m_data.data();
const double *d4 = k4.m_data.data();
int n = size();
for ( int i = 0; i < n; ++i )
d[i] += (dx/6) * (d1[i] + 2*d2[i] + 2*d3[i] + d4[i]);
}
//END class Vector
......@@ -41,18 +41,29 @@ class Vector
int size() const { return m_data.size(); }
void resize( int s ) { if ( size() != s ) m_data.resize( s ); }
double *data() { return m_data.data(); }
const double *data() const { return m_data.data(); }
Vector operator * ( double x ) const;
Vector & operator *= ( double x );
Vector operator + ( const Vector & other ) const;
Vector & operator += ( const Vector & other );
Vector operator - ( const Vector & other ) const;
Vector & operator -= ( const Vector & other );
Vector & operator = ( const Vector & other ) { m_data = other.m_data; return *this; }
Vector & operator = ( const Vector & other );
Vector & operator = ( const QVector<Value> & other );
bool operator==( const Vector & other ) const { return m_data == other.m_data; }
bool operator!=( const Vector & other ) const { return m_data != other.m_data; }
double & operator[] ( int i ) { return m_data[i]; }
double operator[] ( int i ) const { return m_data[i]; }
/**
* Optimization for use in solving differential equations. Sets the
* contents of this vector to a+k*b.
*/
void combine( const Vector & a, double k, const Vector & b );
/**
* Another optimization for use in solving differential equaions.
*/
void addRK4( double dx, const Vector & k1, const Vector & k2, const Vector & k3, const Vector & k4 );
protected:
QVector<double> m_data;
......
......@@ -288,12 +288,8 @@ Vector XParser::rk4_f( int order, Equation * eq, double x, const Vector & y )
if ( useParameter )
m_arg[1] = eq->parent()->k;
for ( int i = 0; i < order; ++i )
{
m_arg[i+1 + (useParameter ? 1 : 0) ] = y[i];
if ( i+1 < order )
m_result[i] = y[i+1];
}
memcpy( m_arg.data() + 1 + (useParameter ? 1 : 0), y.data(), order*sizeof(double) );
memcpy( m_result.data(), y.data() + 1, (order-1)*sizeof(double) );
m_result[order-1] = XParser::fkt( eq, m_arg );
......@@ -329,6 +325,7 @@ double XParser::differential( Equation * eq, DifferentialState * state, double x
m_k2.resize( order );
m_k3.resize( order );
m_k4.resize( order );
m_y_temp.resize( order );
double x = state->x;
m_y = state->y;
......@@ -345,12 +342,19 @@ double XParser::differential( Equation * eq, DifferentialState * state, double x
x = state->x + i*dx;
m_k1 = rk4_f( order, eq, x, m_y );
m_k2 = rk4_f( order, eq, x + dx/2, m_y + (dx/2)*m_k1 );
m_k3 = rk4_f( order, eq, x + dx/2, m_y + (dx/2)*m_k2 );
m_k4 = rk4_f( order, eq, x + dx, m_y + dx*m_k3 );
m_y += (dx/6)*(m_k1 + 2*m_k2 + 2*m_k3 + m_k4);
m_k1 = rk4_f( order, eq, x, m_y );
m_y_temp.combine( m_y, dx/2, m_k1 );
m_k2 = rk4_f( order, eq, x + dx/2, m_y_temp);
m_y_temp.combine( m_y, dx/2, m_k2 );
m_k3 = rk4_f( order, eq, x + dx/2, m_y_temp );
m_y_temp.combine( m_y, dx, m_k3 );
m_k4 = rk4_f( order, eq, x + dx, m_y_temp );
m_y.addRK4( dx, m_k1, m_k2, m_k3, m_k4 );
if ( !std::isfinite(m_y[0]) )
{
......
......@@ -149,7 +149,7 @@ class XParser : public Parser
private:
/// for use in differential
Vector rk4_f( int order, Equation * eq, double x, const Vector & y );
Vector m_k1, m_k2, m_k3, m_k4, m_y, m_result, m_arg;
Vector m_k1, m_k2, m_k3, m_k4, m_y_temp, m_y, m_result, m_arg;
/// indicates if the widget is changed
bool & m_modified;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment