var/da/da_biascorr_airmass/tred2.f90

   1 SUBROUTINE tred2(a,n,np,d,e)
   2 !USE nrutil, ONLY : assert_eq,outerprod
   3 IMPLICIT NONE
   4 INTEGER, PARAMETER :: sp = 8
   5 INTEGER n,np
   6 REAL(SP), DIMENSION(np,np), INTENT(INOUT) :: a
   7 REAL(SP), DIMENSION(np), INTENT(OUT) :: d,e
   8 !LOGICAL(LGT), OPTIONAL, INTENT(IN) :: novectors
   9 INTEGER :: i,j,l,k
  10 REAL(SP) :: f,g,h,hh,scale
  11 REAL(SP), DIMENSION(size(a,1)) :: gg
  12 !LOGICAL(LGT) :: yesvec
  13 !n=assert_eq(size(a,1),size(a,2),size(d),size(e),’tred2’)
  14 !if (present(novectors)) then
  15 !yesvec=.not. novectors
  16 !else
  17 !yesvec=.true.
  18 !end if
  19 do i=n,2,-1
  20 l=i-1
  21 h=0.0
  22 if (l > 1) then
  23 scale=sum(abs(a(i,1:l)))
  24 if (scale == 0.0) then
  25 e(i)=a(i,l)
  26 else
  27 a(i,1:l)=a(i,1:l)/scale
  28 h=sum(a(i,1:l)**2)
  29 f=a(i,l)
  30 g=-sign(sqrt(h),f)
  31 e(i)=scale*g
  32 h=h-f*g
  33 a(i,l)=f-g
  34 a(1:l,i)=a(i,1:l)/h
  35 do j=1,l
  36 e(j)=(dot_product(a(j,1:j),a(i,1:j)) &
  37 +dot_product(a(j+1:l,j),a(i,j+1:l)))/h
  38 end do
  39 f=dot_product(e(1:l),a(i,1:l))
  40 hh=f/(h+h)
  41 e(1:l)=e(1:l)-hh*a(i,1:l)
  42
  43 do j=1,l
  44 a(j,1:j)=a(j,1:j)-a(i,j)*e(1:j)-e(j)*a(i,1:j)
  45 end do
  46 end if
  47 else
  48 e(i)=a(i,l)
  49 end if
  50 d(i)=h
  51 end do
  52 d(1)=0.0
  53 e(1)=0.0
  54 do i=1,n
  55
  56 l=i-1
  57 if (d(i) /= 0.0) then
  58 gg(1:l)=matmul(a(i,1:l),a(1:l,1:l))
  59 !a(1:l,1:l)=a(1:l,1:l)-outerprod(a(1:l,i),gg(1:l))
  60 do j=1,l
  61 do k=1,l
  62    a(k,j)=a(k,j)-gg(j)*a(k,i)
  63 end do
  64 end do
  65 end if
  66 d(i)=a(i,i)
  67 a(i,i)=1.0
  68 a(i,1:l)=0.0
  69 a(1:l,i)=0.0
  70 !else
  71 !d(i)=a(i,i)
  72 !end if
  73 end do
  74 END SUBROUTINE tred2