|
|
Blitz Support : |
From: Hendrik Belitz (h.belitz_at_[hidden])
Date: 2004-02-27 03:50:18
Am Donnerstag, 26. Februar 2004 19:16 schrieb Todd Veldhuizen:
> Hi Hendrik,
>
> It might help if you gave specifics. Can you show the code you use
> to create the blitz arrays vs. the non-blitz version arrays?
// THis is my array initialisation, AIPSArray is my own storage format
// (which is currently identical to ColumnMajor I think)
Array<double,3> bg( w, h, d, AIPSArray<3>() );
Array<double,3> bhf( w, h, d, AIPSArray<3>() );
Array<TVector3D, 3> bgradient( ptrToCArray, shape( w, h, d ),
neverDeleteData, AIPSArray<3>() );
Array<TVector3D, 3> bflowt = bgradient;
Array<TVector3D, 3> bflowdt( w, h, d, AIPSArray<3>() );
/* TVector3D is just a macro abbreviation for TinyVector<double,3>
// And this is the initialisation with my own classes, which just wrap
// around std::vector for easier access
TField3D flowt( 3 /* Dimension */ , grad.getExtends() /* Extends */ );
TField3D flowdt( 3, grad.getExtends() );
// TField3D are arrays of TinyVector<double,3>, so I don't think TinyVector
// is problem as Todd stated
CTypedData<double> g ( 3, grad.getExtends() );
CTypedData<double> hf ( 3, grad.getExtends() );
// PDE loop (This is the same for both versions)
while ( time <= maxtime )
{
for ( int z = 0; z < d; ++z )
for ( int y = 0; y < h; ++y )
for ( int x = 0; x < w; ++x )
{
ushort times = 0;
if ( x < w-1) { neigh += bflowt( x+1,y,z ); times++; }
if ( x > 0) { neigh += bflowt( x-1,y,z ); times++; }
if ( y < h-1) { neigh += bflowt( x ,y+1,z); times++; }
if ( y > 0) { neigh += bflowt( x ,y-1,z); times++; }
if ( z < d-1) { neigh += bflowt( x ,y,z+1); times++; }
if ( z > 0) { neigh += bflowt( x ,y,z-1); times++; }
neigh -= static_cast<double>( times ) * bflowt(x,y,z);
bflowdt(x,y,z)[0] = ( 1.0 - parametersVec[2].paramVal.dVal * bhf(x,y,z) )
* bflowt(x,y,z)[0] + parametersVec[3].paramVal.dVal
* parametersVec[2].paramVal.dVal * neigh[0] * bg(x,y,z)
+ bgradient(x,y,z)[0] * bhf(x,y,z) * parametersVec[2].paramVal.dVal;
bflowdt(x,y,z)[1] = ( 1.0-parametersVec[2].paramVal.dVal*bhf(x,y,z) )
* bflowt(x,y,z)[1] + parametersVec[3].paramVal.dVal
* parametersVec[2].paramVal.dVal * neigh[1] * bg(x,y,z)
+ bgradient(x,y,z)[1] * bhf(x,y,z) * parametersVec[2].paramVal.dVal;
bflowdt(x,y,z)[2] = ( 1.0 - parametersVec[2].paramVal.dVal*bhf(x,y,z) )
* bflowt(x,y,z)[2] + parametersVec[3].paramVal.dVal
* parametersVec[2].paramVal.dVal * neigh[2] * bg(x,y,z)
+ bgradient(x,y,z)[2] * bhf(x,y,z) * parametersVec[2].paramVal.dVal;
/* Using the following for the blitz version is much slower than the above
(why?)
bflowdt(x,y,z) = ( 1.0 - parametersVec[2].paramVal.dVal*bhf(x,y,z) )
* bflowt(x,y,z) + parametersVec[3].paramVal.dVal
* parametersVec[2].paramVal.dVal * neigh * bg(x,y,z)
+ bgradient(x,y,z) * bhf(x,y,z) * parametersVec[2].paramVal.dVal;
*/
}
cycleArrays( bflowt, bflowdt );
time++;
}
----------------------------------------------------------------
Dipl.-Inform. Hendrik Belitz
Zentralinstitut für Elektronik
Forschungszentrum Jülich GmbH
D-52428 Jülich, Germany
Tel.: (++49)2461 61 4578
Fax: (++49)2461 61 3990
email: h.belitz_at_[hidden]