Fortran プログラミングガイド

第 6 章

浮動小数点演算

この章では、浮動小数点演算について説明し、数値計算の誤差を回避および検出するための方針を示します。

SPARC および x86 プロセッサ上の浮動小数点計算に関する詳細は、『数値計算ガイ
ド』を参照してください。

はじめに

SPARC および x86 上のサンの浮動小数点環境は、『IEEE Standard 754 for Binary Floating Point Arithmetic』で指定された演算モデルを実装しています。この環境では、頑丈な、高性能の、移植可能な数値アプリケーションを開発できます。また、数値プログラムによる異常な動作を調査するためのツールも提供します。

数値プログラムでは、計算誤差を引き起こす次のような潜在的な原因がたくさんあります。

• 計算モデルが間違っている
• 使用されているアルゴリズムが数値的に不安定である
• データが正確でない
• ハードウェアが予測できない結果を生み出す

間違った数値計算の誤差の原因を見つけることは大変困難です。市販の検査済みライブラリパッケージを使用することで、コーディングの誤りの可能性を減らすことはできます。優れたアルゴリズムを選択することも重要な問題です。また、コンピュータの特性を考慮した優れた算術演算法を使用することも重要です。

この章では、数値誤差の解析については説明していません。ここでは、Sun WorkShop Fortran コンパイラによって実装された IEEE 浮動小数点モデルを紹介します。

IEEE 浮動小数点演算

IEEE 演算は、無効演算、ゼロ除算、オーバーフロー、アンダーフロー、結果不正確などの問題を発生させる算術演算を取り扱う、比較的新しい方式です。丸め、ゼロに近い数の扱い方、その機種で取り扱える最大数に近い数の扱い方に違いがあります。

IEEE 規格は、例外、丸め、精度のユーザー処理をサポートします。その結果、IEEE 規格は、区間演算と変則性の診断をサポートします。IEEE 規格 754 は、exp や cos などの基本関数の標準化、高精度の演算、数値計算と記号代数計算の結合を実現します。

IEEE 演算では、他の浮動小数点演算と比べると、ユーザーが計算を大きく制御できます。IEEE 規格は、数値的に精練された移植性のあるプログラムを作成する作業を簡単にします。浮動小数点演算についての多くの質問は、基本的な数の演算に関連します。たとえば、次のようなものがあります。

• 無限に正確な結果をコンピュータのハードウェア内で表現できない場合、演算結果はどうなるか
• 乗算や加算のような基本演算は可換か

他のクラスの質問は、浮動小数点例外や例外処理に関連するもので、たとえば、次のようなものがあります。

• 同符号の 2 つの巨大な数を乗算するとどうなるか
• ゼロ以外の値をゼロで除算するとどうなるか
• ゼロをゼロで除算するとどうなるか

旧式の演算モデルでは、最初のクラスの質問は期待された答えにならず、2 番目のクラスの例外ケースはすべて同じ結果になります。つまり、プログラムは当該箇所で異常終了するか、演算は続行されるが結果は意味のないものとなります。

IEEE 規格は、演算が数学的に期待される結果を、期待される特性で出すことを保証します。また、ユーザーが特に他を選択しない限り、例外ケースが指定された結果を出すことも保証します。

たとえば、例外値 +Inf、-Inf、NaN は次のように直感的に理解できます。

big*big = +Inf 正の無限大 big*(-big) = -Inf 負の無限大 num/0.0 = +Inf num > 0.0 のとき num/0.0 = -Inf num < 0.0 のとき 0.0/0.0 = NaN 非数

また、次のような 5 つの種類の浮動小数点例外も発生します。

• 無効演算 − 数学的に定義できない演算。0.0/0.0、sqrt(-1.0)、log(-37.8) など。
• ゼロ除算 − 除数がゼロで、被除数が有限かつゼロでない数。9.9/0.0 など。
• オーバーフロー − 指数の範囲を超える結果を出す演算。
  MAXDOUBLE+0.0000000000001e308 など。
• アンダーフロー − 通常の数として表現できないほど小さな結果を出す演算。
  MINDOUBLE * MINDOUBLE など。
• 結果不正確 − 無限に正確に表現できない結果を出す演算。2.0/3.0、log(1.1)、0.1 の入力など。

IEEE 規格の実装については、『数値計算ガイド』を参照してください。

-ftrap=mode コンパイラオプション

-ftrap=mode オプションは、浮動小数点例外用のトラップを有効にします。ieee_handler() 呼び出しによってシグナルハンドラが設定されていない場合、例外はプログラムを終了し、メモリーダンプコアファイルを作成します。このコンパイラオプションに関する詳細は、『Fortran ユーザーズガイド』を参照してください。たとえば、オーバーフロー、ゼロ除算、無効演算のトラップを有効にするには、-ftrap=common を付けてコンパイルします。

---

注 - トラップを有効にするには、アプリケーションの主プログラムを -ftrap= を付けてコンパイルする必要があります。

---

浮動小数点演算の例外と Fortran

f77 によってコンパイルされたプログラムは、プログラムの終了時、発生した浮動小数点演算の例外リストを自動的に表示します。一般的には、無効演算、ゼロ除算、オーバーフローのいずれかの例外が発生すると、メッセージが表示されます。実際のプログラムでは結果不正確の例外は多く発生するため、結果不正確の例外のメッセージは表示されません。

f95 プログラムは、プログラムの終了時、例外を自動的に報告しません。
ieee_retrospective(3M) への明示的な呼び出しが必要です。

ieee_flags() で例外ステータスフラグをクリアーすると、メッセージの一部また全部を表示させないようにできます。これは、通常はプログラムの最後で行います。

例外処理

IEEE 規格準拠の例外処理は、SPARC および x86 プロセッサ上ではデフォルトで行われます。ただし、浮動小数点例外の検出と、浮動小数点例外に対するシグナル (SIGFPE) の生成には、違いがあります。

浮動小数点演算中にトラップしていない例外が発生すると、IEEE 規格に従って、次の 2 つの処理が行われます。

• システムはデフォルトの結果を返します。たとえば、0/0 (演算不可能) の場合は NaN を結果として返します。
• 例外が発生したことを示すフラグが設定されます。たとえば、0/0 (演算不可能) の場合は「無効演算」のフラグが設定されます。

浮動小数点演算の例外をトラップする

浮動小数点演算の例外を処理する方法は、f77 と f95 では大きく異なります。f77 の場合、SPARC および x86 システムのデフォルトでは、浮動小数点演算例外において実行プログラムを中断するためのシグナルを自動的に生成しません。これは、シグナルはパフォーマンスを低下させるということと、期待された値が戻ってくれば、ほとんどの例外は重要でないという前提に基づいています。

f95 のデフォルトでは、ゼロ除算、オーバーフロー、無効演算の場合に自動的にトラップします。

f77 および f95 のコマンド行オプション -ftrap を使用してデフォルトを変更できます。-ftrap を使用した場合、f77 のデフォルトは -ftrap=%none で、f95 のデフォルトは -ftrap=common です。

例外トラップを有効に設定するためには、-ftrap オプションの 1 つを付けて (たとえば、-ftrap=common)、主プログラムをコンパイルします。

SPARC: 非標準の算術演算

段階的なアンダーフローと呼ばれる、標準の IEEE 算術演算の 1 つは手作業で無効にできます。無効にしたとき、プログラムは非標準の算術演算で実行していると考えられます。

算術演算に関する IEEE 規格では、アンダーフローとなった結果は、有効桁の小数部の基点を動的に調整することにより、段階的に扱うように指定しています。IEEE の浮動小数点の形式では、有効桁の前に小数部の基点が現れ、暗黙的な 1 の先行ビットがあります。浮動小数点の演算結果がアンダーフローとなったときに、段階的なアンダーフローでは、暗黙的な先行ビットをゼロにクリアーし、小数部の基点を有効桁方向にシフトさせるようになっています。これは、SPARC プロセッサではハードウェアではなく、ソフトウェアで行われます。このため、SPARC プロセッサ上で実行しているときに、プログラムにアンダーフローが多数発生すると (アルゴリズムに問題があることを示す可能性があります)、パフォーマンスが低下することになります。

段階的なアンダーフローを無効にするには、-fns オプションを付けてコンパイルします。または、プログラムの中からライブラリルーチン nonstandard_arithmetic() を呼び出して、段階的なアンダーフローを無効にします。standard_arithmetic() を呼び出して、段階的なアンダーフローを有効に戻します。

---

注 - 有効にするためには、アプリケーションの主プログラムを -fns を使用してコンパイルする必要があります。詳細は『Fortran ユーザーズガイド』を参照してください。

---

古いアプリケーションの場合、次のことに注意してください。

• standard_arithmetic() サブルーチンは、以前の gradual_underflow() という名前のルーチンに対応しています。
• nonstandard_arithmetic() サブルーチンは、以前の abrupt_underflow() という名前のルーチンに対応しています。

  ---

  注 - -fns オプションと nonstandard_arithmetic() ライブラリルーチンは、いくつかの SPARC システム上だけで有効です。x86 プラットフォーム上では、段階的なアンダーフローはハードウェアによって行われます。

  ---

IEEE ルーチン

次のインタフェースは、IEEE 演算を使用するユーザーの役に立ちます。これらのほとんどは、数学ライブラリ libsunmath と、いくつかの .h ファイルに置かれています。

• ieee_flags(3m) - 丸めの方向と丸めの精度を制御し、例外ステータスの問い合わせと例外ステータスのクリアーを行います。
• ieee_handler(3m) - 例外ハンドラルーチンを設定します。
• ieee_functions(3m) - 個々の IEEE 関数の名前と目的をリストします。
• ieee_values(3m) - 特別な値を返す関数をリストします。
• その他の libm 関数 (本節で説明)
  
  
  • ieee_retrospective
    
  • nonstandard_arithmetic
    
  • standard_arithmetic

SPARC プロセッサは、さまざまな側面において、ハードウェアとソフトウェアサポートを組み合わせることによって、IEEE 規格に準拠しています。x86 プロセッサは、ハードウェアサポートによって、IEEE 規格に完全に準拠しています。

最新の SPARC プロセッサには浮動小数点ユニットが含まれており、整数の乗算と除算の命令とハードウェアによる平方根演算機能を備えています。

コンパイルしたコードが実行時の浮動小数点ハードウェアに適切に一致したとき、最高のパフォーマンスが得られます。コンパイラの -xtarget= オプションは、実行時ハードウェアの指定を許可します。たとえば、-xtarget=ultra は、UltraSPARC プロセッサ上で最高のパフォーマンスを得られるオブジェクトコードを生成することをコンパイラに伝えます。

SPARC プラットフォーム: fpversion ユーティリティは、浮動小数点のどのハードウェアがインストールされているかを表示し、指定すべき適切な -xtarget 値を示します。このユーティリティは、すべての Sun SPARC アーキテクチャ上で動作します。詳細は、fpversion(1) のマニュアルページ、『Fortran ユーザーズガイド』の -xtarget に関する箇所、『数値計算ガイド』を参照してください。

フラグと ieee_flags()

ieee_flags 関数は、例外ステータスフラグの問い合わせやクリアーに使用します。この関数は、Sun コンパイラとともに出荷される libsunmath ライブラリに組み込まれていて、次の作業を行うことができます。

• 丸めの方向と丸めの精度の制御
• 例外ステータスフラグの検査
• 例外ステータスフラグのクリアー

ieee_flags の一般的な呼び出し方法は次のとおりです。

flags = ieee_flags( action, mode, in, out )

4 つの引数はすべて文字列です。入力は、action、mode、in です。出力は、out と flags です。ieee_flags は、整数値の関数です。flags は 1 ビットフラグの集合で、ここには有用な情報が返されます。詳細は、ieee_flags(3m) のマニュアルページを参照してください。

パラメータの取り得る値は次のとおりです。

表 6-1   ieee_flags ( action, mode, in, out ) の引数の値

action	mode	In, out
get set clear clearall	direction precision exception	nearest tozero negative positive extended double single inexact division underflow overflow invalid all common

precision モードは x86 プラットフォームだけで使用可能です。

これらはリテラル文字列であること、出力パラメータ out は最低でも CHARACTER*9 でなければならないことに注目してください。in と out の取り得る値の意味は、使用時の動作とモードによって異なります。これらの関係を次の表に要約します。

表 6-2   ieee_flags の引数の意味  

in と out の値	意味
nearest, tozero, negative, positive	丸めの方向
extended, double, single	丸めの精度
inexact, division, underflow, overflow, invalid	例外
all	5 種類すべての例外
common	3 種類の一般的な例外。 invalid (無効演算)、 division (ゼロ除算)、 overflow (オーバーフロー)

たとえば、フラグが立てられている例外のうちで何が最も優先順位が高いかを判別するときは、引数の in にヌル文字列を渡します。

CHARACTER *9, out ieeer = ieee_flags( 'get', 'exception', '', out ) PRINT *, out, ' flag raised'

また、overflow 例外フラグが立てられているかどうかを判別するときは、引数の in に overflow を設定します。復帰時に、out が overflow になっていれば、overflow 例外のフラグが立てられているということです。そうでない場合は、その例外は発生していません。

ieeer = ieee_flags( 'get', 'exception', 'overflow', out ) IF ( out.eq. 'overflow') PRINT *,'overflow flag raised'

例 : invalid 例外をクリアーします。

ieeer = ieee_flags( 'clear', 'exception', 'invalid', out )

例 : すべての例外をクリアーします。

ieeer = ieee_flags( 'clear', 'exception', 'all', out )

例 : ゼロに近づけるように丸めます。

ieeer = ieee_flags( 'set', 'direction', 'tozero', out )

例 : 丸めの精度を double に設定します。

ieeer = ieee_flags( 'set', 'precision', 'double', out )

ieee_flags を使用して警告メッセージを抑制する

次の例のように、動作 clear で ieee_flags を呼び出すと、クリアーされてない例外すべてがリセットされます。プログラムが終了する前にこの呼び出しを置くと、プログラム終了時の浮動小数点例外に関するシステム警告メッセージを抑制します。

例 : ieee_flags() を使用して、発生していたすべての例外をクリアーします。

i = ieee_flags('clear', 'exception', 'all', out )

ieee_flags を使用して例外を検出する

次の例は、前の計算によってどの浮動小数点例外フラグが立ったかを決定する方法を示しています。システム include ファイル floatingpoint.h に定義されたビットマスクは、ieee_flags により返された値に適用されます。

---

注 - Fortran 95 (f95) プログラムは floatingpoint.h ファイルをインクルードします。また Fortran 77 (f77) プログラムは f77_floatingpoint.h ファイルをインクルードします。

---

この例 DetExcFlg.F では、インクルードファイルは #include 前処理部指令を使用して導入されます。この指令には、.F 接尾辞のソースファイルを指定しなければなりません。アンダーフローの原因は、最小の倍精度数を 2 で除算しているためです。

例 : ieee_flags を使用して例外を検出し、デコードします。

#include "f77 floatingpoint.h" CHARACTER*16 out DOUBLE PRECISION d_max_subnormal, x INTEGER div, flgs, inv, inx, over, under x = d_max_subnormal() / 2.0 ! アンダーフローを発生 flgs=ieee_flags('get','exception','',out) ! どの例外が発生し たのか ? inx = and(rshift(flgs, fp_inexact) , 1) ! ieee_flags div = and(rshift(flgs, fp_division) , 1) ! によって under = and(rshift(flgs, fp_underflow), 1) ! 返された over = and(rshift(flgs, fp_overflow) , 1) ! 値を inv = and(rshift(flgs, fp_invalid) , 1) ! デコードする PRINT *, "Highest priority exception is: ", out PRINT *, ' invalid divide overflo underflo inexact' PRINT '(5i8)', inv, div, over, under, inx PRINT *, '(1 = exception is raised; 0 = it is not)' i = ieee_flags('clear', 'exception', 'all', out) ! すべてクリアー END

例 : 前出の例をコンパイルし、実行する (DetExcFlg.F)

demo% f95 DetExcFlg.F demo% a.out Highest priority exception is: underflow invalid divide overflo underflo inexact 0 0 0 1 1 (1 = exception is raised; 0 = it is not) demo%

IEEE 極値関数

コンパイラは、特別な IEEE 極値を返すために呼び出すことができる関数のセットを提供します。無限大や最小の正規数などの値は、アプリケーションプログラムの中で直接使用できます。

例 : 収束のテストはハードウェアによってサポートされる最小数に基づき、次のようになります。

IF ( delta .LE. r_min_normal() ) RETURN

利用できる値を次の表にリストします。

表 6-3   IEEE の値を使用する関数

IEEE の値	倍精度	単精度
infinity (無限大)	d_infinity()	r_infinity()
quiet NaN (シグナルを発しない非数)	d_quiet_nan()	r_quiet_nan()
signaling NaN (シグナルを発する非数)	d_signaling_nan()	r_signaling_nan()
min normal (最小の正規数)	d_min_normal()	r_min_normal()
min subnormal (最小の非正規数)	d_min_subnormal()	r_min_subnormal()
max subnormal (最大の非正規数)	d_max_subnormal()	r_max_subnormal()
max normal (最大の正規数)	d_max_normal()	r_max_normal()

2 つの NaN 値 (quiet と signaling) は順序がないものなので、
IF(X.ne.r_quiet_nan())THEN... のように比較の中で使用してはなりません。値が NaN であるかどうかを判別するときは、ir_isnan(r) か id_isnan(d) 関数を使用します。

これらの関数の Fortran 名は、次のマニュアルページにリストされています。

• libm_double(3f)
• libm_single(3f)
• ieee_functions(3m)

また、次も参照してください。

• ieee_values(3m)
• floatingpoint.h および f77_floatingpoint.h ヘッダーファイル

例外ハンドラと ieee_handler()

通常、IEEE 例外について、次のことを知っておく必要があります。

• 例外が発生するときのシステムの動作。
• ieee_handler() を使用して、ユーザー関数を例外ハンドラとして設定する方法。
• 例外ハンドラとして使用できる関数を作成する方法。
• 例外の発生場所を調べる方法。

ユーザールーチンへの例外トラップは、システムの浮動小数点例外に対するシグナルの生成から始まります。「浮動小数点例外のシグナル」を表す UNIX の公式名は SIGFPE です。SPARC および x86 プラットフォームは、デフォルトでは、例外が発生しても SIGFPE を生成しません。システムが SIGFPE を生成するためには、まず、例外トラップを有効にしなければなりません (これは通常は、ieee_handler() への呼び出しによって行います)。

例外ハンドラ関数を設定する

例外ハンドラとして関数を設定するときは、監視したい例外や対応動作といっしょに、関数の名前を ieee_handler() に渡します。ハンドラを一度設定すると、特定の浮動小数点の例外が発生するたびに、SIGFPE シグナルが生成され、指定した関数が呼び出されます。

ieee_handler() を起動する形式を次の表に示します。

表 6-4   ieee_handler (action, exception, handler) の変数

引数	種類	可能な値
action	文字列	get、set、clear のいずれか
exception	文字列	invalid、division、overflow、underflow、inexact のいずれか
handler	関数名	ユーザーハンドラ関数の名前、または、SIGFPE_DEFAULT、SIGFPE_IGNORE、SIGFPE_ABORT のいずれか
戻り値	整数	0 =OK

ieee_handler() を呼び出し、f77 によりコンパイルする Fortran 77 ルーチンでは、次の宣言も必要です。

#include 'f77_floatingpoint.h'

f95 の場合は、次のように宣言します。

#include 'f95/floatingpoint.h'

特別な引数 SIGFPE_DEFAULT、SIGFPE_IGNORE および SIGFPE_ABORT は、これらのインクルードファイルで定義され、特定の例外に対するプログラムの動作を変更するのに使用できます。

SIGFPE_DEFAULT または SIGFPE_IGNORE	指定した例外が発生しても何も動作しない。
SIGFPE_ABORT	例外発生時にはプログラムが異常終了し、恐らくダンプファイルを生成する。

ユーザー例外ハンドラ関数を作成する

ユーザーの例外ハンドラが行う動作はユーザーが自由に設定できます。しかし、ルーチンは、次に示す 3 つの引数を取る、整数型関数でなければなりません。

• handler_name( sig, sip, uap )
  
  
  • handler_name は整数関数の名前です。
    
  • sig は整数です。
    
  • sip は構造体 siginfo を持つ記録です。
    
  • uap は使用されません。

例 : 例外ハンドラ関数

INTEGER FUNCTION hand( sig, sip, uap ) INTEGER sig, location STRUCTURE /fault/ INTEGER address INTEGER trapno END STRUCTURE STRUCTURE /siginfo/ INTEGER si_signo INTEGER si_code INTEGER si_errno RECORD /fault/ fault END STRUCTURE RECORD /siginfo/ sip location = sip.fault.address ... ユーザーの行う処理 ... END

STRUCTURE 内のすべての INTEGER 宣言を INTEGER*8 で置き換えて、この f77 の例を SPARC V9 アーキテクチャ (-xarch=v9 または v9a) で実行できるよう変更する必要があります。

ieee_handler() によって有効にされるハンドラルーチンが例で示すように Fortran で書かれている場合、ハンドラルーチンは 1 番目の引数 (sig) に対しては、一切の参照を行ってはなりません。1 番目の引数は値でルーチンに渡され、loc(sig) としてのみ参照できます。この値はシグナル番号です。

ハンドラを使用して例外を検出する

次の例では、浮動小数点の例外を検出するためのハンドラルーチンを作成する方法を示します。

例 : 例外を検出して、異常終了します。

demo% cat DetExcHan.f EXTERNAL myhandler REAL r / 14.2 /, s / 0.0 / i = ieee_handler ('set', 'division', myhandler ) t = r/s END INTEGER FUNCTION myhandler(sig,code,context) INTEGER sig, code, context(5) CALL abort() END demo% f77 -silent DetExcHan.f demo% a.out 異常終了: abort が呼び出されました 異常終了 demo%

SIGFPE は、浮動小数点演算の例外が発生すると必ず生成されます。SIGFPE が検出されると、制御が myhandler 関数に渡され、即座に異常終了します。-g を付けてコンパイルし、dbx を使用すると、例外箇所を突き止めることができます。

ハンドラを使用して例外箇所を突き止める

例 : 例外箇所を突き止め (アドレスを出力して)、異常終了します。

demo% cat LocExcHan.F #include "f77_floatingpoint.h" EXTERNAL Exhandler INTEGER Exhandler, i, ieee_handler REAL r / 14.2 /, s / 0.0 /, t C ゼロ除算の検出 i = ieee_handler( 'set', 'division', Exhandler ) t = r/s END INTEGER FUNCTION Exhandler( sig, sip, uap) INTEGER sig STRUCTURE /fault/ INTEGER address INTEGER trapno END STRUCTURE STRUCTURE /siginfo/ INTEGER si_signo INTEGER si_code INTEGER si_errno RECORD /fault/ fault END STRUCTURE RECORD /siginfo/ sip WRITE (*,10) sip.si_signo, sip.si_code, sip.fault.address 10 FORMAT('Signal ',i4,' code ',i4,' at hex address ', z8 ) CALL abort() END demo% f77 -silent -g LocExcHan.F demo% a.out Signal 8 code 3 at hex address 11230 異常終了: abort が呼び出されました 異常終了 demo%

SPARC V9 環境では、各 STRUCTURE 内の INTEGER 宣言を INTEGER*8で置き換えて、書式中の i4 を i8 で置き換えます。

ほとんどの場合、例外の実際のアドレスを知るということは、dbx だけに意味があります。

demo% dbx a.out (dbx) stopi at 0x10d00 ブレークポイントをアドレスに設定する (2) stopi at &MAIN+0x68 (dbx) run プログラムを実行する 実行中: a.out (プロセス id 18803) MAIN で停止しました 0x10d00 で 0x10d00: MAIN+0x0068:fdivs %f3, %f2, %f2 (dbx) where 例外の行番号を表示する =>[1] MAIN()、"LocExcHan.F" の 7 行目 (dbx) list 7 ソースコード行を表示する 7 t = r/s (dbx) cont ブレークポイント後、継続してハンドラルーチンに入る Signal 8 code 3 at hex address 11230 異常終了: abort が呼び出されました シグナル ABRT (異常終了) 関数 _kill 0xef6e18a4 で 0xef6e18a4: _kill+0x0008:bgeu _kill+0x30 現関数 :exhandler 24 CALL abort() (dbx) quit demo%

もちろん、エラーの原因となるソース行を決定するためのより簡単な方法があります。しかし、この例は、例外処理の基本を示すのが目的です。

すべてのシグナルハンドラを無効にする

f77 の場合、バスエラー、セグメンテーション違反、不当命令をトラップするためのシステムシグナルハンドラは、デフォルトにより自動的に有効になります。

通常は、デフォルトの動作を停止したいとは思わないでしょう。しかし、大域的 C 変数 f77_no_handlers を 1 に設定する C のプログラムをコンパイルし、ユーザーの実行可能プログラムにリンクすれば、デフォルトの動作を無効にできます。

demo% cat NoHandlers.c int f77_no_handlers=1 ; demo% cc -c NoHandlers.c demo% f77 NoHandlers.o MyProgram.f

そうでない場合 (デフォルトでは)、f77_no_handlers は 0 です。この設定は、制御がユーザープログラムに移行される直前に効力を生じます。

この変数は、プログラムの大域的な名前空間にあります。したがって、f77_no_handlers をプログラムの他の場所で変数名として使用しないでください。

f95 の場合、デフォルトではシグナルハンドラは有効になりません。

発生種類の追求

ieee_retrospective 関数は、浮動小数点のステータスレジスタを調べて、どの例外が発生したのかを突き止め、標準エラーにメッセージを表示し、どの例外が発生してクリアーされていないのかをプログラマに知らせます。この関数は、プログラムが正常に終了するとき (CALL EXIT)、自動的に FORTRAN 77 プログラムによって呼び出されます。メッセージは通常、次のようになります (リリースによって若干異なります)。

注：IEEE 浮動小数点演算例外が発生： ゼロによる除算； IEEE 浮動小数点例外トラップが有効： 不正確；アンダーフロー；オーバーフロー；無効なオペランド； 『数値演算ガイド』、ieee_flags(3M)、ieee_handler(3M) を参照。

Fortran 95 プログラムは ieee_retrospective を自動的に呼び出さないので、明示的に呼び出して -lf77compat をつけてリンクする必要があります。

IEEE の例外のデバッグ

ほとんどの場合、オーバーフロー、アンダーフロー、無効演算などの浮動小数点の例外が発生したことを示すものは、プログラム終了時の要約メッセージだけです。どこで例外が発生したかを突き止めるためには、トラップを有効にした例外が必要です。これは、
-ftrap=common オプションを付けてコンパイルするか、ieee_handler() により例外ハンドラルーチンを呼び出すかによって行うことができます。例外トラップを有効にして、dbx か debugger からプログラムを実行し、dbx catch FPE コマンドを使用すれば、どこでエラーが発生するかを見つけることができます。

-ftrap=common を付けてコンパイルし直すことの利点は、例外をトラップするようにソースコードを変更する必要がないことです。しかし、ieee_handler() を呼び出すことによって、探すべき例外をより限定できます。

例 : -ftrap=common を付けてコンパイルし直して、dbx を使用します。

demo% f77 -g -ftrap=common -silent myprogram.f demo% dbx a.out a.out の読み込み中 rtld /usr/lib/ld.so.1 の読み込み中 libF77.so.3 の読み込み中 libc.so.1 の読み込み中 libdl.so.1 の読み込み中 (dbx) catch FPE (dbx) run 実行中: a.out (プロセス id 19739) シグナル FPE (ゼロによる浮動小数点除算) 関数 MAIN 行番号 212 ファイル "myprogram.f" 212 Z = X/Y (dbx) print Y y = 0.0 (dbx)

プログラムが終了したときオーバーフローと他の例外が発生していた場合、
ieee_handler() を呼び出してオーバーフローだけをトラップすることによって、最初のオーバーフローを突き止めることができます。このためには、次の例に示すように、主プログラムのソースコードを必要最小限だけ修正する必要があります。

例 : 他の例外も発生しているときにオーバーフローを突き止めます。

demo% cat myprog.F #include "f77_floatingpoint.h" program myprogram ... ier = ieee_handler('set','overflow',SIGFPE_ABORT) ... demo% f77 -g -silent myprog.F demo% dbx a.out a.out の読み込み中 rtld /usr/lib/ld.so.1 の読み込み中 libF77.so.3 の読み込み中 libc.so.1 の読み込み中 libdl.so.1 の読み込み中 (dbx) catch FPE (dbx) run 実行中: a.out (プロセス id 19793) シグナル FPE (浮動小数点オーバーフロー) 関数 MAIN 行番号 55 ファイル "myprog.F" 55 w = rmax * 200 ! オーバーフローの原因 (dbx) cont ! 実行を継続して、終了する Note: IEEE floating-point exception flags raised: Inexact; Division by Zero; Underflow; ! 他の例外もある IEEE floating-point exception traps enabled: overflow; See the Numerical Computation Guide ... 実行完了。終了コードは、0 です (dbx)

例外を選択するため、この例では、#include を導入しています。それで、ソースファイルの名前を .F 接尾辞に変更し、ieee_handler() を呼び出す必要があります。さらに一歩進んで、オーバーフロー例外時に呼び出されるユーザー独自のハンドラ関数を作成し、アプリケーション特有な解析を行い、異常終了する前に中間結果またはデバッグ結果を出力することもできます。

その他の例

この節では、無効演算、ゼロ除算、オーバーフロー、アンダーフロー、結果不正確の例外を予想外に生成する算術演算に関連する、より実際的な問題について説明します。

たとえば、IEEE 規格より前には、2 つの極めて小さな数をコンピュータで乗算すると、結果はゼロになっていました。メインフレームやミニコンピュータのほとんどが、この方式でした。これに対し、IEEE の算術演算には、段階的なアンダーフローが演算上の範囲を動的に拡張します。

たとえば、表現可能な最小値が 1.0E-38 である 32 ビットプロッサ を想定し、2 つの小さな数値を乗算してみます。

a = 1.0E-30 b = 1.0E-15 x = a * b

旧式の算術演算では結果が 0.0 になりますが、IEEE の算術演算では (同じワード長で) 結果は 1.40130E-45 となります。アンダーフローは、マシンが本来表せる値よりも小さい結果になったということを示します。この結果は、いくつかのビットが仮数部から「盗まれ」、指数部へシフトされたことによってできました。結果 (非正規化数) は、ある場合には精度が低くなりますが、逆に高くなる場合もあります。詳細は、このマニュアルで扱う範囲を超えているため、興味のある方は、1980 年 1 月 に発行された『Computer』誌 (第 13 巻 No.1) 、特に J. Coonen 氏の論文「Underflow and the Denormalized Numbers」を参照してください。

科学技術プログラムのほとんどは、方程式を解いたり、行列の因数分解など、丸めに敏感に反応するコードセクションがあります。段階的なアンダーフローがなければ、プログラマは不正確なしきい値への接近を検出する独自の方法を実装します。実装できなければ、独自のアルゴリズムの安定した実装はあきらめるかしかありません。

これらの話題に関する詳細は、『数値計算ガイド』を参照してください。

簡単なアンダーフローを防ぐ

アプリケーションの中には、実際に、非常にゼロに近いところで多くの処理をしているものがあります。これは、誤差や差分補正のアルゴリズムでよく発生します。数値的に安全で最大のパフォーマンスを得るには、重要な演算は拡張精度で行うべきです。アプリケーションが単精度の場合は、重要な演算を倍精度にすればかまいません。

例 : 単精度の、簡単なドット積の演算

sum = 0 DO i = 1, n sum = sum + a(i) * b(i) END DO

a(i) と b(i) が極めて小さい数であれば、多数のアンダーフローが発生することになります。演算を倍精度に移行すると、ドット積をより正確に計算でき、アンダーフローもなくなります。

DOUBLE PRECISION sum DO i = 1, n sum = sum + dble(a(i)) * dble(b(i)) END DO result = sum

SPARC プラットフォーム: アンダーフロー に関して、SPARC プロセッサを旧式のシステムのように動作 (Store Zero) させることができます。このためには、ライブラリルーチン nonstandard_arithmetic() への呼び出しを追加するか、アプリケーションの主プログラムを -fns オプションを付けてコンパイルします。

間違った答えのまま継続する

結果が明らかに間違っている場合でも、処理が継続されることに疑問を思われるかもしれません。IEEE の算術演算では、NaN や Inf など、どの種別の間違った答えを無視できるかをユーザーが指定できます。そして、そのような区別に基づいて決定が行われます。

たとえば、回路シミュレーションを想定してみましょう。特有の 50 行の演算の中で、引数の目的となる重要な変数は電圧量だけであり、この変数のとり得る値は +5v、0、
-5v だけであると仮定します。

途中の演算結果が正当な範囲にくるように、演算の各部分を詳細に調整できます。

計算された値が 4.0 よりも大きい場合は、5.0 を返します。 計算された値が -4.0 以上 +4.0以下の場合は、0 を返します。 計算された値が -4.0 よりも小さい場合は、-5.0 を返します。

さらに、Inf は許可されない値であるので、大きな数を乗算しないような特別なロジックが必要です。

IEEE 算術演算ではロジックはかなり単純にできます。演算を明確な形式で記述し、最終結果を正しい値に導くだけでかまいません。なぜなら、±Inf の発生する可能性があり、それは簡単にテストできるからです。

さらに、0/0 の特殊なケースも検出でき、ユーザーの望む形で処理できます。不必要な比較を行わなくてもすむため、結果は読みやすく、実行も高速です。

SPARC: アンダーフローの頻発

2 つの極めて小さな数を乗算すると、結果はアンダーフローとなります。

あらかじめ、乗算 (または減算) の被演算子が小さくなり、アンダーフローしそうであることがわかっていれば、計算を倍精度で行い、その後で結果を単精度に変換します。

たとえば、ドット積ループは次のようになります。

real sum, a(maxn), b(maxn) ... do i =1, n sum = sum + a(i)*b(i) enddo

a(*) と b(*) は、小さな要素が入ることがわかっているので、倍精度で実行し、数値の正確性を保持します。

real a(maxn), b(maxn) double sum ... do i =1, n sum = sum + a(i)*dble(b(i)) enddo

このようにすることによって、オリジナルのループが原因であるアンダーフロー頻発をソフトウェア的に解決し、パフォーマンスを上げることができます。しかし、これに関しては、確実で手間のかからない方法はありません。計算が多いコードを使用するユーザー自身の経験のほうが、より適切な解決策を決定することでしょう。

区間演算

Sun WorkShop 6 Fortran 95 コンパイラの f95 では、内的データタイプとして intervals をサポートしています。区間は、[a,b] = {z|a≦z≦b} で表され、1 組の数次は a≦b となります。区間は次の目的で使用できます。

• 非線形問題を解決します
• 厳密なエラー分析を実行します
• 数値の不安定性のソースを検出します

区間を内的データタイプとして Fortran 95 に導入したので、開発者は Fortran 95 の適用可能な構文や記号をすべてすぐに使用できるようになります。INTERVAL データタイプのほかに、f95 により次の区間拡張機能が Fortran 95 に取り込まれます。

• 3 クラスの INTERVAL 関係演算子
  
  
  • Certainly (断定的な関係)
    
  • possibly (可能性のある関係)
    
  • Set (集合の関係)
• INF、SUP、WID、HULL などの内的 INTERVAL 固有の演算子
• 1 数字入出力などの INTERVAL 入出力編集記述子
• 算術関数、三角関数、およびそのほかの数学関数に対する区間拡張機能
• 式コンテキストに依存した INTERVAL 定数
• 混合モードの区間式処理

Fortran 95 の区間固有の機能を使用するには、f95 コマンド行で -xia または -xinterval を指定します。

Fortran 95 の区間演算の詳細については、『Fortran 95 区間演算プログラミングリファレンス』を参照してください。