fortran 固有値を求めるプログラム

質問

２x２の実行列の固有値を求めるモジュール関数の中で、d > 0 のときに出てくる　sign(squt(d), -b)の値が示す意味と　eval(2) = cmplx(c/e,0.0d0) でなぜ c/eの値が入るのかがわかりません。教えて下さい。よろしくお願いします。
プログラム
module subprogs
implicit none
contains
function eval2x2mat(a) result(eval)
real(8), intent(in) :: a(:,:)
complex(8) eval(2)
real(8) b, c, d, e
if (size(a,1) /= size(a,2)) stop ' not square '
if (size(a,1) /= 2) stop ' not 2x2 matrix '
b = -0.5d0*(a(1,1)+a(2,2))
c = a(1,1)*a(2,2)-a(1,2)*a(2,1)
d = b**2-c
if ( d < 0.0d0 ) then
eval(1) = cmplx(-b,sqrt(-d))
eval(2) = conjg(eval(1))
else if ( d > 0.0d0 ) then
e = -b+sign(sqrt(d),-b)　←ここの部分
eval(1) = cmplx(e,0.0d0)
eval(2) = cmplx(c/e,0.0d0)
else　　　　　　　　　　↑ここの部分
eval(1) = cmplx(-b,0.0d0)
eval(2) = eval(1)
endif
end function eval2x2mat
end module subprogs

program main
use subprogs
implicit none
real(8), allocatable :: a(:,:)
integer :: n = 2
allocate(a(n,n))
a(1,1:2) = (/-1,1/)
a(2,1:2) = (/-1,-1/)
write(*,*) a(:,:), eval2x2mat(a)
end program main

回答No.1

f272 · ベストアンサー

> sign(sqrt(d), -b)の値が示す意味

これは2次方程式の判別式が正のときです。
このとき-b+sqrt(d)と-b-sqrt(d)のどちらも求める解ですが，
-b+sqrt(d)と-b-sqrt(d)のどちらが絶対値が大きいかを考えてください。
-bが正なら-b+sqrt(d)の絶対値が大きく，-bが負なら-b-sqrt(d)の絶対値が大きくことが分かるでしょう。
だから絶対値の大きい方を第1の解としています。

> eval(2) = cmplx(c/e,0.0d0) でなぜ c/eの値が入るのか

第2の解は，第1の解eが分かっていれば簡単です。2つの解の積が定数項のcですね。だからc/eが第2の解になります。
なぜこのように解の絶対値によって解の順序を決めているかと言えば，計算精度を高めるためです。
例えば
a(1,1:2) = (/1.0d15,1.0d15+1.0d0/)
a(2,1:2) = (/1.0d0,1.0d0/)
でこのような処理をした場合としない場合を比較してください。

お礼

なるほど、sign(sqrt(d),-b)の値の示す意味理解できました。絶対値の大きい方を第1解にしてるんですね。
それと、c/eの計算の意味も納得しました。普段あまりそのような解の求め方に馴染みがなかったのでわかってませんでした。ありがとうございました。