qrdcmp Subroutine

public subroutine qrdcmp(a, c, d, sing)

Arguments

Type IntentOptional Attributes Name
real(kind=8), intent(inout), dimension(:,:) :: a
real(kind=8), intent(out), dimension(:) :: c
real(kind=8), intent(out), dimension(:) :: d
logical, intent(out) :: sing

Calls

proc~~qrdcmp~~CallsGraph proc~qrdcmp qrdcmp proc~assert_eq4 assert_eq4 proc~qrdcmp->proc~assert_eq4 proc~outerprod outerprod proc~qrdcmp->proc~outerprod proc~vabs vabs proc~qrdcmp->proc~vabs

Called by

proc~~qrdcmp~~CalledByGraph proc~qrdcmp qrdcmp proc~broyden1~2 broyden1 proc~broyden1~2->proc~qrdcmp

Source Code

  subroutine qrdcmp(a,c,d,sing)
    real(8), dimension(:,:), intent(inout) :: a
    real(8), dimension(:), intent(out) :: c,d
    logical, intent(out) :: sing
    integer :: k,n
    real(8) :: scale,sigma
    n=assert_eq4(size(a,1),size(a,2),size(c),size(d),'qrdcmp')
    sing=.false.
    do k=1,n-1
       scale=maxval(abs(a(k:n,k)))
       if (scale == 0.0) then
          sing=.true.
          c(k)=0.0
          d(k)=0.0
       else
          a(k:n,k)=a(k:n,k)/scale
          sigma=sign(vabs(a(k:n,k)),a(k,k))
          a(k,k)=a(k,k)+sigma
          c(k)=sigma*a(k,k)
          d(k)=-scale*sigma
          a(k:n,k+1:n)=a(k:n,k+1:n)-outerprod(a(k:n,k),&
               matmul(a(k:n,k),a(k:n,k+1:n)))/c(k)
       end if
    end do
    d(n)=a(n,n)
    if (d(n) == 0.0) sing=.true.
  end subroutine qrdcmp