imslc2016 1 math IMSL Cライブラリ v2016.1 ユーザーズガイド : Math

IMSL

®

_{C イブ}

_{ユーザーズガイド}

Math

Version 2016.1

本資料 Rogue Wave Software, Inc. が作成し ローグウ ー ソフ ウ ア ジャパン株式会社が日本語訳したも

す 英語原文 日本語訳 間 内容に相違があ 場合に 英語原文が優先さ ます

目次

イントロダクション

IMSL C ライブラリ··· イントロダクション-2

本マニュアルの構成 ··· イントロダクション-2

適切な関数を見つける方法 ··· イントロダクション-3

命名規約 ··· イントロダクション-3

使用はじめに。 imsl.h ファイル ··· イントロダクション-4

エラー処理、アンダーフロー、オーバーフロー、ドキュメント例 ··· イントロダクション-4

出力配列のためのメモリ割り当て ··· イントロダクション-5

結果を印刷する方法 ··· イントロダクション-5

複素数演算 ··· イントロダクション-5

欠測値 ··· イントロダクション-5

ユーザ定義関数にデータを渡す方法 ··· イントロダクション-6

ユーザ定義関数から値を受け取る方法 ··· イントロダクション-7

スレッドセーフ使用法 ··· イントロダクション-8 OpenMPの使用法 ··· イントロダクション-8

ベンダー提供ライブラリの使用法 ··· イントロダクション-9 C++ からの使用法 ··· イントロダクション-10

行列格納モード ··· イントロダクション-13

第１章：線形方程式

フル行列の線形方程式

実行列 ··· lin_sol_gen 1-4

複素行列 ··· lin_sol_gen (複素数) 1-11

正定値行列の分解、解法、逆行列

実行列 ··· lin_sol_posdef 1-16

複素行列 ··· lin_sol_posdef (複素数) 1-20

帯行列の線形方程式

実行列 ··· lin_sol_gen_band 1-24

複素行列 ··· lin_sol_gen_band (複素数) 1-28

正定値対称行列の分解と解法

実行列 ··· lin_sol_posdef_band 1-32

複素行列 ··· lin_sol_posdef_band (複素数) 1-36

一般疎行列の線形方程式

実行列 ··· lin_sol_gen_coordinate 1-40

複素行列 ··· lin_sol_gen_coordinate (複素数) 1-49

疎行列の分解と解法 II

実行列 ··· superlu 1-55

複素行列 ··· superlu (複素数) 1-65

疎行列のOpenMP並列での分解と解法

複素行列 ··· superlu_smp (複素数) 1-85

正定値行列の分解と解法

実行列 ··· lin_sol_posdef_coordinate 1-95

複素行列 ··· lin_sol_posdef_coordinate (複素数) 1-101

正定値行列のOpenMP並列での分解と解法

実行列 ··· sparse_col_smp 1-107

複素行列 ··· sparse_col_smp (複素数) 1-115

反復法

一般化最小残差法(GMRES) ··· lin_sol_gen_min_residual 1-123

共役勾配法 ··· lin_sol_def_cg 1-127

フル行列の最小二乗法

最小二乗解法, QR 分解 ··· lin_least_squares_gen 1-133

非負の最小二乗解法 ··· nonneg_least_squares 1-140

線形拘束 ··· lin_lsq_lin_constraints 1-145

非負の行列の分解 (NNMF)

非負行列の分解と解 ··· nonneg_matrix_factorization 1-149

特異値分解 (SVD) と一般化逆行列

実行列 ··· lin_svd_gen 1-153

複素行列 ··· lin_svd_gen (複素数) 1-158

半正定値行列の分解、解法、一般化逆行列

実行列 ··· lin_sol_nonnegdef 1-163

第２章：固有値解析

線形固有値問題

一般行列

固有値と固有ベクトル ··· eig_gen 2-4

固有値と固有ベクトル ··· eig_gen (複素数) 2-7

実対称行列

固有値と固有ベクトル ··· eig_sym 2-10

複素エルミート行列

固有値と固有ベクトル ··· eig_herm (複素数) 2-13

一般化固有値問題

実対称行列とB正定値

固有値と固有ベクトル ··· eig_symgen 2-16

実行列 ··· geneig 2-19

複素行列 ··· geneig(複素数) 2-23

第３章

:

補間と近似

3

次スプライン補間

端点の条件によるスプライン補間 ... cub_spline_interp_e_cnd 3-9

形状保存のスプライン補間... cub_spline_interp_shape 3-15

張力連続性バイアス(TCB)スプライン補間 ... cub_spline_tcb 3-19

3

次スプライン計算と積分

計算と微分... cub_spline_value 3-25

スプライン補間

1次元補間 ... spline_interp 3-30

補間データが与えられたノット列 ... spline_knots 3-35

2次元テンソル積補間 ... spline_2d_interp 3-39

スプライン計算と積分

1次元計算と微分 ... spline_value 3-44

1次元積分 ... spline_integral 3-47

2次元計算と微分 ... spline_2d_value 3-49

2次元積分 ... spline_2d_integral 3-52

多次元補間

多次元補間と微分 ... spline_nd_interp 3-55

最小二乗近似と平滑化

ユーザ提供の関数 ...user_fcn_least_squares 3-59

固定ノットを持つスプライン ... spline_least_squares 3-65

固定ノットを持つテンソル積スプライン ... spline_2d_least_squares 3-70

3次元平滑化スプライン ... cub_spline_smooth 3-75

制約のあるスプライン ... spline_lsq_constrained 3-79

エラー検出による1次元データの平滑化 ... smooth_1d_data 3-85

散布データ補間

秋間Akima の曲面当てはめ法 ... scattered_2d_interp 3-89

散布データ最小二乗法

動径基底関数を使用した当てはめ ... radial_scattered_fit 3-93

動径基底当てはめの計算 ... radial_evaluate 3-99

第４章：積分と微分

単変量求積

アダプティブ汎用、端点特異 ··· int_fcn_sing 4-4 1Dアダプティブ汎用、内点あるいは端点に特異 ··· int_fnc_sing_1d 4-8

アダプティブ汎用 ··· int_fcn 4-13

アダプティブ汎用、特異点 ··· int_fcn_sing_pts 4-17

アダプティブ、加重付き代数的特異 ··· int_fcn_alg_log 4-21

アダプティブ、無限区間 ··· int_fcn_inf 4-25

アダプティブ、加重付き振動(三角関数) ··· int_fcn_trig 4-29

アダプティブ、加重付きフーリエ(三角関数) ··· int_fcn_fourier 4-33

コーシーの主値 ··· int_fcn_cauchy 4-37

非アダプティブ汎用 ··· int_fcn_smooth 4-41

多変量求積

2D反復積分 ··· int_fcn_2d 4-44 2D積分、内点あるいは端点に特異 ··· int_fcn_sing_2d 4-48 3D積分、内点あるいは端点に特異 ··· int_fcn_sing_3d 4-54

超矩形上の反復積分. ··· int_fcn_hyper_rect 4-61

擬似モンテカルロ法を使用する反復積分 ··· int_fcn_qmc 4-64

ガウス求積

ガウス求積公式 ··· gauss_quad_rule 4-67

微分

第５章

:

微分方程式

一階常微分方程式

常微分方程式の初期値問題の解

Runge-Kutta法 ··· ode_runge_kutta 5-4

常微分方程式の初期値問題の解

有限差分法 ··· bvp_finite_difference 5-10

微分代数方程式の解 ··· differential_algebraic_eqs 5-20

一階と二階常微分方程式

常微分方程式の初期値問題の解、可変次数Adams法 ··· ode_adams_krogh 5－34

偏微分方程式

一次元連立偏微分方程式の解

pde_1d_mgの紹介 ··· 5-41

移動グリッドで

一次元時間依存連立偏微分方程式を解く ··· pde_1d_mg 5-43

エルミート3次多項式での直線法 ··· modified_method_of_lines 5-74

エルミート5次スプラインを用いて

一般化Feyman-kac方程式を解く ··· feynman_kac 5-88

エルミート５次スプラインの値

またはその微係数を計算する。 ··· feynman_kac_evaluate 5-117

二次元の偏微分法程式の解

高速Poissonソルバー ··· fast_poisson_2d 5-120

第６章：変換

実高速フーリエ変換

実高速フーリエ変換 ···fft_real 6-3

実高速フーリエ変換初期化 ··· fft_real_init 6-7

複素高速フーリエ変換

複素高速フーリエ変換 ··· fft_complex 6-9

複素高速フーリエ変換初期化 ··· fft_complex_init 6-12

実正弦と実余弦高速フーリエ変換

フーリエ余弦変換 ··· fft_cosine 6-14

フーリエ余弦変換初期化 ··· fft_cosine_init 6-16

フーリエ正弦変換 ··· fft_sine 6-18

フーリエ正弦変換初期化 ··· .fft_sine_init 6-20

2

次元高速フーリエ変換

複素2次元高速フーリエ変換 ··· fft_2d_complex 6-22

たたみ込みと相関

実たたみ込みと相関 ··· convolution 6-26

複素たたみ込みと相関 ··· convolution (複素数) 6-31

ラプラス変換

複素関数の逆ラプラス変換 ··· inverse_laplace 6-36

多項式のゼロ点

Jenkins-Traub 法を使用した実係数の多項式のゼロ点 ··· zeros_poly 7-3 Jenkins-Traub 法を使用した複素係数の多項式のゼロ点 ··· zeros_poly (複素数) 7-5

関数のゼロ点

実単変量関数のゼロ点 ··· zeros_univariate 7-7

実関数のゼロ点 ··· zeros_function 7-10

非線形連立方程式の根

第８章：最適化

非拘束最小化

単変量関数

関数値だけを使用 ··· min_uncon 8-4

関数値と1次導関数を使用 ··· min_uncon_deriv 8-7

単変量非平滑関数の最小点を見つける ··· min_uncon_golden 8-11

多変量関数

準Newton 法を使用 ··· min_uncon_multivar 8-14

直接探査ポリトープ法を用いて非平滑関数の最小点を求める ··· min_uncon_polytope 8-20

非線形最小二乗法

Levenberg-Marquardt 法を使用 ··· nonlin_least_squares 8-24

線形拘束最小化

線形計画問題あるいは二次計画問題のMPSファイルの読み込み ··· read_mps 8-32

線形計画問題を解く ··· linear_programing 8-39

密線形計画法 ··· lin_prog 8-44

輸送問題を解く ··· transport 8-48

二次計画法 ··· quadratic_prog 8-51

疎の線形計画法 ··· sparse_lin_prog 8-55

疎の二次計画法 ··· sparse_quadratic_prog 8-65

一般目的関数の最小化 ··· min_con_gen_lin 8-76

単純境界の非線形最小二乗法 ··· bounded_least_squares 8-82

直接探査コンプレックス法による非平滑関数の最小化 ··· min_con_polytope 8-89

信頼領域法による線形制約のある関数の最小化 ··· min_con_lin_trust_region 8-97

非線形拘束最小化

逐次等式制約のある二次計画法 ··· constrained_nlp 8-101

サービスルーチン

有限差分ヤコビアン ··· jacobian 8-108

第９章：特殊関数

誤差関数とガンマ関数

誤差関数 ... erf 9-6

補誤差関数 ... erfc 9-8

スケールされた関数 ... erfe 9-11

逆誤差関数 ... erf_inverse 9-12

逆補誤差関数 ... erfc_inverse 9-14

ベータ関数 ... beta 9-16

対数ベータ関数 ... log_beta 9-18

不完全ベータ関数 ... beta_incomplete 9-19

ガンマ関数

ガンマ関数 ... gamma 9-21

対数ガンマ関数 ... log_gamma 9-23

不完全ガンマ関数 ... gamma_incomplete 9-25

Psi関数

対数ガンマ関数の導関数... psi 9-27

実 psi1 関数ψ1(x) ... psi1 9-28

ベッセル関数

ベッセル関数 J0 ... bessel_J0 9-29 ベッセル関数 J1 ... bessel_J1 9-31 ベッセル関数 Jn ... bessel_Jx 9-32 ベッセル関数 Y0 ... bessel_Y0 9-34 ベッセル関数 Y1 ... bessel_Y1 9-36 ベッセル関数 Yv ... bessel_Yx 9-38 ベッセル関数 I0 ... bessel_I0 9-40 ベッセル関数 e-|x| _I

0(x) ... bessel_exp_I0 9-42

ベッセル関数 I1 ... bessel_I1 9-43

ベッセル関数 e-|x| _I

1(x) ... bessel_exp_I1 9-44

ベッセル関数 Iv ... bessel_Ix 9-45 ベッセル関数 K0 ... bessel_K0 9-47 ベッセル関数 ex _{K0(x) ... bessel_exp_K0 9-49}

ベッセル関数 K1 ... bessel_K1 9-50 ベッセル関数 ex _{K1(x) ... bessel_exp_K1 9-51}

ベッセル関数 Kv . ... bessel_Kx 9-52

楕円積分

第１種完全楕円積分 ... elliptic_integral_K 9-54

第２種完全楕円積分 ... elliptic_integral_E 9-55

第１種カールソンの楕円積分 ... elliptic_integral_RF 9-56

第２種カールソンの楕円積分 ... elliptic_integral_RD 9-57

第３種カールソンの楕円積分 ... elliptic_integral_RJ 9-58

カールソンの楕円積分の特例 ... elliptic_integral_RC 9-59

フレネル積分

余弦フレネル積分 ... fresnel_integral_C 9-60

正弦フレネル積分 ... fresnel_integral_S 9-61

エアリー関数

エアリー関数 ... airy_Ai 9-62

第２種エアリー関数 ... airy_Bi 9-63

エアリー関数の導関数 ...airy_Ai_derivative 9-64

第２種エアリー関数の導関数 ...airy_Bi_derivative 9-65

ケルビン関数

第１種０次ケルビン関数 ber ... kelvin_ber0 9-66

第１種０次ケルビン関数 bei ... kelvin_bei0 9-67

第２種０次ケルビン関数 ker ... kelvin_ker0 9-68

第２種０次ケルビン関数 kei ... kelvin_kei0 9-69

ケルビン関数 bei の導関数 ... kelvin_bei0_derivative 9-71

ケルビン関数 ker の導関数 ... kelvin_ker0_derivative 9-72

ケルビン関数 kei の導関数 ... kelvin_kei0_derivative 9-73

統計関数

正規 (ガウス) 分布関数 ··· normal_cdf 9-74

逆正規分布関数 ··· normal_inverse_cdf 9-76

χ2 _{分布関数 ··· chi_squared_cdf 9-77}

逆_χ2 _{分布関数 ··· chi_squared_inverse_cdf 9-79}

F 分布関数 ··· F_cdf 9-81

逆 F 分布関数 ··· F_inverse_cdf 9-83

スチューデント t 分布関数 ··· t_cdf 9-85

逆スチューデント t 分布関数 ··· t_inverse_cdf 9-87

ガンマ分布関数 ··· gamma_cdf 9-88

二項分布関数 ··· binomial_cdf 9-90

超幾何分布関数 ··· hypergeometric_cdf 9-92

ポアソン分布関数 ··· poisson_cdf 9-94

ベータ分布関数 ··· beta_cdf 9-96

逆ベータ分布関数 ··· beta_inverse_cdf 9-98

二変量正規分布関数 ··· bivariate_normal_cdf 9-99

基本金融関数

累積利息の計算 ··· cumulative_interest 9-101

元本累計の計算 ··· cumulative_principal 9-103

減価償却の計算（定率法） ··· depreciation_db 9-105

減価償却の計算（倍額定率法） ··· depreciation_ddb 9-107

減価償却の計算（定額法） ··· depreciation_sln 9-109

減価償却の計算（算術級数法） ··· depreciation_syd 9-110

減価償却の計算（倍額定率法） ··· depreciation_vdb 9-111

分数価格から小数価格への変換 ··· dollar_decimal 9-113

小数価格から分数価格への変換 ··· dollar_fraction 9-114

実効年利率の計算 ···effective_rate 9-115

将来価値の計算 ··· future_value 9-116

複利予定表に基づく初期元金の将来価値の計算 ··· future_value_schedule 9-118

利息の計算 ··· interest_payment 9-119

年金の期間当たりの利率計算 ··· interest_rate_annuity 9-121

キャッシュフローに対する内部収益率の計算 ··· internal_rate_of_return 9-123

キャッシュフローに対する内部収益率の計算

（キャッシュフローは定期的である必要はない） ··· internal_rate_schedule 9-125

定期キャッシュフローに対する

修正内部収益率の計算 ··· modified_internal_rate 9-127

正味現在価値の計算 ··· net_present_value 9-129

名目年利率の計算 ··· nominal_rate 9-131

特定の目的の投資に必要な支払いの期間数の計算 ··· number_of_periods 9-132

投資の定期的支払い額の計算 ··· payment 9-134

現在価値の計算 ··· present_value 9-136

キャッシュフローに対する現在価値の計算

（キャッシュフローは定期的である必要はない） ··· present_value_schedule 9-138

特定期間の元金支払い額の計算 ··· principal_payment 9-140

債券関数

満期日に利息を払う債券の累積利息の計算 ··· accr_interest_maturity 9-142

定期的に利息を払う債券の累積利息の計算 ··· accr_interest_periodic 9-144

米国財務省短期債券の債券換算利回りの計算 ··· bond_equivalent_yield 9-146

証券のコンベクシティの計算 ··· convexity 9-148

決済日を含む利払い期間の日数の計算 ··· coupon_days 9-150

決済日と満期日の間の支払うクーポン数の計算 ··· coupon_number 9-152

クーポン期間最初の日から決済日までの

決済日から次の利払い日までの日数の計算 ··· days_to_next_coupon 9-156

資産耐用年数に基づく減価償却係数を

適用した各会計期間の減価償却費の計算 ··· depreciation_amordegrc 9-158

資産耐用年数に基づく減価償却係数を

適用しない各会計期間の減価償却費の計算 ··· depreciation_amorlinc 9-160

割引債の価格の計算 ··· discount_price 9-162

債券の割引率の計算 ··· discount_rate 9-164

割引債の年利回りの計算 ··· discount_yield 9-166

マコーレー・デュレーションの計算 ··· duration 9-168

債券の利率の計算 ··· interest_rate_security 9-170

修正マコーレー・デュレーションの計算 ··· modified_duration 9-172

決済日以降のクーポン期日の計算 ··· next_coupon_date 9-174

決済日直前のクーポン期日の計算 ··· previous_coupon_date 9-176

定期的に利息を払う額面100ドルあたりの債券価格の計算 ··· price 9-178

満期日に利息を払う額面100ドルあたりの債券価格の計算 ··· price_maturity 9-180

満期日の受取額の計算 ··· received_maturity 9-182

額面100ドルあたりの米国財務省短期債券価格の計算 ··· treasury_bill_price 9-184

米国財務省短期債券利回りの計算 ··· treasury_bill_yield 9-186

２つの日付間の総日数で表される１年の端数の計算 ··· year_fraction 9-188

満期日に利息を払う債券の年利回りの計算 ··· yield_maturity 9-190

定期的に利息を払う債券の年利回りの計算 ··· yield_periodic 9-192

第１０章：統計と乱数生成

統計

単変量要約統計量 simple_statistics 10-4

一元度数表 table_oneway 10-8

χ2 _{１標本適合度検定 chi_squared_test 10-12}

相関 covariances 10-19

多重線形回帰 regression 10-24

多項式回帰 poly_regression 10-31

数値順位解析 ranks 10-37

乱数

シードの現在値の検索 ... random_seed_get 10-42

乱数シードの初期化 ... random_seed_set 10-44

一様 (0, 1) 乱数生成器の選択 ... random_option 10-45

疑似一様乱数を生成 ... random_uniform 10-46

疑似正規乱数を生成 ... random_normal 10-48

疑似ポアソン乱数を生成... random_poisson 10-50

疑似ガンマ乱数を生成 ... random_gamma 10-52

疑似ベータ乱数を生成 ... random_beta 10-54

疑似標準指数乱数を生成... random_exponential 10-56

低不適合度乱数列

(low-discrepancy sequence)

シャッフルされたFaure列を生成 ...faure_next_point 10-58

第１１章：プリント関数

行列、又は、ベクトルをプリントする ··· write_matrix 11-2

ページの幅と高さを設定する ···page 11-7

印刷オプションを設定する ··· write_options 11-9

第１２章：ユーティリティ

出力ファイルの設定 ··· output_file 12-3

ライブラリバージョンとライセンス番号の取得 ··· version 12-6

時間と日付

CPU 使用時間 ··· ctime 12-7

今世紀の日付から日数 ··· date_to_days 12-8

今世紀の日数から日付 ··· days_to_date 12-9

エラー処理

エラーメッセージオプション ··· error_options 12-10

エラータイプを取得 ··· error_type 12-15

エラーメッセージを取得 ··· error_message 12-16

エラーコードを取得 ··· error_code 12-17

エラー処理システムを初期化する ··· initialize_error_hander 12-18

ユーザ提供の関数から戻る場合の設定 ··· set_user_fnc_return_flag 12-20

C

ランタイム

ライブラリ

メモリを解放する ··· free 12-25

ファイルをOpenする ··· fopen 12-26

ファイルをCloseする ··· fclose 12-28

OpenMP

OpenMP オプション ··· omp_options 12-29

定数

自然定数と数学定数 ··· constant 12-31

整数マシン定数 ··· machine (整数) 12-35

浮動小数点マシン定数 ··· machine (浮動小数点) 12-37

ソート

浮動小数点ベクトルのソート ··· sort 12-39

整数ベクトルのソート ··· sort (整数) 12-41

ベクトルノルムを計算する方法

種々のノルムの計算 ··· vector_norm 12-43

種々のノルムの計算 ··· vector_norm (複素数) 12-45

線形代数支援ルーチン

ベクトル-ベクトル、行列-ベクトル、行列-行列 乗算

実行列 ··· mat_mul_rect 12-48

複素行列 ··· mat_mul_rect (複素数) 12-51

実帯行列 ··· mat_mul_rect_band 12-54

複素帯行列 ··· mat_mul_rect_band (複素数) 12-58

実座標対応行列 ··· mat_mul_rect_coordinate 12-62

複素座標対応行列 ··· mat_mul_rect_coordinate (複素数) 12-66

ベクトル-ベクトル、行列-ベクトル、行列-行列 加算

実帯行列 ··· mat_add_band 12-71

複素帯行列 ··· mat_add_band (複素数) 12-74

実座標対応行列 ··· mat_add_coordinate 12-78

複素座標対応行列 ··· mat_add_coordinate (複素数) 12-81

行列のノルム

実行列 ··· matrix_norm 12-84

実帯行列 ··· matrix_norm_band 12-86

実座標対応行列 ··· matrix_norm_coordinate 12-89

実行列 ··· generate_test_band 12-92

複素行列 ··· generate_test_band (複素数) 12-94

実座標対応行列 ··· generate_test_coordinate 12-96

複素座標対応行列 ··· generate_test_coordinate (複素数) 12-100

GPU

サポート

NVIDIA CUDA Toolkit を使うためのプログラミングノード ··· 12-104

閾値を得る ··· cuda_get 12-106

閾値を設定する ··· cuda_set 12-108 NVIDIAメモリを解放する ··· cuda_free 12-109

参照資料

イントロダクション

目次

IMSL C ライブラリ··· イントロダクション-2

本マニュアルの構成 ··· イントロダクション-2

適切な関数を見つける方法 ··· イントロダクション-3

命名規約 ··· イントロダクション-3

使用はじめに。 imsl.h ファイル ··· イントロダクション-4

エラー処理、アンダーフロー、オーバーフロー、ドキュメント例 ··· イントロダクション-4

出力配列のためのメモリ割り当て ··· イントロダクション-5

結果をプリントする方法 ··· イントロダクション-5

複素数演算 ··· イントロダクション-5

欠測値 ··· イントロダクション-5

ユーザ定義関数にデータを渡す方法 ··· イントロダクション-6

ユーザ定義関数から値を受け取る方法 ··· イントロダクション-7

スレッドセーフ使用法 ··· イントロダクション-8 OpenMPの使用法 ··· イントロダクション-8

ベンダー提供ライブラリの使用法 ··· イントロダクション-9 C++ からの使用法 ··· イントロダクション-10

IMSL C

ライブラリ

IMSL C ライブラリは科学技術計算のプログラミングに役立つC関数のライブラリです。各関数は科学

技術専門家の研究活動に使用されるように設計され文書化されています。いくつかの例題プログラムは、 出力結果をグラフで表示しています。

本マニュアルの構成

本マニュアルは各関数の簡潔な記述を含みます。各関数には少なくとも1つの例題があり、サンプルの

入力と結果があります。

各章はその章で紹介されている関数を一覧した目次から始まり、イントロダクションへと続きます。各 関数のドキュメントは、以下の情報から構成されています。

･ Section Name（セクション名）: 通常、その関数のfloatタイプとdoubleタイプの共通のルート名 が記述される。

･ Purpose（目的）: その関数の目的。

･ Synopsis（梗概）： リストされた必須な引数でその関数を呼ぶ形式。

･ Required Arguments（必須な引数）: 必須な引数。次の項目が順に記されます。

Input（入力）:引数は初期化される必要があります。この値は関数によって変更されません。 Input/Output（入力/出力）: 引数は初期化される必要があります。関数はこの引数に出力を

返します。この引数は、定数あるいは数式であってはいけません。

Output（出力）:初期化する必要はありません。この引数は、定数あるいは数式であってはい けません。関数は、この引数に出力を返します。

･ Return Value（戻り値）: 関数によって返される値。

･ Synopsis with Optional Arguments（オプション引数の梗概）:リストされた必須の引数とオプション 引数で関数を呼び出す場合の形式。

･ Optional Arguments（オプション引数）: 各オプション引数の説明。

･ Description（説明）:アルゴリズムの説明と詳細な情報のための参考文献。多くの場合、同様な機能

もしくは、補足的な機能をもつ他のIMSL関数についても言及されます。

･ Examples（例題）: 入力とオプション引数の使い方を示した1つ以上の例題。

･ Errors（エラー）:この関数で起こり得るエラーのリスト。エラータイプに関する説明は、参照資料 の「ユーザエラー」にあります。エラーは、以下のタイプに分類されます。

適切な関数を見つける方法

IMSL C ライブラリは章ごとに分類され、各章は、同じ系統の計算・分析をする関数を集めています。

与えられた問題に対する適切な関数を見つけるには、各章の始めにある目次か、このマニュアルの最後 にあるアルファベット順の関数リストのいずれかを使用することができます。

IMSL C ライブラリを迅速に使用するには、使用者の問題に似ている例題を見つけて、その例題をまね

る方法がしばしば有効です。各関数の例題は1つ以上あります。

命名規約

ほとんどの関数には、floatタイプとdoubleタイプがあり、これらは共通のルート名を持っています。

いくつかの関数はintタイプ、あるいは IMSLが定義するf_complexタイプとd_complexタイプがありま す。複数のタイプがある場合の、各タイプの関数名の接頭辞は次の通りです。

タイプ 接頭辞

float double int f_complex d_complex

imsl_f_ imsl_d_ imsl_i_ imsl_c_ imsl_z_

関数のセクション名は、その関数をより簡単に見つけられるように共通のルート名のみが記されます。 例えば、関数 imsl_f_lin_sol_gen と imsl_d_lin_sol_gen は、第１章「線形方程式」の lin_sol_gen にあります。

このマニュアルで使用される変数名は、本マニュアルの各章を通して一貫して使用されます。 例えば、

固有値解析の関数では、evalは固有値のベクトルで、n_evalは計算された、あるいは計算される固

有値の数を示します。

IMSL C ライブラリを呼び出すプログラムを書く場合、ユーザは、IMSLの外部関数名と競合しないC

の名前を選ばなければなりません。名前を選ぶ際に、以下の規則に従うことで、IMSL名との競合を避

けることができます。

使用はじめに。

_imsl.h

ファイル

使用はじめに

IMSL C ライブラリの関数を使用するには、最初にCでその関数を呼び出すプログラムを書かなければ

ならなりません。各関数は、プログラミングとドキュメントにおいて決まった規定に従っています。こ の製品は開発における優先順位を、効率の良いアルゴリズム、明瞭なドキュメント、正確な計算結果と しています。これらのルーチンの均一な設計は作成するアプリケーションの中で１つ以上のルーチンを

使用することを容易にします。そして又この設計の一貫性により、１つの IMSL ライブラリルーチン

に対する経験を、その後使用する全ての IMSL ルーチンに適用することができます。

imsl.h

ファイル

マニュアル内の全ての例題は、インクルードファイル <imsl.h> を使用します。このファイルは、次

の全てのIMSL定義の関数のプロトタイプを含みます：

スプライン構造体 Imsl_f_ppoly、 Imsl_d_ppoly、 Imsl_f_spline、 Imsl_d_spline； 列挙データ型 Imsl_quad、 Imsl_write_options、Imsl_page_options, Imsl_ode,Imsl_error； IMSL定義のデータタイプ f_complex (float complex型) d_complex (double complex型)。

エラー処理、アンダーフロー、オーバーフロー、ドキュメント例

IMSL C ライブラリの関数は、エラーと不当な入力を検出し、レポートします。このエラー処理は、エ

ラー条件の処置のためにユーザが特定の用意をすることなく自動的に行なわれます。 エラーは重要度

で分類され、コード番号が割り当てられます。デフォルトでは、重要度が中程度以上のエラーでは自動

的にメッセージが出力されます。 更に、もっとも重要度の高いエラーではプログラムの実行を停止し

ます。 それらのエラーの一般的な性質と重要性のレベルは、シンボリック名 IMSL_FATAL、

IMSL_WARNING、などを持つ「エラータイプ」によって指定されます。詳しいことは、「参照資料」の

「ユーザエラー」節を参照して下さい。

一般に、システム（ハードウェア、又は、ソフトウェア）がアンダーフローをゼロに置き換える場合、

IMSL C ライブラリは、計算がアンダーフローによって影響されないように書かれています。通常、ア

ンダーフローを示すシステム・エラーメッセージは無視することができます。

また、IMSLコードは、オーバーフローを回避するように書かれています。オーバーフローのシステム・

エラーメッセージが出力された場合は、不正確な入力データ、引数タイプの間違い、不適当な次元数の ようなプログラミングミスを調査すべきです。

多くの場合、計算がうまくいかなくなる注意点を各関数の説明で指摘してあります。

出力配列のためのメモリ割り当て

多くの関数は、計算結果を含んだ配列へのポインタを返します。 デフォルトではIMSL C ライブラリ

関数の値として返された配列はその関数によって割り当てられたメモリに格納されます。この領域を解

放するには、imsl_freeを使用します。 呼び出しプログラムにより割り当てられた領域の配列を返す

には次のオプション引数を使用します。 IMSL_RETURN_USER, float a[]

このように、計算結果のための領域の割り当てはユーザによるか、あるいは、関数により内部的に行な われます。

同様に、他のオプション引数は、追加の計算結果の配列がユーザによって割り当てられるか、又は内部 的にその関数が割り当てるかを指定します。例えば「線形システム」の多くの関数で、オプション引数

IMSL_INVERSE_USER, float inva[] (出力) IMSL_INVERSE, float **p_inva (出力)

は2つの互に排他的なオプション引数を指定します。最初のオプションが選ばれた場合、その行列の逆

行列は、ユーザ提供の配列inva に格納されます。

2番目のオプションでは、 float **p_inva は逆行列へのポインタのアドレスを参照します。呼ばれた関

数はその配列のためのメモリを割り当て、このメモリ領域へのポインタを *p_inva とセットします。

一般的に float *p_inva が宣言されると、この関数への引数として &p_inva を使用します。この領域

を解放するにはimsl_free(p_inva) を使用します。

結果をプリントする方法

IMSL C ライブラリのほとんどの関数は、結果をプリントしません。その結果は、C変数で返されます。

IMSL C/ ライブラリには、配列のプリント用の幾つかの特殊な関数があります。 例えば、

imsl_f_write_matrix はタイプ float の行列をプリントするのに便利な関数です。これらの関数の詳

細な説明は、第11章「プリント関数」を参照してください。

複素数演算

ユーザは、IMSLの定義済みデータタイプを使用することによって複素数演算を行なうことができます。

これらのタイプは次の2つの浮動小数点演算精度で使用できます：

◍ 単精度複素数値に対して f_complex ◍ 倍精度複素数値に対して d_complex

複素データタイプと関数の説明は、「参照資料」にあります。

欠測値

IMSL C ライブラリの幾つかの関数は、そのデータが欠測値を含むことを許します。これらの関数は欠

測値を、"not a number" 又は NaN として参照される特殊な値として認識します。実際の値はコンピ ュータシステムにより異なりますが、それは第12章「ユーティリティ」のIMSL関数imsl_f_machine で求められます。

ユーザ定義関数にデータを渡す方法

問題特有のデータを、IMSL C ライブラリのインターフェイスを通して、ユーザ定義関数に渡すことが

必要な場合があります。 この機能は、ユーザ定義関数がユーザの呼び出し関数にローカルなデータを

必要とする場合に、グローバルデータを用いてユーザ定義関数がそのデータをアクセスすることを避け

たい場合有用です。 IMSL C ライブラリのいくつかの関数は、ユーザ指定のデータをその関数に渡せ

るように、オプション引数でデータへのポインタをユーザ定義関数が持てるようになっています。 次

の例題は、IMSL C ライブラリの関数 imsl_f_min_uncon でオプション引数 IMSL_FCN_W_DATA を 使用したこの機能の例です。

例

#include <imsl.h> #include <math.h> #include <stdio.h>

float fcn_w_data(float x, void *data);

int main() {

float a = -100.0; float b = 100.0; float fx, x;

float usr_data[] = {5.0, 10.0}; x = imsl_f_min_uncon (NULL, a, b,

IMSL_FCN_W_DATA, fcn_w_data, usr_data, 0);

fx = fcn_w_data(x, (void*)usr_data);

printf ("The solution is: %8.4f\n", x);

printf ("The function evaluated at the solution is: %8.4f\n", fx);

}

/*

* User function that accepts additional data in a (void*) pointer. * This (void*) pointer can be cast to any type and dereferenced to * get at any sort of data-type or structure that is needed.

* For example, to get at the data in this example

* *((float*)data) and usr_data[0] contains the value 5.0 * *((float*)data+1) and usr_data[1] contains the value 10.0 */

float fcn_w_data(float x, void *data) {

float *usr_data = (float*)data;

ユーザ定義関数から値を受け取る方法

ユーザ定義の関数から戻される全ての値は正当な実数でなければなりません。 ユーザ定義関数により

戻される値が、NaN、無限大、又は，負の無限値を含まないことを確かめることはユーザの責任です。

下の例題に示した方法のほかに、ユーザ定義関数の中で回復不可能なエラーが発生する場合、IMSL C

ライブラリが呼び出しプログラムに制御を返すようにすることもできます。この機能については関数

imsl_set_user_fcn_return_flag を参照してください。

例

#include <imsl.h>

#include <math.h>

void fcn(int, int, float[], float[]);

int main()

{

int m=3, n=1;

float *result, fx[3];

float xguess[]={1.0};

result = imsl_f_nonlin_least_squares(fcn, m, n, IMSL_XGUESS,

xguess, 0);

fcn(m, n, result, fx);

/* Print results */

imsl_f_write_matrix("The solution is", 1, 1, result, 0);

imsl_f_write_matrix("The function values are", 1, 3, fx, 0);

}

void fcn(int m, int n, float x[], float f[])

{

int i;

float y[3] = {2.0, 4.0, 3.0};

float t[3] = {1.0, 2.0, 3.0};

for (i=0; i<m; i++)

{

/* check for x=0

do not want to return infinity to nonlin_least_squares */

if (x[0] == 0.0) {

f[i] = 10000.;

} else {

f[i] = t[i]/x[0] - y[i];

}

}

スレッドセーフ使用法

IMSL C ライブラリはOpenMPに基づいたスレッドセーフになっています。これが意味するのは、呼び

出しプログラムがガイドラインに従っていればマルチスレッド・アプリケーションから安全に呼び出す ことができることです。 特に、IMSL C ライブラリのエラー処理と I/O が理解されていなければなり ません。

エラー処理

IMSL C ライブラリのマルチスレッドアプリケーションでのエラー処理は、シングルスレッドの場合と

同じように振舞います。 主な違いは、IMSL C ライブラリのルーチンを呼ぶ各スレッドにエラースタ

ックがあることです。各スレッドに別々なエラースタックがあることにより、各スレッドがエラー処理

オプションをより適切に制御できます。各スレッドはIMSLのエラー処理ルーチン

imsl_error_options を用いて独自のオプションを設定することができます。スレッドで別々にエ

ラー処理オプションを設定する例は第12章「ユーティリティ」の imsl_error_options の例題3 を参照してください。

出力を生成するルーチン

IMSL C ライブラリのいくつかのルーチンは結果を出力することができます。 関数

imsl_output_file は出力するファイルを制御するのに使用されます。 シングルスレッドで実行さ

れるアプリケーションでは、imsl_output_file を1回呼び出して、出力するファイルを設定するこ

とができます。 マルチスレッドのアプリケーションでは、各スレッドは出力先のファイルのデフォル

ト設定を変更するのに imsl_output_file をそれぞれ呼ぶ必要があります。 詳細は ｢ユーティリ

ティ｣章の imsl_output_file の例題2を参照してください。

OpenMP

の使用法

IMSL C ライブラリのスレッドセーフは OpenMP に基づきます。ユーザはまたこのライブラリの関数

の OpenMP API 仕様のサポートにより共有メモリ並列処理を行なうことができます。これらの関数は

次のOpenMP アイコンで示されます。

OpenMP 並列はスレッドにより行なわれます。OpenMP プログラミングモデルでは、記憶域はマルチ

コアマシンのようにスレッド間で共有されます。 これらのスレッドはソースコードのOpenMPディレ

クティブにより生成されます。

IMSLライブラリのOpenMP の使用はユーザに透過的です。OpenMPディレクティブで強化されたコー ドは逐次型の実行環境でも適切に実行されます。エラー処理ルーチンは、並列実行での重大なエラーは ユーザに返されるように拡張されています。

OpenMP は次の主な目的で使用されます：

1. IMSL C ライブラリの中でスレッドセーフとするために。

2. データ並列化による計算量の多い関数の処理速度の向上のために。

3. 数値積分などのユーザ定義関数の処理を並列化するために。

この３番目の場合では、ユーザはユーザ定義関数がスレッドセーフであることを明示しなければなりま

せん。デフォルトではユーザ定義関数は並列に処理されません。 ユーザはユーティリティ

imsl_omp_options を使ってIMSLライブラリに渡されるルーチンがスレッドセーフであることを

明示することができます。

スレッドセーフであれば、関数は複数のスレッドにより同時に、かつ正しく実行されます。ユーザ定義

関数がスレッドセーフであることは非常に重要ですが、その理由は、OpenMP により生成される異なる

ないからです。従って並列化されたアルゴリズムはその逐次の「原型」と同じ様に動作することが必須 です。計算結果が実行される順序に依存する関数はスレッドセーフではなく、並列化の良い候補ではあ りません。広域データをアクセスしたり、修正したりする関数も同様です。

OpenMP 標準の仕様は (http://www.openmp.org/specifications/) にあります。

ベンダー提供のライブラリの使用法

IMSL C ライブラリは、IntelのMath Kernel Library (MKL) 、あるいは Sunの High Performance Library などのベンダー提供のライブラリの関数を利用する関数を含みます。ベンダー提供のライブラリの関

数はそれらの動作環境に十分に最適化されています。 これらの関数に対して、IMSL C ライブラリの

ユーザはIMSL 自身の関数かあるいはベンダー提供のライブラリの関数のどちらかにリンクするオプ

ションを持ちます。ベンダー提供のライブラリがある場合は、関数の説明の中で次のアイコンにより

示されます：

C++

からの使用法

IMSL C ライブラリの関数は C と C++ アプリケーションの両方から使用されます。 またIMSLライ

ブラリ関数を C++ クラスの中にラップすることもできます。

関数 imsl_f_int_fcn_sing はユーザ定義の関数を積分します。C++ では、ユーザ定義の関数は次

のように定義されるアブストラクトクラス IntFcnSingFunction のメンバー関数として定義されま

す。

#include <imsl.h> #include <math.h> #include <stdio.h> class IntFcnSingFunction {

public:

virtual float f(float x) = 0; };

関数 imsl_f_int_fcn_sing はC++ クラスIntFcnSing としてラップされます。この方法はオプ

ション引数 IMSL_FCN_W_DATA を使用して local_function を呼び、それがユーザ定義関数を計

算するメソッド f を呼びます。 簡単にするために、ここでは一つのオプション引数

IMSL_MAX_SUBINTER（部分区間の最大数）だけをラップします。同様の方法でもっと多くを含める ことができます。

#include <imsl.h>

class IntFcnSing

{

public:

int max_subinter;

IntFcnSing();

float integrate(IntFcnSingFunction *F, float a, float b);

};

static float local_function(float x, void *data)

{

IntFcnSingFunction *F = (IntFcnSingFunction*)data;

return F->f(x);

}

IntFcnSing::IntFcnSing()

{

max_subinter = 500;

}

float IntFcnSing::integrate(IntFcnSingFunction *F, float a, float b)

{

float result;

result = imsl_f_int_fcn_sing(NULL, a, b,

IMSL_FCN_W_DATA, local_function, F,

0);

if (imsl_error_type() >= 3)

{

throw imsl_error_message();

}

return result;

}

この IntFcnSing を使用するために、ユーザ定義関数は IntFcnSingFunction を拡張するクラス

のメソッド f として定義されなければなりません。 次のクラス MyClass は関数

_f

_x

_

_e

_

_ax

)

(

ここで_a はパラメータ、を定義します。

class MyClass : public IntFcnSingFunction

{

public:

MyClass();

float f(double x);

private:

float my_parameter;

};

MyClass::MyClass()

{

my_parameter = 5.0;

}

float MyClass::f(float x)

{

return exp(x) - my_parameter*x;

}

これらのクラスの使用例を以下に示します。C++ は重大あるいは停止の IMSL エラーの例外を出すの

で、imsl_error_options を呼ぶことによってこれらのエラーのプリントと処理の停止はオフにさ

れます。また、ユーザ定義関数はスレッドセーフであるので、imsl_omp_options を呼びこれを宣告

します。この設定により求積コードは OpenMP を使用して並列に計算されます。これらの両方の呼び

出しは実行毎に一回行うことが必要です。

この例題の２番目の部分は、エラー処理の方法を示すために部分区間の最大値を非現実的な小さい数5

にしています。 int main()

{

imsl_error_options(

IMSL_SET_PRINT, IMSL_FATAL, 0,

IMSL_SET_PRINT, IMSL_TERMINAL, 0,

IMSL_SET_STOP, IMSL_FATAL, 0,

IMSL_SET_STOP, IMSL_TERMINAL, 0,

0);

imsl_omp_options(IMSL_SET_FUNCTIONS_THREAD_SAFE, 1, 0);

IntFcnSing *intFcnSing = new IntFcnSing();

MyClass *myClass = new MyClass();

float x = intFcnSing->integrate(myClass, -1.0, 1.0);

printf("Solution in [-1,+1]: %g\n", x);

try {

intFcnSing->max_subinter = 5;

x = intFcnSing -> integrate (myClass, -100.0, 1000.0);

printf("Solution in [-100,1000]: %g\n", x);

} catch(char * exception) {

printf("Exception raised: %s\n", exception);

}

}

出力

Integral over [-1,+1] = 2.3504

行列格納モード

本節では、行列 (matrix) という単語は数学的オブジェクトを参照するために使用し、配列 (array) と

いう単語はCのデータ構造としての表現に使用します。 以下の配列タイプのリストで、IMSL C ライ

ブラリの関数は行列の次元数とエントリの全ての値からなる入力を必要します。これらの値はその配列 の中で行順で格納されます。

各関数は入力配列を処理し、一般的には、「結果」へのポインタを返します。 例えば、線形方程式を

解く際には、ポインタはその解に対するものです。 一般実固有値問題では、ポインタは固有値に対す

るものです。 通常、入力配列値は関数によって変更されません。

IMSL C ライブラリでは、配列はデータの連続するブロックへのポインタです。それらは、行列の行へ

のポインタのポインタではありません。宣言の例を示します： float *a = {1, 2, 3, 4};

float b[2][2] = {1, 2, 3, 4}; float c[] = {1, 2, 3, 4};

一般モード

一般行列は、正方 n ×n 行列です。 一般配列のデータタイプは、float、double、f_complex又は d_complex です。

矩形モード

矩形行列は m × n 行列です。矩形配列のデータタイプは、float、double、f_complex又はd_complex で す。

対称モード

対称行列は、AT _{= A のような正方 n ×n 行列です（行列A}T _{は Aの転置行列です）。} 対称配列のデータタイプは、float あるいは doubleです。

エルミートモード

エルミート行列は、

_A



_A



_A

行列

_A

A



A

A

_A



_A

疎行列の座標対応格納モード

疎行列では非ゼロの要素だけを関数に知らせる必要があります。疎行列の座標対応格納フォーマット

は、入力の行と列のインデックスと共に各エントリを格納します。次の4つのデータ構造体が定義さ

れます：

typedef struct { int row; int col; float val; } Imsl_f_sparse_elem;

typedef struct { int row; int col; double val; } Imsl_d_sparse_elem;

} Imsl_c_sparse_elem;

typedef struct { int row; int col; d_complex val; } Imsl_z_sparse_elem;

複素データタイプf_complex とd_complex の説明は本マニュアルの末尾の参照資料を参照して下さい。

これらの構造体ではデータタイプだけが違います。疎行列は、これらのデータタイプの1つで、疎行列

座標対応フォーマットを受け取る関数に渡されます。この配列の要素の数は疎行列中の非ゼロの数に等 しくなります。

例として、6 × 6の行列を考えてみます：

                               6 0 0 0 2 1 3 1 5 0 0 1 0 1 7 0 0 2 0 0 0 5 0 0 0 0 1 3 9 0 0 0 0 0 0 2 A

この行列Aは15の非ゼロ要素を持ち、その座標対応の表現は次の様になります。

行 0 1 1 1 2 3 3 3 4 4 4 4 5 5 5 列 0 1 2 3 2 0 3 4 0 3 4 5 0 1 5 値 2 9 3 1 5 2 7 1 1 5 1 3 1 2 6 この表現は順序に依存しないので、次の形式でも同じです。

行 5 4 3 0 5 1 2 1 4 3 1 4 3 5 4

列 0 0 0 0 1 1 2 2 3 3 3 4 4 5 5

値 1 1 2 2 2 9 5 3 5 7 1 1 1 6 3

このデータを、例えば、Imsl_f_sparse_elemのタイプの配列を初期化するのに使用するにはいくつかの方

法があります。 以下のプログラムを考えてみます：

#include <imsl.h> int main()

{

Imsl_f_sparse_elem a[] = { {0, 0, 2.0},

{1, 1, 9.0}, {1, 2, -3.0}, {1, 3, -1.0}, {2, 2, 5.0}, {3, 0, -2.0}, {3, 3, -7.0}, {3, 4, -1.0}, {4, 0, -1.0}, {4, 3, -5.0}, {4, 4, 1.0}, {4, 5, -3.0}, {5, 0, -1.0}, {5, 1, -2.0}, {5, 5, 6.0} };

b[0].row = b[0].col = 0; b[0].val = 2.0; b[1].row = b[1].col = 1; b[1].val = 9.0; b[2].row = 1; b[2].col = 2; b[2].val = -3.0; b[3].row = 1; b[3].col = 3; b[3].val = -1.0; b[4].row = b[4].col = 2; b[4].val = 5.0; b[5].row = 3; b[5].col = 0; b[5].val = -2.0; b[6].row = b[6].col = 3; b[6].val = -7.0; b[7].row = 3; b[7].col = 4; b[7].val = -1; b[8].row = 4; b[8].col = 0; b[8].val = -1.0; b[9].row = 4; b[9].col = 3; b[9].val = -5.0; b[10].row = b[10].col = 4; b[10].val = 1.0; b[11].row = 4; b[11].col = 5; b[11].val = -3.0; b[12].row = 5; b[12].col = 0; b[12].val = -1.0; b[13].row = 5; b[13] = 1; b[13].val = -2.0; b[14].row = b[14].col = 5; b[14].val = 6.0; }

aとbは共に疎行列Aを表現し、このモジュールの関数は、どの識別子が引数リストから送られたかに

は関係なく全く同じ結果となります。

疎の対称、あるいは、エルミート行列は特別なケースで、対角項と、上、又は、下三角だけを格納しま す。 例として、5 × 5の線形連立方程式を考えてみます：

































)

0

,

4

(

)

1

,

1

(

0

0

)

1

,

1

(

)

0

,

4

(

)

1

,

1

(

0

0

)

1

,

1

(

)

0

,

4

(

)

1

,

1

(

0

0

)

1

,

1

(

)

0

,

4

(

H

IMSLライブラリのエルミートと対称正定値用の連立方程式解法ルーチンは、対角項と下3角が指定さ

れることを期待しています。下3角の疎の座標対応形式は以下で与えられます。

行 0 1 2 3 1 2 3

列 0 1 2 3 0 1 2

値 (4,0) (4,0) (4,0) (4,0) (1,1) (1,1) (1,1) また、前と同様、以下の形式でも表現できます。

行 0 1 1 2 2 3 3

列 0 0 1 1 2 2 3

値 (4,0) (1,1) (4,0) (1,1) (4,0) (1,1) (4,0) 以下のプログラムは、aとbの両方をHに初期化します。

#include <imsl.h> int main()

{

Imsl_c_sparse_elem a[] = { {0, 0, {4.0, 0.0}}, {1, 1, {4.0, 0.0}}, {2, 2, {4.0, 0.0}}, {3, 3, {4.0, 0.0}}, {1, 0, {1.0, 1.0}}, {2, 1, {1.0, 1.0}}, {3, 2, {1.0, 1.0}} }

Imsl_c_sparse_elem b[7];

b[0].row = b[0].col = 0;

b[1].row = 1; b[1].col = 0;

b[1].val = imsl_cf_convert (1.0, 1.0); b[2].row = b[2].col = 1;

b[2].val = imsl_cf_convert (4.0, 0.0); b[3].row = 2; b[3].col = 1;

b[3].val = imsl_cf_convert (1.0, 1.0); b[4].row = b[4].col = 2;

b[4].val = imsl_cf_convert (4.0, 0.0); b[5].row = 3; b[5].col = 2;

b[5].val = imsl_cf_convert (1.0, 1.0); b[6].row = b[6].col = 3;

b[6].val = imsl_cf_convert (4.0, 0.0); }

ここで注意すべき幾つかの重要な点があります。 Hは対称ではなくエルミート形式です。 エルミー

トデータを受け取る関数はこれを理解して、次を仮定して処理します：

h



h

IMSL C ライブラリは、正定値ではない行列の対称性を有効に使用することはできません。ここに意味

することは、不定である対称行列はこの簡潔な対称モードで格納することができないということです。 上と下の両方の三角行列を指定して、疎の一般の解法ルーチンを呼ばねばなりません。

帯格納モード

帯行列は、その非ゼロ要素の全てが主対角項に ”近い” M xN 行列です。つまり、i - j > nlca又は j - i > nuca であれば Aij= 0 です。整数m = nlca + nuca + １ は、全バンド幅です。 主対角以外

の対角は共対角と呼ばれます。任意の M_xN 行列は帯行列ですが、帯格納フォーマットは、非ゼロの

共対角の数がNより非常に小さいときだけ有用です。

帯格納フォーマットでは、nlca 個の下共対角と、nuca 個の上共対角は大きさM × N の配列の行に 格納されます。その要素は、それらがその行列の中にあるように、配列の同じ列に格納されます。バン ド幅の内側の値Aij は、位置[(i - j+ nuca +１) * n + j ]に線形の配列で格納されます。これは行並

びの、その行列の2次元記述からの一次元写像となります。

例えば、1つの下と2つの上共対角を持つ5×5の行列Aを考えます。

                 4 , 4 3 , 4 4 , 3 3 , 3 2 , 3 4 , 2 3 , 2 2 , 2 1 , 2 3 , 1 2 , 1 1 , 1 0 , 1 2 , 0 1 , 0 0 , 0 0 0 0 0 0 0 0 0 0 A A A A A A A A A A A A A A A A A

帯格納フォーマットでは、このデータは次のように並べられます。

























0

0

0

0

A

A

A

A

A

A

A

A

A

A

A

A

A

A

A

A

このデータは長さ20の配列に連続して行並びで保存されます。

                 50 8 0 0 0 4 40 7 0 0 0 3 30 6 0 0 0 2 20 5 0 0 0 1 10 A

次の宜言は、この行列を帯格納フォーマットで保存します。 float a[] = {

0.0, 1.0, 2.0, 3.0, 4.0, 10.0, 20.0, 30.0, 40.0, 50.0, 5.0, 6.0, 7.0, 8.0, 0.0}；

疎行列の座標対応フォーマットのように、スペースを節約する帯格納の対称バージョンがあります。 例として以下の5 × 5の対称行列を見て下さい。



































0

0

0

0

0

0

A

A

A

A

A

A

A

A

A

A

A

A

A

A

A

A

A

A

A

A

帯対称格納フォーマットでは、データは次のように配列されます。





















0

0

0

A

A

A

A

A

A

A

A

A

A

A

A

次のエルミート形式の例は、この手順を示します。

                        ) 0 , 8 ( ) 1 , 1 ( ) 1 , 1 ( 0 0 ) 1 , 1 ( ) 0 , 8 ( ) 1 , 1 ( ) 1 , 1 ( 0 ) 1 , 1 ( ) 1 , 1 ( ) 0 , 8 ( ) 1 , 1 ( ) 1 , 1 ( 0 ) 1 , 1 ( ) 1 , 1 ( ) 0 , 8 ( ) 1 , 1 ( 0 0 ) 1 , 1 ( ) 1 , 1 ( ) 0 , 8 ( H

以下では、帯対称の格納フォーマットを使用して h に H を格納します。

f_complex h[ ] = {

{0.0, 0.0}, {0.0, 0.0}, {1.0, 1.0}, {1.0, 1.0}, {1.0, 1.0}, {0.0, 0.0}, {1.0, 1.O}, {1.0, 1.0}, {1.0, 1.0}, {1.0, 1.0}, {8.0, 0.0}, {8.0, 0.0}, {8.0, 0.0}, {8.0, 0.0}, {8.0, 0.0}

あるいは、

f_complex h[15]；

h[0] = h[l] = h[5] = imsl_cf_convert（0.0, 0.0）； h[2] = h[3] = h[4] = h[6] = h[7] = h[8] = h[9] =

imsl_cf_convert（1.0, 1.0）；

帯形式と座標対応形式の選択方法

どのような行列も疎の座標対応形式か帯形式で格納することができます。 選択は行列の疎の分布パタ ーンによります。主対角近くの帯の中に全ての非ゼロのデータを持つ行列は、帯形式が良いでしょう。 非ゼロの要素がほぼ一様に行列に散在するならば、疎の座標対応形式が最良の選択です。極端な例とし

て、次の2つのケースを考えます。

(1) 主対角と (0、n - 1) と (n - 1, 0) の要素が非ゼロの n × n 行列。疎の座標対応ベクトルは、n + 2 の 長さです。帯格納では長さ n (2n −１) の配列が必要とされ、密な解法ルーチンが要求するほぼ2倍の 記憶域となります。

(2) 全ての対角項と上対角項と下対角項が非ゼロである三重対角行列。 帯形式では長さ 3n の配列が

必要です。疎の座標対応形式では長さ 3n − 2 のベクトルが必要になります。しかし、例えば、単精度

浮動小数点では、座標対応形式の 3n −2 個のそれぞれが、帯形式に必要な 3n の3倍もの記憶域が必 要とされます。これは座標対応形式では、行と列のインデックスを伴うためです。一方帯形式では、順 序付けられたリストとリストの中の位置によって元の行列の中の位置を定義します。

圧縮された疎の列（

CSC : Compressed Sparse Column

）フォーマット

座標対応形式でのデータを受け取る関数は、Harwell-Boeing疎行列のユーザ・ガイドに記述されたフォ

ーマットで格納されたデータを受けることができます。 これは列指向で、各列は疎ベクトルとして保

持され、整数配列のエントリの行インデックスのリスト(下の "rowind")と別の float(double, f_complex, d_complex)配列の対応する値のリスト(下の "values")によって表現されます。各列のデー タは連続的に格納され、列は順に格納されます。3番目の配列（下の "colptr"）は元の疎行列の連続 する各列の最初の非ゼロ値を置く配列 "array" の位置を示します。 従って、 colptr[i] は元の 疎行列の i 番目の列からの値を格納する配列 "values" の最初の自由な位置のインデックスを含み

ます。 言い換えれば、values[colptr[i]] は元の疎行列のi番目の列の最初の非ゼロ値を持ちま

す。 対称行列とエルミート行列では下三角と対角のエントリだけが格納されます。 全ての配列は、

Harwell-Boeingテストセットの、1基準の配列と異なり、ゼロを基準にします。

Harwell-Boeingユーザ・ガイドにあるように、以下の例でこの格納方式を説明します。 5 × 5行列











































6

0

5

0

5

0

4

0

4

0

0

0

0

0

2

3

0

2

0

0

0

1

0

3

1

は次のように、配列colptr（最初のエントリの位置）、rowind（行インデックス）とvalues（非ゼロ のエントリ）に格納されます。

添字 0 1 2 3 4 5 6 7 8 9 10

colptr 0 3 5 7 9 11

rowind 0 4 2 3 0 1 4 0 3 4 1

values 1 5 2 4 -3 -2 -5 -1 -4 6 3

以下のプログラムは、CSCフォーマットと座標対応フォーマットの関係を示します。

void main() {

int i, j, k, n=5, nz, start, stop; int colptr[] = {0, 3, 5, 7, 9, 11};

Imsl_d_sparse_elem a[11]; k = 0;

for (i=0; i<n; i++) { start = colptr[i]; stop = colptr[i+1];

for (j=start; j<stop; j++) { a[k].row = rowind[j]; a[k].col = i;

a[k++].val = values[j]; }