Последний раз редактировалось aitap 18.05.2016, 00:42, всего редактировалось 1 раз.
Попытки рассчитать  для больших по модулю аргументов  (куда входит  ) при помощи библиотеки PFQ не увенчались успехом ("ERROR - VALUE OF EXPONENT REQUIRED FOR SUMMATION WAS LARGER THAN THE MAXIMUM MACHINE EXPONENT" и "ERROR IN IPREMAX--DID NOT FIND MAXIMUM EXPONENT" в разных ситуациях); я пошёл на сделку с совестью и проигнорировал возвращаемый из внутреннего интеграла  код ошибки 2 (даже при  , как ни странно), поскольку несмотря на это, график полученного интеграла от  выглядит точь-в-точь, как экспериментальные спектры. При помощи Maxima взял частные производные от подинтегрального выражения по параметрам и написал для библиотеки MINPACK минимизируемую функцию (  по всем точкам, в которых снят спектр) и якобиан по параметрам:
module holtsmark_errf_param
implicit none
! minpack fits real*8
real (kind=8) :: lambda_0
real (kind=8), dimension(:,:), pointer :: spectrum
end module
module j_deriv_d
use holtsmark_errf_param
use constants
implicit none
real (kind=8) :: z_param, alpha_param
!$omp threadprivate(z_param, alpha_param)
contains
function d_j_d_d(lambda, d, omega, alpha, A_0)
! this is nonsense because integrands are real*4
! but minpack fits real*8
real (kind=8) :: d_j_d_d, omega, alpha, lambda, d, A_0
! quadpack integrates real*4
real (kind=4) :: integral, abserr
integer (kind=4) :: neval, ier
z_param = (lambda - d - lambda_0)/omega
alpha_param = alpha**(4./3.)
call qagi ( &
d_j_d_d_outer_int, &
0.0_4, 1_4, &
epsabs, epsrel, &
integral, &
abserr, neval, ier &
)
d_j_d_d = A_0 * 4 / (pi**2 * omega) * integral
return
end function
function d_j_d_d_outer_int(beta)
use holtsmark, only: beta_inner => beta_param, &
holtsmark_inner => inner_integrand
! quadpack integrates real*4
real (kind=4) :: integral, abserr, d_j_d_d_outer_int, beta
integer (kind=4) :: neval, ier
beta_inner = beta
call qagi( &
holtsmark_inner, &
0.0_4, 1_4, &
epsabs, epsrel, &
integral, &
abserr, neval, ier &
)
d_j_d_d_outer_int = integral &
* beta * (real(z_param,4)-beta**2.*real(alpha_param,4)) &
/ ( ((real(z_param,4)-beta**2.*real(alpha_param,4))**2.+1)**2. )
return
end function
end module
module j_deriv_omega
use holtsmark_errf_param
use constants
implicit none
real (kind=8) :: z_param, alpha_param
!$omp threadprivate(z_param, alpha_param)
contains
function d_j_d_omega(lambda, d, omega, alpha, A_0)
! this is nonsense because integrands are real*4
! but minpack fits real*8
real (kind=8) :: d_j_d_omega, omega, alpha, lambda, d, A_0
! quadpack integrates real*4
real (kind=4) :: integral, abserr
integer (kind=4) :: neval, ier
z_param = (lambda - d - lambda_0)/omega
alpha_param = alpha**(4./3.)
call qagi( &
d_j_d_omega_outer_int, &
0.0_4, 1_4, &
epsabs, epsrel, &
integral, abserr, &
neval, ier &
)
d_j_d_omega = A_0 * 4 * z_param * integral / (pi**2 * omega)
return
end function
function d_j_d_omega_outer_int(beta)
use holtsmark, only: beta_inner => beta_param, &
holtsmark_inner => inner_integrand
! quadpack integrates real*4
real (kind=4) :: integral, abserr, d_j_d_omega_outer_int, beta
integer (kind=4) :: neval, ier
beta_inner = beta
call qagi( &
holtsmark_inner, &
0.0_4, 1_4, &
epsabs, epsrel, &
integral, &
abserr, neval, ier &
)
d_j_d_omega_outer_int = integral * beta &
* (real(z_param,4) - beta**2 * real(alpha_param,4)) &
/ ( 1 + (real(z_param,4) - beta**2 * real(alpha_param,4))**2 )**2
return
end function
end module
module j_deriv_alpha
use holtsmark_errf_param
use constants
implicit none
real (kind=8) :: z_param, alpha_param
!$omp threadprivate(z_param,alpha_param)
contains
function d_j_d_alpha(lambda, d, omega, alpha, A_0)
! this is nonsense because integrands are real*4
! but minpack fits real*8
real (kind=8) :: d_j_d_alpha, omega, alpha, lambda, d, A_0
! quadpack integrates real*4
real (kind=4) :: integral, abserr
integer (kind=4) :: neval, ier
z_param = (lambda - d - lambda_0)/omega
alpha_param = alpha**(4./3.)
call qagi( &
d_j_d_alpha_outer_int, &
0.0_4, 1_4, &
epsabs, epsrel, &
integral, abserr, &
neval, ier &
)
d_j_d_alpha = A_0 * 16./(3. * PI**2) * alpha**(1./3.) * integral
return
end function
function d_j_d_alpha_outer_int(beta)
use holtsmark, only: beta_inner => beta_param, &
holtsmark_inner => inner_integrand
! quadpack integrates real*4
real (kind=4) :: integral, abserr, d_j_d_alpha_outer_int, beta
integer (kind=4) :: neval, ier
beta_inner = beta
call qagi( &
holtsmark_inner, &
0.0_4, 1_4, &
epsabs, epsrel, &
integral, &
abserr, neval, ier &
)
d_j_d_alpha_outer_int = integral * beta**3 &
* (real(z_param,4) - beta**2 * real(alpha_param,4)) &
/ ( 1 + (real(z_param,4) - beta**2 * real(alpha_param,4))**2 )**2
return
end function
end module
module holtsmark_errf
!use, intrinsic :: ieee_arithmetic, only: ieee_value !missing until gfortran 5, sigh
use holtsmark_errf_param
implicit none
contains
subroutine error_function(m, n, x, fvec, fjac, ldfjac, iflag)
use j_deriv_d
use j_deriv_omega
use j_deriv_alpha
use holtsmark
implicit none
integer (kind = 4) ldfjac
integer (kind = 4) n, m ! parameters, points
real (kind=8) fjac(ldfjac,n)
real (kind=8) fvec(m)
integer (kind = 4) iflag ! 1 => functions->fvec; 2 => jacobian->fjac
real (kind=8) x(n)
integer (kind=4) index
print*, 'x:', x(1:4)
if (iflag == 1) then
if ((x(2)<0).or.(x(3)<0).or.(x(4)<0)) then
! omega, alpha, A0 should not be <0; let's not waste time on that
fvec(1:m) = sqrt(-real(n))!ieee_value(fvec,ieee_quiet_nan)
return
end if
!$omp parallel do shared(fvec,spectrum) private(index) schedule(guided)
do index=1,m
fvec(index) = j_holtsmark(spectrum(index,1),lambda_0,x(1),x(2),x(3),x(4))-spectrum(index,2)
!print*, 'err(',index,',',spectrum(index,1),'):', fvec(index)
end do
!$omp end parallel do
print*, 'sumsq(err):', sum(fvec**2)
else if (iflag == 2) then
if ((x(2)<0).or.(x(3)<0).or.(x(4)<0)) then
! omega, alpha, A0 should not be <0; let's not waste time on that
fjac(1:ldfjac,1:n) = sqrt(-real(n))!ieee_value(fvec,ieee_quiet_nan)
return
end if
!$omp parallel do shared(fjac,spectrum,x) private(index) schedule(dynamic)
do index=1,m
!$omp parallel sections
!$omp section
fjac(index,1) = d_j_d_d(spectrum(index,1),x(1),x(2),x(3),x(4)) ! d
!$omp section
fjac(index,2) = d_j_d_omega(spectrum(index,1),x(1),x(2),x(3),x(4)) ! omega
!$omp section
fjac(index,3) = d_j_d_alpha(spectrum(index,1),x(1),x(2),x(3),x(4)) ! alpha
!$omp section
fjac(index,4) = j_holtsmark(spectrum(index,1),lambda_0,x(1),x(2),x(3),x(4))/x(4) ! A_0
!$omp end parallel sections
!print*, 'J(',index,',',spectrum(index,1),'):', fjac(index,:)
end do
!$omp end parallel do
print*, 'sumsq(J):', sum(fjac**2)
end if
return
end subroutine
end module
(процедуру error_function я передаю в LMDER1 библиотеки MINPACK после считывания спектра и начальных значений в память) После добавления OpenMP и правильной расстановки атрибута threadprivate по переменным модулей фит даже стал иногда заканчиваться за ночь. Больше всего мне помог совет доброго человека равномерно отобрать из 1390 экспериметнальных точек  штук и предварительно фитить только по ним; это позволило закончить за ночь все расчёты. Зато оказалось, что базовая линия имеет вид наклонной, а не горизонтальной прямой, так что её придётся вернуть в аппроксимацию в виде слагаемого  , благо, частные производные по новым параметрам будут проще некуда,  и  соответственно. Параметр  пришлось тоже вернуть, потому что моя теоретическая функция нормирована на  , а с экспериментальными данными такое с достаточной точностью сделать не удаётся. На данный момент хороший фит (единственный, окончания которого по 1390 точкам я дождался) выглядит вот так:  а плохой (нуждающийся в нормальной поправке на базовую линию) - так:  Ещё один путь, которым я попробую ускорить расчёты - табулировать  через определённый интервал  (пока не решил, какой); вне пределов табулирования возвращать 0, а внутри - линейно интерполировать.
|