ssht_test.f90

Go to the documentation of this file.
00001 ! SSHT package to perform spin spherical harmonic transforms
00002 ! Copyright (C) 2011  Jason McEwen
00003 ! See LICENSE.txt for license details
00004 
00005 
00006 !------------------------------------------------------------------------------
00007 ! ssht_test
00008 !
00021 
00022 program ssht_test
00023 
00024   use ssht_types_mod
00025   use ssht_error_mod
00026   use ssht_sampling_mod
00027   use ssht_core_mod
00028 #ifdef NAGFOR
00029   use F90_UNIX_ENV
00030 #endif
00031 
00032   implicit none
00033 
00034   interface
00035      subroutine ssht_test_gen_flm_real(L, flm, seed)
00036        use ssht_types_mod, only: dpc
00037        implicit none
00038        integer, intent(in) :: L
00039        complex(dpc), intent(out) :: flm(0:L**2-1)
00040        integer, intent(in) :: seed
00041      end subroutine ssht_test_gen_flm_real
00042 
00043      subroutine ssht_test_gen_flm_complex(L, spin, flm, seed)
00044        use ssht_types_mod, only: dpc
00045        implicit none
00046        integer, intent(in) :: L
00047        integer, intent(in) :: spin
00048        complex(dpc), intent(out) :: flm(0:L**2-1)
00049        integer, intent(in) :: seed
00050      end subroutine ssht_test_gen_flm_complex
00051 
00052   end interface
00053 
00054   character(len=64) :: arg
00055   integer, parameter :: N_repeat = 3
00056   integer :: verbosity = 0
00057   integer :: fail = 0, seed, i_repeat
00058   real :: time_start, time_end
00059 
00060   real(dp) :: errors_dh(0:N_repeat-1)
00061   real :: durations_forward_dh(0:N_repeat-1)
00062   real :: durations_inverse_dh(0:N_repeat-1)
00063   real(dp) :: errors_gl(0:N_repeat-1)
00064   real :: durations_forward_gl(0:N_repeat-1)
00065   real :: durations_inverse_gl(0:N_repeat-1)
00066   real(dp) :: errors_mweo(0:N_repeat-1)
00067   real :: durations_forward_mweo(0:N_repeat-1)
00068   real :: durations_inverse_mweo(0:N_repeat-1)
00069   real(dp) :: errors_mw(0:N_repeat-1)
00070   real :: durations_forward_mw(0:N_repeat-1)
00071   real :: durations_inverse_mw(0:N_repeat-1)
00072 
00073   real(dp) :: errors_dh_real(0:N_repeat-1)
00074   real :: durations_forward_dh_real(0:N_repeat-1)
00075   real :: durations_inverse_dh_real(0:N_repeat-1)
00076   real(dp) :: errors_gl_real(0:N_repeat-1)
00077   real :: durations_forward_gl_real(0:N_repeat-1)
00078   real :: durations_inverse_gl_real(0:N_repeat-1)
00079   real(dp) :: errors_mweo_real(0:N_repeat-1)
00080   real :: durations_forward_mweo_real(0:N_repeat-1)
00081   real :: durations_inverse_mweo_real(0:N_repeat-1)
00082   real(dp) :: errors_mw_real(0:N_repeat-1)
00083   real :: durations_forward_mw_real(0:N_repeat-1)
00084   real :: durations_inverse_mw_real(0:N_repeat-1)
00085 
00086   integer :: L, ind, ind_check, el, el_check, m, m_check
00087   integer :: spin
00088   complex(dpc), allocatable :: flm_orig(:), flm_syn(:)
00089   complex(dpc), allocatable :: f_dh(:,:), f_gl(:,:), f_mweo(:,:), f_mw(:,:)
00090   real(dp), allocatable :: f_dh_real(:,:), f_gl_real(:,:), f_mweo_real(:,:), f_mw_real(:,:)
00091   real(dp) :: phi_sp_mweo, phi_sp_mw
00092   complex(dpc) :: f_sp_mweo, f_sp_mw
00093   real(dp) :: f_real_sp_dh, f_real_sp_gl, f_real_sp_mweo, f_real_sp_mw
00094 
00095   ! Initialise parameters.
00096   call getarg(1, arg)
00097   read(arg,*) L
00098   call getarg(2, arg)
00099   read(arg,*) spin
00100   seed = 1
00101 
00102   ! Allocate memory.
00103   allocate(flm_orig(0:L**2-1), stat=fail)
00104   allocate(flm_syn(0:L**2-1), stat=fail)  
00105   allocate(f_dh(0:2*L-1, 0:2*L-2), stat=fail)
00106   allocate(f_gl(0:L-1, 0:2*L-2), stat=fail)
00107   allocate(f_mweo(0:L-2, 0:2*L-2), stat=fail)
00108   allocate(f_mw(0:L-2, 0:2*L-2), stat=fail)
00109   allocate(f_dh_real(0:2*L-1, 0:2*L-2), stat=fail)
00110   allocate(f_gl_real(0:L-1, 0:2*L-2), stat=fail)
00111   allocate(f_mweo_real(0:L-2, 0:2*L-2), stat=fail)
00112   allocate(f_mw_real(0:L-2, 0:2*L-2), stat=fail)
00113   if(fail /= 0) then
00114      call ssht_error(SSHT_ERROR_MEM_ALLOC_FAIL, 'ssht_test')
00115   end if
00116   flm_orig(0:L**2-1) = cmplx(0d0, 0d0)
00117   flm_syn(0:L**2-1) = cmplx(0d0, 0d0)
00118   f_dh(0:2*L-1, 0:2*L-2) = cmplx(0d0, 0d0)
00119   f_gl(0:L-1, 0:2*L-2) = cmplx(0d0, 0d0)
00120   f_mweo(0:L-2, 0:2*L-2) = cmplx(0d0, 0d0)
00121   f_mw(0:L-2, 0:2*L-2) = cmplx(0d0, 0d0)
00122   f_dh_real(0:2*L-1, 0:2*L-2) = 0d0
00123   f_gl_real(0:L-1, 0:2*L-2) = 0d0
00124   f_mweo_real(0:L-2, 0:2*L-2) = 0d0
00125   f_mw_real(0:L-2, 0:2*L-2) = 0d0
00126 
00127   ! Write program name.
00128   write(*,*)
00129   write(*,'(a)') 'SSHT test program'
00130   write(*,'(a)') '==============================================================='
00131   write(*,*)
00132 
00133   ! Test index conversion between (el,m) and ind.
00134   ind_check = 0
00135   do el = 0, L-1     
00136      do m = -el, el       
00137         call ssht_sampling_elm2ind(ind, el, m)
00138         if (ind /= ind_check) then
00139            call ssht_error(SSHT_ERROR_INDEX_INVALID, 'ssht_test')
00140         end if
00141         call ssht_sampling_ind2elm(el_check, m_check, ind)
00142         if (el /= el_check .and. m /= m_check) then
00143            call ssht_error(SSHT_ERROR_INDEX_INVALID, 'ssht_test')
00144         end if
00145         ind_check = ind_check + 1
00146      end do
00147   end do
00148   if (ind_check /= L**2) then
00149      call ssht_error(SSHT_ERROR_INDEX_INVALID, 'ssht_test')
00150   end if
00151 
00152   ! Run algorithm error and timing tests.
00153   do i_repeat = 0,N_repeat-1
00154 
00155      ! If spin=0 run tests on algorithms optimised for real spin=0 signal.
00156      if (spin == 0) then 
00157 
00158         !=========================================================================
00159         ! DH real spin=0
00160         write(*,'(a,i2)') 'DH real spin=0 test no.', i_repeat
00161         flm_orig(0:L**2-1) = cmplx(0d0, 0d0)
00162         flm_syn(0:L**2-1) = cmplx(0d0, 0d0)
00163         call ssht_test_gen_flm_real(L, flm_orig, seed)
00164         call cpu_time(time_start)
00165         !-------------------------------------------------------------------------
00166         call ssht_core_dh_inverse_real(f_dh_real, flm_orig, L, verbosity)
00167         !-------------------------------------------------------------------------
00168         call cpu_time(time_end)
00169         durations_inverse_dh_real(i_repeat) = time_end - time_start
00170         call cpu_time(time_start)
00171         !-------------------------------------------------------------------------
00172         call ssht_core_dh_forward_real(flm_syn, f_dh_real, L, verbosity)
00173         !-------------------------------------------------------------------------
00174         call cpu_time(time_end)
00175         durations_forward_dh_real(i_repeat) = time_end - time_start
00176         errors_dh_real(i_repeat) = maxval(abs(flm_orig(0:L**2-1) - flm_syn(0:L**2-1)))
00177         write(*,'(a,f40.4)') ' duration_inverse (s) =', &
00178              durations_inverse_dh_real(i_repeat)
00179         write(*,'(a,f40.4)') ' duration_forward (s) =', &
00180              durations_forward_dh_real(i_repeat)
00181         write(*,'(a,e40.5)') ' error                =', &
00182              errors_dh_real(i_repeat)
00183         write(*,*)
00184 
00185         !=========================================================================
00186         ! GL real spin=0
00187         write(*,'(a,i2)') 'GL real spin=0 test no.', i_repeat
00188         flm_orig(0:L**2-1) = cmplx(0d0, 0d0)
00189         flm_syn(0:L**2-1) = cmplx(0d0, 0d0)
00190         call ssht_test_gen_flm_real(L, flm_orig, seed)
00191         call cpu_time(time_start)
00192         !-------------------------------------------------------------------------
00193         call ssht_core_gl_inverse_real(f_gl_real, flm_orig, L, verbosity)
00194         !-------------------------------------------------------------------------
00195         call cpu_time(time_end)
00196         durations_inverse_gl_real(i_repeat) = time_end - time_start
00197         call cpu_time(time_start)
00198         !-------------------------------------------------------------------------
00199         call ssht_core_gl_forward_real(flm_syn, f_gl_real, L, verbosity)
00200         !-------------------------------------------------------------------------
00201         call cpu_time(time_end)
00202         durations_forward_gl_real(i_repeat) = time_end - time_start
00203         errors_gl_real(i_repeat) = maxval(abs(flm_orig(0:L**2-1) - flm_syn(0:L**2-1)))
00204         write(*,'(a,f40.4)') ' duration_inverse (s) =', &
00205              durations_inverse_gl_real(i_repeat)
00206         write(*,'(a,f40.4)') ' duration_forward (s) =', &
00207              durations_forward_gl_real(i_repeat)
00208         write(*,'(a,e40.5)') ' error                =', &
00209              errors_gl_real(i_repeat)
00210         write(*,*)
00211 
00212         !=========================================================================
00213         ! MWEO real spin=0
00214         write(*,'(a,i2)') 'MWEO real spin=0 test no.', i_repeat
00215         flm_orig(0:L**2-1) = cmplx(0d0, 0d0)
00216         flm_syn(0:L**2-1) = cmplx(0d0, 0d0)
00217         call ssht_test_gen_flm_real(L, flm_orig, seed)
00218         call cpu_time(time_start)
00219         !-------------------------------------------------------------------------
00220         call ssht_core_mweo_inverse_real(f_mweo_real, &
00221              f_real_sp_mweo, flm_orig, L, verbosity)
00222         !-------------------------------------------------------------------------
00223         call cpu_time(time_end)
00224         durations_inverse_mweo_real(i_repeat) = time_end - time_start
00225         call cpu_time(time_start)
00226         !-------------------------------------------------------------------------
00227         call ssht_core_mweo_forward_real(flm_syn, &
00228              f_mweo_real, f_real_sp_mweo, L, verbosity)
00229         !-------------------------------------------------------------------------
00230         call cpu_time(time_end)
00231         durations_forward_mweo_real(i_repeat) = time_end - time_start
00232         errors_mweo_real(i_repeat) = maxval(abs(flm_orig(0:L**2-1) - flm_syn(0:L**2-1)))
00233         write(*,'(a,f40.4)') ' duration_inverse (s) =', &
00234              durations_inverse_mweo_real(i_repeat)
00235         write(*,'(a,f40.4)') ' duration_forward (s) =', &
00236              durations_forward_mweo_real(i_repeat)
00237         write(*,'(a,e40.5)') ' error                =', &
00238              errors_mweo_real(i_repeat)
00239         write(*,*)
00240 
00241         !=========================================================================
00242         ! MW real spin=0
00243         write(*,'(a,i2)') 'MW real spin=0 test no.', i_repeat
00244         flm_orig(0:L**2-1) = cmplx(0d0, 0d0)
00245         flm_syn(0:L**2-1) = cmplx(0d0, 0d0)
00246         call ssht_test_gen_flm_real(L, flm_orig, seed)
00247         call cpu_time(time_start)
00248         !-------------------------------------------------------------------------
00249         call ssht_core_mw_inverse_real(f_mw_real, &
00250              f_real_sp_mw, flm_orig, L, verbosity)
00251         !-------------------------------------------------------------------------
00252         call cpu_time(time_end)
00253         durations_inverse_mw_real(i_repeat) = time_end - time_start
00254         call cpu_time(time_start)
00255         !-------------------------------------------------------------------------
00256         call ssht_core_mw_forward_real(flm_syn, &
00257              f_mw_real, f_real_sp_mw, L, verbosity)
00258         !-------------------------------------------------------------------------
00259         call cpu_time(time_end)
00260         durations_forward_mw_real(i_repeat) = time_end - time_start
00261         errors_mw_real(i_repeat) = maxval(abs(flm_orig(0:L**2-1) - flm_syn(0:L**2-1)))
00262         write(*,'(a,f40.4)') ' duration_inverse (s) =', &
00263              durations_inverse_mw_real(i_repeat)
00264         write(*,'(a,f40.4)') ' duration_forward (s) =', &
00265              durations_forward_mw_real(i_repeat)
00266         write(*,'(a,e40.5)') ' error                =', &
00267              errors_mw_real(i_repeat)
00268         write(*,*)
00269 
00270      end if
00271 
00272      !=========================================================================
00273      ! DH
00274      write(*,'(a,i2)') 'DH test no.', i_repeat
00275      flm_orig(0:L**2-1) = cmplx(0d0, 0d0)
00276      flm_syn(0:L**2-1) = cmplx(0d0, 0d0)
00277      call ssht_test_gen_flm_complex(L, spin, flm_orig, seed)
00278      call cpu_time(time_start)
00279      !-------------------------------------------------------------------------
00280      call ssht_core_dh_inverse(f_dh, flm_orig, L, spin, verbosity)
00281      !-------------------------------------------------------------------------
00282      call cpu_time(time_end)
00283      durations_inverse_dh(i_repeat) = time_end - time_start
00284      call cpu_time(time_start)
00285      !-------------------------------------------------------------------------
00286      call ssht_core_dh_forward(flm_syn, f_dh, L, spin, verbosity)
00287      !-------------------------------------------------------------------------
00288      call cpu_time(time_end)
00289      durations_forward_dh(i_repeat) = time_end - time_start
00290      errors_dh(i_repeat) = maxval(abs(flm_orig(0:L**2-1) - flm_syn(0:L**2-1)))
00291      write(*,'(a,f40.4)') ' duration_inverse (s) =', &
00292           durations_inverse_dh(i_repeat)
00293      write(*,'(a,f40.4)') ' duration_forward (s) =', &
00294           durations_forward_dh(i_repeat)
00295      write(*,'(a,e40.5)') ' error                =', &
00296           errors_dh(i_repeat)
00297      write(*,*)
00298 
00299      !=========================================================================
00300      ! GL
00301      write(*,'(a,i2)') 'GL test no.', i_repeat
00302      flm_orig(0:L**2-1) = cmplx(0d0, 0d0)
00303      flm_syn(0:L**2-1) = cmplx(0d0, 0d0)
00304      call ssht_test_gen_flm_complex(L, spin, flm_orig, seed)
00305      call cpu_time(time_start)
00306      !-------------------------------------------------------------------------
00307      call ssht_core_gl_inverse(f_gl, flm_orig, L, spin, verbosity)
00308      !-------------------------------------------------------------------------
00309      call cpu_time(time_end)
00310      durations_inverse_gl(i_repeat) = time_end - time_start
00311      call cpu_time(time_start)
00312      !-------------------------------------------------------------------------
00313      call ssht_core_gl_forward(flm_syn, f_gl, L, spin, verbosity)
00314      !-------------------------------------------------------------------------
00315      call cpu_time(time_end)
00316      durations_forward_gl(i_repeat) = time_end - time_start
00317      errors_gl(i_repeat) = maxval(abs(flm_orig(0:L**2-1) - flm_syn(0:L**2-1)))
00318      write(*,'(a,f40.4)') ' duration_inverse (s) =', &
00319           durations_inverse_gl(i_repeat)
00320      write(*,'(a,f40.4)') ' duration_forward (s) =', &
00321           durations_forward_gl(i_repeat)
00322      write(*,'(a,e40.5)') ' error                =', &
00323           errors_gl(i_repeat)
00324      write(*,*)
00325 
00326      !=========================================================================
00327      ! MWEO
00328      write(*,'(a,i2)') 'MWEO test no.', i_repeat
00329      flm_orig(0:L**2-1) = cmplx(0d0, 0d0)
00330      flm_syn(0:L**2-1) = cmplx(0d0, 0d0)
00331      call ssht_test_gen_flm_complex(L, spin, flm_orig, seed)
00332      call cpu_time(time_start)
00333      !-------------------------------------------------------------------------
00334      call ssht_core_mweo_inverse(f_mweo, f_sp_mweo, phi_sp_mweo, &
00335           flm_orig, L, spin, verbosity)
00336      !-------------------------------------------------------------------------
00337      call cpu_time(time_end)
00338      durations_inverse_mweo(i_repeat) = time_end - time_start
00339      call cpu_time(time_start)
00340      !-------------------------------------------------------------------------
00341      call ssht_core_mweo_forward(flm_syn, f_mweo, &
00342           f_sp_mweo, phi_sp_mweo, L, spin, verbosity)
00343      !-------------------------------------------------------------------------
00344      call cpu_time(time_end)
00345      durations_forward_mweo(i_repeat) = time_end - time_start
00346      errors_mweo(i_repeat) = maxval(abs(flm_orig(0:L**2-1) - flm_syn(0:L**2-1)))
00347      write(*,'(a,f40.4)') ' duration_inverse (s) =', &
00348           durations_inverse_mweo(i_repeat)
00349      write(*,'(a,f40.4)') ' duration_forward (s) =', &
00350           durations_forward_mweo(i_repeat)
00351      write(*,'(a,e40.5)') ' error                =', &
00352           errors_mweo(i_repeat)
00353      write(*,*)
00354 
00355      !=========================================================================
00356      ! MW
00357      write(*,'(a,i2)') 'MW test no.', i_repeat
00358      flm_orig(0:L**2-1) = cmplx(0d0, 0d0)
00359      flm_syn(0:L**2-1) = cmplx(0d0, 0d0)
00360      call ssht_test_gen_flm_complex(L, spin, flm_orig, seed)
00361      call cpu_time(time_start)
00362      !-------------------------------------------------------------------------
00363      call ssht_core_mw_inverse(f_mw, f_sp_mw, phi_sp_mw, &
00364           flm_orig, L, spin, verbosity)
00365      !-------------------------------------------------------------------------
00366      call cpu_time(time_end)
00367      durations_inverse_mw(i_repeat) = time_end - time_start
00368      call cpu_time(time_start)
00369      !-------------------------------------------------------------------------
00370      call ssht_core_mw_forward(flm_syn, f_mw, &
00371           f_sp_mw, phi_sp_mw, L, spin, verbosity)
00372      !-------------------------------------------------------------------------
00373      call cpu_time(time_end)
00374      durations_forward_mw(i_repeat) = time_end - time_start
00375      errors_mw(i_repeat) = maxval(abs(flm_orig(0:L**2-1) - flm_syn(0:L**2-1)))
00376      write(*,'(a,f40.4)') ' duration_inverse (s) =', &
00377           durations_inverse_mw(i_repeat)
00378      write(*,'(a,f40.4)') ' duration_forward (s) =', &
00379           durations_forward_mw(i_repeat)
00380      write(*,'(a,e40.5)') ' error                =', &
00381           errors_mw(i_repeat)
00382      write(*,*)
00383 
00384   end do
00385 
00386   ! Print summary.
00387   write(*,'(a)') '==============================================================='
00388   write(*,'(a)') 'Summary'
00389   write(*,*)
00390   write(*,'(a,i40)') 'N_repeat              =', N_repeat
00391   write(*,'(a,i40)') 'L                     =', L
00392   write(*,'(a,i40)') 'spin                  =', spin
00393   write(*,*)
00394 
00395   if (spin == 0) then
00396      write(*,'(a,i2)') 'DH real'
00397      write(*,'(a,f30.5)') ' Average forward transform time =', &
00398           sum(durations_forward_dh_real(0:N_repeat-1)) / real(N_repeat)
00399      write(*,'(a,f30.5)') ' Average inverse transform time =', &
00400           sum(durations_inverse_dh_real(0:N_repeat-1)) / real(N_repeat)
00401      write(*,'(a,e30.5)') ' Average max error              =', &
00402           sum(errors_dh_real(0:N_repeat-1)) / real(N_repeat)
00403      write(*,*)
00404 
00405      write(*,'(a,i2)') 'GL real'
00406      write(*,'(a,f30.5)') ' Average forward transform time =', &
00407           sum(durations_forward_gl_real(0:N_repeat-1)) / real(N_repeat)
00408      write(*,'(a,f30.5)') ' Average inverse transform time =', &
00409           sum(durations_inverse_gl_real(0:N_repeat-1)) / real(N_repeat)
00410      write(*,'(a,e30.5)') ' Average max error              =', &
00411           sum(errors_gl_real(0:N_repeat-1)) / real(N_repeat)
00412      write(*,*)
00413 
00414      write(*,'(a,i2)') 'MWEO real'
00415      write(*,'(a,f30.5)') ' Average forward transform time =', &
00416           sum(durations_forward_mweo_real(0:N_repeat-1)) / real(N_repeat)
00417      write(*,'(a,f30.5)') ' Average inverse transform time =', &
00418           sum(durations_inverse_mweo_real(0:N_repeat-1)) / real(N_repeat)
00419      write(*,'(a,e30.5)') ' Average max error              =', &
00420           sum(errors_mweo_real(0:N_repeat-1)) / real(N_repeat)
00421      write(*,*)
00422 
00423      write(*,'(a,i2)') 'MW real'
00424      write(*,'(a,f30.5)') ' Average forward transform time =', &
00425           sum(durations_forward_mw_real(0:N_repeat-1)) / real(N_repeat)
00426      write(*,'(a,f30.5)') ' Average inverse transform time =', &
00427           sum(durations_inverse_mw_real(0:N_repeat-1)) / real(N_repeat)
00428      write(*,'(a,e30.5)') ' Average max error              =', &
00429           sum(errors_mw_real(0:N_repeat-1)) / real(N_repeat)
00430      write(*,*)
00431   end if
00432 
00433   write(*,'(a,i2)') 'DH'
00434   write(*,'(a,f30.5)') ' Average forward transform time =', &
00435        sum(durations_forward_dh(0:N_repeat-1)) / real(N_repeat)
00436   write(*,'(a,f30.5)') ' Average inverse transform time =', &
00437        sum(durations_inverse_dh(0:N_repeat-1)) / real(N_repeat)
00438   write(*,'(a,e30.5)') ' Average max error              =', &
00439        sum(errors_dh(0:N_repeat-1)) / real(N_repeat)
00440   write(*,*)
00441 
00442   write(*,'(a,i2)') 'GL'
00443   write(*,'(a,f30.5)') ' Average forward transform time =', &
00444        sum(durations_forward_gl(0:N_repeat-1)) / real(N_repeat)
00445   write(*,'(a,f30.5)') ' Average inverse transform time =', &
00446        sum(durations_inverse_gl(0:N_repeat-1)) / real(N_repeat)
00447   write(*,'(a,e30.5)') ' Average max error              =', &
00448        sum(errors_gl(0:N_repeat-1)) / real(N_repeat)
00449   write(*,*)
00450 
00451   write(*,'(a,i2)') 'MWEO'
00452   write(*,'(a,f30.5)') ' Average forward transform time =', &
00453        sum(durations_forward_mweo(0:N_repeat-1)) / real(N_repeat)
00454   write(*,'(a,f30.5)') ' Average inverse transform time =', &
00455        sum(durations_inverse_mweo(0:N_repeat-1)) / real(N_repeat)
00456   write(*,'(a,e30.5)') ' Average max error              =', &
00457        sum(errors_mweo(0:N_repeat-1)) / real(N_repeat)
00458   write(*,*)
00459 
00460   write(*,'(a,i2)') 'MW'
00461   write(*,'(a,f30.5)') ' Average forward transform time =', &
00462        sum(durations_forward_mw(0:N_repeat-1)) / real(N_repeat)
00463   write(*,'(a,f30.5)') ' Average inverse transform time =', &
00464        sum(durations_inverse_mw(0:N_repeat-1)) / real(N_repeat)
00465   write(*,'(a,e30.5)') ' Average max error              =', &
00466        sum(errors_mw(0:N_repeat-1)) / real(N_repeat)
00467   write(*,*)
00468 
00469   ! Deallocate memory.
00470   deallocate(flm_orig, flm_syn)
00471   deallocate(f_dh, f_gl, f_mweo, f_mw)
00472   deallocate(f_dh_real, f_gl_real, f_mweo_real, f_mw_real)
00473 
00474 end program ssht_test
00475 
00476 
00477 !--------------------------------------------------------------------------
00478 ! ssht_test_gen_flm_real
00479 !
00480 !! Generate random spherical harmonic coefficients of a real spin=0
00481 !! signal.
00482 !!
00483 !! Variables:
00484 !!   - L: Harmonic band-limit [input].
00485 !!   - flm(0:L**2-1): Random spherical harmonic coefficients generated 
00486 !!     [output].
00487 !!   - seed: Integer seed required for random number generator [input].
00488 !
00489 !! @author J. D. McEwen
00490 !! @version 0.1 October 2010
00491 ! 
00492 ! Revisions:
00493 !   October 2010 - Jason McEwen
00494 !--------------------------------------------------------------------------
00495 
00496 subroutine ssht_test_gen_flm_real(L, flm, seed)
00497 
00498   use ssht_types_mod, only: dpc
00499   use ssht_sampling_mod
00500 
00501   implicit none
00502 
00503   interface 
00504      function ran2_dp(idum)
00505        use ssht_types_mod, only: dp
00506        real(dp) :: ran2_dp
00507        integer :: idum
00508      end function ran2_dp
00509   end interface
00510 
00511   integer, intent(in) :: L
00512   complex(dpc), intent(out) :: flm(0:L**2-1)
00513   integer, intent(in) :: seed
00514 
00515   integer :: el, m, ind
00516   complex(dpc) :: tmp
00517 
00518   flm(0:L**2-1) = cmplx(0d0, 0d0)
00519   do el = 0,L-1
00520      m = 0
00521      call ssht_sampling_elm2ind(ind, el, m)  
00522      tmp = cmplx(2d0*ran2_dp(seed)-1d0, 0d0)
00523      flm(ind) = tmp
00524 
00525      do m = 1,el
00526         call ssht_sampling_elm2ind(ind, el, m)  
00527         tmp = cmplx(2d0*ran2_dp(seed)-1d0, 2d0*ran2_dp(seed)-1d0)
00528         flm(ind) = tmp
00529         call ssht_sampling_elm2ind(ind, el, -m)  
00530         flm(ind) = (-1)**m * conjg(tmp)
00531      end do
00532 
00533   end do
00534 
00535 end subroutine ssht_test_gen_flm_real
00536 
00537 
00538 !--------------------------------------------------------------------------
00539 ! ssht_test_gen_flm_complex
00540 !
00541 !! Generate random spherical harmonic coefficients of a complex signal.
00542 !!
00543 !! Variables:
00544 !!   - L: Harmonic band-limit [input].
00545 !!   - spin: Spin number [input].
00546 !!   - flm(0:L**2-1): Random spherical harmonic coefficients generated 
00547 !!     [output].
00548 !!   - seed: Integer seed required for random number generator [input].
00549 !
00550 !! @author J. D. McEwen
00551 !! @version 0.1 October 2010
00552 ! 
00553 ! Revisions:
00554 !   October 2010 - Jason McEwen
00555 !--------------------------------------------------------------------------
00556 
00557 subroutine ssht_test_gen_flm_complex(L, spin, flm, seed)
00558 
00559   use ssht_types_mod, only: dpc
00560   use ssht_sampling_mod
00561 
00562   implicit none
00563 
00564   interface 
00565      function ran2_dp(idum)
00566        use ssht_types_mod, only: dp
00567        real(dp) :: ran2_dp
00568        integer :: idum
00569      end function ran2_dp
00570   end interface
00571 
00572   integer, intent(in) :: L
00573   integer, intent(in) :: spin
00574   complex(dpc), intent(out) :: flm(0:L**2-1)
00575   integer, intent(in) :: seed
00576 
00577   integer :: ind, ind_lo, el, m
00578 
00579   flm(0:L**2-1) = cmplx(0d0, 0d0)
00580 
00581   call ssht_sampling_elm2ind(ind_lo, abs(spin), 0)
00582   do ind = ind_lo,L**2 - 1
00583 !     call ssht_core_ind2elm(el, m, ind)  
00584 !     if(el >= abs(spin)) then
00585         flm(ind) = cmplx(2d0*ran2_dp(seed)-1d0, 2d0*ran2_dp(seed)-1d0)
00586 !     end if
00587   end do
00588 
00589 end subroutine ssht_test_gen_flm_complex
00590 
00591 
00592 !--------------------------------------------------------------------------
00593 ! ran2_dp
00594 !
00595 !! Generate uniform deviate in range [0,1) given seed.
00596 !! (Using double precision.)
00597 !!
00598 !! Notes: 
00599 !!   - Uniform deviate (Num rec 1992, chap 7.1), original routine
00600 !!     said to be 'perfect'.
00601 !!
00602 !! Variables:
00603 !!   - idum: Seed.
00604 !
00605 !! @author J. D. McEwen
00606 !! @version 0.1 May 2007
00607 ! 
00608 ! Revisions:
00609 !   May 2007 - Integrated by Jason McEwen
00610 !--------------------------------------------------------------------------
00611 
00612 function ran2_dp(idum)
00613 
00614   use ssht_types_mod, only: dp
00615 
00616   implicit none
00617 
00618   INTEGER :: idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
00619   REAL(dp) :: ran2_dp,AM,EPS,RNMX
00620   PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1, &
00621        & IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791, &
00622        & NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
00623   INTEGER :: idum2,j,k,iv(NTAB),iy
00624   SAVE iv,iy,idum2
00625   DATA idum2/123456789/, iv/NTAB*0/, iy/0/
00626   if (idum.le.0) then
00627      idum=max(-idum,1)
00628      idum2=idum
00629      do j=NTAB+8,1,-1
00630         k=idum/IQ1
00631         idum=IA1*(idum-k*IQ1)-k*IR1
00632         if (idum.lt.0) idum=idum+IM1
00633         if (j.le.NTAB) iv(j)=idum
00634      end do
00635      iy=iv(1)
00636   end if
00637   k=idum/IQ1
00638   idum=IA1*(idum-k*IQ1)-k*IR1
00639   if (idum.lt.0) idum=idum+IM1
00640   k=idum2/IQ2
00641   idum2=IA2*(idum2-k*IQ2)-k*IR2
00642   if (idum2.lt.0) idum2=idum2+IM2
00643   j=1+iy/NDIV
00644   iy=iv(j)-idum2
00645   iv(j)=idum
00646   if(iy.lt.1)iy=iy+IMM1
00647   ran2_dp=min(AM*iy,RNMX)
00648   return
00649 
00650 end function ran2_dp
00651