Changeset 2958 in Sophya


Ignore:
Timestamp:
Jun 1, 2006, 1:34:50 PM (19 years ago)
Author:
ansari
Message:

1/ passage en sa_size_t (au lieu de int) dans Alm<T> et Bm<T>
2/ Ajout des methodes optimisees (statiques) pour calcul transforme Ylm
ds LambdaLMBuilder et utilisation ds SphericalTransformServer

Reza , 1/06/2006

Location:
trunk/SophyaLib/Samba
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/Samba/alm.h

    r2885 r2958  
    2121Alm() : TriangularMatrix<complex<T> >()  {;};
    2222    /* instanciate the class from the maximum value of the index \f$l\f$ in the sequence of the \f$a_{lm}\f$'s */
    23 Alm(int nlmax) : TriangularMatrix<complex<T> >(nlmax+1) {;}
     23Alm(sa_size_t nlmax) : TriangularMatrix<complex<T> >(nlmax+1) {;}
    2424Alm(const Alm<T>& a,  bool share=false)  : TriangularMatrix<complex<T> >(a, share)  {;}
    2525
    2626Alm(const TVector<T>& clin, const r_8 fwhm) ;
    2727/*! resize with a new lmax */
    28 inline void ReSizeToLmax(int_4 nlmax) {this->ReSizeRow(nlmax+1);}
    29 inline int_4 Lmax() const {return this->rowNumber()-1;}
     28inline void ReSizeToLmax(sa_size_t nlmax) {this->ReSizeRow(nlmax+1);}
     29inline sa_size_t Lmax() const {return this->rowNumber()-1;}
    3030TVector<T> powerSpectrum() const;
    3131
    … …  
    4848Bm(): nmmax_(0) {;};
    4949    /* instanciate from the maximum absolute value of the index */
    50 Bm(int mmax) : nmmax_((uint_4)mmax)   { bm_.ReSize( (uint_4)(2*mmax+1) );}
     50Bm(sa_size_t mmax) : nmmax_(mmax)   { bm_.ReSize( (2*mmax+1) );}
    5151 Bm(const Bm<T>& b,  bool share=false) : nmmax_(b.nmmax_), bm_(b.bm_, share) {;}
    5252 /*! resize with a new value of mmax */
    53 inline void ReSizeToMmax(int_4 mmax)
     53inline void ReSizeToMmax(sa_size_t mmax)
    5454  {
    55     nmmax_=(uint_4)mmax;
     55    nmmax_= mmax;
    5656    bm_.ReSize(2* nmmax_+1);
    5757  }
    58 //inline int_4 Size() const {return 2*nmmax_+1;};
    59 inline T& operator()(int m) {return bm_( adr_i(m));};
    60 inline T const& operator()(int m) const {return *(bm_.Begin()+ adr_i(m));};
     58//inline sa_size_t Size() const {return 2*nmmax_+1;};
     59inline T& operator()(sa_size_t m) {return bm_( adr_i(m));};
     60inline T const& operator()(sa_size_t m) const {return *(bm_.Begin()+ adr_i(m));};
    6161/*! return the current value of the maximum absolute value of the index */
    62 inline int_4 Mmax() const
     62inline sa_size_t Mmax() const
    6363   {
    64      return (int_4)nmmax_;
     64     return nmmax_;
    6565   }
    6666inline Bm<T>& operator = (const Bm<T>& b)
    … …  
    7676  private:
    7777/*! return the address in the array representing the vector */
    78 inline uint_4 adr_i(int i) const
     78inline sa_size_t adr_i(sa_size_t i) const
    7979{
    80   return(uint_4)(i+nmmax_);
     80  return (i+nmmax_);
    8181}
    8282
    8383
    8484
    85  uint_4 nmmax_;
     85 sa_size_t nmmax_;
    8686 NDataBlock<T> bm_;
    8787

  • trunk/SophyaLib/Samba/lambdaBuilder.cc

    r2615 r2958  
    8282void LambdaLMBuilder::array_init(int lmax, int mmax)
    8383   {
    84      if (a_recurrence_.Size() == 0)
    85        {
    86          //      a_recurrence_ = new TriangularMatrix<r_8>;
    87          updateArrayRecurrence(lmax);
    88        }
    89      else
    90        if ( lmax > (a_recurrence_).rowNumber()-1     )
    91          {
    92            cout << " WARNING : The classes LambdaXXBuilder will be more efficient if instanciated with parameter lmax = maximum value of l index which will be needed in the whole application (arrays not recomputed) " << endl;
    93            cout << "lmax= " << lmax << " previous instanciation with lmax= " <<  (a_recurrence_).rowNumber() << endl;
    94          updateArrayRecurrence(lmax);
    95          }
     84     updateArrayRecurrence(lmax);
     85
    9686     lmax_=lmax;
    9787     mmax_=mmax;
    … …  
    129119   }
    130120
     121/*!
     122  \brief : Specialized/optimized static function for fast spherical harmonic transform.
     123  Computes bm(m) = Sum_l>=m [ lambda(l,m) * alm(l,m) ]
     124  See SphericalTransformServer<T>::GenerateFromAlm(map, pixsize, alm)
     125*/
     126void LambdaLMBuilder::ComputeBmFrAlm(r_8 theta,int_4 lmax, int_4 mmax,
     127                                     const Alm<r_8>& alm, Bm< complex<r_8> >& bm)
     128{
     129  updateArrayRecurrence(lmax);
     130  r_8 cth = cos(theta);
     131  r_8 sth = sin(theta);
     132 
     133  r_8 bignorm2 = 1.e268; // = 1e-20*1.d288
     134 
     135  r_8 lam_mm = 1. / sqrt(4.*Pi) *bignorm2;
     136  r_8 lam_0, lam_1, lam_2;
     137  register r_8 a_rec, b_rec;
     138
     139  const complex<r_8>* almp = alm.columnData(0);
     140  r_8* arecp = a_recurrence_.columnData(0);
     141
     142  int_4 l, m;
     143
     144  for (m=0; m<=mmax;m++)   {
     145    //MOV^     const complex<r_8>* almp = alm.columnData(m);
     146    complex<r_8>* bmp = &(bm(m));
     147   
     148    *bmp = (lam_mm / bignorm2)*(*almp);  almp++;
     149
     150    lam_0=0.;
     151    lam_1=1. /bignorm2 ;
     152   
     153    a_rec = *arecp; arecp++; // a_recurrence_(m,m);
     154    b_rec = 0.;
     155    for (l=m+1; l<=lmax; l++) {
     156      lam_2 = (cth*lam_1-b_rec*lam_0)*a_rec;
     157     
     158      //DEL        lambda_(l,m) = lam_2*lam_mm;
     159      *bmp += lam_2*lam_mm*(*almp);   almp++;
     160      b_rec=1./a_rec;
     161      //           a_rec= LWK->a_recurr(l,m);
     162      a_rec= *arecp; arecp++; // a_recurrence_(l,m);
     163      lam_0 = lam_1;
     164      lam_1 = lam_2;
     165    }
     166    lam_mm = -lam_mm*sth* sqrt( (2.*m+3.)/ (2.*m+2.) );
     167  }
     168}
     169
     170/*!
     171  \brief : Specialized/optimized static function for fast spherical harmonic transform.
     172*/
     173void LambdaLMBuilder::ComputeBmFrAlm(r_8 theta,int_4 lmax, int_4 mmax, 
     174                                     const Alm<r_4>& alm, Bm< complex<r_4> >& bm)
     175{
     176  Alm<r_8> alm8(alm.Lmax());
     177  for(sa_size_t k=0; k<alm8.Size(); k++)
     178    alm8(k) = complex<r_8>((r_8)alm(k).real() , (r_8)alm(k).imag());
     179  Bm< complex<r_8> > bm8(bm.Mmax());
     180  ComputeBmFrAlm(theta, lmax, mmax, alm8, bm8);
     181  for(sa_size_t kk=-bm.Mmax(); kk<=bm.Mmax(); kk++)
     182    bm(kk)= complex<r_4>((r_4)bm8(kk).real() , (r_4)bm8(kk).imag());
     183  return;
     184}
     185
     186
     187/*!
     188  \brief : Specialized/optimized static function for fast spherical harmonic transform.
     189  Computes alm(l,m) = Sum_l>=m [ lambda(l,m) * phase(l,m) ]
     190  See SphericalTransformServer<T>::carteVersAlm(SphericalMap<T>& map, nlmax, ctcut, alm)
     191*/
     192void LambdaLMBuilder::ComputeAlmFrPhase(r_8 theta,int_4 lmax, int_4 mmax,
     193                                        TVector< complex<r_8> >& phase, Alm<r_8> & alm)
     194{
     195  updateArrayRecurrence(lmax);
     196  r_8 cth = cos(theta);
     197  r_8 sth = sin(theta);
     198 
     199  r_8 bignorm2 = 1.e268; // = 1e-20*1.d288
     200 
     201  r_8 lam_mm = 1. / sqrt(4.*Pi) *bignorm2;
     202  r_8 lam_0, lam_1, lam_2;
     203  register r_8 a_rec, b_rec;
     204 
     205  complex<r_8>* almp = alm.columnData(0);
     206  r_8* arecp = a_recurrence_.columnData(0);
     207
     208  int_4 l, m;
     209  for (m=0; m<=mmax;m++)  {
     210    //MOV^ complex<r_8>* almp = alm.columnData(m);
     211    complex<r_8>  phi = phase(m);
     212   
     213    *almp += (lam_mm / bignorm2)*phi;    almp++;
     214
     215    lam_0=0.;
     216    lam_1=1. /bignorm2 ;
     217   
     218    //MOV^ r_8* arecp = a_recurrence_.columnData(m);
     219    a_rec = *arecp; arecp++; // a_recurrence_(m,m);
     220    b_rec = 0.;
     221    for (l=m+1; l<=lmax; l++)   {
     222      lam_2 = (cth*lam_1-b_rec*lam_0)*a_rec;
     223     
     224      //DEL        lambda_(l,m) = lam_2*lam_mm;
     225      *almp += (lam_2*lam_mm) * phi;
     226      almp++;
     227     
     228      b_rec=1./a_rec;
     229      //           a_rec= LWK->a_recurr(l,m);
     230      a_rec= *arecp; arecp++; // a_recurrence_(l,m);
     231      lam_0 = lam_1;
     232      lam_1 = lam_2;
     233    }
     234    lam_mm = -lam_mm*sth* sqrt( (2.*m+3.)/ (2.*m+2.) );
     235  }
     236}
     237
     238/*!
     239  \brief : Specialized/optimized static function for fast spherical harmonic transform.
     240  Computes alm(l,m) = Sum_l>=m [ lambda(l,m) * phase(l,m) ]
     241  See SphericalTransformServer<T>::carteVersAlm(SphericalMap<T>& map, nlmax, ctcut, alm)
     242*/
     243void LambdaLMBuilder::ComputeAlmFrPhase(r_8 theta,int_4 lmax, int_4 mmax,
     244                                        TVector< complex<r_4> >& phase, Alm<r_4> & alm)
     245{
     246  TVector< complex<r_8> > phase8;
     247  phase8 = phase;
     248  Alm<r_8> alm8(alm.Lmax());
     249  ComputeAlmFrPhase(theta, lmax, mmax, phase8, alm8);
     250  for(sa_size_t k=0; k<alm.Size(); k++)
     251    alm(k) = complex<r_4>((r_4)alm8(k).real() , (r_4)alm8(k).imag());
     252  return;
     253}
     254
    131255
    132256/*! \fn void  SOPHYA::LambdaLMBuilder::updateArrayRecurrence(int_4 lmax)
    … …  
    136260void  LambdaLMBuilder::updateArrayRecurrence(int_4 lmax)
    137261   {
     262     if ( (a_recurrence_.Size() > 0) && (lmax < a_recurrence_.rowNumber()) )
     263       return;  // Pas besoin de recalculer le tableau de recurrence
     264     if (a_recurrence_.Size() > 0)  {
     265       cout << " WARNING : The classes LambdaXXBuilder will be more efficient "
     266            << "if instanciated with parameter lmax = maximum value of l index "
     267            << "which will be needed in the whole application (arrays not "
     268            << "recomputed) " << endl;
     269       cout << "   lmax= " << lmax << " previous instanciation with lmax= "
     270            <<  a_recurrence_.rowNumber()-1 << endl;
     271         }
     272
    138273     a_recurrence_.ReSizeRow(lmax+1);
    139274     for (int m=0; m<=lmax;m++)

  • trunk/SophyaLib/Samba/lambdaBuilder.h

    r2291 r2958  
    4646LambdaLMBuilder(r_8 theta,int_4 lmax, int_4 mmax);
    4747LambdaLMBuilder(r_8 ct, r_8 st,int_4 lmax, int_4 mmax);
     48
    4849virtual ~LambdaLMBuilder() {};
     50
     51
     52//--- Optimisation Reza mai 2006 Serie de fonctions statiques
     53// specialisees/optimisees pour SphericalTransform
     54static void ComputeBmFrAlm(r_8 theta,int_4 lmax, int_4 mmax,
     55                           const Alm<r_8>& alm, Bm< complex<r_8> >& bm);
     56static void ComputeBmFrAlm(r_8 theta,int_4 lmax, int_4 mmax,
     57                           const Alm<r_4>& alm, Bm< complex<r_4> >& bm);
     58static void ComputeAlmFrPhase(r_8 theta,int_4 lmax, int_4 mmax,
     59                              TVector< complex<r_8> >& phase, Alm<r_8> & alm);
     60static void ComputeAlmFrPhase(r_8 theta,int_4 lmax, int_4 mmax,
     61                              TVector< complex<r_4> >& phase, Alm<r_4> & alm);
     62//---- Fin foctions specialisees
    4963
    5064 /*! return the value of the coefficient \f$  \lambda_l^m \f$ */
    … …  
    6579
    6680 private:
    67  void updateArrayRecurrence(int_4 lmax);
     81 static void updateArrayRecurrence(int_4 lmax);
    6882
    6983 void array_init(int lmax, int mmax);

  • trunk/SophyaLib/Samba/sphericaltransformserver.cc

    r2872 r2958  
    232232      //              b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m)
    233233      //              ------------------------------------------------------
     234      // ===> Optimisation Reza, Mai 2006
     235      /*---  Le bout de code suivant est remplace par l'appel a la nouvelle fonction
     236        qui calcule la somme au vol
    234237      LambdaLMBuilder lb(theta,nlmax,nmmax);
    235238      //  somme sur m de 0 a l'infini
    236       int m;
    237       for (m = 0; m <= nmmax; m++)
     239      for (int_4 m = 0; m <= nmmax; m++)
    238240        {
    239241          b_m_theta(m) = (T)( lb.lamlm(m,m) ) * alm(m,m);
    … …  
    243245            }
    244246        }
     247      ------- Fin version PRE-Mai2006 */
     248      LambdaLMBuilder::ComputeBmFrAlm(theta,nlmax,nmmax, alm, b_m_theta);
     249      //Fin Optimisation Reza, Mai 2006 <====
     250
    245251    //        obtains the negative m of b(m,theta) (= complex conjugate)
    246252
    247       for (m=1;m<=nmmax;m++)
     253      for (int_4 m=1;m<=nmmax;m++)
    248254        {
    249255          b_m_theta(-m) = conj(b_m_theta(m));
    … …  
    555561//      for each m and l
    556562//      -----------------------------------------------------------------------
     563
     564      // ===> Optimisation Reza, Mai 2006
     565      /*---  Le bout de code suivant est remplace par l'appel a la nouvelle fonction
     566        qui calcule la somme au vol
    557567      //        PrtTim("avant instanciation LM ");
    558568      LambdaLMBuilder lb(theta,nlmax,nmmax);
    … …  
    576586            }
    577587        }
     588      ------- Fin version PRE-Mai2006 */
     589      r_8 domega=map.PixSolAngle(map.PixIndexSph(theta,phi0));
     590      phase *= complex<T>((T)domega, 0.);
     591      LambdaLMBuilder::ComputeAlmFrPhase(theta,nlmax,nmmax, phase, alm);
     592      //Fin Optimisation Reza, Mai 2006 <====
    578593     
    579594
Note: See TracChangeset for help on using the changeset viewer.