算法识别(规范)



我正在处理一些旧代码,我偶然发现了一个"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)。这是一个欧氏范数例程。

相关内容

最新更新