Блог

Я работаю с некоторым старым кодом и наткнулся на подпрограмму «CalcNorm». Я попытался определить, какой тип нормы вычисляет эта функция, но я не могу точно указать на это. Это выглядит как 2-норма с некоторыми исключениями. Кто-нибудь узнает этот алгоритм? Комментарии в коде были сделаны мной и не обязательно могут быть правильными.

Спасибо,

Дэйв

 subroutine CalcNorm(n,x,norm)

  implicit none
  !Number of dimensions in the vector x
  integer(kind = int4),intent(in) :: n
  !Input vector
  real(kind = real8), dimension(n), intent(in) :: x
  !Norm of x
  real(kind = real8), intent(out) :: norm
  !Counter
  integer(kind = int4) i
  !Absolute value of x
  real(kind = real8) :: xabs
  !Largest x > rgiant
  real(kind = real8) :: x1max
  !Largest x > 0 but < rdwarf
  real(kind = real8) :: x3max
  !Sample Quantities
  real(kind=real8) :: s1
  !Variance of x where rdwarf < xabs < rgiant
  real(kind=real8) :: s2
  real(kind=real8) :: s3
  !Limits for the three methods for collecting samples
  real(kind=real8), parameter :: rdwarf = 3.834d-20, rgiant = dsqrt(1.304d19)
  real(kind=real8), parameter :: adwarf = dsqrt(rdwarf)

  s1 = 0.d0
  s2 = 0.d0
  s3 = 0.d0
  x1max = 0.d0
  x3max = 0.d0
  do i = 1, n
     !Determine relative size of x(i)
     xabs = dabs(x(i))
     !The most likely case
     if (xabs > rdwarf .and. xabs < rgiant) then
        s2 = s2   xabs**2
     !x is tiny, quite likely
     else if (xabs <= rdwarf) then
        !Both of these happen quite often
        if (xabs <= x3max) then
           if(x3max /= 0.d0) s3 = s3   (xabs/x3max)**2
        else
           s3 = 1.d0   s3*(x3max/xabs)**2
           x3max = xabs
        end if
     !These last two cases are the rarest and aren't typically envoked
     !xabs >= rgiant .and. xabs < x1max
        !This won't happen the first time xabs >= rgiant
     else if (xabs <= x1max) then
        s1 = s1   (xabs/x1max)**2
     else !xabs is huge: xabs >= rgiant .and. xabs > x1max
        s1 = 1.d0   s1*(x1max/xabs)**2
        x1max = xabs
     end if
  end do

  !all xabs < rgiant
  if (s1 == 0.d0) then

     !All xabs < rdwarf
     if (s2 == 0.d0) then
        norm = x3max*dsqrt(s3)

     !At least one xabs >= rdwarf and sum of xabs**2 >= adwarf
     elseif (s2 >= adwarf) then
        norm = dsqrt(s2 x3max**2*s3)

     !Any xabs > rdwarf and sum of xabs**2 < adwarf
     else !s2 > 0 .and. < adwarf
        !selective 2-norm
        norm = dsqrt(s2 x3max**2)

     end if
  !At least one xabs > rgiant
  else
     norm = dsqrt(s1*x1max**2 s2)
  end if 

  return
end subroutine CalcNorm
  

Спасибо @Jim Mischel за помощь. Я смог найти исходный код f77, просмотрев enorm (http://www.netlib.no/netlib/minpack/enorm.f). Это процедура евклидовой нормы.