34 throw ProSHADE_exception (
"Attempted allocating Wigner D matrices before\n : allocating E matrices memory.",
"EW00024", __FILE__, __LINE__, __func__,
"The E matrices and Wigner matrices both require to know\n : the bandwidth of the comparison (which may differ from the\n : object bandwidth). This is set when allocating E matrices\n : and therefore if it is 0 now, E matrices were not yet\n : allocated." );
42 for ( proshade_unsign bandIter = 0; bandIter < this->
maxCompBand; bandIter++ )
45 this->
wignerMatrices[bandIter] =
new proshade_complex* [(bandIter * 2) + 1];
49 for ( proshade_unsign order1Iter = 0; order1Iter < ( (bandIter * 2) + 1 ); order1Iter++ )
51 this->
wignerMatrices[bandIter][order1Iter] =
new proshade_complex [(bandIter * 2) + 1];
77 void ProSHADE_internal_wigner::allocateWignerWorkspace ( proshade_double*& matIn, proshade_double*& matOut, proshade_double*& sqrts, proshade_double*& workspace, proshade_double*& alphaExponentReal, proshade_double*& alphaExponentImag, proshade_double*& gammaExponentReal, proshade_double*& gammaExponentImag, proshade_double*& trigs, proshade_unsign compBand )
80 matIn =
new proshade_double[
static_cast<proshade_unsign
> ( 4 * pow( compBand, 2.0 ) - 4 * compBand + 1 )];
82 matOut =
new proshade_double[
static_cast<proshade_unsign
> ( 4 * pow( compBand, 2.0 ) - 4 * compBand + 1 )];
83 sqrts =
new proshade_double[
static_cast<proshade_unsign
> ( 2 * compBand )];
84 workspace =
new proshade_double[
static_cast<proshade_unsign
> ( 4 * pow( compBand, 2.0 ) )];
85 alphaExponentReal =
new proshade_double[
static_cast<proshade_unsign
> ( 2 * compBand - 1 )];
86 alphaExponentImag =
new proshade_double[
static_cast<proshade_unsign
> ( 2 * compBand - 1 )];
87 gammaExponentReal =
new proshade_double[
static_cast<proshade_unsign
> ( 2 * compBand - 1 )];
88 gammaExponentImag =
new proshade_double[
static_cast<proshade_unsign
> ( 2 * compBand - 1 )];
89 trigs =
new proshade_double[2];
118 void ProSHADE_internal_wigner::releaseWignerWorkspace ( proshade_double*& matIn, proshade_double*& matOut, proshade_double*& sqrts, proshade_double*& workspace, proshade_double*& alphaExponentReal, proshade_double*& alphaExponentImag, proshade_double*& gammaExponentReal, proshade_double*& gammaExponentImag, proshade_double*& trigs )
121 if ( matIn != NULL ) {
delete[] matIn; }
122 if ( matOut != NULL ) {
delete[] matOut; }
123 if ( sqrts != NULL ) {
delete[] sqrts; }
124 if ( workspace != NULL ) {
delete[] workspace; }
125 if ( trigs != NULL ) {
delete[] trigs; }
126 if ( alphaExponentReal != NULL ) {
delete[] alphaExponentReal; }
127 if ( alphaExponentImag != NULL ) {
delete[] alphaExponentImag; }
128 if ( gammaExponentReal != NULL ) {
delete[] gammaExponentReal; }
129 if ( gammaExponentImag != NULL ) {
delete[] gammaExponentImag; }
149 void ProSHADE_internal_wigner::prepareTrigsSqrtsAndExponents ( proshade_double* sqrts, proshade_double* alphaExponentReal, proshade_double* alphaExponentImag, proshade_double* gammaExponentReal, proshade_double* gammaExponentImag, proshade_double* trigs, proshade_unsign compBand, proshade_double angAlpha, proshade_double angBeta, proshade_double angGamma )
152 for ( proshade_unsign iter = 0; iter < ( 2 * compBand ); iter++ )
154 sqrts[iter] =
static_cast<proshade_double
> ( sqrt (
static_cast<proshade_double
> ( iter ) ) );
158 trigs[0] =
static_cast<proshade_double
> ( cos ( 0.5 * -angBeta ) );
159 trigs[1] =
static_cast<proshade_double
> ( sin ( 0.5 * -angBeta ) );
162 genExp ( compBand, angAlpha, alphaExponentReal, alphaExponentImag );
163 genExp ( compBand, angGamma, gammaExponentReal, gammaExponentImag );
189 void ProSHADE_internal_wigner::computeWignerMatrices (
ProSHADE_settings* settings,
ProSHADE_internal_data::ProSHADE_data* obj, proshade_double* alphaExponentReal, proshade_double* alphaExponentImag, proshade_double* gammaExponentReal, proshade_double* gammaExponentImag, proshade_double* matIn, proshade_double* matOut, proshade_double* trigs, proshade_double* sqrts, proshade_double* workspace )
195 proshade_double *expARStart, *expAIStart, *expGRStart, *expGIStart;
196 proshade_double Dij, eARi, eAIi, eGRj, eGIj, iSign, rSign;
197 proshade_complex hlpVal;
198 proshade_unsign noOrders, arrConvIter;
199 for ( proshade_unsign bandIter = 0; bandIter < obj->
getComparisonBand(); bandIter++ )
202 noOrders = 2 * bandIter + 1;
212 wignerdmat ( bandIter, matIn, matOut, trigs, sqrts, workspace );
215 for ( proshade_unsign d1Iter = 0; d1Iter < noOrders; d1Iter++ )
217 eARi = expARStart[d1Iter];
218 eAIi = expAIStart[d1Iter];
220 for ( proshade_unsign d2Iter = 0; d2Iter < noOrders; d2Iter++ )
222 Dij = matOut[arrConvIter];
223 eGRj = expGRStart[d2Iter];
224 eGIj = expGIStart[d2Iter];
226 hlpVal[0] = ( Dij * eGRj * eARi - Dij * eGIj * eAIi ) * rSign;
227 hlpVal[1] = ( Dij * eGRj * eAIi + Dij * eGIj * eARi ) * iSign;
237 memcpy ( matIn, matOut,
sizeof ( proshade_double ) * ( noOrders * noOrders ) );
263 proshade_double *matIn, *matOut, *sqrts, *workspace, *alphaExponentReal, *alphaExponentImag, *gammaExponentReal, *gammaExponentImag, *trigs;
277 computeWignerMatrices ( settings, obj, alphaExponentReal, alphaExponentImag, gammaExponentReal, gammaExponentImag,
278 matIn, matOut, trigs, sqrts, workspace );
282 gammaExponentReal, gammaExponentImag, trigs );