- メインページ
- モジュール
- ネームスペース
- データ構造
- ファイル

- ファイル一覧
- グローバル

numeric.h
説明を見る。
1 //
2 // Copyright (c) 2003-2011, MIST Project, Nagoya University
3 // All rights reserved.
4 //
5 // Redistribution and use in source and binary forms, with or without modification,
6 // are permitted provided that the following conditions are met:
7 //
8 // 1. Redistributions of source code must retain the above copyright notice,
9 // this list of conditions and the following disclaimer.
10 //
11 // 2. Redistributions in binary form must reproduce the above copyright notice,
12 // this list of conditions and the following disclaimer in the documentation
13 // and/or other materials provided with the distribution.
14 //
15 // 3. Neither the name of the Nagoya University nor the names of its contributors
16 // may be used to endorse or promote products derived from this software
17 // without specific prior written permission.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
20 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
21 // FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
22 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
24 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
25 // IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
26 // THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 //
28
33 #ifndef __INCLUDE_NUMERIC__
34 #define __INCLUDE_NUMERIC__
35
36
37 #ifndef __INCLUDE_MIST_TYPE_TRAIT_H__
38 #include "config/type_trait.h"
39 #endif
40
41 // 行列クラスがインクルードされていない場合はインポートする．
42 #ifndef __INCLUDE_MIST_MATRIX__
43 #include "matrix.h"
44 #endif
45
46 #include <vector>
47 #include <algorithm>
48
49
50 // mist名前空間の始まり
51 _MIST_BEGIN
52
53
54
62
63
68 struct matrix_style
69 {
71 enum style
72 {
73 ge,
74 gb,
75 gt,
76 sy,
77 sb,
78 st,
79 he,
80 hb,
81 ht,
82 pb,
83 };
84 };
85
87 // 行列演算グループの終わり
88
89
90 struct numerical_exception
91 {
92 int error;
93 std::string message;
94
95 numerical_exception( ) : error( 0 ), message( )
96 {
97 }
98
99 numerical_exception( int err, const std::string &msg = "" ) : error( err ), message( msg )
100 {
101 }
102
103 numerical_exception( const std::string &msg ) : error( 0 ), message( msg )
104 {
105 }
106
107 numerical_exception( const numerical_exception &e ) : error( e.error ), message( e.message )
108 {
109 }
110 };
111
112
113 namespace __clapack__
114 {
115
116 // インテルのMKLとの互換性を保つための，関数名の変換マクロ
117 #if defined(_USE_INTEL_MATH_KERNEL_LIBRARY_) && _USE_INTEL_MATH_KERNEL_LIBRARY_ != 0
118 #define LPFNAME( name ) name // LAPACK用
119 #define BLFNAME( name ) name // BLAS用
120 #else
121 #define LPFNAME( name ) name ## _ // LAPACK用
122 #define BLFNAME( name ) f2c_ ## name // BLAS用
123 #endif
124
125
126 extern "C"
127 {
128 typedef long int integer;
129 typedef float real;
130 typedef double doublereal;
131 struct complex { real r, i; };
132 struct doublecomplex { doublereal r, i; };
133 //typedef struct { real r, i; } complex;
134 //typedef struct { doublereal r, i; } doublecomplex;
135
136 // 一般行列同士の掛け算を行う
137 int BLFNAME( sgemm ) ( char *transa, char *transb, integer *m, integer *n, integer *k, real *alpha, real *a, integer *lda, real *b, integer *ldb, real *beta, real *c__, integer *ldc );
138 int BLFNAME( dgemm ) ( char *transa, char *transb, integer *m, integer *n, integer *k, doublereal *alpha, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *beta, doublereal *c__, integer *ldc );
139 int BLFNAME( cgemm ) ( char *transa, char *transb, integer *m, integer *n, integer *k, complex *alpha, complex *a, integer *lda, complex *b, integer *ldb, complex *beta, complex *c__, integer *ldc );
140 int BLFNAME( zgemm ) ( char *transa, char *transb, integer *m, integer *n, integer *k, doublecomplex *alpha, doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, doublecomplex *beta, doublecomplex *c__, integer *ldc );
141
142 // 一般行列とベクトルの掛け算を行う
143 int BLFNAME( sgemv ) ( char *trans, integer *m, integer *n, real *alpha, real *a, integer *lda, real *x, integer *incx, real *beta, real *y, integer *incy );
144 int BLFNAME( dgemv ) ( char *trans, integer *m, integer *n, doublereal *alpha, doublereal *a, integer *lda, doublereal *x, integer *incx, doublereal *beta, doublereal *y, integer *incy );
145 int BLFNAME( cgemv ) ( char *trans, integer *m, integer *n, complex *alpha, complex *a, integer *lda, complex *x, integer *incx, complex *beta, complex *y, integer *incy );
146 int BLFNAME( zgemv ) ( char *trans, integer *m, integer *n, doublecomplex *alpha, doublecomplex *a, integer *lda, doublecomplex *x, integer *incx, doublecomplex *beta, doublecomplex *y, integer *incy );
147
148 // 一般正方行列の連立方程式を解く関数
149 int LPFNAME( sgesv ) ( integer *n, integer *nrhs, real *a, integer *lda, integer *ipiv, real *b, integer *ldb, integer *info );
150 int LPFNAME( dgesv ) ( integer *n, integer *nrhs, doublereal *a, integer *lda, integer *ipiv, doublereal *b, integer *ldb, integer *info );
151 int LPFNAME( cgesv ) ( integer *n, integer *nrhs, complex *a, integer *lda, integer *ipiv, complex *b, integer *ldb, integer *info );
152 int LPFNAME( zgesv ) ( integer *n, integer *nrhs, doublecomplex *a, integer *lda, integer *ipiv, doublecomplex *b, integer *ldb, integer *info );
153 // 一般帯行列の連立方程式を解く関数
154 int LPFNAME( sgbsv ) ( integer *n, integer *kl, integer *ku, integer *nrhs, real *ab, integer *ldab, integer *ipiv, real *b, integer *ldb, integer *info );
155 int LPFNAME( dgbsv ) ( integer *n, integer *kl, integer *ku, integer *nrhs, doublereal *ab, integer *ldab, integer *ipiv, doublereal *b, integer *ldb, integer *info );
156 int LPFNAME( cgbsv ) ( integer *n, integer *kl, integer *ku, integer *nrhs, complex *ab, integer *ldab, integer *ipiv, complex *b, integer *ldb, integer *info );
157 int LPFNAME( zgbsv ) ( integer *n, integer *kl, integer *ku, integer *nrhs, doublecomplex *ab, integer *ldab, integer *ipiv, doublecomplex *b, integer *ldb, integer *info );
158 // 一般三重対角行列の連立方程式を解く関数
159 int LPFNAME( sgtsv ) ( integer *n, integer *nrhs, real *dl, real *d__, real *du, real *b, integer *ldb, integer *info );
160 int LPFNAME( dgtsv ) ( integer *n, integer *nrhs, doublereal *dl, doublereal *d__, doublereal *du, doublereal *b, integer *ldb, integer *info );
161 int LPFNAME( cgtsv ) ( integer *n, integer *nrhs, complex *dl, complex *d__, complex *du, complex *b, integer *ldb, integer *info );
162 int LPFNAME( zgtsv ) ( integer *n, integer *nrhs, doublecomplex *dl, doublecomplex *d__, doublecomplex *du, doublecomplex *b, integer *ldb, integer *info );
163 // 対称正方行列の連立方程式を解く関数
164 int LPFNAME( ssysv ) ( char *uplo, integer *n, integer *nrhs, real *a, integer *lda, integer *ipiv, real *b, integer *ldb, real *work, integer *lwork, integer *info );
165 int LPFNAME( dsysv ) ( char *uplo, integer *n, integer *nrhs, doublereal *a, integer *lda, integer *ipiv, doublereal *b, integer *ldb, doublereal *work, integer *lwork, integer *info );
166 int LPFNAME( csysv ) ( char *uplo, integer *n, integer *nrhs, complex *a, integer *lda, integer *ipiv, complex *b, integer *ldb, complex *work, integer *lwork, integer *info );
167 int LPFNAME( zsysv ) ( char *uplo, integer *n, integer *nrhs, doublecomplex *a, integer *lda, integer *ipiv, doublecomplex *b, integer *ldb, doublecomplex *work, integer *lwork, integer *info );
168 // エルミート行列の連立方程式を解く関数
169 int LPFNAME( chesv ) ( char *uplo, integer *n, integer *nrhs, complex *a, integer *lda, integer *ipiv, complex *b, integer *ldb, complex *work, integer *lwork, integer *info );
170 int LPFNAME( zhesv ) ( char *uplo, integer *n, integer *nrhs, doublecomplex *a, integer *lda, integer *ipiv, doublecomplex *b, integer *ldb, doublecomplex *work, integer *lwork, integer *info );
171 // 対称帯行列の連立方程式を解く関数
172 int LPFNAME( spbsv ) ( char *uplo, integer *n, integer *kd, integer *nrhs, real *ab, integer *ldab, real *b, integer *ldb, integer *info );
173 int LPFNAME( dpbsv ) ( char *uplo, integer *n, integer *kd, integer *nrhs, doublereal *ab, integer *ldab, doublereal *b, integer *ldb, integer *info );
174 int LPFNAME( cpbsv ) ( char *uplo, integer *n, integer *kd, integer *nrhs, complex *ab, integer *ldab, complex *b, integer *ldb, integer *info );
175 int LPFNAME( zpbsv ) ( char *uplo, integer *n, integer *kd, integer *nrhs, doublecomplex *ab, integer *ldab, doublecomplex *b, integer *ldb, integer *info );
176
177 // 一般行列のLU分解
178 int LPFNAME( sgetrf ) ( integer *m, integer *n, real *a, integer *lda, integer *ipiv, integer *info );
179 int LPFNAME( dgetrf ) ( integer *m, integer *n, doublereal *a, integer *lda, integer *ipiv, integer *info );
180 int LPFNAME( cgetrf ) ( integer *m, integer *n, complex *a, integer *lda, integer *ipiv, integer *info );
181 int LPFNAME( zgetrf ) ( integer *m, integer *n, doublecomplex *a, integer *lda, integer *ipiv, integer *info );
182 // 一般帯行列の連立方程式を解く関数
183 int LPFNAME( sgbtrf ) ( integer *m, integer *n, integer *kl, integer *ku, real *ab, integer *ldab, integer *ipiv, integer *info );
184 int LPFNAME( dgbtrf ) ( integer *m, integer *n, integer *kl, integer *ku, doublereal *ab, integer *ldab, integer *ipiv, integer *info );
185 int LPFNAME( cgbtrf ) ( integer *m, integer *n, integer *kl, integer *ku, complex *ab, integer *ldab, integer *ipiv, integer *info );
186 int LPFNAME( zgbtrf ) ( integer *m, integer *n, integer *kl, integer *ku, doublecomplex *ab, integer *ldab, integer *ipiv, integer *info );
187 // 一般三重対角行列の連立方程式を解く関数
188 int LPFNAME( sgttrf ) ( integer *n, real *dl, real *d__, real *du, real *du2, integer *ipiv, integer *info );
189 int LPFNAME( dgttrf ) ( integer *n, doublereal *dl, doublereal *d__, doublereal *du, doublereal *du2, integer *ipiv, integer *info );
190 int LPFNAME( cgttrf ) ( integer *n, complex *dl, complex *d__, complex *du, complex *du2, integer *ipiv, integer *info );
191 int LPFNAME( zgttrf ) ( integer *n, doublecomplex *dl, doublecomplex *d__, doublecomplex *du, doublecomplex *du2, integer *ipiv, integer *info );
192 // 対称行列のLU分解
193 int LPFNAME( ssytrf ) ( char *uplo, integer *n, real *a, integer *lda, integer *ipiv, real *work, integer *lwork, integer *info );
194 int LPFNAME( dsytrf ) ( char *uplo, integer *n, doublereal *a, integer *lda, integer *ipiv, doublereal *work, integer *lwork, integer *info );
195 int LPFNAME( csytrf ) ( char *uplo, integer *n, complex *a, integer *lda, integer *ipiv, complex *work, integer *lwork, integer *info );
196 int LPFNAME( zsytrf ) ( char *uplo, integer *n, doublecomplex *a, integer *lda, integer *ipiv, doublecomplex *work, integer *lwork, integer *info );
197 // エルミート行列のLU分解
198 int LPFNAME( chetrf ) ( char *uplo, integer *n, complex *a, integer *lda, integer *ipiv, complex *work, integer *lwork, integer *info );
199 int LPFNAME( zhetrf ) ( char *uplo, integer *n, doublecomplex *a, integer *lda, integer *ipiv, doublecomplex *work, integer *lwork, integer *info );
200
201 // 実対称行列のコレスキー分解
202 int LPFNAME( spotrf ) ( char *uplo, integer *n, real *a, integer *lda, integer *info );
203 int LPFNAME( dpotrf ) ( char *uplo, integer *n, doublereal *a, integer *lda, integer *info );
204 int LPFNAME( cpotrf ) ( char *uplo, integer *n, complex *a, integer *lda, integer *info );
205 int LPFNAME( zpotrf ) ( char *uplo, integer *n, doublecomplex *a, integer *lda, integer *info );
206
207 // 一般行列のQR分解
208 int LPFNAME( sgeqrf ) ( integer *m, integer *n, real *a, integer *lda, real *tau, real *work, integer *lwork, integer *info );
209 int LPFNAME( dgeqrf ) ( integer *m, integer *n, doublereal *a, integer *lda, doublereal *tau, doublereal *work, integer *lwork, integer *info );
210 int LPFNAME( cgeqrf ) ( integer *m, integer *n, complex *a, integer *lda, complex *tau, complex *work, integer *lwork, integer *info );
211 int LPFNAME( zgeqrf ) ( integer *m, integer *n, doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *work, integer *lwork, integer *info );
212 // 一般行列のQR分解(Level3 Blas)
213 int LPFNAME( sgeqp3 ) ( integer *m, integer *n, real *a, integer *lda, integer *jpvt, real *tau, real *work, integer *lwork, integer *info );
214 int LPFNAME( dgeqp3 ) ( integer *m, integer *n, doublereal *a, integer *lda, integer *jpvt, doublereal *tau, doublereal *work, integer *lwork, integer *info );
215 int LPFNAME( cgeqp3 ) ( integer *m, integer *n, complex *a, integer *lda, integer *jpvt, complex *tau, complex *work, integer *lwork, real *rwork, integer *info );
216 int LPFNAME( zgeqp3 ) ( integer *m, integer *n, doublecomplex *a, integer *lda, integer *jpvt, doublecomplex *tau, doublecomplex *work, integer *lwork, doublereal *rwork, integer *info );
217 // 一般行列のQR分解後のQ行列を作成するルーチン
218 int LPFNAME( sorgqr ) ( integer *m, integer *n, integer *k, real *a, integer *lda, real *tau, real *work, integer *lwork, integer *info );
219 int LPFNAME( dorgqr ) ( integer *m, integer *n, integer *k, doublereal *a, integer *lda, doublereal *tau, doublereal *work, integer *lwork, integer *info );
220 int LPFNAME( cungqr ) ( integer *m, integer *n, integer *k, complex *a, integer *lda, complex *tau, complex *work, integer *lwork, integer *info );
221 int LPFNAME( zungqr ) ( integer *m, integer *n, integer *k, doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *work, integer *lwork, integer *info );
222
223 // LU分解の結果を用いた一般行列の逆行列の計算
224 int LPFNAME( sgetri ) ( integer *n, real *a, integer *lda, integer *ipiv, real *work, integer *lwork, integer *info );
225 int LPFNAME( dgetri ) ( integer *n, doublereal *a, integer *lda, integer *ipiv, doublereal *work, integer *lwork, integer *info );
226 int LPFNAME( cgetri ) ( integer *n, complex *a, integer *lda, integer *ipiv, complex *work, integer *lwork, integer *info );
227 int LPFNAME( zgetri ) ( integer *n, doublecomplex *a, integer *lda, integer *ipiv, doublecomplex *work, integer *lwork, integer *info );
228 // LU分解の結果を用いた対称行列の逆行列の計算
229 int LPFNAME( ssytri ) ( char *uplo, integer *n, real *a, integer *lda, integer *ipiv, real *work, integer *info );
230 int LPFNAME( dsytri ) ( char *uplo, integer *n, doublereal *a, integer *lda, integer *ipiv, doublereal *work, integer *info );
231 int LPFNAME( csytri ) ( char *uplo, integer *n, complex *a, integer *lda, integer *ipiv, complex *work, integer *info );
232 int LPFNAME( zsytri ) ( char *uplo, integer *n, doublecomplex *a, integer *lda, integer *ipiv, doublecomplex *work, integer *info );
233 // LU分解の結果を用いたエルミート行列の逆行列の計算
234 int LPFNAME( chetri ) ( char *uplo, integer *n, complex *a, integer *lda, integer *ipiv, complex *work, integer *info );
235 int LPFNAME( zhetri ) ( char *uplo, integer *n, doublecomplex *a, integer *lda, integer *ipiv, doublecomplex *work, integer *info );
236
237
238 // 一般行列に対する固有値・固有ベクトルを計算
239 int LPFNAME( sgeev ) ( char *jobvl, char *jobvr, integer *n, real *a, integer *lda, real *wr, real *wi, real *vl, integer *ldvl, real *vr, integer *ldvr, real *work, integer *lwork, integer *info );
240 int LPFNAME( dgeev ) ( char *jobvl, char *jobvr, integer *n, doublereal *a, integer *lda, doublereal *wr, doublereal *wi, doublereal *vl, integer *ldvl, doublereal *vr, integer *ldvr, doublereal *work, integer *lwork, integer *info );
241 int LPFNAME( cgeev ) ( char *jobvl, char *jobvr, integer *n, complex *a, integer *lda, complex *w, complex *vl, integer *ldvl, complex *vr, integer *ldvr, complex *work, integer *lwork, real *rwork, integer *info );
242 int LPFNAME( zgeev ) ( char *jobvl, char *jobvr, integer *n, doublecomplex *a, integer *lda, doublecomplex *w, doublecomplex *vl, integer *ldvl, doublecomplex *vr, integer *ldvr, doublecomplex *work, integer *lwork, doublereal *rwork, integer *info );
243 // 対称行列に対する固有値・固有ベクトルを計算
244 int LPFNAME( ssyev ) ( char *jobz, char *uplo, integer *n, real *a, integer *lda, real *w, real *work, integer *lwork, integer *info );
245 int LPFNAME( dsyev ) ( char *jobz, char *uplo, integer *n, doublereal *a, integer *lda, doublereal *w, doublereal *work, integer *lwork, integer *info );
246 // 対称帯行列に対する固有値・固有ベクトルを計算
247 int LPFNAME( ssbev ) ( char *jobz, char *uplo, integer *n, integer *kd, real *ab, integer *ldab, real *w, real *z__, integer *ldz, real *work, integer *info );
248 int LPFNAME( dsbev ) ( char *jobz, char *uplo, integer *n, integer *kd, doublereal *ab, integer *ldab, doublereal *w, doublereal *z__, integer *ldz, doublereal *work, integer *info );
249 // 対称3重対角行列に対する固有値・固有ベクトルを計算
250 int LPFNAME( sstev ) ( char *jobz, integer *n, real *d__, real *e, real *z__, integer *ldz, real *work, integer *info );
251 int LPFNAME( dstev ) ( char *jobz, integer *n, doublereal *d__, doublereal *e, doublereal *z__, integer *ldz, doublereal *work, integer *info );
252 // エルミート行列に対する固有値・固有ベクトルを計算
253 int LPFNAME( cheev ) ( char *jobz, char *uplo, integer *n, complex *a, integer *lda, real *w, complex *work, integer *lwork, real *rwork, integer *info );
254 int LPFNAME( zheev ) ( char *jobz, char *uplo, integer *n, doublecomplex *a, integer *lda, doublereal *w, doublecomplex *work, integer *lwork, doublereal *rwork, integer *info );
255 // エルミート帯行列に対する固有値・固有ベクトルを計算
256 int LPFNAME( chbev ) ( char *jobz, char *uplo, integer *n, integer *kd, complex *ab, integer *ldab, real *w, complex *z__, integer *ldz, complex *work, real *rwork, integer *info );
257 int LPFNAME( zhbev ) ( char *jobz, char *uplo, integer *n, integer *kd, doublecomplex *ab, integer *ldab, doublereal *w, doublecomplex *z__, integer *ldz, doublecomplex *work, doublereal *rwork, integer *info );
258
259
260 // 一般行列に対する固有値・固有ベクトルを計算．行列を対角化して精度を上げるバージョン
261 int LPFNAME( sgeevx ) ( char *balanc, char *jobvl, char *jobvr, char *sense, integer *n, real *a, integer *lda, real *wr, real *wi, real *vl, integer *ldvl, real *vr, integer *ldvr, integer *ilo, integer *ihi, real *scale, real *abnrm, real *rconde, real *rcondv, real *work, integer *lwork, integer *iwork, integer *info );
262 int LPFNAME( dgeevx ) ( char *balanc, char *jobvl, char *jobvr, char *sense, integer *n, doublereal *a, integer *lda, doublereal *wr, doublereal *wi, doublereal *vl, integer *ldvl, doublereal *vr, integer *ldvr, integer *ilo, integer *ihi, doublereal *scale, doublereal *abnrm, doublereal *rconde, doublereal *rcondv, doublereal *work, integer *lwork, integer *iwork, integer *info );
263 int LPFNAME( cgeevx ) ( char *balanc, char *jobvl, char *jobvr, char *sense, integer *n, complex *a, integer *lda, complex *w, complex *vl, integer *ldvl, complex *vr, integer *ldvr, integer *ilo, integer *ihi, real *scale, real *abnrm, real *rconde, real *rcondv, complex *work, integer *lwork, real *rwork, integer *info );
264 int LPFNAME( zgeevx ) ( char *balanc, char *jobvl, char *jobvr, char *sense, integer *n, doublecomplex *a, integer *lda, doublecomplex *w, doublecomplex *vl, integer *ldvl, doublecomplex *vr, integer *ldvr, integer *ilo, integer *ihi, doublereal *scale, doublereal *abnrm, doublereal *rconde, doublereal *rcondv, doublecomplex *work, integer *lwork, doublereal *rwork, integer *info );
265 // 対称行列に対する固有値・固有ベクトルを計算．行列を対角化して精度を上げるバージョン
266 int LPFNAME( ssyevx ) ( char *jobz, char *range, char *uplo, integer *n, real *a, integer *lda, real *vl, real *vu, integer *il, integer *iu, real *abstol, integer *m, real *w, real *z__, integer *ldz, real *work, integer *lwork, integer *iwork, integer *ifail, integer *info );
267 int LPFNAME( dsyevx ) ( char *jobz, char *range, char *uplo, integer *n, doublereal *a, integer *lda, doublereal *vl, doublereal *vu, integer *il, integer *iu, doublereal *abstol, integer *m, doublereal *w, doublereal *z__, integer *ldz, doublereal *work, integer *lwork, integer *iwork, integer *ifail, integer *info );
268 // 対称帯行列に対する固有値・固有ベクトルを計算．行列を対角化して精度を上げるバージョン
269 int LPFNAME( ssbevx ) ( char *jobz, char *range, char *uplo, integer *n, integer *kd, real *ab, integer *ldab, real *q, integer *ldq, real *vl, real *vu, integer *il, integer *iu, real *abstol, integer *m, real *w, real *z__, integer *ldz, real *work, integer *iwork, integer *ifail, integer *info );
270 int LPFNAME( dsbevx ) ( char *jobz, char *range, char *uplo, integer *n, integer *kd, doublereal *ab, integer *ldab, doublereal *q, integer *ldq, doublereal *vl, doublereal *vu, integer *il, integer *iu, doublereal *abstol, integer *m, doublereal *w, doublereal *z__, integer *ldz, doublereal *work, integer *iwork, integer *ifail, integer *info );
271 // 対称3重対角行列に対する固有値・固有ベクトルを計算．行列を対角化して精度を上げるバージョン
272 int LPFNAME( sstevx ) ( char *jobz, char *range, integer *n, real *d__, real *e, real *vl, real *vu, integer *il, integer *iu, real *abstol, integer *m, real *w, real *z__, integer *ldz, real *work, integer *iwork, integer *ifail, integer *info );
273 int LPFNAME( dstevx ) ( char *jobz, char *range, integer *n, doublereal *d__, doublereal *e, doublereal *vl, doublereal *vu, integer *il, integer *iu, doublereal *abstol, integer *m, doublereal *w, doublereal *z__, integer *ldz, doublereal *work, integer *iwork, integer *ifail, integer *info );
274 // エルミート行列に対する固有値・固有ベクトルを計算．行列を対角化して精度を上げるバージョン
275 int LPFNAME( cheevx ) ( char *jobz, char *range, char *uplo, integer *n, complex *a, integer *lda, real *vl, real *vu, integer *il, integer *iu, real *abstol, integer *m, real *w, complex *z__, integer *ldz, complex *work, integer *lwork, real *rwork, integer *iwork, integer *ifail, integer *info );
276 int LPFNAME( zheevx ) ( char *jobz, char *range, char *uplo, integer *n, doublecomplex *a, integer *lda, doublereal *vl, doublereal *vu, integer *il, integer *iu, doublereal *abstol, integer *m, doublereal *w, doublecomplex *z__, integer *ldz, doublecomplex *work, integer *lwork, doublereal *rwork, integer *iwork, integer *ifail, integer *info );
277 // エルミート帯行列に対する固有値・固有ベクトルを計算．行列を対角化して精度を上げるバージョン
278 int LPFNAME( chbevx ) ( char *jobz, char *range, char *uplo, integer *n, integer *kd, complex *ab, integer *ldab, complex *q, integer *ldq, real *vl, real *vu, integer *il, integer *iu, real *abstol, integer *m, real *w, complex *z__, integer *ldz, complex *work, real *rwork, integer *iwork, integer *ifail, integer *info );
279 int LPFNAME( zhbevx ) ( char *jobz, char *range, char *uplo, integer *n, integer *kd, doublecomplex *ab, integer *ldab, doublecomplex *q, integer *ldq, doublereal *vl, doublereal *vu, integer *il, integer *iu, doublereal *abstol, integer *m, doublereal *w, doublecomplex *z__, integer *ldz, doublecomplex *work, doublereal *rwork, integer *iwork, integer *ifail, integer *info );
280 // エルミート3重対角行列に対する固有値・固有ベクトルを計算．行列を対角化して精度を上げるバージョン
281
282
283 // 一般行列に対する特異値分解を計算
284 int LPFNAME( sgesvd ) ( char *jobu, char *jobvt, integer *m, integer *n, real *a, integer *lda, real *s, real *u, integer *ldu, real *vt, integer *ldvt, real *work, integer *lwork, integer *info );
285 int LPFNAME( dgesvd ) ( char *jobu, char *jobvt, integer *m, integer *n, doublereal *a, integer *lda, doublereal *s, doublereal *u, integer *ldu, doublereal *vt, integer *ldvt, doublereal *work, integer *lwork, integer *info );
286 int LPFNAME( cgesvd ) ( char *jobu, char *jobvt, integer *m, integer *n, complex *a, integer *lda, real *s, complex *u, integer *ldu, complex *vt, integer *ldvt, complex *work, integer *lwork, real *rwork, integer *info );
287 int LPFNAME( zgesvd ) ( char *jobu, char *jobvt, integer *m, integer *n, doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u, integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work, integer *lwork, doublereal *rwork, integer *info );
288
289
290 // 一般行列に対する特異値分解を計算．分割統治法を用いた高速化バージョン
291 int LPFNAME( sgesdd ) ( char *jobz, integer *m, integer *n, real *a, integer *lda, real *s, real *u, integer *ldu, real *vt, integer *ldvt, real *work, integer *lwork, integer *iwork, integer *info );
292 int LPFNAME( dgesdd ) ( char *jobz, integer *m, integer *n, doublereal *a, integer *lda, doublereal *s, doublereal *u, integer *ldu, doublereal *vt, integer *ldvt, doublereal *work, integer *lwork, integer *iwork, integer *info );
293 int LPFNAME( cgesdd ) ( char *jobz, integer *m, integer *n, complex *a, integer *lda, real *s, complex *u, integer *ldu, complex *vt, integer *ldvt, complex *work, integer *lwork, real *rwork, integer *iwork, integer *info );
294 int LPFNAME( zgesdd ) ( char *jobz, integer *m, integer *n, doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u, integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work, integer *lwork, doublereal *rwork, integer *iwork, integer *info );
295
296
297
298 // 対称行列に対する固有値・固有ベクトルを計算を計算．分割統治法を用いた高速バージョン
299 int LPFNAME( ssyevd ) ( char *jobz, char *uplo, integer *n, real *a, integer *lda, real *w, real *work, integer *lwork, integer *iwork, integer *liwork, integer *info );
300 int LPFNAME( dsyevd ) ( char *jobz, char *uplo, integer *n, doublereal *a, integer *lda, doublereal *w, doublereal *work, integer *lwork, integer *iwork, integer *liwork, integer *info );
301 // 対称帯行列に対する固有値・固有ベクトルを計算を計算．分割統治法を用いた高速バージョン
302 int LPFNAME( ssbevd ) ( char *jobz, char *uplo, integer *n, integer *kd, real *ab, integer *ldab, real *w, real *z__, integer *ldz, real *work, integer *lwork, integer *iwork, integer *liwork, integer *info );
303 int LPFNAME( dsbevd ) ( char *jobz, char *uplo, integer *n, integer *kd, doublereal *ab, integer *ldab, doublereal *w, doublereal *z__, integer *ldz, doublereal *work, integer *lwork, integer *iwork, integer *liwork, integer *info );
304 // 対称3重対角行列に対する固有値・固有ベクトルを計算を計算．分割統治法を用いた高速バージョン
305 int LPFNAME( sstevd ) ( char *jobz, integer *n, real *d__, real *e, real *z__, integer *ldz, real *work, integer *lwork, integer *iwork, integer *liwork, integer *info );
306 int LPFNAME( dstevd ) ( char *jobz, integer *n, doublereal *d__, doublereal *e, doublereal *z__, integer *ldz, doublereal *work, integer *lwork, integer *iwork, integer *liwork, integer *info );
307 // エルミート行列に対する固有値・固有ベクトルを計算を計算．分割統治法を用いた高速バージョン
308 int LPFNAME( cheevd ) ( char *jobz, char *uplo, integer *n, complex *a, integer *lda, real *w, complex *work, integer *lwork, real *rwork, integer *lrwork, integer *iwork, integer *liwork, integer *info );
309 int LPFNAME( zheevd ) ( char *jobz, char *uplo, integer *n, doublecomplex *a, integer *lda, doublereal *w, doublecomplex *work, integer *lwork, doublereal *rwork, integer *lrwork, integer *iwork, integer *liwork, integer *info );
310 // エルミート帯行列に対する固有値・固有ベクトルを計算を計算．分割統治法を用いた高速バージョン
311 int LPFNAME( chbevd ) ( char *jobz, char *uplo, integer *n, integer *kd, complex *ab, integer *ldab, real *w, complex *z__, integer *ldz, complex *work, integer *lwork, real *rwork, integer *lrwork, integer *iwork, integer *liwork, integer *info );
312 int LPFNAME( zhbevd ) ( char *jobz, char *uplo, integer *n, integer *kd, doublecomplex *ab, integer *ldab, doublereal *w, doublecomplex *z__, integer *ldz, doublecomplex *work, integer *lwork, doublereal *rwork, integer *lrwork, integer *iwork, integer *liwork, integer *info );
313 }
314
315
316 // 複素数及び実数の両方を区別することなく，実数部の値を取り出す関数
317 inline real get_real( const real &r ){ return( r ); }
318 inline doublereal get_real( const doublereal &r ){ return( r ); }
319 inline real get_real( const std::complex< real > &r ){ return( r.real( ) ); }
320 inline doublereal get_real( const std::complex< doublereal > &r ){ return( r.real( ) ); }
321
322 // 以降，BLASルーチン
323
324 // 一般行列同士の掛け算を行う
325 inline int gemm( char *transa, char *transb, integer &m, integer &n, integer &k, real &alpha, real *a, integer &lda, real *b, integer &ldb, real &beta, real *c__, integer &ldc )
326 {
327 return( BLFNAME( sgemm ) ( transa, transb, &m, &n, &k, &alpha, a, &lda, b, &ldb, &beta, c__, &ldc ) );
328 }
329 inline int gemm( char *transa, char *transb, integer &m, integer &n, integer &k, doublereal &alpha, doublereal *a, integer &lda, doublereal *b, integer &ldb, doublereal &beta, doublereal *c__, integer &ldc )
330 {
331 return( BLFNAME( dgemm ) ( transa, transb, &m, &n, &k, &alpha, a, &lda, b, &ldb, &beta, c__, &ldc ) );
332 }
333 inline int gemm( char *transa, char *transb, integer &m, integer &n, integer &k, std::complex< real > &alpha, std::complex< real > *a, integer &lda, std::complex< real > *b, integer &ldb, std::complex< real > &beta, std::complex< real > *c__, integer &ldc )
334 {
335 return( BLFNAME( cgemm ) ( transa, transb, &m, &n, &k, reinterpret_cast< complex* >( &alpha ), reinterpret_cast< complex* >( a ), &lda, reinterpret_cast< complex* >( b ), &ldb, reinterpret_cast< complex* >( &beta ), reinterpret_cast< complex* >( c__ ), &ldc ) );
336 }
337 inline int gemm( char *transa, char *transb, integer &m, integer &n, integer &k, std::complex< doublereal > &alpha, std::complex< doublereal > *a, integer &lda, std::complex< doublereal > *b, integer &ldb, std::complex< doublereal > &beta, std::complex< doublereal > *c__, integer &ldc )
338 {
339 return( BLFNAME( zgemm ) ( transa, transb, &m, &n, &k, reinterpret_cast< doublecomplex* >( &alpha ), reinterpret_cast< doublecomplex* >( a ), &lda, reinterpret_cast< doublecomplex* >( b ), &ldb, reinterpret_cast< doublecomplex* >( &beta ), reinterpret_cast< doublecomplex* >( c__ ), &ldc ) );
340 }
341
342 // 一般行列とベクトルの掛け算を行う
343 inline int gemv( char *trans, integer &m, integer &n, real &alpha, real *a, integer &lda, real *x, integer &incx, real &beta, real *y, integer &incy )
344 {
345 return( BLFNAME( sgemv ) ( trans, &m, &n, &alpha, a, &lda, x, &incx, &beta, y, &incy ) );
346 }
347 inline int gemv( char *trans, integer &m, integer &n, doublereal &alpha, doublereal *a, integer &lda, doublereal *x, integer &incx, doublereal &beta, doublereal *y, integer &incy )
348 {
349 return( BLFNAME( dgemv ) ( trans, &m, &n, &alpha, a, &lda, x, &incx, &beta, y, &incy ) );
350 }
351 inline int gemv( char *trans, integer &m, integer &n, std::complex< real > &alpha, std::complex< real > *a, integer &lda, std::complex< real > *x, integer &incx, std::complex< real > &beta, std::complex< real > *y, integer &incy )
352 {
353 return( BLFNAME( cgemv ) ( trans, &m, &n, reinterpret_cast< complex* >( &alpha ), reinterpret_cast< complex* >( a ), &lda, reinterpret_cast< complex* >( x ), &incx, reinterpret_cast< complex* >( &beta ), reinterpret_cast< complex* >( y ), &incy ) );
354 }
355 inline int gemv( char *trans, integer &m, integer &n, std::complex< doublereal > &alpha, std::complex< doublereal > *a, integer &lda, std::complex< doublereal > *x, integer &incx, std::complex< doublereal > &beta, std::complex< doublereal > *y, integer &incy )
356 {
357 return( BLFNAME( zgemv ) ( trans, &m, &n, reinterpret_cast< doublecomplex* >( &alpha ), reinterpret_cast< doublecomplex* >( a ), &lda, reinterpret_cast< doublecomplex* >( x ), &incx, reinterpret_cast< doublecomplex* >( &beta ), reinterpret_cast< doublecomplex* >( y ), &incy ) );
358 }
359
360
361 // 以降，LAPACKルーチン
362
363 // 一般正方行列の連立方程式を解く関数
364 inline int gesv( integer &n, integer &nrhs, real *a, integer &lda, integer *ipiv, real *b, integer &ldb, integer &info )
365 {
366 return( LPFNAME( sgesv ) ( &n, &nrhs, a, &lda, ipiv, b, &ldb, &info ) );
367 }
368 inline int gesv( integer &n, integer &nrhs, doublereal *a, integer &lda, integer *ipiv, doublereal *b, integer &ldb, integer &info )
369 {
370 return( LPFNAME( dgesv ) ( &n, &nrhs, a, &lda, ipiv, b, &ldb, &info ) );
371 }
372 inline int gesv( integer &n, integer &nrhs, std::complex< real > *a, integer &lda, integer *ipiv, std::complex< real > *b, integer &ldb, integer &info )
373 {
374 return( LPFNAME( cgesv ) ( &n, &nrhs, reinterpret_cast< complex* >( a ), &lda, ipiv, reinterpret_cast< complex* >( b ), &ldb, &info ) );
375 }
376 inline int gesv( integer &n, integer &nrhs, std::complex< doublereal > *a, integer &lda, integer *ipiv, std::complex< doublereal > *b, integer &ldb, integer &info )
377 {
378 return( LPFNAME( zgesv ) ( &n, &nrhs, reinterpret_cast< doublecomplex* >( a ), &lda, ipiv, reinterpret_cast< doublecomplex* >( b ), &ldb, &info ) );
379 }
380 // 一般帯行列の連立方程式を解く関数
381 inline int gbsv( integer &n, integer &kl, integer &ku, integer &nrhs, real *ab, integer &ldab, integer *ipiv, real *b, integer &ldb, integer &info )
382 {
383 return( LPFNAME( sgbsv ) ( &n, &kl, &ku, &nrhs, ab, &ldab, ipiv, b, &ldb, &info ) );
384 }
385 inline int gbsv( integer &n, integer &kl, integer &ku, integer &nrhs, doublereal *ab, integer &ldab, integer *ipiv, doublereal *b, integer &ldb, integer &info )
386 {
387 return( LPFNAME( dgbsv ) ( &n, &kl, &ku, &nrhs, ab, &ldab, ipiv, b, &ldb, &info ) );
388 }
389 inline int gbsv( integer &n, integer &kl, integer &ku, integer &nrhs, complex *ab, integer &ldab, integer *ipiv, complex *b, integer &ldb, integer &info )
390 {
391 return( LPFNAME( cgbsv ) ( &n, &kl, &ku, &nrhs, reinterpret_cast< complex* >( ab ), &ldab, ipiv, reinterpret_cast< complex* >( b ), &ldb, &info ) );
392 }
393 inline int gbsv( integer &n, integer &kl, integer &ku, integer &nrhs, doublecomplex *ab, integer &ldab, integer *ipiv, doublecomplex *b, integer &ldb, integer &info )
394 {
395 return( LPFNAME( zgbsv ) ( &n, &kl, &ku, &nrhs, reinterpret_cast< doublecomplex* >( ab ), &ldab, ipiv, reinterpret_cast< doublecomplex* >( b ), &ldb, &info ) );
396 }
397 // 一般三重対角行列の連立方程式を解く関数
398 inline int gtsv( integer &n, integer &nrhs, real *dl, real *d__, real *du, real *b, integer &ldb, integer &info )
399 {
400 return( LPFNAME( sgtsv ) ( &n, &nrhs, dl, d__, du, b, &ldb, &info ) );
401 }
402 inline int gtsv( integer &n, integer &nrhs, doublereal *dl, doublereal *d__, doublereal *du, doublereal *b, integer &ldb, integer &info )
403 {
404 return( LPFNAME( dgtsv ) ( &n, &nrhs, dl, d__, du, b, &ldb, &info ) );
405 }
406 inline int gtsv( integer &n, integer &nrhs, complex *dl, complex *d__, complex *du, complex *b, integer &ldb, integer &info )
407 {
408 return( LPFNAME( cgtsv ) ( &n, &nrhs, reinterpret_cast< complex* >( dl ), reinterpret_cast< complex* >( d__ ), reinterpret_cast< complex* >( du ), reinterpret_cast< complex* >( b ), &ldb, &info ) );
409 }
410 inline int gtsv( integer &n, integer &nrhs, doublecomplex *dl, doublecomplex *d__, doublecomplex *du, doublecomplex *b, integer &ldb, integer &info )
411 {
412 return( LPFNAME( zgtsv ) ( &n, &nrhs, reinterpret_cast< doublecomplex* >( dl ), reinterpret_cast< doublecomplex* >( d__ ), reinterpret_cast< doublecomplex* >( du ), reinterpret_cast< doublecomplex* >( b ), &ldb, &info ) );
413 }
414 // 対称正方行列の連立方程式を解く関数
415 inline int sysv( char *uplo, integer &n, integer &nrhs, real *a, integer &lda, integer *ipiv, real *b, integer &ldb, real *work, integer &lwork, integer &info )
416 {
417 return( LPFNAME( ssysv ) ( uplo, &n, &nrhs, a, &lda, ipiv, b, &ldb, work, &lwork, &info ) );
418 }
419 inline int sysv( char *uplo, integer &n, integer &nrhs, doublereal *a, integer &lda, integer *ipiv, doublereal *b, integer &ldb, doublereal *work, integer &lwork, integer &info )
420 {
421 return( LPFNAME( dsysv ) ( uplo, &n, &nrhs, a, &lda, ipiv, b, &ldb, work, &lwork, &info ) );
422 }
423 inline int sysv( char *uplo, integer &n, integer &nrhs, std::complex< real > *a, integer &lda, integer *ipiv, std::complex< real > *b, integer &ldb, std::complex< real > *work, integer &lwork, integer &info )
424 {
425 return( LPFNAME( csysv ) ( uplo, &n, &nrhs, reinterpret_cast< complex* >( a ), &lda, ipiv, reinterpret_cast< complex* >( b ), &ldb, reinterpret_cast< complex* >( work ), &lwork, &info ) );
426 }
427 inline int sysv( char *uplo, integer &n, integer &nrhs, std::complex< doublereal > *a, integer &lda, integer *ipiv, std::complex< doublereal > *b, integer &ldb, std::complex< doublereal > *work, integer &lwork, integer &info )
428 {
429 return( LPFNAME( zsysv ) ( uplo, &n, &nrhs, reinterpret_cast< doublecomplex* >( a ), &lda, ipiv, reinterpret_cast< doublecomplex* >( b ), &ldb, reinterpret_cast< doublecomplex* >( work ), &lwork, &info ) );
430 }
431 // エルミート行列の連立方程式を解く関数
432 inline int hesv( char *uplo, integer &n, integer &nrhs, std::complex< real > *a, integer &lda, integer *ipiv, std::complex< real > *b, integer &ldb, std::complex< real > *work, integer &lwork, integer&info )
433 {
434 return( LPFNAME( chesv ) ( uplo, &n, &nrhs, reinterpret_cast< complex* >( a ), &lda, ipiv, reinterpret_cast< complex* >( b ), &ldb, reinterpret_cast< complex* >( work ), &lwork, &info ) );
435 }
436 inline int hesv( char *uplo, integer &n, integer &nrhs, std::complex< doublereal > *a, integer &lda, integer *ipiv, std::complex< doublereal > *b, integer &ldb, std::complex< doublereal > *work, integer &lwork, integer &info )
437 {
438 return( LPFNAME( zsysv ) ( uplo, &n, &nrhs, reinterpret_cast< doublecomplex* >( a ), &lda, ipiv, reinterpret_cast< doublecomplex* >( b ), &ldb, reinterpret_cast< doublecomplex* >( work ), &lwork, &info ) );
439 }
440 // 対称帯行列の連立方程式を解く関数
441 inline int pbsv( char *uplo, integer &n, integer &kd, integer &nrhs, real *ab, integer &ldab, real *b, integer &ldb, integer &info )
442 {
443 return( LPFNAME( spbsv ) ( uplo, &n, &kd, &nrhs, ab, &ldab, b, &ldb, &info ) );
444 }
445 inline int pbsv( char *uplo, integer &n, integer &kd, integer &nrhs, doublereal *ab, integer &ldab, doublereal *b, integer &ldb, integer &info )
446 {
447 return( LPFNAME( dpbsv ) ( uplo, &n, &kd, &nrhs, ab, &ldab, b, &ldb, &info ) );
448 }
449 inline int pbsv( char *uplo, integer &n, integer &kd, integer &nrhs, std::complex< real > *ab, integer &ldab, std::complex< real > *b, integer &ldb, integer &info )
450 {
451 return( LPFNAME( cpbsv ) ( uplo, &n, &kd, &nrhs, reinterpret_cast< complex* >( ab ), &ldab, reinterpret_cast< complex* >( b ), &ldb, &info ) );
452 }
453 inline int pbsv( char *uplo, integer &n, integer &kd, integer &nrhs, std::complex< doublereal > *ab, integer &ldab, std::complex< doublereal > *b, integer &ldb, integer &info )
454 {
455 return( LPFNAME( zpbsv ) ( uplo, &n, &kd, &nrhs, reinterpret_cast< doublecomplex* >( ab ), &ldab, reinterpret_cast< doublecomplex* >( b ), &ldb, &info ) );
456 }
457
458 // 一般行列のLU分解
459 inline int getrf( integer &m, integer &n, real *a, integer &lda, integer *ipiv, integer &info )
460 {
461 return( LPFNAME( sgetrf ) ( &m, &n, a, &lda, ipiv, &info ) );
462 }
463 inline int getrf( integer &m, integer &n, doublereal *a, integer &lda, integer *ipiv, integer &info )
464 {
465 return( LPFNAME( dgetrf ) ( &m, &n, a, &lda, ipiv, &info ) );
466 }
467 inline int getrf( integer &m, integer &n, std::complex< real > *a, integer &lda, integer *ipiv, integer &info )
468 {
469 return( LPFNAME( cgetrf ) ( &m, &n, reinterpret_cast< complex* >( a ), &lda, ipiv, &info ) );
470 }
471 inline int getrf( integer &m, integer &n, std::complex< doublereal > *a, integer &lda, integer *ipiv, integer &info )
472 {
473 return( LPFNAME( zgetrf ) ( &m, &n, reinterpret_cast< doublecomplex* >( a ), &lda, ipiv, &info ) );
474 }
475 // 一般帯行列の連立方程式を解く関数
476 inline int gbtrf( integer &m, integer &n, integer &kl, integer &ku, real *ab, integer &ldab, integer *ipiv, integer &info )
477 {
478 return( LPFNAME( sgbtrf ) ( &m, &n, &kl, &ku, ab, &ldab, ipiv, &info ) );
479 }
480 inline int gbtrf( integer &m, integer &n, integer &kl, integer &ku, doublereal *ab, integer &ldab, integer *ipiv, integer &info )
481 {
482 return( LPFNAME( dgbtrf ) ( &m, &n, &kl, &ku, ab, &ldab, ipiv, &info ) );
483 }
484 inline int gbtrf( integer &m, integer &n, integer &kl, integer &ku, std::complex< real > *ab, integer &ldab, integer *ipiv, integer &info )
485 {
486 return( LPFNAME( cgbtrf ) ( &m, &n, &kl, &ku, reinterpret_cast< complex* >( ab ), &ldab, ipiv, &info ) );
487 }
488 inline int gbtrf( integer &m, integer &n, integer &kl, integer &ku, std::complex< doublereal > *ab, integer &ldab, integer *ipiv, integer &info )
489 {
490 return( LPFNAME( zgbtrf ) ( &m, &n, &kl, &ku, reinterpret_cast< doublecomplex* >( ab ), &ldab, ipiv, &info ) );
491 }
492 // 一般三重対角行列の連立方程式を解く関数
493 inline int gttrf( integer &n, real *dl, real *d__, real *du, real *du2, integer *ipiv, integer &info )
494 {
495 return( LPFNAME( sgttrf ) ( &n, dl, d__, du, du2, ipiv, &info ) );
496 }
497 inline int gttrf( integer &n, doublereal *dl, doublereal *d__, doublereal *du, doublereal *du2, integer *ipiv, integer &info )
498 {
499 return( LPFNAME( dgttrf ) ( &n, dl, d__, du, du2, ipiv, &info ) );
500 }
501 inline int gttrf( integer &n, std::complex< real > *dl, std::complex< real > *d__, std::complex< real > *du, std::complex< real > *du2, integer *ipiv, integer &info )
502 {
503 return( LPFNAME( cgttrf ) ( &n, reinterpret_cast< complex* >( dl ), reinterpret_cast< complex* >( d__ ), reinterpret_cast< complex* >( du ), reinterpret_cast< complex* >( du2 ), ipiv, &info ) );
504 }
505 inline int gttrf( integer &n, std::complex< doublereal > *dl, std::complex< doublereal > *d__, std::complex< doublereal > *du, std::complex< doublereal > *du2, integer *ipiv, integer &info )
506 {
507 return( LPFNAME( zgttrf ) ( &n, reinterpret_cast< doublecomplex* >( dl ), reinterpret_cast< doublecomplex* >( d__ ), reinterpret_cast< doublecomplex* >( du ), reinterpret_cast< doublecomplex* >( du2 ), ipiv, &info ) );
508 }
509 // 対称行列のLU分解
510 inline int sytrf( char *uplo, integer &n, real *a, integer &lda, integer *ipiv, real *work, integer &lwork, integer &info )
511 {
512 return( LPFNAME( ssytrf ) ( uplo, &n, a, &lda, ipiv, work, &lwork, &info ) );
513 }
514 inline int sytrf( char *uplo, integer &n, doublereal *a, integer &lda, integer *ipiv, doublereal *work, integer &lwork, integer &info )
515 {
516 return( LPFNAME( dsytrf ) ( uplo, &n, a, &lda, ipiv, work, &lwork, &info ) );
517 }
518 inline int sytrf( char *uplo, integer &n, std::complex< real > *a, integer &lda, integer *ipiv, std::complex< real > *work, integer &lwork, integer &info )
519 {
520 return( LPFNAME( csytrf ) ( uplo, &n, reinterpret_cast< complex* >( a ), &lda, ipiv, reinterpret_cast< complex* >( work ), &lwork, &info ) );
521 }
522 inline int sytrf( char *uplo, integer &n, std::complex< doublereal > *a, integer &lda, integer *ipiv, std::complex< doublereal > *work, integer &lwork, integer &info )
523 {
524 return( LPFNAME( zsytrf ) ( uplo, &n, reinterpret_cast< doublecomplex* >( a ), &lda, ipiv, reinterpret_cast< doublecomplex* >( work ), &lwork, &info ) );
525 }
526 // エルミート行列のLU分解
527 inline int hetrf( char *uplo, integer &n, std::complex< real > *a, integer &lda, integer *ipiv, std::complex< real > *work, integer &lwork, integer &info )
528 {
529 return( LPFNAME( chetrf ) ( uplo, &n, reinterpret_cast< complex* >( a ), &lda, ipiv, reinterpret_cast< complex* >( work ), &lwork, &info ) );
530 }
531 inline int hetrf( char *uplo, integer &n, std::complex< doublereal > *a, integer &lda, integer *ipiv, std::complex< doublereal > *work, integer &lwork, integer &info )
532 {
533 return( LPFNAME( zhetrf ) ( uplo, &n, reinterpret_cast< doublecomplex* >( a ), &lda, ipiv, reinterpret_cast< doublecomplex* >( work ), &lwork, &info ) );
534 }
535
536 // 対称行列のコレスキー分解
537 inline int potrf( char *uplo, integer &n, real *a, integer &lda, integer &info )
538 {
539 return( LPFNAME( spotrf ) ( uplo, &n, a, &lda, &info ) );
540 }
541 inline int potrf( char *uplo, integer &n, doublereal *a, integer &lda, integer &info )
542 {
543 return( LPFNAME( dpotrf ) ( uplo, &n, a, &lda, &info ) );
544 }
545 inline int potrf( char *uplo, integer &n, std::complex< real > *a, integer &lda, integer &info )
546 {
547 return( LPFNAME( cpotrf ) ( uplo, &n, reinterpret_cast< complex* >( a ), &lda, &info ) );
548 }
549 inline int potrf( char *uplo, integer &n, std::complex< doublereal > *a, integer &lda, integer &info )
550 {
551 return( LPFNAME( zpotrf ) ( uplo, &n, reinterpret_cast< doublecomplex* >( a ), &lda, &info ) );
552 }
553
554 // QR分解
555 inline int geqrf( integer &m, integer &n, real *a, integer &lda, real *tau, real *work, integer &lwork, integer &info )
556 {
557 return( LPFNAME( sgeqrf ) ( &m, &n, a, &lda, tau, work, &lwork, &info ) );
558 }
559 inline int geqrf( integer &m, integer &n, doublereal *a, integer &lda, doublereal *tau, doublereal *work, integer &lwork, integer &info )
560 {
561 return( LPFNAME( dgeqrf ) ( &m, &n, a, &lda, tau, work, &lwork, &info ) );
562 }
563 inline int geqrf( integer &m, integer &n, std::complex< real > *a, integer &lda, std::complex< real > *tau, std::complex< real > *work, integer &lwork, integer &info )
564 {
565 return( LPFNAME( cgeqrf ) ( &m, &n, reinterpret_cast< complex* >( a ), &lda, reinterpret_cast< complex* >( tau ), reinterpret_cast< complex* >( work ), &lwork, &info ) );
566 }
567 inline int geqrf( integer &m, integer &n, std::complex< doublereal > *a, integer &lda, std::complex< doublereal > *tau, std::complex< doublereal > *work, integer &lwork, integer &info )
568 {
569 return( LPFNAME( zgeqrf ) ( &m, &n, reinterpret_cast< doublecomplex* >( a ), &lda, reinterpret_cast< doublecomplex* >( tau ), reinterpret_cast< doublecomplex* >( work ), &lwork, &info ) );
570 }
571
572 // QR分解(Level3 Blas)
573 inline int geqp3( integer &m, integer &n, real *a, integer &lda, integer *jpvt, real *tau, real *work, integer &lwork, integer &info )
574 {
575 return( LPFNAME( sgeqp3 ) ( &m, &n, a, &lda, jpvt, tau, work, &lwork, &info ) );
576 }
577 inline int geqp3( integer &m, integer &n, doublereal *a, integer &lda, integer *jpvt, doublereal *tau, doublereal *work, integer &lwork, integer &info )
578 {
579 return( LPFNAME( dgeqp3 ) ( &m, &n, a, &lda, jpvt, tau, work, &lwork, &info ) );
580 }
581 inline int geqp3( integer &m, integer &n, std::complex< real > *a, integer &lda, integer *jpvt, std::complex< real > *tau, std::complex< real > *work, integer &lwork, real *rwork, integer &info )
582 {
583 return( LPFNAME( cgeqp3 ) ( &m, &n, reinterpret_cast< complex* >( a ), &lda, jpvt, reinterpret_cast< complex* >( tau ), reinterpret_cast< complex* >( work ), &lwork, rwork, &info ) );
584 }
585 inline int geqp3( integer &m, integer &n, std::complex< doublereal > *a, integer &lda, integer *jpvt, std::complex< doublereal > *tau, std::complex< doublereal > *work, integer &lwork, doublereal *rwork, integer &info )
586 {
587 return( LPFNAME( zgeqp3 ) ( &m, &n, reinterpret_cast< doublecomplex* >( a ), &lda, jpvt, reinterpret_cast< doublecomplex* >( tau ), reinterpret_cast< doublecomplex* >( work ), &lwork, rwork, &info ) );
588 }
589
590 // 一般行列のQR分解後のQ行列を作成するルーチン
591 inline int orgqr( integer &m, integer &n, integer &k, real *a, integer &lda, real *tau, real *work, integer &lwork, integer &info )
592 {
593 return( LPFNAME( sorgqr ) ( &m, &n, &k, a, &lda, tau, work, &lwork, &info ) );
594 }
595 inline int orgqr( integer &m, integer &n, integer &k, doublereal *a, integer &lda, doublereal *tau, doublereal *work, integer &lwork, integer &info )
596 {
597 return( LPFNAME( dorgqr ) ( &m, &n, &k, a, &lda, tau, work, &lwork, &info ) );
598 }
599 inline int ungqr( integer &m, integer &n, integer &k, std::complex< real > *a, integer &lda, std::complex< real > *tau, std::complex< real > *work, integer &lwork, integer &info )
600 {
601 return( LPFNAME( cungqr ) ( &m, &n, &k, reinterpret_cast< complex* >( a ), &lda, reinterpret_cast< complex* >( tau ), reinterpret_cast< complex* >( work ), &lwork, &info ) );
602 }
603 inline int ungqr( integer &m, integer &n, integer &k, std::complex< doublereal > *a, integer &lda, std::complex< doublereal > *tau, std::complex< doublereal > *work, integer &lwork, integer &info )
604 {
605 return( LPFNAME( zungqr ) ( &m, &n, &k, reinterpret_cast< doublecomplex* >( a ), &lda, reinterpret_cast< doublecomplex* >( tau ), reinterpret_cast< doublecomplex* >( work ), &lwork, &info ) );
606 }
607
608 // LU分解の結果を用いた一般行列の逆行列の計算
609 inline int getri( integer &n, real *a, integer &lda, integer *ipiv, real *work, integer &lwork, integer &info )
610 {
611 return( LPFNAME( sgetri ) ( &n, a, &lda, ipiv, work, &lwork, &info ) );
612 }
613 inline int getri( integer &n, doublereal *a, integer &lda, integer *ipiv, doublereal *work, integer &lwork, integer &info )
614 {
615 return( LPFNAME( dgetri ) ( &n, a, &lda, ipiv, work, &lwork, &info ) );
616 }
617 inline int getri( integer &n, std::complex< real > *a, integer &lda, integer *ipiv, std::complex< real > *work, integer &lwork, integer &info )
618 {
619 return( LPFNAME( cgetri ) ( &n, reinterpret_cast< complex* >( a ), &lda, ipiv, reinterpret_cast< complex* >( work ), &lwork, &info ) );
620 }
621 inline int getri( integer &n, std::complex< doublereal > *a, integer &lda, integer *ipiv, std::complex< doublereal > *work, integer &lwork, integer &info )
622 {
623 return( LPFNAME( zgetri ) ( &n, reinterpret_cast< doublecomplex* >( a ), &lda, ipiv, reinterpret_cast< doublecomplex* >( work ), &lwork, &info ) );
624 }
625 // LU分解の結果を用いた対称行列の逆行列の計算
626 inline int sytri( char *uplo, integer &n, real *a, integer &lda, integer *ipiv, real *work, integer &info )
627 {
628 return( LPFNAME( ssytri ) ( uplo, &n, a, &lda, ipiv, work, &info ) );
629 }
630 inline int sytri( char *uplo, integer &n, doublereal *a, integer &lda, integer *ipiv, doublereal *work, integer &info )
631 {
632 return( LPFNAME( dsytri ) ( uplo, &n, a, &lda, ipiv, work, &info ) );
633 }
634 inline int sytri( char *uplo, integer &n, std::complex< real > *a, integer &lda, integer *ipiv, std::complex< real > *work, integer &info )
635 {
636 return( LPFNAME( csytri ) ( uplo, &n, reinterpret_cast< complex* >( a ), &lda, ipiv, reinterpret_cast< complex* >( work ), &info ) );
637 }
638 inline int sytri( char *uplo, integer &n, std::complex< doublereal > *a, integer &lda, integer *ipiv, std::complex< doublereal > *work, integer &info )
639 {
640 return( LPFNAME( zsytri ) ( uplo, &n, reinterpret_cast< doublecomplex* >( a ), &lda, ipiv, reinterpret_cast< doublecomplex* >( work ), &info ) );
641 }
642 // LU分解の結果を用いたエルミート行列の逆行列の計算
643 inline int hetri( char *uplo, integer &n, std::complex< real > *a, integer &lda, integer *ipiv, std::complex< real > *work, integer &info )
644 {
645 return( LPFNAME( chetri ) ( uplo, &n, reinterpret_cast< complex* >( a ), &lda, ipiv, reinterpret_cast< complex* >( work ), &info ) );
646 }
647 inline int hetri( char *uplo, integer &n, std::complex< doublereal > *a, integer &lda, integer *ipiv, std::complex< doublereal > *work, integer &info )
648 {
649 return( LPFNAME( zhetri ) ( uplo, &n, reinterpret_cast< doublecomplex* >( a ), &lda, ipiv, reinterpret_cast< doublecomplex* >( work ), &info ) );
650 }
651
652
653 // 一般行列に対する固有値・固有ベクトルを計算
654 inline int geev( char *jobvl, char *jobvr, integer &n, real *a, integer &lda, real *wr, real *wi,
655 real *vl, integer &ldvl, real *vr, integer &ldvr, real *work, integer &lwork, integer &info )
656 {
657 return( LPFNAME( sgeev ) ( jobvl, jobvr, &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr, work, &lwork, &info ) );
658 }
659 inline int geev( char *jobvl, char *jobvr, integer &n, doublereal *a, integer &lda, doublereal *wr, doublereal *wi,
660 doublereal *vl, integer &ldvl, doublereal *vr, integer &ldvr, doublereal *work, integer &lwork, integer &info )
661 {
662 return( LPFNAME( dgeev ) ( jobvl, jobvr, &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr, work, &lwork, &info ) );
663 }
664 inline int geev( char *jobvl, char *jobvr, integer &n, std::complex< real > *a, integer &lda, std::complex< real > *w, std::complex< real > *vl, integer &ldvl,
665 std::complex< real > *vr, integer &ldvr, std::complex< real > *work, integer &lwork, real *rwork, integer &info )
666 {
667 return( LPFNAME( cgeev ) ( jobvl, jobvr, &n, reinterpret_cast< complex* >( a ), &lda, reinterpret_cast< complex* >( w ), reinterpret_cast< complex* >( vl ), &ldvl,
668 reinterpret_cast< complex* >( vr ), &ldvr, reinterpret_cast< complex* >( work ), &lwork, rwork, &info ) );
669 }
670 inline int geev( char *jobvl, char *jobvr, integer &n, std::complex< doublereal > *a, integer &lda, std::complex< doublereal > *w, std::complex< doublereal > *vl, integer &ldvl,
671 std::complex< doublereal > *vr, integer &ldvr, std::complex< doublereal > *work, integer &lwork, doublereal *rwork, integer &info )
672 {
673 return( LPFNAME( zgeev ) ( jobvl, jobvr, &n, reinterpret_cast< doublecomplex* >( a ), &lda, reinterpret_cast< doublecomplex* >( w ), reinterpret_cast< doublecomplex* >( vl ), &ldvl,
674 reinterpret_cast< doublecomplex* >( vr ), &ldvr, reinterpret_cast< doublecomplex* >( work ), &lwork, rwork, &info ) );
675 }
676 // 対称行列に対する固有値・固有ベクトルを計算
677 inline int syev( char *jobz, char *uplo, integer &n, real *a, integer &lda, real *w, real *work, integer &lwork, integer &info )
678 {
679 return( LPFNAME( ssyev ) ( jobz, uplo, &n, a, &lda, w, work, &lwork, &info ) );
680 }
681 inline int syev( char *jobz, char *uplo, integer &n, doublereal *a, integer &lda, doublereal *w, doublereal *work, integer &lwork, integer &info )
682 {
683 return( LPFNAME( dsyev ) ( jobz, uplo, &n, a, &lda, w, work, &lwork, &info ) );
684 }
685 // 対称帯行列に対する固有値・固有ベクトルを計算
686 inline int sbev( char *jobz, char *uplo, integer &n, integer &kd, real *ab, integer &ldab, real *w, real *z__, integer &ldz, real *work, integer &info )
687 {
688 return( LPFNAME( ssbev ) ( jobz, uplo, &n, &kd, ab, &ldab, w, z__, &ldz, work, &info ) );
689 }
690 inline int sbev( char *jobz, char *uplo, integer &n, integer &kd, doublereal *ab, integer &ldab, doublereal *w, doublereal *z__, integer &ldz, doublereal *work, integer &info )
691 {
692 return( LPFNAME( dsbev ) ( jobz, uplo, &n, &kd, ab, &ldab, w, z__, &ldz, work, &info ) );
693 }
694 // 対称3重対角行列に対する固有値・固有ベクトルを計算
695 inline int stev( char *jobz, integer &n, real *d__, real *e, real *z__, integer &ldz, real *work, integer &info )
696 {
697 return( LPFNAME( sstev ) ( jobz, &n, d__, e, z__, &ldz, work, &info ) );
698 }
699 inline int stev( char *jobz, integer &n, doublereal *d__, doublereal *e, doublereal *z__, integer &ldz, doublereal *work, integer &info )
700 {
701 return( LPFNAME( dstev ) ( jobz, &n, d__, e, z__, &ldz, work, &info ) );
702 }
703 // エルミート行列に対する固有値・固有ベクトルを計算
704 inline int heev( char *jobz, char *uplo, integer &n, std::complex< real > *a, integer &lda, real *w, std::complex< real > *work, integer &lwork, real *rwork, integer &info )
705 {
706 return( LPFNAME( cheev ) ( jobz, uplo, &n, reinterpret_cast< complex* >( a ), &lda, w, reinterpret_cast< complex* >( work ), &lwork, rwork, &info ) );
707 }
708 inline int heev( char *jobz, char *uplo, integer &n, std::complex< doublereal > *a, integer &lda, doublereal *w, std::complex< doublereal > *work, integer &lwork, doublereal *rwork, integer &info )
709 {
710 return( LPFNAME( zheev ) ( jobz, uplo, &n, reinterpret_cast< doublecomplex* >( a ), &lda, w, reinterpret_cast< doublecomplex* >( work ), &lwork, rwork, &info ) );
711 }
712 // エルミート帯行列に対する固有値・固有ベクトルを計算
713 inline int hbev( char *jobz, char *uplo, integer &n, integer &kd, std::complex< real > *ab, integer &ldab, real *w, std::complex< real > *z__, integer &ldz, std::complex< real > *work, real *rwork, integer &info )
714 {
715 return( LPFNAME( chbev ) ( jobz, uplo, &n, &kd, reinterpret_cast< complex* >( ab ), &ldab, w, reinterpret_cast< complex* >( z__ ), &ldz, reinterpret_cast< complex* >( work ), rwork, &info ) );
716 }
717 inline int hbev( char *jobz, char *uplo, integer &n, integer &kd, std::complex< doublereal > *ab, integer &ldab, doublereal *w, std::complex< doublereal > *z__, integer &ldz, std::complex< doublereal > *work, doublereal *rwork, integer &info )
718 {
719 return( LPFNAME( zhbev ) ( jobz, uplo, &n, &kd, reinterpret_cast< doublecomplex* >( ab ), &ldab, w, reinterpret_cast< doublecomplex* >( z__ ), &ldz, reinterpret_cast< doublecomplex* >( work ), rwork, &info ) );
720 }
721
722
723
724 // 一般行列に対する固有値・固有ベクトルを計算．行列を対角化して精度を上げるバージョン
725 inline int geevx( char *balanc, char *jobvl, char *jobvr, char *sense, integer &n, real *a, integer &lda, real *wr, real *wi,
726 real *vl, integer &ldvl, real *vr, integer &ldvr, integer &ilo, integer &ihi, real *scale,
727 real &abnrm, real *rconde, real *rcondv, real *work, integer &lwork, integer *iwork, integer &info )
728 {
729 return( LPFNAME( sgeevx ) ( balanc, jobvl, jobvr, sense, &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr, &ilo, &ihi, scale, &abnrm, rconde, rcondv, work, &lwork, iwork, &info ) );
730 }
731 inline int geevx( char *balanc, char *jobvl, char *jobvr, char *sense, integer &n, doublereal *a, integer &lda, doublereal *wr, doublereal *wi,
732 doublereal *vl, integer &ldvl, doublereal *vr, integer &ldvr, integer &ilo, integer &ihi, doublereal *scale,
733 doublereal &abnrm, doublereal *rconde, doublereal *rcondv, doublereal *work, integer &lwork, integer *iwork, integer &info )
734 {
735 return( LPFNAME( dgeevx ) ( balanc, jobvl, jobvr, sense, &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr, &ilo, &ihi, scale, &abnrm, rconde, rcondv, work, &lwork, iwork, &info ) );
736 }
737 inline int geevx( char *balanc, char *jobvl, char *jobvr, char *sense, integer &n, std::complex< real > *a, integer &lda, std::complex< real > *w,
738 std::complex< real > *vl, integer &ldvl, std::complex< real > *vr, integer &ldvr, integer &ilo, integer &ihi, real *scale,
739 real &abnrm, real *rconde, real *rcondv, std::complex< real > *work, integer &lwork, real *rwork, integer &info )
740 {
741 return( LPFNAME( cgeevx ) ( balanc, jobvl, jobvr, sense, &n, reinterpret_cast< complex* >( a ), &lda, reinterpret_cast< complex* >( w ), reinterpret_cast< complex* >( vl ), &ldvl,
742 reinterpret_cast< complex* >( vr ), &ldvr, &ilo, &ihi, scale, &abnrm, rconde, rcondv, reinterpret_cast< complex* >( work ), &lwork, rwork, &info ) );
743 }
744 inline int geevx( char *balanc, char *jobvl, char *jobvr, char *sense, integer &n, std::complex< doublereal > *a, integer &lda, std::complex< doublereal > *w,
745 std::complex< doublereal > *vl, integer &ldvl, std::complex< doublereal > *vr, integer &ldvr, integer &ilo, integer &ihi, doublereal *scale,
746 doublereal &abnrm, doublereal *rconde, doublereal *rcondv, std::complex< doublereal > *work, integer &lwork, doublereal *rwork, integer &info )
747 {
748 return( LPFNAME( zgeevx ) ( balanc, jobvl, jobvr, sense, &n, reinterpret_cast< doublecomplex* >( a ), &lda, reinterpret_cast< doublecomplex* >( w ), reinterpret_cast< doublecomplex* >( vl ), &ldvl,
749 reinterpret_cast< doublecomplex* >( vr ), &ldvr, &ilo, &ihi, scale, &abnrm, rconde, rcondv, reinterpret_cast< doublecomplex* >( work ), &lwork, rwork, &info ) );
750 }
751 // 対称行列に対する固有値・固有ベクトルを計算
752 inline int syevx( char *jobz, char *range, char *uplo, integer &n, real *a, integer &lda,
753 real &vl, real &vu, integer &il, integer &iu, real &abstol, integer &m, real *w, real *z__,
754 integer &ldz, real *work, integer &lwork, integer *iwork, integer *ifail, integer &info )
755 {
756 return( LPFNAME( ssyevx ) ( jobz, range, uplo, &n, a, &lda, &vl, &vu, &il, &iu, &abstol, &m, w, z__, &ldz, work, &lwork, iwork, ifail, &info ) );
757 }
758 inline int syevx( char *jobz, char *range, char *uplo, integer &n, doublereal *a, integer &lda,
759 doublereal &vl, doublereal &vu, integer &il, integer &iu, doublereal &abstol, integer &m, doublereal *w, doublereal *z__,
760 integer &ldz, doublereal *work, integer &lwork, integer *iwork, integer *ifail, integer &info )
761 {
762 return( LPFNAME( dsyevx ) ( jobz, range, uplo, &n, a, &lda, &vl, &vu, &il, &iu, &abstol, &m, w, z__, &ldz, work, &lwork, iwork, ifail, &info ) );
763 }
764 // 対称帯行列に対する固有値・固有ベクトルを計算．行列を対角化して精度を上げるバージョン
765 inline int sbevx( char *jobz, char *range, char *uplo, integer &n, integer &kd, real *ab, integer &ldab, real *q, integer &ldq, real &vl, real &vu, integer &il, integer &iu, real &abstol, integer &m, real *w, real *z__, integer &ldz, real *work, integer *iwork, integer *ifail, integer &info )
766 {
767 return( LPFNAME( ssbevx ) ( jobz, range, uplo, &n, &kd, ab, &ldab, q, &ldq, &vl, &vu, &il, &iu, &abstol, &m, w, z__, &ldz, work, iwork, ifail, &info ) );
768 }
769 inline int dsbevx( char *jobz, char *range, char *uplo, integer &n, integer &kd, doublereal *ab, integer &ldab, doublereal *q, integer &ldq, doublereal &vl, doublereal &vu, integer &il, integer &iu, doublereal &abstol, integer &m, doublereal *w, doublereal *z__, integer &ldz, doublereal *work, integer *iwork, integer *ifail, integer &info )
770 {
771 return( LPFNAME( dsbevx ) ( jobz, range, uplo, &n, &kd, ab, &ldab, q, &ldq, &vl, &vu, &il, &iu, &abstol, &m, w, z__, &ldz, work, iwork, ifail, &info ) );
772 }
773 // 対称3重対角行列に対する固有値・固有ベクトルを計算．行列を対角化して精度を上げるバージョン
774 inline int stevx( char *jobz, char *range, integer &n, real *d__, real *e, real &vl, real &vu, integer &il, integer &iu, real &abstol, integer &m, real *w, real *z__, integer &ldz, real *work, integer *iwork, integer *ifail, integer &info )
775 {
776 return( LPFNAME( sstevx ) ( jobz, range, &n, d__, e, &vl, &vu, &il, &iu, &abstol, &m, w, z__, &ldz, work, iwork, ifail, &info ) );
777 }
778 inline int stevx( char *jobz, char *range, integer &n, doublereal *d__, doublereal *e, doublereal &vl, doublereal &vu, integer &il, integer &iu, doublereal &abstol, integer &m, doublereal *w, doublereal *z__, integer &ldz, doublereal *work, integer *iwork, integer *ifail, integer &info )
779 {
780 return( LPFNAME( dstevx ) ( jobz, range, &n, d__, e, &vl, &vu, &il, &iu, &abstol, &m, w, z__, &ldz, work, iwork, ifail, &info ) );
781 }
782 // エルミート行列に対する固有値・固有ベクトルを計算．行列を対角化して精度を上げるバージョン
783 inline int heevx( char *jobz, char *range, char *uplo, integer &n, std::complex< real > *a, integer &lda, real &vl, real &vu, integer &il, integer &iu, real &abstol, integer &m, real *w, std::complex< real > *z__, integer &ldz, std::complex< real > *work, integer &lwork, real *rwork, integer *iwork, integer *ifail, integer &info )
784 {
785 return( LPFNAME( cheevx ) ( jobz, range, uplo, &n, reinterpret_cast< complex* >( a ), &lda, &vl, &vu, &il, &iu, &abstol, &m, w, reinterpret_cast< complex* >( z__ ), &ldz, reinterpret_cast< complex* >( work ), &lwork, rwork, iwork, ifail, &info ) );
786 }
787 inline int heevx( char *jobz, char *range, char *uplo, integer &n, std::complex< doublereal > *a, integer &lda, doublereal &vl, doublereal &vu, integer &il, integer &iu, doublereal &abstol, integer &m, doublereal *w, std::complex< doublereal > *z__, integer &ldz, std::complex< doublereal > *work, integer &lwork, doublereal *rwork, integer *iwork, integer *ifail, integer &info )
788 {
789 return( LPFNAME( zheevx ) ( jobz, range, uplo, &n, reinterpret_cast< doublecomplex* >( a ), &lda, &vl, &vu, &il, &iu, &abstol, &m, w, reinterpret_cast< doublecomplex* >( z__ ), &ldz, reinterpret_cast< doublecomplex* >( work ), &lwork, rwork, iwork, ifail, &info ) );
790 }
791 // エルミート帯行列に対する固有値・固有ベクトルを計算．行列を対角化して精度を上げるバージョン
792 inline int hbevx( char *jobz, char *range, char *uplo, integer &n, integer &kd, complex *ab, integer &ldab, complex *q, integer &ldq, real &vl, real &vu, integer &il, integer &iu, real &abstol, integer &m, real *w, complex *z__, integer &ldz, complex *work, real *rwork, integer *iwork, integer *ifail, integer &info )
793 {
794 return( LPFNAME( chbevx ) ( jobz, range, uplo, &n, &kd, reinterpret_cast< complex* >( ab ), &ldab, reinterpret_cast< complex* >( q ), &ldq, &vl, &vu, &il, &iu, &abstol, &m, w, reinterpret_cast< complex* >( z__ ), &ldz, reinterpret_cast< complex* >( work ), rwork, iwork, ifail, &info ) );
795 }
796 inline int hbevx( char *jobz, char *range, char *uplo, integer &n, integer &kd, complex *ab, integer &ldab, complex *q, integer &ldq, doublereal &vl, doublereal &vu, integer &il, integer &iu, doublereal &abstol, integer &m, doublereal *w, complex *z__, integer &ldz, complex *work, doublereal *rwork, integer *iwork, integer *ifail, integer &info )
797 {
798 return( LPFNAME( zhbevx ) ( jobz, range, uplo, &n, &kd, reinterpret_cast< doublecomplex* >( ab ), &ldab, reinterpret_cast< doublecomplex* >( q ), &ldq, &vl, &vu, &il, &iu, &abstol, &m, w, reinterpret_cast< doublecomplex* >( z__ ), &ldz, reinterpret_cast< doublecomplex* >( work ), rwork, iwork, ifail, &info ) );
799 }
800
801
802
803
804 // 一般行列に対する特異値分解を計算
805 inline int gesvd( char *jobu, char *jobvt, integer &m, integer &n, real *a, integer &lda, real *s, real *u,
806 integer &ldu, real *vt, integer &ldvt, real *work, integer &lwork, integer &info )
807 {
808 return( LPFNAME( sgesvd ) ( jobu, jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, &info ) );
809 }
810 inline int gesvd( char *jobu, char *jobvt, integer &m, integer &n, doublereal *a, integer &lda, doublereal *s, doublereal *u,
811 integer &ldu, doublereal *vt, integer &ldvt, doublereal *work, integer &lwork, integer &info )
812 {
813 return( LPFNAME( dgesvd ) ( jobu, jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, &info ) );
814 }
815 inline int gesvd( char *jobu, char *jobvt, integer &m, integer &n, std::complex< real > *a, integer &lda, real *s, std::complex< real > *u,
816 integer &ldu, std::complex< real > *vt, integer &ldvt, std::complex< real > *work, integer &lwork, real *rwork, integer &info )
817 {
818 return( LPFNAME( cgesvd ) ( jobu, jobvt, &m, &n, reinterpret_cast< complex* >( a ), &lda, s, reinterpret_cast< complex* >( u ),
819 &ldu, reinterpret_cast< complex* >( vt ), &ldvt, reinterpret_cast< complex* >( work ), &lwork, rwork, &info ) );
820 }
821 inline int gesvd( char *jobu, char *jobvt, integer &m, integer &n, std::complex< doublereal > *a, integer &lda, doublereal *s, std::complex< doublereal > *u,
822 integer &ldu, std::complex< doublereal > *vt, integer &ldvt, std::complex< doublereal > *work, integer &lwork, doublereal *rwork, integer &info )
823 {
824 return( LPFNAME( zgesvd ) ( jobu, jobvt, &m, &n, reinterpret_cast< doublecomplex* >( a ), &lda, s, reinterpret_cast< doublecomplex* >( u ),
825 &ldu, reinterpret_cast< doublecomplex* >( vt ), &ldvt, reinterpret_cast< doublecomplex* >( work ), &lwork, rwork, &info ) );
826 }
827
828 // 一般行列に対する特異値分解を計算．分割統治法を用いた高速バージョン
829 inline int gesdd( char *jobz, integer &m, integer &n, real *a, integer &lda, real *s, real *u, integer &ldu,
830 real *vt, integer &ldvt, real *work, integer &lwork, integer *iwork, integer &info )
831 {
832 return( LPFNAME( sgesdd ) ( jobz, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, iwork, &info ) );
833 }
834 inline int gesdd( char *jobz, integer &m, integer &n, doublereal *a, integer &lda, doublereal *s, doublereal *u, integer &ldu,
835 doublereal *vt, integer &ldvt, doublereal *work, integer &lwork, integer *iwork, integer &info )
836 {
837 return( LPFNAME( dgesdd ) ( jobz, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, iwork, &info ) );
838 }
839 inline int gesdd( char *jobz, integer &m, integer &n, std::complex< real > *a, integer &lda, real *s, std::complex< real > *u, integer &ldu,
840 std::complex< real > *vt, integer &ldvt, std::complex< real > *work, integer &lwork, real *rwork, integer *iwork, integer &info )
841 {
842 return( LPFNAME( cgesdd ) ( jobz, &m, &n, reinterpret_cast< complex* >( a ), &lda, s, reinterpret_cast< complex* >( u ), &ldu,
843 reinterpret_cast< complex* >( vt ), &ldvt, reinterpret_cast< complex* >( work ), &lwork, rwork, iwork, &info ) );
844 }
845 inline int gesdd( char *jobz, integer &m, integer &n, std::complex< doublereal > *a, integer &lda, doublereal *s, std::complex< doublereal > *u, integer &ldu,
846 std::complex< doublereal > *vt, integer &ldvt, std::complex< doublereal > *work, integer &lwork, doublereal *rwork, integer *iwork, integer &info )
847 {
848 return( LPFNAME( zgesdd ) ( jobz, &m, &n, reinterpret_cast< doublecomplex* >( a ), &lda, s, reinterpret_cast< doublecomplex* >( u ), &ldu,
849 reinterpret_cast< doublecomplex* >( vt ), &ldvt, reinterpret_cast< doublecomplex* >( work ), &lwork, rwork, iwork, &info ) );
850 }
851
852
853 // 対称行列に対する固有値・固有ベクトルを計算を計算．分割統治法を用いた高速バージョン
854 inline int syevd( char *jobz, char *uplo, integer &n, real *a, integer &lda, real *w, real *work, integer &lwork, integer *iwork, integer &liwork, integer &info )
855 {
856 return( LPFNAME( ssyevd ) ( jobz, uplo, &n, a, &lda, w, work, &lwork, iwork, &liwork, &info ) );
857 }
858 inline int syevd( char *jobz, char *uplo, integer &n, doublereal *a, integer &lda, doublereal *w, doublereal *work, integer &lwork, integer *iwork, integer &liwork, integer &info )
859 {
860 return( LPFNAME( dsyevd ) ( jobz, uplo, &n, a, &lda, w, work, &lwork, iwork, &liwork, &info ) );
861 }
862 // 対称帯行列に対する固有値・固有ベクトルを計算を計算．分割統治法を用いた高速バージョン
863 inline int sbevd( char *jobz, char *uplo, integer &n, integer &kd, real *ab, integer &ldab, real *w, real *z__, integer &ldz, real *work, integer &lwork, integer *iwork, integer &liwork, integer &info )
864 {
865 return( LPFNAME( ssbevd ) ( jobz, uplo, &n, &kd, ab, &ldab, w, z__, &ldz, work, &lwork, iwork, &liwork, &info ) );
866 }
867 inline int sbevd( char *jobz, char *uplo, integer &n, integer &kd, doublereal *ab, integer &ldab, doublereal *w, doublereal *z__, integer &ldz, doublereal *work, integer &lwork, integer *iwork, integer &liwork, integer &info )
868 {
869 return( LPFNAME( dsbevd ) ( jobz, uplo, &n, &kd, ab, &ldab, w, z__, &ldz, work, &lwork, iwork, &liwork, &info ) );
870 }
871 // 対称3重対角行列に対する固有値・固有ベクトルを計算を計算．分割統治法を用いた高速バージョン
872 inline int stevd( char *jobz, integer &n, real *d__, real *e, real *z__, integer &ldz, real *work, integer &lwork, integer *iwork, integer &liwork, integer &info )
873 {
874 return( LPFNAME( sstevd ) ( jobz, &n, d__, e, z__, &ldz, work, &lwork, iwork, &liwork, &info ) );
875 }
876 inline int stevd( char *jobz, integer &n, doublereal *d__, doublereal *e, doublereal *z__, integer &ldz, doublereal *work, integer &lwork, integer *iwork, integer &liwork, integer &info )
877 {
878 return( LPFNAME( dstevd ) ( jobz, &n, d__, e, z__, &ldz, work, &lwork, iwork, &liwork, &info ) );
879 }
880 // エルミート行列に対する固有値・固有ベクトルを計算を計算．分割統治法を用いた高速バージョン
881 inline int heevd( char *jobz, char *uplo, integer &n, std::complex< real > *a, integer &lda, real *w, std::complex< real > *work, integer &lwork, real *rwork, integer &lrwork, integer *iwork, integer &liwork, integer &info )
882 {
883 return( LPFNAME( cheevd ) ( jobz, uplo, &n, reinterpret_cast< complex* >( a ), &lda, w, reinterpret_cast< complex* >( work ), &lwork, rwork, &lrwork, iwork, &liwork, &info ) );
884 }
885 inline int heevd( char *jobz, char *uplo, integer &n, std::complex< doublereal > *a, integer &lda, doublereal *w, std::complex< doublereal > *work, integer &lwork, doublereal *rwork, integer &lrwork, integer *iwork, integer &liwork, integer &info )
886 {
887 return( LPFNAME( zheevd ) ( jobz, uplo, &n, reinterpret_cast< doublecomplex* >( a ), &lda, w, reinterpret_cast< doublecomplex* >( work ), &lwork, rwork, &lrwork, iwork, &liwork, &info ) );
888 }
889 // エルミート帯行列に対する固有値・固有ベクトルを計算を計算．分割統治法を用いた高速バージョン
890 inline int hbevd( char *jobz, char *uplo, integer &n, integer &kd, std::complex< real > *ab, integer &ldab, real *w, std::complex< real > *z__, integer &ldz, std::complex< real > *work, integer &lwork, real *rwork, integer &lrwork, integer *iwork, integer &liwork, integer &info )
891 {
892 return( LPFNAME( chbevd ) ( jobz, uplo, &n, &kd, reinterpret_cast< complex* >( ab ), &ldab, w, reinterpret_cast< complex* >( z__ ), &ldz, reinterpret_cast< complex* >( work ), &lwork, rwork, &lrwork, iwork, &liwork, &info ) );
893 }
894 inline int hbevd( char *jobz, char *uplo, integer &n, integer &kd, std::complex< doublereal > *ab, integer &ldab, doublereal *w, std::complex< doublereal > *z__, integer &ldz, std::complex< doublereal > *work, integer &lwork, doublereal *rwork, integer &lrwork, integer *iwork, integer &liwork, integer &info )
895 {
896 return( LPFNAME( zhbevd ) ( jobz, uplo, &n, &kd, reinterpret_cast< doublecomplex* >( ab ), &ldab, w, reinterpret_cast< doublecomplex* >( z__ ), &ldz, reinterpret_cast< doublecomplex* >( work ), &lwork, rwork, &lrwork, iwork, &liwork, &info ) );
897 }
898
899
900 // インテルのMKLとの互換性を保つための，関数名の変換マクロ
901 #undef LPFNAME
902
903
904 template < class T >
905 class __value_index_pair__
906 {
907 public:
908 typedef size_t size_type;
909 typedef T value_type;
910
911 value_type value;
912 size_type index;
913
914 __value_index_pair__( ){ }
915
916 __value_index_pair__( value_type v, size_type indx ) : value( v ), index( indx ){ }
917
918 static bool ascending( const __value_index_pair__ &a, const __value_index_pair__ &b )
919 {
920 return( a.value < b.value );
921 }
922
923 static bool descending( const __value_index_pair__ &a, const __value_index_pair__ &b )
924 {
925 return( a.value > b.value );
926 }
927 };
928
929 template < class T >
930 class __value_index_pair__< std::complex< T > >
931 {
932 public:
933 typedef size_t size_type;
934 typedef std::complex< T > value_type;
935
936 value_type value;
937 size_type index;
938
939 __value_index_pair__( ){ }
940
941 __value_index_pair__( value_type v, size_type indx ) : value( v ), index( indx ){ }
942
943 static bool ascending( const __value_index_pair__ &a, const __value_index_pair__ &b )
944 {
945 if( a.value.real( ) < b.value.real( ) )
946 {
947 return( true );
948 }
949 else if( a.value.real( ) > b.value.real( ) )
950 {
951 return( false );
952 }
953 else
954 {
955 return( a.value.imag( ) < b.value.imag( ) );
956 }
957 }
958
959 static bool descending( const __value_index_pair__ &a, const __value_index_pair__ &b )
960 {
961 if( a.value.real( ) > b.value.real( ) )
962 {
963 return( true );
964 }
965 else if( a.value.real( ) < b.value.real( ) )
966 {
967 return( false );
968 }
969 else
970 {
971 return( a.value.imag( ) > b.value.imag( ) );
972 }
973 }
974 };
975
976 inline const std::string to_string( integer value )
977 {
978 char buff[ 100 ];
979 sprintf( buff, "%ld", value );
980 return( std::string( buff ) );
981 }
982 }
983
984
985 namespace __solve__
986 {
987 // 行列の連立一次方程式を解く関数
988 template < bool >
989 struct __solve__
990 {
991 // 実数バージョン
992 template < class T, class Allocator >
993 static matrix< T, Allocator >& solve( matrix< T, Allocator > &a, matrix< T, Allocator > &b, matrix_style::style style )
994 {
995 typedef __clapack__::integer integer;
996 typedef typename matrix< T, Allocator >::value_type value_type;
997
998 integer info = 0;
999
1000 switch( style )
1001 {
1002 case matrix_style::sy:
1003 {
1004 if( a.rows( ) != b.rows( ) || a.empty( ) || b.empty( ) )
1005 {
1006 // 行列のサイズが正しくないので例外をスローする
1007 throw( numerical_exception( 1, "Incorrect matrix size is specified." ) );
1008 }
1009
1010 // LAPACK関数の引数
1011 integer n = static_cast< integer >( a.cols( ) );
1012 integer nrhs = static_cast< integer >( b.cols( ) );
1013 integer lda = static_cast< integer >( a.rows( ) );
1014 integer ldb = static_cast< integer >( b.rows( ) );
1015 value_type dmy;
1016 integer *ipiv = new integer[ n ];
1017 integer lwork = -1;
1018 char uplo[] = "U";
1019
1020 // まず最適な作業用配列のサイズを取得する
1021 __clapack__::sysv( uplo, n, nrhs, NULL, lda, NULL, NULL, ldb, &dmy, lwork, info );
1022 if( info == 0 )
1023 {
1024 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
1025 matrix< T, Allocator > work( lwork, 1 );
1026 __clapack__::sysv( uplo, n, nrhs, &( a[0] ), lda, ipiv, &( b[0] ), ldb, &( work[0] ), lwork, info );
1027 }
1028 delete [] ipiv;
1029
1030 if( info < 0 )
1031 {
1032 // 行列計算が正しく終了しなかったので例外をスローする
1033 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
1034 }
1035 else if( info > 0 )
1036 {
1037 // 行列計算が正しく終了しなかったので例外をスローする
1038 std::string val = __clapack__::to_string( info );
1039 std::string msg = "D( " + val + ", " + val + " ) is exactly zero. The factorization has been completed, but the block diagonal matrix D is exactly singular, so the solution could not be computed.";
1040 throw( numerical_exception( info, msg ) );
1041 }
1042 }
1043 break;
1044
1045 case matrix_style::pb:
1046 {
1047 if( a.cols( ) != b.rows( ) || a.empty( ) || b.empty( ) )
1048 {
1049 // 行列のサイズが正しくないので例外をスローする
1050 throw( numerical_exception( 1, "Incorrect matrix size is specified." ) );
1051 }
1052
1053 // LAPACK関数の引数
1054 integer n = static_cast< integer >( a.cols( ) );
1055 integer kd = static_cast< integer >( a.rows( ) ) - 1;
1056 integer nrhs = static_cast< integer >( b.cols( ) );
1057 integer ldab = static_cast< integer >( a.rows( ) );
1058 integer ldb = static_cast< integer >( b.rows( ) );
1059 char uplo[] = "U";
1060
1061 __clapack__::pbsv( uplo, n, kd, nrhs, &( a[0] ), ldab, &( b[0] ), ldb, info );
1062
1063 if( info < 0 )
1064 {
1065 // 行列計算が正しく終了しなかったので例外をスローする
1066 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
1067 }
1068 else if( info > 0 )
1069 {
1070 // 行列計算が正しく終了しなかったので例外をスローする
1071 std::string val = __clapack__::to_string( info );
1072 std::string msg = "The leading minor of order " + val + " of A is not positive definite, so the factorization could not be completed, and the solution has not been computed.";
1073 throw( numerical_exception( info, msg ) );
1074 }
1075 }
1076 break;
1077
1078 case matrix_style::ge:
1079 default:
1080 {
1081 if( a.rows( ) != b.rows( ) || a.empty( ) || b.empty( ) )
1082 {
1083 // 行列のサイズが正しくないので例外をスローする
1084 throw( numerical_exception( 1, "Incorrect matrix size is specified." ) );
1085 }
1086
1087 // LAPACK関数の引数
1088 integer n = static_cast< integer >( a.cols( ) );
1089 integer nrhs = static_cast< integer >( b.cols( ) );
1090 integer lda = static_cast< integer >( a.rows( ) );
1091 integer ldb = static_cast< integer >( b.rows( ) );
1092 integer *ipiv = new integer[ n ];
1093
1094 // integer n: マトリックス a の列数
1095 // integer nrhs: マトリックス b の列数
1096 // double* a: 連立方程式を解く係数の配列（ lda×n 要素 ）
1097 // integer lda: マトリックス a の行数
1098 // integer* ipiv: ピボット選択に用いる配列（ lda要素 ）
1099 // double* b: 連立方程式を解く係数の配列（ ldb×n 要素 ）
1100 // integer ldb: マトリックス b の行数
1101 // integer info: 正常終了したか否かを知らせてくれる変数
1102 //
1103 // この関数を呼ぶ場合，入力となる a の内容は変更される
1104 // 最終的な結果は b に代入される
1105 __clapack__::gesv( n, nrhs, &( a[0] ), lda, ipiv, &( b[0] ), ldb, info );
1106 delete [] ipiv;
1107
1108 if( info < 0 )
1109 {
1110 // 行列計算が正しく終了しなかったので例外をスローする
1111 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
1112 }
1113 else if( info > 0 )
1114 {
1115 // 行列計算が正しく終了しなかったので例外をスローする
1116 std::string val = __clapack__::to_string( info );
1117 std::string msg = "U( " + val + ", " + val + " ) is exactly zero. The factorization has been completed, but the factor U is exactly singular, so the solution could not be computed";
1118 throw( numerical_exception( info, msg ) );
1119 }
1120 }
1121 break;
1122 }
1123
1124 return( b );
1125 }
1126 };
1127
1128 template < >
1129 struct __solve__< true >
1130 {
1131 // 複素数バージョン
1132 template < class T, class Allocator >
1133 static matrix< T, Allocator >& solve( matrix< T, Allocator > a, matrix< T, Allocator > &b, matrix_style::style style )
1134 {
1135 typedef __clapack__::integer integer;
1136 typedef typename T::value_type value_type;
1137
1138 integer info = 0;
1139
1140 switch( style )
1141 {
1142 case matrix_style::sy:
1143 {
1144 if( a.rows( ) != b.rows( ) || a.empty( ) || b.empty( ) )
1145 {
1146 // 行列のサイズが正しくないので例外をスローする
1147 throw( numerical_exception( 1, "Incorrect matrix size is specified." ) );
1148 }
1149
1150 // LAPACK関数の引数
1151 integer n = static_cast< integer >( a.cols( ) );
1152 integer nrhs = static_cast< integer >( b.cols( ) );
1153 integer lda = static_cast< integer >( a.rows( ) );
1154 integer ldb = static_cast< integer >( b.rows( ) );
1155 value_type dmy;
1156 integer *ipiv = new integer[ n ];
1157 integer lwork = -1;
1158 char uplo[] = "U";
1159
1160 // まず最適な作業用配列のサイズを取得する
1161 __clapack__::sysv( uplo, n, nrhs, NULL, lda, NULL, NULL, ldb, &dmy, lwork, info );
1162 if( info == 0 )
1163 {
1164 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
1165 matrix< T, Allocator > work( lwork, 1 );
1166 __clapack__::sysv( uplo, n, nrhs, &( a[0] ), lda, ipiv, &( b[0] ), ldb, &( work[0] ), lwork, info );
1167 }
1168 delete [] ipiv;
1169
1170 if( info < 0 )
1171 {
1172 // 行列計算が正しく終了しなかったので例外をスローする
1173 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
1174 }
1175 else if( info > 0 )
1176 {
1177 // 行列計算が正しく終了しなかったので例外をスローする
1178 std::string val = __clapack__::to_string( info );
1179 std::string msg = "D( " + val + ", " + val + " ) is exactly zero. The factorization has been completed, but the block diagonal matrix D is exactly singular, so the solution could not be computed.";
1180 throw( numerical_exception( info, msg ) );
1181 }
1182 }
1183 break;
1184
1185 case matrix_style::pb:
1186 {
1187 if( a.cols( ) != b.rows( ) || a.empty( ) || b.empty( ) )
1188 {
1189 // 行列のサイズが正しくないので例外をスローする
1190 throw( numerical_exception( 1, "Incorrect matrix size is specified." ) );
1191 }
1192
1193 // LAPACK関数の引数
1194 integer n = static_cast< integer >( a.cols( ) );
1195 integer kd = static_cast< integer >( a.rows( ) ) - 1;
1196 integer nrhs = static_cast< integer >( b.cols( ) );
1197 integer ldab = static_cast< integer >( a.rows( ) );
1198 integer ldb = static_cast< integer >( b.rows( ) );
1199 char uplo[] = "U";
1200
1201 __clapack__::pbsv( uplo, n, kd, nrhs, &( a[0] ), ldab, &( b[0] ), ldb, info );
1202
1203 if( info < 0 )
1204 {
1205 // 行列計算が正しく終了しなかったので例外をスローする
1206 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
1207 }
1208 else if( info > 0 )
1209 {
1210 // 行列計算が正しく終了しなかったので例外をスローする
1211 std::string val = __clapack__::to_string( info );
1212 std::string msg = "The leading minor of order " + val + " of A is not positive definite, so the factorization could not be completed, and the solution has not been computed.";
1213 throw( numerical_exception( info, msg ) );
1214 }
1215 }
1216 break;
1217
1218 case matrix_style::ge:
1219 default:
1220 {
1221 if( a.rows( ) != b.rows( ) || a.empty( ) || b.empty( ) )
1222 {
1223 // 行列のサイズが正しくないので例外をスローする
1224 throw( numerical_exception( 1, "Incorrect matrix size is specified." ) );
1225 }
1226
1227 // LAPACK関数の引数
1228 integer n = static_cast< integer >( a.cols( ) );
1229 integer nrhs = static_cast< integer >( b.cols( ) );
1230 integer lda = static_cast< integer >( a.rows( ) );
1231 integer ldb = static_cast< integer >( b.rows( ) );
1232 integer *ipiv = new integer[ n ];
1233
1234 // integer n: マトリックス a の列数
1235 // integer nrhs: マトリックス b の列数
1236 // double* a: 連立方程式を解く係数の配列（ lda×n 要素 ）
1237 // integer lda: マトリックス a の行数
1238 // integer* ipiv: ピボット選択に用いる配列（ lda要素 ）
1239 // double* b: 連立方程式を解く係数の配列（ ldb×n 要素 ）
1240 // integer ldb: マトリックス b の行数
1241 // integer info: 正常終了したか否かを知らせてくれる変数
1242 //
1243 // この関数を呼ぶ場合，入力となる a の内容は変更される
1244 // 最終的な結果は b に代入される
1245 __clapack__::gesv( n, nrhs, &( a[0] ), lda, ipiv, &( b[0] ), ldb, info );
1246 delete [] ipiv;
1247
1248 if( info < 0 )
1249 {
1250 // 行列計算が正しく終了しなかったので例外をスローする
1251 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
1252 }
1253 else if( info > 0 )
1254 {
1255 // 行列計算が正しく終了しなかったので例外をスローする
1256 std::string val = __clapack__::to_string( info );
1257 std::string msg = "U( " + val + ", " + val + " ) is exactly zero. The factorization has been completed, but the factor U is exactly singular, so the solution could not be computed";
1258 throw( numerical_exception( info, msg ) );
1259 }
1260 }
1261 break;
1262 }
1263
1264 return( b );
1265 }
1266 };
1267 }
1268
1269
1270
1271
1272 namespace __lu__
1273 {
1274 // 行列の連立一次方程式を解く関数
1275 template < bool b >
1276 struct __lu__
1277 {
1278 // 実数バージョン
1279 template < class T, class Allocator1, class Allocator2 >
1280 static matrix< T, Allocator1 >& lu_factorization( matrix< T, Allocator1 > &a, matrix< __clapack__::integer, Allocator2 > &pivot, matrix_style::style style )
1281 {
1282 typedef __clapack__::integer integer;
1283 typedef typename matrix< T, Allocator1 >::value_type value_type;
1284
1285 if( a.empty( ) )
1286 {
1287 // 行列のサイズが正しくないので例外をスローする
1288 throw( numerical_exception( 1, "Empty matrix." ) );
1289 }
1290
1291 integer info = 0;
1292
1293 switch( style )
1294 {
1295 case matrix_style::sy:
1296 case matrix_style::ge:
1297 default:
1298 {
1299 // LAPACK関数の引数
1300 integer m = static_cast< integer >( a.rows( ) );
1301 integer n = static_cast< integer >( a.cols( ) );
1302 integer lda = m;
1303
1304 pivot.resize( n, 1 );
1305
1306 // LU分解を行う
1307 __clapack__::getrf( m, n, &( a[0] ), lda, &( pivot[0] ), info );
1308 }
1309 break;
1310 }
1311
1312 if( info < 0 )
1313 {
1314 // 行列計算が正しく終了しなかったので例外をスローする
1315 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
1316 }
1317 else if( info > 0 )
1318 {
1319 // 行列計算が正しく終了しなかったので例外をスローする
1320 std::string val = __clapack__::to_string( info );
1321 std::string msg = "U( " + val + ", " + val + " ) is exactly zero. The factorization has been completed, but the factor U is exactly singular, and division by zero will occur if it is used to solve a system of equations.";
1322 throw( numerical_exception( info, msg ) );
1323 }
1324
1325 return( a );
1326 }
1327 };
1328
1329 template < >
1330 struct __lu__< true >
1331 {
1332 // 複素数バージョン
1333 template < class T, class Allocator1, class Allocator2 >
1334 static matrix< T, Allocator1 >& lu_factorization( matrix< T, Allocator1 > &a, matrix< __clapack__::integer, Allocator2 > &pivot, matrix_style::style style )
1335 {
1336 typedef __clapack__::integer integer;
1337 typedef typename T::value_type value_type;
1338
1339 if( a.empty( ) )
1340 {
1341 // 行列のサイズが正しくないので例外をスローする
1342 throw( numerical_exception( 1, "Empty matrix." ) );
1343 }
1344
1345 integer info = 0;
1346
1347 switch( style )
1348 {
1349 case matrix_style::sy:
1350 case matrix_style::ge:
1351 default:
1352 {
1353 // LAPACK関数の引数
1354 integer m = static_cast< integer >( a.rows( ) );
1355 integer n = static_cast< integer >( a.cols( ) );
1356 integer lda = m;
1357
1358 pivot.resize( n, 1 );
1359
1360 // LU分解を行う
1361 __clapack__::getrf( m, n, &( a[0] ), lda, &( pivot[0] ), info );
1362 }
1363 break;
1364 }
1365
1366 if( info < 0 )
1367 {
1368 // 行列計算が正しく終了しなかったので例外をスローする
1369 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
1370 }
1371 else if( info > 0 )
1372 {
1373 // 行列計算が正しく終了しなかったので例外をスローする
1374 std::string val = __clapack__::to_string( info );
1375 std::string msg = "U( " + val + ", " + val + " ) is exactly zero. The factorization has been completed, but the factor U is exactly singular, and division by zero will occur if it is used to solve a system of equations.";
1376 throw( numerical_exception( info, msg ) );
1377 }
1378
1379 return( a );
1380 }
1381 };
1382 }
1383
1384
1385
1386
1387 namespace __qr__
1388 {
1389 // 行列の連立一次方程式を解く関数
1390 template < bool b >
1391 struct __qr__
1392 {
1393 // 実数バージョン
1394 template < class T, class Allocator >
1395 static void qr_factorization( matrix< T, Allocator > &A, matrix< T, Allocator > &Q, matrix< T, Allocator > &R, matrix_style::style style )
1396 {
1397 typedef __clapack__::integer integer;
1398 typedef typename matrix< T, Allocator >::value_type value_type;
1399 typedef typename matrix< T, Allocator >::size_type size_type;
1400
1401 if( A.empty( ) )
1402 {
1403 // 行列のサイズが正しくないので例外をスローする
1404 throw( numerical_exception( 1, "Empty matrix." ) );
1405 }
1406 else if( A.rows( ) < A.cols( ) )
1407 {
1408 // 行列のサイズが正しくないので例外をスローする
1409 throw( numerical_exception( 1, "The number of rows should be larger than the number of columns." ) );
1410 }
1411
1412 integer info = 0;
1413
1414 switch( style )
1415 {
1416 case matrix_style::sy:
1417 case matrix_style::ge:
1418 default:
1419 {
1420 // LAPACK関数の引数
1421 matrix< T, Allocator > &a = A;
1422 integer m = static_cast< integer >( a.rows( ) );
1423 integer n = static_cast< integer >( a.cols( ) );
1424 integer dim = m < n ? m : n;
1425 integer lda = m;
1426 typename matrix< T, Allocator >::value_type dmy;
1427 integer lwork = -1;
1428
1429 // QR分解を行う前に，必要な作業用配列のサイズを取得する
1430 __clapack__::geqrf( m, n, NULL, lda, NULL, &dmy, lwork, info );
1431 if( info == 0 )
1432 {
1433 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
1434 matrix< T, Allocator > work( lwork, 1 ), tau( dim, 1 );
1435 __clapack__::geqrf( m, n, &( a[0] ), lda, &( tau[0] ), &( work[0] ), lwork, info );
1436
1437 if( info == 0 )
1438 {
1439 // QR分解に成功したので、結果からQとRの行列を求める
1440 R.resize( m, n );
1441 Q.resize( m, m );
1442 for( size_type c = 0 ; c < R.cols( ) ; c++ )
1443 {
1444 for( size_type r = 0 ; r <= c ; r++ )
1445 {
1446 R( r, c ) = a( r, c );
1447 }
1448 for( size_type r = 0 ; r < R.rows( ) ; r++ )
1449 {
1450 Q( r, c ) = a( r, c );
1451 }
1452 }
1453
1454 __clapack__::orgqr( m, m, dim, &( Q[0] ), lda, &( tau[0] ), &( work[0] ), lwork, info );
1455 }
1456 }
1457 }
1458 break;
1459 }
1460
1461 if( info != 0 )
1462 {
1463 // 行列計算が正しく終了しなかったので例外をスローする
1464 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
1465 }
1466 }
1467 };
1468
1469 template < >
1470 struct __qr__< true >
1471 {
1472 // 複素数バージョン
1473 template < class T, class Allocator >
1474 static void qr_factorization( matrix< T, Allocator > &A, matrix< T, Allocator > &Q, matrix< T, Allocator > &R, matrix_style::style style )
1475 {
1476 typedef __clapack__::integer integer;
1477 typedef typename matrix< T, Allocator >::value_type value_type;
1478 typedef typename matrix< T, Allocator >::size_type size_type;
1479
1480 if( A.empty( ) )
1481 {
1482 // 行列のサイズが正しくないので例外をスローする
1483 throw( numerical_exception( 1, "Empty matrix." ) );
1484 }
1485 else if( A.rows( ) < A.cols( ) )
1486 {
1487 // 行列のサイズが正しくないので例外をスローする
1488 throw( numerical_exception( 1, "The number of rows should be larger than the number of columns." ) );
1489 }
1490
1491 integer info = 0;
1492
1493 switch( style )
1494 {
1495 case matrix_style::ge:
1496 default:
1497 {
1498 // LAPACK関数の引数
1499 matrix< T, Allocator > &a = A;
1500 integer m = static_cast< integer >( a.rows( ) );
1501 integer n = static_cast< integer >( a.cols( ) );
1502 integer dim = m < n ? m : n;
1503 integer lda = m;
1504 typename matrix< T, Allocator >::value_type dmy;
1505 integer lwork = -1;
1506
1507 // QR分解を行う前に，必要な作業用配列のサイズを取得する
1508 __clapack__::geqrf( m, n, NULL, lda, NULL, &dmy, lwork, info );
1509 if( info == 0 )
1510 {
1511 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
1512 matrix< T, Allocator > work( lwork, 1 ), tau( dim, 1 );
1513 __clapack__::geqrf( m, n, &( a[0] ), lda, &( tau[0] ), &( work[0] ), lwork, info );
1514
1515 if( info == 0 )
1516 {
1517 // QR分解に成功したので、結果からQとRの行列を求める
1518 R.resize( m, n );
1519 Q.resize( m, m );
1520 for( size_type c = 0 ; c < R.cols( ) ; c++ )
1521 {
1522 for( size_type r = 0 ; r <= c ; r++ )
1523 {
1524 R( r, c ) = a( r, c );
1525 }
1526 for( size_type r = 0 ; r < R.rows( ) ; r++ )
1527 {
1528 Q( r, c ) = a( r, c );
1529 }
1530 }
1531
1532 __clapack__::ungqr( m, n, dim, &( Q[0] ), lda, &( tau[0] ), &( work[0] ), lwork, info );
1533 }
1534 }
1535 }
1536 break;
1537 }
1538
1539 if( info != 0 )
1540 {
1541 // 行列計算が正しく終了しなかったので例外をスローする
1542 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
1543 }
1544 }
1545 };
1546 }
1547
1548
1549 namespace __cholesky__
1550 {
1551 // コレスキー分解を行う
1552 template < class T, class Allocator >
1553 static matrix< T, Allocator >& cholesky_factorization( matrix< T, Allocator > &a, matrix_style::style style )
1554 {
1555 typedef matrix< T, Allocator > matrix_type;
1556 typedef typename matrix_type::value_type value_type;
1557 typedef typename matrix_type::size_type size_type;
1558 typedef __clapack__::integer integer;
1559
1560 if( a.empty( ) )
1561 {
1562 // 行列のサイズが正しくないので例外をスローする
1563 throw( numerical_exception( 1, "Empty matrix." ) );
1564 }
1565 else if( a.rows( ) != a.cols( ) )
1566 {
1567 // 行列のサイズが正しくないので例外をスローする
1568 throw( numerical_exception( 1, "The number of rows should be equal to the number of columns." ) );
1569 }
1570
1571 integer info = 0;
1572
1573 // コレスキー分解で上三角は参照しないので，下三角行列にしておく
1574 for( size_type r = 0 ; r < a.rows( ) ; r++ )
1575 {
1576 for( size_type c = r + 1 ; c < a.cols( ) ; c++ )
1577 {
1578 a( r, c ) = static_cast< value_type >( 0 );
1579 }
1580 }
1581
1582 switch( style )
1583 {
1584 case matrix_style::sy:
1585 //case matrix_style::ge:
1586 default:
1587 {
1588 // LAPACK関数の引数
1589 integer n = static_cast< integer >( a.cols( ) );
1590 integer lda = static_cast< integer >( a.rows( ) );
1591 char uplo[] = "L";
1592
1593 // まず最適な作業用配列のサイズを取得する
1594 __clapack__::potrf( uplo, n, &( a[0] ), lda, info );
1595 }
1596 break;
1597 }
1598
1599 if( info < 0 )
1600 {
1601 // 行列計算が正しく終了しなかったので例外をスローする
1602 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
1603 }
1604 else if( info > 0 )
1605 {
1606 // 行列計算が正しく終了しなかったので例外をスローする
1607 std::string val = __clapack__::to_string( info );
1608 std::string msg = "The leading minor of order " + val + " is not positive definite, and the factorization could not be completed.";
1609 throw( numerical_exception( info, msg ) );
1610 }
1611
1612 return( a );
1613 }
1614 }
1615
1616
1617 namespace __inverse__
1618 {
1619 // 行列の連立一次方程式を解く関数
1620 template < bool b >
1621 struct __inverse__
1622 {
1623 // 実数バージョン
1624 template < class T, class Allocator >
1625 static matrix< T, Allocator >& inverse( matrix< T, Allocator > &a, matrix_style::style style )
1626 {
1627 typedef matrix< T, Allocator > matrix_type;
1628 typedef typename matrix< T, Allocator >::value_type value_type;
1629 typedef __clapack__::integer integer;
1630
1631 if( a.empty( ) )
1632 {
1633 // 行列のサイズが正しくないので例外をスローする
1634 throw( numerical_exception( 1, "Empty matrix." ) );
1635 }
1636
1637 integer info = 0;
1638
1639 if( a.rows( ) == a.cols( ) )
1640 {
1641 // ２行２列と３行３列用の特殊バージョン
1642 switch( a.rows( ) )
1643 {
1644 case 2:
1645 {
1646 double a11 = a( 0, 0 );
1647 double a12 = a( 0, 1 );
1648 double a21 = a( 1, 0 );
1649 double a22 = a( 1, 1 );
1650 double detA = a11 * a22 - a12 * a21;
1651 if( detA != 0 )
1652 {
1653 a( 0, 0 ) = static_cast< value_type >( a22 / detA );
1654 a( 0, 1 ) = static_cast< value_type >( -a12 / detA );
1655 a( 1, 0 ) = static_cast< value_type >( -a21 / detA );
1656 a( 1, 1 ) = static_cast< value_type >( a11 / detA );
1657 return( a );
1658 }
1659 else
1660 {
1661 // 行列式がゼロ
1662 throw( numerical_exception( 1, "Determinant of the matrix is zero." ) );
1663 }
1664 }
1665 break;
1666
1667 case 3:
1668 {
1669 double a11 = a( 0, 0 ), a12 = a( 0, 1 ), a13 = a( 0, 2 );
1670 double a21 = a( 1, 0 ), a22 = a( 1, 1 ), a23 = a( 1, 2 );
1671 double a31 = a( 2, 0 ), a32 = a( 2, 1 ), a33 = a( 2, 2 );
1672
1673 // 共通で利用する係数の計算
1674 double _22x33_23x32_ = a22 * a33 - a23 * a32;
1675 double _21x32_22x31_ = a21 * a32 - a22 * a31;
1676 double _23x31_21x33_ = a23 * a31 - a21 * a33;
1677
1678 // 行列式を計算
1679 double detA = a11 * _22x33_23x32_ + a13 * _21x32_22x31_ + a12 * _23x31_21x33_;
1680
1681 // 逆行列が存在する場合のも逆行列を計算する
1682 if( detA != 0 )
1683 {
1684 // 各要素の値を計算
1685 double A11 = _22x33_23x32_;
1686 double A12 = a13 * a32 - a12 * a33;
1687 double A13 = a12 * a23 - a13 * a22;
1688 double A21 = _23x31_21x33_;
1689 double A22 = a11 * a33 - a13 * a31;
1690 double A23 = a13 * a21 - a11 * a23;
1691 double A31 = _21x32_22x31_;
1692 double A32 = a12 * a31 - a11 * a32;
1693 double A33 = a11 * a22 - a12 * a21;
1694
1695 double _1_detA = 1.0 / detA;
1696 a( 0, 0 ) = static_cast< value_type >( A11 * _1_detA );
1697 a( 0, 1 ) = static_cast< value_type >( A12 * _1_detA );
1698 a( 0, 2 ) = static_cast< value_type >( A13 * _1_detA );
1699 a( 1, 0 ) = static_cast< value_type >( A21 * _1_detA );
1700 a( 1, 1 ) = static_cast< value_type >( A22 * _1_detA );
1701 a( 1, 2 ) = static_cast< value_type >( A23 * _1_detA );
1702 a( 2, 0 ) = static_cast< value_type >( A31 * _1_detA );
1703 a( 2, 1 ) = static_cast< value_type >( A32 * _1_detA );
1704 a( 2, 2 ) = static_cast< value_type >( A33 * _1_detA );
1705
1706 return( a );
1707 }
1708 else
1709 {
1710 // 行列式がゼロ
1711 throw( numerical_exception( 1, "Determinant of the matrix is zero." ) );
1712 }
1713 }
1714 break;
1715
1716 case 4:
1717 {
1718 double a11 = a( 0, 0 ), a12 = a( 0, 1 ), a13 = a( 0, 2 ), a14 = a( 0, 3 );
1719 double a21 = a( 1, 0 ), a22 = a( 1, 1 ), a23 = a( 1, 2 ), a24 = a( 1, 3 );
1720 double a31 = a( 2, 0 ), a32 = a( 2, 1 ), a33 = a( 2, 2 ), a34 = a( 2, 3 );
1721 double a41 = a( 3, 0 ), a42 = a( 3, 1 ), a43 = a( 3, 2 ), a44 = a( 3, 3 );
1722
1723 // 共通で利用する係数の計算
1724 double _33x44_34x43_ = a33 * a44 - a34 * a43;
1725 double _34x42_32x44_ = a34 * a42 - a32 * a44;
1726 double _32x43_33x42_ = a32 * a43 - a33 * a42;
1727 double _31x44_34x41_ = a31 * a44 - a34 * a41;
1728 double _33x41_31x43_ = a33 * a41 - a31 * a43;
1729 double _31x42_32x41_ = a31 * a42 - a32 * a41;
1730
1731 // 行列式を計算
1732 double detA = a11 * ( a22 * _33x44_34x43_ + a23 * _34x42_32x44_ + a24 * _32x43_33x42_ )
1733 - a12 * ( a21 * _33x44_34x43_ - a23 * _31x44_34x41_ - a24 * _33x41_31x43_ )
1734 - a13 * ( a21 * _34x42_32x44_ + a22 * _31x44_34x41_ - a24 * _31x42_32x41_ )
1735 - a14 * ( a21 * _32x43_33x42_ + a22 * _33x41_31x43_ + a23 * _31x42_32x41_ );
1736
1737 // 逆行列が存在する場合のも逆行列を計算する
1738 if( detA != 0 )
1739 {
1740 // 共通で利用する係数の計算
1741 double _23x44_24x43_ = a23 * a44 - a24 * a43;
1742 double _24x42_22x44_ = a24 * a42 - a22 * a44;
1743 double _22x43_23x42_ = a22 * a43 - a23 * a42;
1744 double _24x33_23x34_ = a24 * a33 - a23 * a34;
1745 double _22x34_24x32_ = a22 * a34 - a24 * a32;
1746 double _23x32_22x33_ = a23 * a32 - a22 * a33;
1747 double _21x44_24x41_ = a21 * a44 - a24 * a41;
1748 double _23x41_21x43_ = a23 * a41 - a21 * a43;
1749 double _24x31_21x34_ = a24 * a31 - a21 * a34;
1750 double _21x33_23x31_ = a21 * a33 - a23 * a31;
1751 double _21x42_22x41_ = a21 * a42 - a22 * a41;
1752 double _22x31_21x32_ = a22 * a31 - a21 * a32;
1753
1754 // 各要素の値を計算
1755 double A11 = a22 * _33x44_34x43_ + a23 * _34x42_32x44_ + a24 * _32x43_33x42_;
1756 double A12 = -a12 * _33x44_34x43_ - a13 * _34x42_32x44_ - a14 * _32x43_33x42_;
1757 double A13 = a12 * _23x44_24x43_ + a13 * _24x42_22x44_ + a14 * _22x43_23x42_;
1758 double A14 = a12 * _24x33_23x34_ + a13 * _22x34_24x32_ + a14 * _23x32_22x33_;
1759 double A21 = -a21 * _33x44_34x43_ + a23 * _31x44_34x41_ + a24 * _33x41_31x43_;
1760 double A22 = a11 * _33x44_34x43_ - a13 * _31x44_34x41_ - a14 * _33x41_31x43_;
1761 double A23 = -a11 * _23x44_24x43_ + a13 * _21x44_24x41_ + a14 * _23x41_21x43_;
1762 double A24 = -a11 * _24x33_23x34_ + a13 * _24x31_21x34_ + a14 * _21x33_23x31_;
1763 double A31 = -a21 * _34x42_32x44_ - a22 * _31x44_34x41_ + a24 * _31x42_32x41_;
1764 double A32 = a11 * _34x42_32x44_ + a12 * _31x44_34x41_ - a14 * _31x42_32x41_;
1765 double A33 = -a11 * _24x42_22x44_ - a12 * _21x44_24x41_ + a14 * _21x42_22x41_;
1766 double A34 = -a11 * _22x34_24x32_ - a12 * _24x31_21x34_ + a14 * _22x31_21x32_;
1767 double A41 = -a21 * _32x43_33x42_ - a22 * _33x41_31x43_ - a23 * _31x42_32x41_;
1768 double A42 = a11 * _32x43_33x42_ + a12 * _33x41_31x43_ + a13 * _31x42_32x41_;
1769 double A43 = -a11 * _22x43_23x42_ - a12 * _23x41_21x43_ - a13 * _21x42_22x41_;
1770 double A44 = -a11 * _23x32_22x33_ - a12 * _21x33_23x31_ - a13 * _22x31_21x32_;
1771
1772 double _1_detA = 1.0 / detA;
1773 a( 0, 0 ) = static_cast< value_type >( A11 * _1_detA );
1774 a( 0, 1 ) = static_cast< value_type >( A12 * _1_detA );
1775 a( 0, 2 ) = static_cast< value_type >( A13 * _1_detA );
1776 a( 0, 3 ) = static_cast< value_type >( A14 * _1_detA );
1777 a( 1, 0 ) = static_cast< value_type >( A21 * _1_detA );
1778 a( 1, 1 ) = static_cast< value_type >( A22 * _1_detA );
1779 a( 1, 2 ) = static_cast< value_type >( A23 * _1_detA );
1780 a( 1, 3 ) = static_cast< value_type >( A24 * _1_detA );
1781 a( 2, 0 ) = static_cast< value_type >( A31 * _1_detA );
1782 a( 2, 1 ) = static_cast< value_type >( A32 * _1_detA );
1783 a( 2, 2 ) = static_cast< value_type >( A33 * _1_detA );
1784 a( 2, 3 ) = static_cast< value_type >( A34 * _1_detA );
1785 a( 3, 0 ) = static_cast< value_type >( A41 * _1_detA );
1786 a( 3, 1 ) = static_cast< value_type >( A42 * _1_detA );
1787 a( 3, 2 ) = static_cast< value_type >( A43 * _1_detA );
1788 a( 3, 3 ) = static_cast< value_type >( A44 * _1_detA );
1789 return( a );
1790 }
1791 else
1792 {
1793 // 行列式がゼロ
1794 throw( numerical_exception( 1, "Determinant of the matrix is zero." ) );
1795 }
1796 }
1797 break;
1798
1799 default:
1800 break;
1801 }
1802 }
1803
1804 if( info == 0 )
1805 {
1806 switch( style )
1807 {
1808 case matrix_style::sy:
1809 {
1810 // LAPACK関数の引数
1811 integer n = static_cast< integer >( a.cols( ) );
1812 integer lda = static_cast< integer >( a.rows( ) );
1813 integer *ipiv = new integer[ n ];
1814 typename matrix< T, Allocator >::value_type dmy;
1815 integer lwork = -1;
1816 char uplo[] = "U";
1817
1818 // LU分解を行う
1819 // まず最適な作業用配列のサイズを取得する
1820 __clapack__::sytrf( uplo, n, NULL, lda, ipiv, &dmy, lwork, info );
1821 if( info == 0 )
1822 {
1823 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
1824 matrix< T, Allocator > work( lwork, 1 );
1825 __clapack__::sytrf( uplo, n, &( a[0] ), lda, ipiv, &( work[0] ), lwork, info );
1826
1827 if( info == 0 )
1828 {
1829 if( lwork < n )
1830 {
1831 work.resize( n, 1 );
1832 }
1833 __clapack__::sytri( uplo, n, &( a[0] ), lda, ipiv, &( work[0] ), info );
1834
1835 // 結果が上三角のみにしか代入されていないため，下三角にもデータをコピーする
1836 if( info == 0 )
1837 {
1838 typedef typename matrix< T, Allocator >::size_type size_type;
1839 for( size_type c = 0 ; c < a.rows( ) ; c++ )
1840 {
1841 for( size_type r = c + 1 ; r < a.rows( ) ; r++ )
1842 {
1843 a( r, c ) = a( c, r );
1844 }
1845 }
1846 }
1847 }
1848 }
1849 delete [] ipiv;
1850 }
1851 break;
1852
1853 case matrix_style::ge:
1854 default:
1855 {
1856 // LAPACK関数の引数
1857 integer n = static_cast< integer >( a.cols( ) );
1858 integer lda = static_cast< integer >( a.rows( ) );
1859 integer *ipiv = new integer[ n ];
1860 typename matrix< T, Allocator >::value_type dmy;
1861 integer lwork = -1;
1862
1863 // LU分解を行う
1864 __clapack__::getrf( lda, n, &( a[0] ), lda, ipiv, info );
1865 if( info == 0 )
1866 {
1867 // まず最適な作業用配列のサイズを取得する
1868 __clapack__::getri( n, NULL, lda, NULL, &dmy, lwork, info );
1869 if( info == 0 )
1870 {
1871 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
1872 matrix< T, Allocator > work( lwork, 1 );
1873 __clapack__::getri( n, &( a[0] ), lda, ipiv, &( work[0] ), lwork, info );
1874 }
1875 }
1876 delete [] ipiv;
1877 }
1878 break;
1879 }
1880 }
1881
1882 if( info < 0 )
1883 {
1884 // 行列計算が正しく終了しなかったので例外をスローする
1885 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
1886 }
1887 else if( info > 0 )
1888 {
1889 // 行列計算が正しく終了しなかったので例外をスローする
1890 std::string val = __clapack__::to_string( info );
1891 std::string msg = "U( " + val + ", " + val + " ) is exactly zero. the matrix is singular and its inverse could not be computed";
1892 throw( numerical_exception( info, msg ) );
1893 }
1894
1895 return( a );
1896 }
1897 };
1898
1899 template < >
1900 struct __inverse__< true >
1901 {
1902 // 複素数バージョン
1903 template < class T, class Allocator >
1904 static matrix< T, Allocator >& inverse( matrix< T, Allocator > &a, matrix_style::style style )
1905 {
1906 typedef __clapack__::integer integer;
1907 typedef typename T::value_type value_type;
1908
1909 if( a.empty( ) )
1910 {
1911 // 行列のサイズが正しくないので例外をスローする
1912 throw( numerical_exception( 1, "Empty matrix." ) );
1913 }
1914
1915 integer info = 0;
1916
1917 switch( style )
1918 {
1919 case matrix_style::sy:
1920 {
1921 // LAPACK関数の引数
1922 integer n = static_cast< integer >( a.cols( ) );
1923 integer lda = static_cast< integer >( a.rows( ) );
1924 integer *ipiv = new integer[ n ];
1925 typename matrix< T, Allocator >::value_type dmy;
1926 integer lwork = -1;
1927 char uplo[] = "U";
1928
1929 // LU分解を行う
1930 // まず最適な作業用配列のサイズを取得する
1931 __clapack__::sytrf( uplo, n, NULL, lda, ipiv, &dmy, lwork, info );
1932 if( info == 0 )
1933 {
1934 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
1935 matrix< T, Allocator > work( lwork, 1 );
1936 __clapack__::sytrf( uplo, n, &( a[0] ), lda, ipiv, &( work[0] ), lwork, info );
1937
1938 if( info == 0 )
1939 {
1940 if( lwork < 2 * n )
1941 {
1942 work.resize( 2 * n, 1 );
1943 }
1944 __clapack__::sytri( uplo, n, &( a[0] ), lda, ipiv, &( work[0] ), info );
1945
1946 // 結果が上三角のみにしか代入されていないため，下三角にもデータをコピーする
1947 if( info == 0 )
1948 {
1949 typedef typename matrix< T, Allocator >::size_type size_type;
1950 for( size_type c = 0 ; c < a.rows( ) ; c++ )
1951 {
1952 for( size_type r = c + 1 ; r < a.rows( ) ; r++ )
1953 {
1954 a( r, c ) = a( c, r );
1955 }
1956 }
1957 }
1958 }
1959 }
1960 delete [] ipiv;
1961 }
1962 break;
1963
1964 case matrix_style::ge:
1965 default:
1966 {
1967 // LAPACK関数の引数
1968 integer n = static_cast< integer >( a.cols( ) );
1969 integer lda = static_cast< integer >( a.rows( ) );
1970 integer *ipiv = new integer[ n ];
1971 typename matrix< T, Allocator >::value_type dmy;
1972 integer lwork = -1;
1973
1974 // LU分解を行う
1975 __clapack__::getrf( lda, n, &( a[0] ), lda, ipiv, info );
1976 if( info == 0 )
1977 {
1978 // まず最適な作業用配列のサイズを取得する
1979 __clapack__::getri( n, NULL, lda, NULL, &dmy, lwork, info );
1980 if( info == 0 )
1981 {
1982 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
1983 matrix< T, Allocator > work( lwork, 1 );
1984 __clapack__::getri( n, &( a[0] ), lda, ipiv, &( work[0] ), lwork, info );
1985 }
1986 }
1987 delete [] ipiv;
1988 }
1989 break;
1990 }
1991
1992 if( info < 0 )
1993 {
1994 // 行列計算が正しく終了しなかったので例外をスローする
1995 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
1996 }
1997 else if( info > 0 )
1998 {
1999 // 行列計算が正しく終了しなかったので例外をスローする
2000 std::string val = __clapack__::to_string( info );
2001 std::string msg = "U( " + val + ", " + val + " ) is exactly zero. the matrix is singular and its inverse could not be computed";
2002 throw( numerical_exception( info, msg ) );
2003 }
2004
2005 return( a );
2006 }
2007 };
2008 }
2009
2010
2011
2012 namespace __eigen__
2013 {
2014 // 一般行列の固有値・固有ベクトルを計算する
2015 #if defined(_USE_BALANCING_MATRIX_EIGEN_) && _USE_BALANCING_MATRIX_EIGEN_ != 0
2016
2017 // 行列の対角化を行うバージョン
2018 template < bool b >
2019 struct __eigen__
2020 {
2021 // 実数バージョン
2022 template < class T, class Allocator >
2023 static matrix< T, Allocator >& eigen( matrix< T, Allocator > &a, matrix< T, Allocator > &eigen_value, matrix< T, Allocator > &eigen_vector, matrix_style::style style )
2024 {
2025 typedef __clapack__::integer integer;
2026 typedef typename matrix< T, Allocator >::value_type value_type;
2027
2028 if( a.empty( ) )
2029 {
2030 // 行列のサイズが正しくないので例外をスローする
2031 throw( numerical_exception( 1, "Empty matrix." ) );
2032 }
2033
2034 integer info = 0;
2035
2036 switch( style )
2037 {
2038 case matrix_style::sy:
2039 {
2040 // LAPACK関数の引数
2041 integer n = static_cast< integer >( a.cols( ) );
2042 integer m = n;
2043 integer lda = static_cast< integer >( a.rows( ) );
2044 integer ldz = n;
2045 value_type dmy = 0;
2046 integer lwork = -1;
2047 value_type vl = 0;
2048 value_type vu = 0;
2049 integer il = 0;
2050 integer iu = 0;
2051 value_type abstol = 0;
2052 char jobz[] = "V";
2053 char range[] = "A";
2054 char uplo[] = "U";
2055
2056 // まず最適な作業用配列のサイズを取得する
2057 __clapack__::syevx( jobz, range, uplo, n, NULL, lda, vl, vu, il, iu, abstol, m, NULL, NULL, ldz, &dmy, lwork, NULL, NULL, info );
2058 if( info == 0 )
2059 {
2060 eigen_value.resize( n, 1 );
2061 matrix< integer > iwork( 5 * n, 1 );
2062 matrix< integer > ifail( n, 1 );
2063 eigen_vector.resize( n, n );
2064
2065 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
2066 matrix< T, Allocator > work( lwork, 1 );
2067 __clapack__::syevx( jobz, range, uplo, n, &( a[0] ), lda, vl, vu, il, iu, abstol, m,
2068 &( eigen_value[0] ), &( eigen_vector[0] ), ldz, &( work[0] ), lwork, &( iwork[ 0 ] ), &( ifail[ 0 ] ), info );
2069 }
2070
2071 if( info < 0 )
2072 {
2073 // 行列計算が正しく終了しなかったので例外をスローする
2074 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
2075 }
2076 else if( info > 0 )
2077 {
2078 // 行列計算が正しく終了しなかったので例外をスローする
2079 throw( numerical_exception( info, __clapack__::to_string( info ) + " eigenvectors failed to converge." ) );
2080 }
2081 }
2082 break;
2083
2084 case matrix_style::ge:
2085 default:
2086 {
2087 // LAPACK関数の引数
2088 integer n = static_cast< integer >( a.cols( ) );
2089 integer lda = static_cast< integer >( a.rows( ) );
2090 typename matrix< T, Allocator >::value_type dmy, abnrm;
2091 integer ldvl = 1;
2092 integer ldvr = n;
2093 integer lwork = -1;
2094 integer ilo = 0;
2095 integer ihi = 0;
2096 char balanc[] = "B";
2097 char jobvl[] = "N";
2098 char jobvr[] = "V";
2099 char sense[] = "N";
2100
2101 // まず最適な作業用配列のサイズを取得する
2102 __clapack__::geevx( balanc, jobvl, jobvr, sense, n, NULL, lda, NULL, NULL, NULL, ldvl, NULL, ldvr, ilo, ihi, NULL, abnrm, NULL, NULL, &dmy, lwork, NULL, info );
2103 if( info == 0 )
2104 {
2105 eigen_value.resize( n, 1 );
2106 matrix< T, Allocator > tmp( n, 1 );
2107 matrix< T, Allocator > scale( n, 1 );
2108 eigen_vector.resize( n, n );
2109
2110 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
2111 matrix< T, Allocator > work( lwork, 1 );
2112 __clapack__::geevx( balanc, jobvl, jobvr, sense, n, &( a[0] ), lda, &( eigen_value[0] ), &( tmp[0] ),
2113 NULL, ldvl, &( eigen_vector[0] ), ldvr, ilo, ihi, &( scale[0] ), abnrm, NULL, NULL, &( work[0] ), lwork, NULL, info );
2114 }
2115
2116 if( info < 0 )
2117 {
2118 // 行列計算が正しく終了しなかったので例外をスローする
2119 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
2120 }
2121 else if( info > 0 )
2122 {
2123 // 行列計算が正しく終了しなかったので例外をスローする
2124 throw( numerical_exception( info, "QR algorithm failed to compute all the eigenvalues, and no eigenvectors or condition numbers have been computed." ) );
2125 }
2126 }
2127 break;
2128 }
2129
2130 return( eigen_value );
2131 }
2132 };
2133
2134 template < >
2135 struct __eigen__< true >
2136 {
2137 // 複素数バージョン
2138 template < class T, class Allocator >
2139 static matrix< T, Allocator >& eigen( matrix< T, Allocator > &a, matrix< T, Allocator > &eigen_value, matrix< T, Allocator > &eigen_vector, matrix_style::style style )
2140 {
2141 typedef __clapack__::integer integer;
2142 typedef typename T::value_type value_type;
2143
2144 if( a.empty( ) )
2145 {
2146 // 行列のサイズが正しくないので例外をスローする
2147 throw( numerical_exception( 1, "Empty matrix." ) );
2148 }
2149
2150 integer info = 0;
2151
2152 switch( style )
2153 {
2154 case matrix_style::ge:
2155 default:
2156 {
2157 // LAPACK関数の引数
2158 integer n = static_cast< integer >( a.cols( ) );
2159 integer lda = static_cast< integer >( a.rows( ) );
2160 typename matrix< T, Allocator >::value_type dmy;
2161 value_type abnrm;
2162 integer ldvl = 1;
2163 integer ldvr = n;
2164 integer lwork = -1;
2165 integer ilo = 0;
2166 integer ihi = 0;
2167 char balanc[] = "B";
2168 char jobvl[] = "N";
2169 char jobvr[] = "V";
2170 char sense[] = "N";
2171
2172 // まず最適な作業用配列のサイズを取得する
2173 __clapack__::geevx( balanc, jobvl, jobvr, sense, n, NULL, lda, NULL, NULL, ldvl, NULL, ldvr, ilo, ihi, NULL, abnrm, NULL, NULL, &dmy, lwork, NULL, info );
2174 if( info == 0 )
2175 {
2176 eigen_value.resize( n, 1 );
2177 eigen_vector.resize( n, n );
2178 value_type *scale = new value_type[ n ];
2179 value_type *rwork = new value_type[ 2 * n ];
2180
2181 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
2182 matrix< T, Allocator > work( lwork, 1 );
2183 __clapack__::geevx( balanc, jobvl, jobvr, sense, n, &( a[0] ), lda, &( eigen_value[0] ),
2184 NULL, ldvl, &( eigen_vector[0] ), ldvr, ilo, ihi, scale, abnrm, NULL, NULL, &( work[0] ), lwork, rwork, info );
2185
2186 delete [] scale;
2187 delete [] rwork;
2188 }
2189
2190 if( info < 0 )
2191 {
2192 // 行列計算が正しく終了しなかったので例外をスローする
2193 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
2194 }
2195 else if( info > 0 )
2196 {
2197 // 行列計算が正しく終了しなかったので例外をスローする
2198 if( info <= n )
2199 {
2200 std::string val = __clapack__::to_string( info );
2201 std::string msg = "U( " + val + ", " + val + " ) is exactly zero. The factorization has been completed, but the factor U is exactly singular, so the solution and error bounds could not be computed. reciprocal condition number = 0 is returned.";
2202 throw( numerical_exception( info, msg ) );
2203 }
2204 else
2205 {
2206 std::string msg = "U is nonsingular, but reciprocal condition number (RCOND) is less than machine precision, meaning that the matrix is singular to working precision. Nevertheless, the solution and error bounds are computed because there are a number of situations where the computed solution can be more accurate than the value of RCOND would suggest.";
2207 throw( numerical_exception( info, msg ) );
2208 }
2209 }
2210 }
2211 break;
2212 }
2213
2214 return( eigen_value );
2215 }
2216 };
2217
2218
2219 #else
2220
2221 // 行列の対角化を行わないバージョン
2222 template < bool b >
2223 struct __eigen__
2224 {
2225 // 実数バージョン
2226 template < class T, class Allocator >
2227 static matrix< T, Allocator >& eigen( matrix< T, Allocator > &a, matrix< T, Allocator > &eigen_value, matrix< T, Allocator > &eigen_vector, matrix_style::style style )
2228 {
2229 typedef __clapack__::integer integer;
2230
2231 if( a.empty( ) )
2232 {
2233 // 行列のサイズが正しくないので例外をスローする
2234 throw( numerical_exception( 1, "Empty matrix." ) );
2235 }
2236
2237 integer info = 0;
2238
2239 switch( style )
2240 {
2241 case matrix_style::sy:
2242 {
2243 // LAPACK関数の引数
2244 integer n = static_cast< integer >( a.cols( ) );
2245 integer lda = static_cast< integer >( a.rows( ) );
2246 typename matrix< T, Allocator >::value_type dmy;
2247 integer lwork = -1;
2248 char jobz[] = "V";
2249 char uplo[] = "U";
2250
2251 // まず最適な作業用配列のサイズを取得する
2252 __clapack__::syev( jobz, uplo, n, NULL, lda, NULL, &dmy, lwork, info );
2253 if( info == 0 )
2254 {
2255 eigen_value.resize( n, 1 );
2256 eigen_vector = a;
2257
2258 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
2259 matrix< T, Allocator > work( lwork, 1 );
2260 __clapack__::syev( jobz, uplo, n, &( eigen_vector[0] ), lda, &( eigen_value[0] ), &( work[0] ), lwork, info );
2261 }
2262
2263 if( info < 0 )
2264 {
2265 // 行列計算が正しく終了しなかったので例外をスローする
2266 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
2267 }
2268 else if( info > 0 )
2269 {
2270 // 行列計算が正しく終了しなかったので例外をスローする
2271 throw( numerical_exception( info, "The algorithm failed to converge. " + __clapack__::to_string( info ) + "off-diagonal elements of an intermediate tridiagonal form did not converge to zero." ) );
2272 }
2273 }
2274 break;
2275
2276 case matrix_style::ge:
2277 default:
2278 {
2279 // LAPACK関数の引数
2280 integer n = static_cast< integer >( a.cols( ) );
2281 integer lda = static_cast< integer >( a.rows( ) );
2282 typename matrix< T, Allocator >::value_type dmy;
2283 integer ldvl = 1;
2284 integer ldvr = n;
2285 integer lwork = -1;
2286 char jobvl[] = "N";
2287 char jobvr[] = "V";
2288
2289 // まず最適な作業用配列のサイズを取得する
2290 __clapack__::geev( jobvl, jobvr, n, NULL, lda, NULL, NULL, NULL, ldvl, NULL, ldvr, &dmy, lwork, info );
2291 if( info == 0 )
2292 {
2293 eigen_value.resize( n, 1 );
2294 matrix< T, Allocator > tmp( n, 1 );
2295 eigen_vector.resize( n, n );
2296
2297 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
2298 matrix< T, Allocator > work( lwork, 1 );
2299 __clapack__::geev( jobvl, jobvr, n, &( a[0] ), lda, &( eigen_value[0] ), &( tmp[0] ),
2300 NULL, ldvl, &( eigen_vector[0] ), ldvr, &( work[0] ), lwork, info );
2301 }
2302
2303 if( info < 0 )
2304 {
2305 // 行列計算が正しく終了しなかったので例外をスローする
2306 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
2307 }
2308 else if( info > 0 )
2309 {
2310 // 行列計算が正しく終了しなかったので例外をスローする
2311 throw( numerical_exception( info, "QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed." ) );
2312 }
2313 }
2314 break;
2315 }
2316
2317 if( info != 0 )
2318 {
2319 // 行列計算が正しく終了しなかったので例外をスローする
2320 throw( numerical_exception( info, "Failed to compute the eigen values and vectors." ) );
2321 }
2322
2323 return( eigen_value );
2324 }
2325 };
2326
2327 template < >
2328 struct __eigen__< true >
2329 {
2330 // 複素数バージョン
2331 template < class T, class Allocator >
2332 static matrix< T, Allocator >& eigen( matrix< T, Allocator > &a, matrix< T, Allocator > &eigen_value, matrix< T, Allocator > &eigen_vector, matrix_style::style style )
2333 {
2334 typedef __clapack__::integer integer;
2335 typedef typename T::value_type value_type;
2336
2337 if( a.empty( ) )
2338 {
2339 // 行列のサイズが正しくないので例外をスローする
2340 throw( numerical_exception( 1, "Empty matrix." ) );
2341 }
2342
2343 integer info = 0;
2344
2345 switch( style )
2346 {
2347 case matrix_style::ge:
2348 default:
2349 {
2350 // LAPACK関数の引数
2351 integer n = static_cast< integer >( a.cols( ) );
2352 integer lda = static_cast< integer >( a.rows( ) );
2353 typename matrix< T, Allocator >::value_type dmy;
2354 integer ldvl = 1;
2355 integer ldvr = n;
2356 integer lwork = -1;
2357 char jobvl[] = "N";
2358 char jobvr[] = "V";
2359
2360 // まず最適な作業用配列のサイズを取得する
2361 __clapack__::geev( jobvl, jobvr, n, NULL, lda, NULL, NULL, ldvl, NULL, ldvr, &dmy, lwork, NULL, info );
2362 if( info == 0 )
2363 {
2364 eigen_value.resize( n, 1 );
2365 eigen_vector.resize( n, n );
2366 value_type *rwork = new value_type[ 2 * n ];
2367
2368 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
2369 matrix< T, Allocator > work( lwork, 1 );
2370 __clapack__::geev( jobvl, jobvr, n, &( a[0] ), lda, &( eigen_value[0] ),
2371 NULL, ldvl, &( eigen_vector[0] ), ldvr, &( work[0] ), lwork, rwork, info );
2372
2373 delete [] rwork;
2374 }
2375
2376 if( info < 0 )
2377 {
2378 // 行列計算が正しく終了しなかったので例外をスローする
2379 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
2380 }
2381 else if( info > 0 )
2382 {
2383 // 行列計算が正しく終了しなかったので例外をスローする
2384 throw( numerical_exception( info, "QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed." ) );
2385 }
2386 }
2387 break;
2388 }
2389
2390 return( eigen_value );
2391 }
2392 };
2393
2394 #endif
2395 }
2396
2397
2398
2399 namespace __svd__
2400 {
2401 // 一般行列の特異値分解を計算する
2402 #if defined(_USE_DIVIDE_AND_CONQUER_SVD_) && _USE_DIVIDE_AND_CONQUER_SVD_ != 0
2403
2404 // 分割統治法を用いるバージョンの特異値分解
2405 template < bool b >
2406 struct __svd__
2407 {
2408 // 実数バージョン
2409 template < class T, class Allocator >
2410 static matrix< T, Allocator >& svd( matrix< T, Allocator > &a, matrix< T, Allocator > &u, matrix< T, Allocator > &s, matrix< T, Allocator > &vt, matrix_style::style style )
2411 {
2412 typedef __clapack__::integer integer;
2413 typedef typename matrix< T, Allocator >::size_type size_type;
2414
2415 if( a.empty( ) )
2416 {
2417 // 行列のサイズが正しくないので例外をスローする
2418 throw( numerical_exception( 1, "Empty matrix." ) );
2419 }
2420
2421 integer info = 0;
2422
2423 switch( style )
2424 {
2425 case matrix_style::ge:
2426 default:
2427 {
2428 // LAPACK関数の引数
2429 integer m = static_cast< integer >( a.rows( ) );
2430 integer n = static_cast< integer >( a.cols( ) );
2431 integer lda = m;
2432 integer ldu = m;
2433 integer size = m < n ? m : n;
2434 typename matrix< T, Allocator >::value_type dmy;
2435 integer ldvt = n;
2436 integer lwork = -1;
2437 char jobz[] = "A";
2438
2439 // まず最適な作業用配列のサイズを取得する
2440 __clapack__::gesdd( jobz, m, n, NULL, lda, NULL, NULL, ldu, NULL, ldvt, &dmy, lwork, NULL, info );
2441 if( info == 0 )
2442 {
2443 u.resize( ldu, m );
2444 matrix< T, Allocator > ss( size, 1 );
2445 vt.resize( ldvt, n );
2446 integer *iwork = new integer[ 8 * size ];
2447
2448 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
2449 matrix< T, Allocator > work( lwork, 1 );
2450 __clapack__::gesdd( jobz, m, n, &( a[0] ), lda, &( ss[0] ), &( u[0] ), ldu, &( vt[0] ), ldvt, &( work[0] ), lwork, iwork, info );
2451
2452 delete [] iwork;
2453
2454 s.resize( m, n );
2455 for( size_type i = 0 ; i < ss.rows( ) ; i++ )
2456 {
2457 s( i, i ) = ss[ i ];
2458 }
2459 }
2460
2461 if( info < 0 )
2462 {
2463 // 行列計算が正しく終了しなかったので例外をスローする
2464 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
2465 }
2466 else if( info > 0 )
2467 {
2468 // 行列計算が正しく終了しなかったので例外をスローする
2469 throw( numerical_exception( info, "DBDSDC did not converge, updating process failed." ) );
2470 }
2471 }
2472 break;
2473 }
2474
2475 return( s );
2476 }
2477
2478 // 実数バージョン
2479 template < class T, class Allocator >
2480 static matrix< T, Allocator >& svd( matrix< T, Allocator > &a, matrix< T, Allocator > &s, matrix< T, Allocator > &vt, matrix_style::style style )
2481 {
2482 typedef __clapack__::integer integer;
2483 typedef typename matrix< T, Allocator >::size_type size_type;
2484
2485 if( a.empty( ) )
2486 {
2487 // 行列のサイズが正しくないので例外をスローする
2488 throw( numerical_exception( 1, "Empty matrix." ) );
2489 }
2490
2491 integer info = 0;
2492
2493 switch( style )
2494 {
2495 case matrix_style::ge:
2496 default:
2497 {
2498 // LAPACK関数の引数
2499 integer m = static_cast< integer >( a.rows( ) );
2500 integer n = static_cast< integer >( a.cols( ) );
2501 integer lda = m;
2502 integer ldu = m;
2503 integer size = m < n ? m : n;
2504 typename matrix< T, Allocator >::value_type dmy;
2505 integer ldvt = n;
2506 integer lwork = -1;
2507 char jobz[] = "O";
2508
2509 // まず最適な作業用配列のサイズを取得する
2510 __clapack__::gesdd( jobz, m, n, NULL, lda, NULL, NULL, ldu, NULL, ldvt, &dmy, lwork, NULL, info );
2511 if( info == 0 )
2512 {
2513 integer *iwork = new integer[ 8 * size ];
2514 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
2515 matrix< T, Allocator > work( lwork, 1 );
2516 s.resize( size, 1 );
2517
2518 if( m >= n )
2519 {
2520 vt.resize( ldvt, n );
2521 __clapack__::gesdd( jobz, m, n, &( a[0] ), lda, &( s[0] ), NULL, ldu, &( vt[0] ), ldvt, &( work[0] ), lwork, iwork, info );
2522 }
2523 else
2524 {
2525 vt.resize( ldu, m );
2526 __clapack__::gesdd( jobz, m, n, &( a[0] ), lda, &( s[0] ), &( vt[0] ), ldu, NULL, ldvt, &( work[0] ), lwork, iwork, info );
2527
2528 vt.resize( ldvt, n );
2529 vt.fill( );
2530 for( integer r = 0 ; r < m ; r++ )
2531 {
2532 for( size_type c = 0 ; c < vt.cols( ) ; c++ )
2533 {
2534 vt( r, c ) = a( r, c );
2535 }
2536 }
2537 }
2538
2539 delete [] iwork;
2540 }
2541
2542 if( info < 0 )
2543 {
2544 // 行列計算が正しく終了しなかったので例外をスローする
2545 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
2546 }
2547 else if( info > 0 )
2548 {
2549 // 行列計算が正しく終了しなかったので例外をスローする
2550 throw( numerical_exception( info, "DBDSDC did not converge, updating process failed." ) );
2551 }
2552 }
2553 break;
2554 }
2555
2556 return( s );
2557 }
2558 };
2559
2560 template < >
2561 struct __svd__< true >
2562 {
2563 // 複素数バージョン
2564 template < class T1, class T2, class Allocator1, class Allocator2 >
2565 static matrix< T2, Allocator2 >& svd( matrix< T1, Allocator1 > &a, matrix< T1, Allocator1 > &u, matrix< T2, Allocator2 > &s, matrix< T1, Allocator1 > &vt, matrix_style::style style )
2566 {
2567 typedef __clapack__::integer integer;
2568 typedef typename matrix< T1, Allocator1 >::size_type size_type;
2569 typedef typename T1::value_type value_type;
2570
2571 if( a.empty( ) )
2572 {
2573 // 行列のサイズが正しくないので例外をスローする
2574 throw( numerical_exception( 1, "Empty matrix." ) );
2575 }
2576
2577 integer info = 0;
2578
2579 switch( style )
2580 {
2581 case matrix_style::ge:
2582 default:
2583 {
2584 // LAPACK関数の引数
2585 integer m = static_cast< integer >( a.rows( ) );
2586 integer n = static_cast< integer >( a.cols( ) );
2587 integer lda = m;
2588 integer ldu = m;
2589 integer size = m < n ? m : n;
2590 typename matrix< T1, Allocator1 >::value_type dmy;
2591 integer ldvt = n;
2592 integer lwork = -1;
2593 char jobz[] = "A";
2594
2595 // まず最適な作業用配列のサイズを取得する
2596 __clapack__::gesdd( jobz, m, n, NULL, lda, NULL, NULL, ldu, NULL, ldvt, &dmy, lwork, NULL, NULL, info );
2597 if( info == 0 )
2598 {
2599 u.resize( ldu, m );
2600 #if defined( __MIST_MSVC__ ) && __MIST_MSVC__ < 7
2601 // VC6ではSTLのアロケータの定義が、標準に準拠していないので、デフォルトで代用する
2602 matrix< value_type > ss( size, 1 );
2603 #else
2604 matrix< value_type, typename Allocator1::template rebind< value_type >::other > ss( size, 1 );
2605 #endif
2606 vt.resize( ldvt, n );
2607 value_type *rwork = new value_type[ 5 * size * size + 5 * size ];
2608 integer *iwork = new integer[ 8 * size ];
2609
2610 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
2611 matrix< T1, Allocator1 > work( lwork, 1 );
2612 __clapack__::gesdd( jobz, m, n, &( a[0] ), lda, &( ss[0] ), &( u[0] ), ldu, &( vt[0] ), ldvt, &( work[0] ), lwork, rwork, iwork, info );
2613
2614 delete [] rwork;
2615 delete [] iwork;
2616
2617 s.resize( m, n );
2618 for( size_type i = 0 ; i < ss.rows( ) ; i++ )
2619 {
2620 s( i, i ) = ss[ i ];
2621 }
2622 }
2623
2624 if( info < 0 )
2625 {
2626 // 行列計算が正しく終了しなかったので例外をスローする
2627 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
2628 }
2629 else if( info > 0 )
2630 {
2631 // 行列計算が正しく終了しなかったので例外をスローする
2632 throw( numerical_exception( info, "DBDSDC did not converge, updating process failed." ) );
2633 }
2634 }
2635 break;
2636 }
2637
2638 return( s );
2639 }
2640
2641
2642 template < class T1, class T2, class Allocator1, class Allocator2 >
2643 static matrix< T2, Allocator2 >& svd( matrix< T1, Allocator1 > &a, matrix< T2, Allocator2 > &s, matrix< T1, Allocator1 > &vt, matrix_style::style style )
2644 {
2645 typedef __clapack__::integer integer;
2646 typedef typename matrix< T1, Allocator1 >::size_type size_type;
2647 typedef typename matrix< T1, Allocator1 >::value_type complex_type;
2648 typedef typename T1::value_type value_type;
2649
2650 if( a.empty( ) )
2651 {
2652 // 行列のサイズが正しくないので例外をスローする
2653 throw( numerical_exception( 1, "Empty matrix." ) );
2654 }
2655
2656 integer info = 0;
2657
2658 switch( style )
2659 {
2660 case matrix_style::ge:
2661 default:
2662 {
2663 // LAPACK関数の引数
2664 integer m = static_cast< integer >( a.rows( ) );
2665 integer n = static_cast< integer >( a.cols( ) );
2666 integer lda = m;
2667 integer ldu = m;
2668 integer size = m < n ? m : n;
2669 typename matrix< T1, Allocator1 >::value_type dmy;
2670 integer ldvt = n;
2671 integer lwork = -1;
2672 char jobz[] = "O";
2673
2674 // まず最適な作業用配列のサイズを取得する
2675 __clapack__::gesdd( jobz, m, n, NULL, lda, NULL, NULL, ldu, NULL, ldvt, &dmy, lwork, NULL, NULL, info );
2676 if( info == 0 )
2677 {
2678 value_type *rwork = new value_type[ 5 * size * size + 5 * size ];
2679 integer *iwork = new integer[ 8 * size ];
2680
2681 #if defined( __MIST_MSVC__ ) && __MIST_MSVC__ < 7
2682 // VC6ではSTLのアロケータの定義が、標準に準拠していないので、デフォルトで代用する
2683 matrix< value_type > ss( size, 1 );
2684 #else
2685 matrix< value_type, typename Allocator1::template rebind< value_type >::other > ss( size, 1 );
2686 #endif
2687
2688 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
2689 matrix< T1, Allocator1 > work( lwork, 1 );
2690 complex_type *dmy = NULL;
2691
2692 if( m >= n )
2693 {
2694 vt.resize( ldvt, n );
2695 __clapack__::gesdd( jobz, m, n, &( a[0] ), lda, &( ss[0] ), NULL, ldu, &( vt[0] ), ldvt, &( work[0] ), lwork, rwork, iwork, info );
2696 }
2697 else
2698 {
2699 vt.resize( ldu, m );
2700 __clapack__::gesdd( jobz, m, n, &( a[0] ), lda, &( ss[0] ), &( vt[0] ), ldu, NULL, ldvt, &( work[0] ), lwork, rwork, iwork, info );
2701
2702 vt.resize( ldvt, n );
2703 vt.fill( );
2704 for( integer r = 0 ; r < m ; r++ )
2705 {
2706 for( size_type c = 0 ; c < vt.cols( ) ; c++ )
2707 {
2708 vt( r, c ) = a( r, c );
2709 }
2710 }
2711 }
2712
2713 delete [] rwork;
2714 delete [] iwork;
2715
2716 s.resize( size, 1 );
2717 for( integer i = 0 ; i < size ; i++ )
2718 {
2719 s[ i ] = ss[ i ];
2720 }
2721 }
2722
2723 if( info < 0 )
2724 {
2725 // 行列計算が正しく終了しなかったので例外をスローする
2726 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
2727 }
2728 else if( info > 0 )
2729 {
2730 // 行列計算が正しく終了しなかったので例外をスローする
2731 throw( numerical_exception( info, "DBDSDC did not converge, updating process failed." ) );
2732 }
2733 }
2734 break;
2735 }
2736
2737 return( s );
2738 }
2739 };
2740
2741 #else
2742
2743 template < bool b >
2744 struct __svd__
2745 {
2746 // 実数バージョン
2747 template < class T, class Allocator >
2748 static matrix< T, Allocator >& svd( matrix< T, Allocator > &a, matrix< T, Allocator > &u, matrix< T, Allocator > &s, matrix< T, Allocator > &vt, matrix_style::style style )
2749 {
2750 typedef __clapack__::integer integer;
2751 typedef typename matrix< T, Allocator >::size_type size_type;
2752
2753 if( a.empty( ) )
2754 {
2755 // 行列のサイズが正しくないので例外をスローする
2756 throw( numerical_exception( 1, "Empty matrix." ) );
2757 }
2758
2759 integer info = 0;
2760
2761 switch( style )
2762 {
2763 case matrix_style::ge:
2764 default:
2765 {
2766 // LAPACK関数の引数
2767 integer m = static_cast< integer >( a.rows( ) );
2768 integer n = static_cast< integer >( a.cols( ) );
2769 integer lda = m;
2770 integer ldu = m;
2771 typename matrix< T, Allocator >::value_type dmy;
2772 integer ldvt = n;
2773 integer lwork = -1;
2774 char jobu[] = "A";
2775 char jobvt[] = "A";
2776
2777 // まず最適な作業用配列のサイズを取得する
2778 __clapack__::gesvd( jobu, jobvt, m, n, NULL, lda, NULL, NULL, ldu, NULL, ldvt, &dmy, lwork, info );
2779 if( info == 0 )
2780 {
2781 u.resize( ldu, m );
2782 matrix< T, Allocator > ss( m < n ? m : n, 1 );
2783 vt.resize( ldvt, n );
2784
2785 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
2786 matrix< T, Allocator > work( lwork, 1 );
2787 __clapack__::gesvd( jobu, jobvt, m, n, &( a[0] ), lda, &( ss[0] ), &( u[0] ), ldu, &( vt[0] ), ldvt, &( work[0] ), lwork, info );
2788
2789 s.resize( m, n );
2790 for( size_type i = 0 ; i < ss.rows( ) ; i++ )
2791 {
2792 s( i, i ) = ss[ i ];
2793 }
2794 }
2795
2796 if( info < 0 )
2797 {
2798 // 行列計算が正しく終了しなかったので例外をスローする
2799 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
2800 }
2801 else if( info > 0 )
2802 {
2803 // 行列計算が正しく終了しなかったので例外をスローする
2804 throw( numerical_exception( info, "DBDSQR did not converge, INFO specifies how many superdiagonals of an intermediate bidiagonal form B did not converge to zero. See the description of WORK above for details." ) );
2805 }
2806 }
2807 break;
2808 }
2809
2810 return( s );
2811 }
2812
2813
2814 template < class T, class Allocator >
2815 static matrix< T, Allocator >& svd( matrix< T, Allocator > &a, matrix< T, Allocator > &s, matrix< T, Allocator > &vt, matrix_style::style style )
2816 {
2817 typedef __clapack__::integer integer;
2818 typedef typename matrix< T, Allocator >::size_type size_type;
2819
2820 if( a.empty( ) )
2821 {
2822 // 行列のサイズが正しくないので例外をスローする
2823 throw( numerical_exception( 1, "Empty matrix." ) );
2824 }
2825
2826 integer info = 0;
2827
2828 switch( style )
2829 {
2830 case matrix_style::ge:
2831 default:
2832 {
2833 // LAPACK関数の引数
2834 integer m = static_cast< integer >( a.rows( ) );
2835 integer n = static_cast< integer >( a.cols( ) );
2836 integer lda = m;
2837 integer ldu = m;
2838 typename matrix< T, Allocator >::value_type dmy;
2839 integer ldvt = n;
2840 integer lwork = -1;
2841 char jobu[] = "N";
2842 char jobvt[] = "A";
2843
2844 // まず最適な作業用配列のサイズを取得する
2845 __clapack__::gesvd( jobu, jobvt, m, n, NULL, lda, NULL, NULL, ldu, NULL, ldvt, &dmy, lwork, info );
2846 if( info == 0 )
2847 {
2848 s.resize( m < n ? m : n, 1 );
2849 vt.resize( ldvt, n );
2850
2851 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
2852 matrix< T, Allocator > work( lwork, 1 );
2853 __clapack__::gesvd( jobu, jobvt, m, n, &( a[0] ), lda, &( s[0] ), NULL, ldu, &( vt[0] ), ldvt, &( work[0] ), lwork, info );
2854 }
2855
2856 if( info < 0 )
2857 {
2858 // 行列計算が正しく終了しなかったので例外をスローする
2859 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
2860 }
2861 else if( info > 0 )
2862 {
2863 // 行列計算が正しく終了しなかったので例外をスローする
2864 throw( numerical_exception( info, "DBDSQR did not converge, INFO specifies how many superdiagonals of an intermediate bidiagonal form B did not converge to zero. See the description of WORK above for details." ) );
2865 }
2866 }
2867 break;
2868 }
2869
2870 return( s );
2871 }
2872 };
2873
2874 template < >
2875 struct __svd__< true >
2876 {
2877 // 複素数バージョン
2878 template < class T1, class T2, class Allocator1, class Allocator2 >
2879 static matrix< T2, Allocator2 >& svd( matrix< T1, Allocator1 > &a, matrix< T1, Allocator1 > &u, matrix< T2, Allocator2 > &s, matrix< T1, Allocator1 > &vt, matrix_style::style style )
2880 {
2881 typedef __clapack__::integer integer;
2882 typedef typename matrix< T1, Allocator1 >::size_type size_type;
2883 typedef typename T1::value_type value_type;
2884
2885 if( a.empty( ) )
2886 {
2887 // 行列のサイズが正しくないので例外をスローする
2888 throw( numerical_exception( 1, "Empty matrix." ) );
2889 }
2890
2891 integer info = 0;
2892
2893 switch( style )
2894 {
2895 case matrix_style::ge:
2896 default:
2897 {
2898 // LAPACK関数の引数
2899 integer m = static_cast< integer >( a.rows( ) );
2900 integer n = static_cast< integer >( a.cols( ) );
2901 integer lda = m;
2902 integer ldu = m;
2903 integer size = m < n ? m : n;
2904 typename matrix< T1, Allocator1 >::value_type dmy;
2905 integer ldvt = n;
2906 integer lwork = -1;
2907 char jobu[] = "A";
2908 char jobvt[] = "A";
2909
2910 // まず最適な作業用配列のサイズを取得する
2911 __clapack__::gesvd( jobu, jobvt, m, n, NULL, lda, NULL, NULL, ldu, NULL, ldvt, &dmy, lwork, NULL, info );
2912 if( info == 0 )
2913 {
2914 u.resize( ldu, m );
2915 #if defined( __MIST_MSVC__ ) && __MIST_MSVC__ < 7
2916 // VC6ではSTLのアロケータの定義が、標準に準拠していないので、デフォルトで代用する
2917 matrix< value_type > ss( size, 1 );
2918 #else
2919 matrix< value_type, typename Allocator1::template rebind< value_type >::other > ss( size, 1 );
2920 #endif
2921 vt.resize( ldvt, n );
2922 value_type *rwork = new value_type[ 5 * size ];
2923
2924 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
2925 matrix< T1, Allocator1 > work( lwork, 1 );
2926 __clapack__::gesvd( jobu, jobvt, m, n, &( a[0] ), lda, &( ss[0] ), &( u[0] ), ldu, &( vt[0] ), ldvt, &( work[0] ), lwork, rwork, info );
2927
2928 s.resize( m, n );
2929 for( size_type i = 0 ; i < ss.rows( ) ; i++ )
2930 {
2931 s( i, i ) = ss[ i ];
2932 }
2933
2934 delete [] rwork;
2935 }
2936
2937 if( info < 0 )
2938 {
2939 // 行列計算が正しく終了しなかったので例外をスローする
2940 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
2941 }
2942 else if( info > 0 )
2943 {
2944 // 行列計算が正しく終了しなかったので例外をスローする
2945 throw( numerical_exception( info, "ZBDSQR did not converge, INFO specifies how many superdiagonals of an intermediate bidiagonal form B did not converge to zero. See the description of WORK above for details." ) );
2946 }
2947 }
2948 break;
2949 }
2950
2951 return( s );
2952 }
2953
2954
2955 template < class T1, class T2, class Allocator1, class Allocator2 >
2956 static matrix< T2, Allocator2 >& svd( matrix< T1, Allocator1 > &a, matrix< T2, Allocator2 > &s, matrix< T1, Allocator1 > &vt, matrix_style::style style )
2957 {
2958 typedef __clapack__::integer integer;
2959 typedef typename matrix< T1, Allocator1 >::size_type size_type;
2960 typedef typename T1::value_type value_type;
2961
2962 if( a.empty( ) )
2963 {
2964 // 行列のサイズが正しくないので例外をスローする
2965 throw( numerical_exception( 1, "Empty matrix." ) );
2966 }
2967
2968 integer info = 0;
2969
2970 switch( style )
2971 {
2972 case matrix_style::ge:
2973 default:
2974 {
2975 // LAPACK関数の引数
2976 integer m = static_cast< integer >( a.rows( ) );
2977 integer n = static_cast< integer >( a.cols( ) );
2978 integer lda = m;
2979 integer ldu = m;
2980 integer size = m < n ? m : n;
2981 typename matrix< T1, Allocator1 >::value_type dmy;
2982 integer ldvt = n;
2983 integer lwork = -1;
2984 char jobu[] = "N";
2985 char jobvt[] = "A";
2986
2987 // まず最適な作業用配列のサイズを取得する
2988 __clapack__::gesvd( jobu, jobvt, m, n, NULL, lda, NULL, NULL, ldu, NULL, ldvt, &dmy, lwork, NULL, info );
2989 if( info == 0 )
2990 {
2991 #if defined( __MIST_MSVC__ ) && __MIST_MSVC__ < 7
2992 // VC6ではSTLのアロケータの定義が、標準に準拠していないので、デフォルトで代用する
2993 s.resize( size, 1 );
2994 #else
2995 matrix< value_type, typename Allocator1::template rebind< value_type >::other > ss( size, 1 );
2996 #endif
2997 vt.resize( ldvt, n );
2998 value_type *rwork = new value_type[ 5 * size ];
2999
3000 lwork = static_cast< integer >( __clapack__::get_real( dmy ) );
3001 matrix< T1, Allocator1 > work( lwork, 1 );
3002 __clapack__::gesvd( jobu, jobvt, m, n, &( a[0] ), lda, &( s[0] ), NULL, ldu, &( vt[0] ), ldvt, &( work[0] ), lwork, rwork, info );
3003
3004 delete [] rwork;
3005 }
3006
3007 if( info < 0 )
3008 {
3009 // 行列計算が正しく終了しなかったので例外をスローする
3010 throw( numerical_exception( info, __clapack__::to_string( info ) + "-th argument had an illegal value" ) );
3011 }
3012 else if( info > 0 )
3013 {
3014 // 行列計算が正しく終了しなかったので例外をスローする
3015 throw( numerical_exception( info, "ZBDSQR did not converge, INFO specifies how many superdiagonals of an intermediate bidiagonal form B did not converge to zero. See the description of WORK above for details." ) );
3016 }
3017 }
3018 break;
3019 }
3020
3021 return( s );
3022 }
3023 };
3024
3025 #endif
3026 }
3027
3028
3029
3032
3033
3034
3049 template < class T, class Allocator >
3050 inline bool multiply( const matrix< T, Allocator > &a, const matrix< T, Allocator > &b, matrix< T, Allocator > &c, bool a_is_transpose, bool b_is_transpose, typename matrix< T, Allocator >::value_type alpha, typename matrix< T, Allocator >::value_type beta )
3051 {
3052 typedef __clapack__::integer integer;
3053 typedef typename matrix< T, Allocator >::value_type value_type;
3054
3055 if( a.empty( ) || b.empty( ) )
3056 {
3057 // 行列のサイズが正しくないので例外をスローする
3058 return( false );
3059 }
3060
3061 // LAPACK関数の引数
3062 integer m = static_cast< integer >( a_is_transpose ? a.cols( ) : a.rows( ) );
3063 integer n = static_cast< integer >( b_is_transpose ? b.rows( ) : b.cols( ) );
3064 integer k = static_cast< integer >( a_is_transpose ? a.rows( ) : a.cols( ) );
3065 integer k_ = static_cast< integer >( b_is_transpose ? b.cols( ) : b.rows( ) );
3066 integer lda = a_is_transpose ? k : m;
3067 integer ldb = b_is_transpose ? n : k;
3068 integer ldc = m;
3069 char transa[ 2 ];
3070 char transb[ 2 ];
3071
3072 if( k != k_ )
3073 {
3074 // 行列のサイズが正しくないので例外をスローする
3075 return( false );
3076 }
3077
3078 transa[ 0 ] = a_is_transpose ? 'T' : 'N';
3079 transa[ 1 ] = '\0';
3080 transb[ 0 ] = b_is_transpose ? 'T' : 'N';
3081 transb[ 1 ] = '\0';
3082
3083 c.resize( m, n );
3084 c.fill( );
3085
3086 // BLASルーチンでは，入力行列AとBの内容は変化しないが，
3087 // インターフェースは const 修飾を受けていないのでキャストを行う
3088 __clapack__::gemm( transa, transb, m, n, k, alpha, const_cast< value_type * >( &( a[ 0 ] ) ), lda, const_cast< value_type * >( &( b[ 0 ] ) ), ldb, beta, &( c[ 0 ] ), ldc );
3089
3090 return( true );
3091 }
3092
3093
3094
3107 template < class T, class Allocator >
3108 inline bool multiply( const matrix< T, Allocator > &a, const matrix< T, Allocator > &b, matrix< T, Allocator > &c, bool a_is_transpose = false, bool b_is_transpose = false )
3109 {
3110 return( multiply( a, b, c, a_is_transpose, b_is_transpose, 1, 0 ) );
3111 }
3112
3113
3114
3122 template < class T1, class T2, class Allocator1, class Allocator2 >
3123 inline void permutation_matrix( const matrix< T1, Allocator1 > &pivot, matrix< T2, Allocator2 > &out )
3124 {
3125 typedef typename matrix< T1, Allocator1 >::size_type size_type;
3126 typedef typename matrix< T1, Allocator1 >::value_type value_type;
3127
3128 size_type i;
3129 matrix< T1, Allocator1 > p( pivot );
3130 for( i = 0 ; i < pivot.size( ) ; i++ )
3131 {
3132 p[ i ] = static_cast< value_type >( i );
3133 }
3134
3135 // ピボットでデータを入れ替える
3136 for( i = 0 ; i < pivot.size( ) ; i++ )
3137 {
3138 value_type tmp = p[ i ];
3139 p[ i ] = p[ pivot[ i ] - 1 ];
3140 p[ pivot[ i ] - 1 ] = tmp;
3141 }
3142
3143 out.resize( pivot.size( ), pivot.size( ) );
3144 out.fill( );
3145
3146 // ピボット行列を作成する
3147 for( i = 0 ; i < pivot.size( ) ; i++ )
3148 {
3149 out( p[ i ], i ) = 1;
3150 }
3151 }
3152
3153
3162 template < class T, class Allocator >
3163 inline const typename matrix< T, Allocator >::value_type trace( const matrix< T, Allocator > &a )
3164 {
3165 typedef typename matrix< T, Allocator >::size_type size_type;
3166 typedef typename matrix< T, Allocator >::value_type value_type;
3167 value_type v = value_type( 0 );
3168 size_type size = a.rows( ) < a.cols( ) ? a.rows( ) : a.cols( );
3169 for( size_type i = 0 ; i < size ; ++i )
3170 {
3171 v += a( i, i );
3172 }
3173 return( v );
3174 }
3175
3176
3208 template < class T, class Allocator >
3209 inline const typename matrix< T, Allocator >::value_type det( const matrix< T, Allocator > &a, matrix_style::style style = matrix_style::ge )
3210 {
3211 typedef typename matrix< T, Allocator >::size_type size_type;
3212 typedef typename matrix< T, Allocator >::value_type value_type;
3213
3214 if( a.empty( ) || a.rows( ) != a.cols( ) )
3215 {
3216 return( value_type( 0 ) );
3217 }
3218
3219 switch( a.rows( ) )
3220 {
3221 case 1:
3222 return( a( 0, 0 ) );
3223 break;
3224
3225 case 2:
3226 return( a( 0, 0 ) * a( 1, 1 ) - a( 0, 1 ) * a( 1, 0 ) );
3227 break;
3228
3229 case 3:
3230 return( a( 0, 0 ) * a( 1, 1 ) * a( 2, 2 ) + a( 0, 1 ) * a( 1, 2 ) * a( 2, 0 ) + a( 0, 2 ) * a( 1, 0 ) * a( 2, 1 )
3231 - a( 0, 2 ) * a( 1, 1 ) * a( 2, 0 ) - a( 0, 1 ) * a( 1, 0 ) * a( 2, 2 ) - a( 0, 0 ) * a( 1, 2 ) * a( 2, 1 ) );
3232 break;
3233
3234 default:
3235 {
3236 #if defined( __MIST_MSVC__ ) && __MIST_MSVC__ < 7
3237 // VC6ではSTLのアロケータの定義が、標準に準拠していないので、デフォルトで代用する
3238 matrix< __clapack__::integer > pivot;
3239 #else
3240 matrix< __clapack__::integer, typename Allocator::template rebind< __clapack__::integer >::other > pivot;
3241 #endif
3242
3243 try
3244 {
3245 matrix< T, Allocator > m = lu_factorization( a, pivot, style );
3246
3247 value_type v = m( 0, 0 );
3248 size_type size = a.rows( ) < a.cols( ) ? a.rows( ) : a.cols( );
3249 size_type count = 0, i;
3250
3251 // LU分解時に行の入れ替えが行われた回数を計算する
3252 for( i = 0 ; i < pivot.size( ) ; ++i )
3253 {
3254 count += static_cast< size_type >( pivot[ i ] ) != i + 1 ? 1 : 0;
3255 }
3256
3257 // 対角成分の積を計算する
3258 for( i = 1 ; i < size ; ++i )
3259 {
3260 v *= m( i, i );
3261 }
3262
3263 if( count % 2 == 0 )
3264 {
3265 return( v );
3266 }
3267 else
3268 {
3269 return( -v );
3270 }
3271 }
3272 catch( const numerical_exception &e )
3273 {
3274 // error が負の場合は根本的に計算不能．
3275 if( e.error < 0 )
3276 {
3277 throw( e );
3278 }
3279 else
3280 {
3281 return( 0 );
3282 }
3283 }
3284 }
3285 break;
3286 }
3287 }
3288
3289
3290
3301 template < class T, class Allocator >
3302 inline const matrix< T, Allocator >& solve( const matrix< T, Allocator > &a, matrix< T, Allocator > &b, matrix_style::style style = matrix_style::ge )
3303 {
3304 try
3305 {
3306 matrix< T, Allocator > a_( a );
3307 return( __solve__::__solve__< __numeric__::is_complex< T >::value >::solve( a_, b, style ) );
3308 }
3309 catch( const numerical_exception &e )
3310 {
3311 throw( e );
3312 }
3313 }
3314
3326 //template < class T, class Allocator >
3327 //inline const matrix< T, Allocator >& solve_band_matrix( const matrix< T, Allocator > &a, matrix< T, Allocator > &b )
3328 //{
3329 // try
3330 // {
3331 // matrix< T, Allocator > a_( a );
3332 // return( __solve__::__solve_band_matrix__< __numeric__::is_complex< T >::value >::solve_band_matrix( a_, b ) );
3333 // }
3334 // catch( const numerical_exception &e )
3335 // {
3336 // throw( e );
3337 // }
3338 //}
3339 //>>> yosidah
3340
3341
3350 template < class T, class Allocator1, class Allocator2 >
3351 const matrix< T, Allocator1 > lu_factorization( const matrix< T, Allocator1 > &a, matrix< __clapack__::integer, Allocator2 > &pivot, matrix_style::style style = matrix_style::ge )
3352 {
3353 try
3354 {
3355 matrix< T, Allocator1 > a_( a );
3356 return( __lu__::__lu__< __numeric__::is_complex< T >::value >::lu_factorization( a_, pivot, style ) );
3357 }
3358 catch( const numerical_exception &e )
3359 {
3360 throw( e );
3361 }
3362 }
3363
3364
3377 template < class T, class Allocator >
3378 bool lu_factorization( const matrix< T, Allocator > &a, matrix< T, Allocator > &L, matrix< T, Allocator > &U, matrix< T, Allocator > &pivot, matrix_style::style style = matrix_style::ge )
3379 {
3380 #if defined( __MIST_MSVC__ ) && __MIST_MSVC__ < 7
3381 // VC6ではSTLのアロケータの定義が、標準に準拠していないので、デフォルトで代用する
3382 matrix< __clapack__::integer > piv( a.cols( ), 1 );
3383 #else
3384 matrix< __clapack__::integer, typename Allocator::template rebind< __clapack__::integer >::other > piv( a.cols( ), 1 );
3385 #endif
3386
3387 try
3388 {
3389 L = a;
3390 __lu__::__lu__< __numeric__::is_complex< T >::value >::lu_factorization( L, piv, style );
3391
3392 typedef typename matrix< T, Allocator >::size_type size_type;
3393
3394 U.resize( L.rows( ), L.cols( ) );
3395 for( size_type r = 0 ; r < a.rows( ) ; r++ )
3396 {
3397 for( size_type c = r ; c < a.cols( ) ; c++ )
3398 {
3399 U( r, c ) = L( r, c );
3400 L( r, c ) = 0;
3401 }
3402 L( r, r ) = 1;
3403 }
3404
3405 permutation_matrix( piv, pivot );
3406 }
3407 catch( const numerical_exception &e )
3408 {
3409 throw( e );
3410 }
3411
3412 return( true );
3413 }
3414
3415
3423 template < class T, class Allocator >
3424 const matrix< T, Allocator > lu_factorization( const matrix< T, Allocator > &a, matrix_style::style style = matrix_style::ge )
3425 {
3426 #if defined( __MIST_MSVC__ ) && __MIST_MSVC__ < 7
3427 // VC6ではSTLのアロケータの定義が、標準に準拠していないので、デフォルトで代用する
3428 matrix< __clapack__::integer > pivot( a.cols( ), 1 );
3429 #else
3430 matrix< __clapack__::integer, typename Allocator::template rebind< __clapack__::integer >::other > pivot( a.cols( ), 1 );
3431 #endif
3432
3433 try
3434 {
3435 matrix< T, Allocator > a_( a );
3436 return( __lu__::__lu__< __numeric__::is_complex< T >::value >::lu_factorization( a_, pivot, style ) );
3437 }
3438 catch( const numerical_exception &e )
3439 {
3440 throw( e );
3441 }
3442 }
3443
3444
3445
3453 template < class T, class Allocator >
3454 const matrix< T, Allocator > cholesky_factorization( const matrix< T, Allocator > &a, matrix_style::style style = matrix_style::sy )
3455 {
3456 try
3457 {
3458 matrix< T, Allocator > a_( a );
3459 return( __cholesky__::cholesky_factorization( a_, style ) );
3460 }
3461 catch( const numerical_exception &e )
3462 {
3463 throw( e );
3464 }
3465 }
3466
3467
3468
3478 template < class T, class Allocator >
3479 void qr_factorization( const matrix< T, Allocator > &a, matrix< T, Allocator > &Q, matrix< T, Allocator > &R, matrix_style::style style = matrix_style::ge )
3480 {
3481 try
3482 {
3483 matrix< T, Allocator > a_( a );
3484 __qr__::__qr__< __numeric__::is_complex< T >::value >::qr_factorization( a_, Q, R, style );
3485 }
3486 catch( const numerical_exception &e )
3487 {
3488 throw( e );
3489 }
3490 }
3491
3492
3493
3505 template < class T, class Allocator >
3506 matrix< T, Allocator > inverse( const matrix< T, Allocator > &a, matrix_style::style style = matrix_style::ge )
3507 {
3508 typedef matrix< T, Allocator > matrix_type;
3509 typedef typename matrix_type::size_type size_type;
3510
3511 try
3512 {
3513 if( a.rows( ) != a.cols( ) )
3514 {
3515 // 特異値分解を用いて一般化逆行列を計算する
3516 matrix_type u, s, vt;
3517 svd( a, u, s, vt, style );
3518 size_type num = s.rows( ) < s.cols( ) ? s.rows( ) : s.cols( );
3519 for( size_type i = 0 ; i < num ; i++ )
3520 {
3521 if( __clapack__::get_real( s( i, i ) ) != 0 )
3522 {
3523 s( i, i ) = 1 / __clapack__::get_real( s( i, i ) );
3524 }
3525 }
3526 return( ( u * s * vt ).t( ) );
3527 }
3528 else
3529 {
3530 matrix_type a_( a );
3531 return( __inverse__::__inverse__< __numeric__::is_complex< T >::value >::inverse( a_, style ) );
3532 }
3533 }
3534 catch( const numerical_exception &e )
3535 {
3536 throw( e );
3537 }
3538 }
3539
3540
3541
3555 template < class T, class Allocator >
3556 const matrix< T, Allocator >& eigen( const matrix< T, Allocator > &a, matrix< T, Allocator > &eigen_value, matrix< T, Allocator > &eigen_vector, matrix_style::style style = matrix_style::ge )
3557 {
3558 typedef typename matrix< T, Allocator >::size_type size_type;
3559
3560 matrix< T, Allocator > a_( a );
3561
3562 try
3563 {
3564 __eigen__::__eigen__< __numeric__::is_complex< T >::value >::eigen( a_, eigen_value, eigen_vector, style );
3565 }
3566 catch( const numerical_exception &e )
3567 {
3568 throw( e );
3569 }
3570
3571 typedef __clapack__::__value_index_pair__< T > value_index_pair;
3572
3573 std::vector< value_index_pair > vips( eigen_value.size( ) );
3574
3575 for( size_type i = 0 ; i < eigen_value.size( ) ; i++ )
3576 {
3577 value_index_pair &vi = vips[ i ];
3578 vi.value = eigen_value[ i ];
3579 vi.index = i;
3580 }
3581
3582
3583 #if defined( _DESCENDING_ORDER_EIGEN_VALUE_ ) && _DESCENDING_ORDER_EIGEN_VALUE_ == 1
3584 // 固有値と固有ベクトルを降順に並び替える
3585 std::sort( vips.begin( ), vips.end( ), value_index_pair::descending );
3586 #else
3587 // 固有値と固有ベクトルを昇順に並び替える
3588 std::sort( vips.begin( ), vips.end( ), value_index_pair::ascending );
3589 #endif
3590
3591 a_ = eigen_vector;
3592
3593 for( size_type i = 0 ; i < vips.size( ) ; i++ )
3594 {
3595 const value_index_pair &vi = vips[ i ];
3596 eigen_value[ i ] = vi.value;
3597
3598 for( size_type j = 0 ; j < eigen_vector.rows( ) ; j++ )
3599 {
3600 eigen_vector( j, i ) = a_( j, vi.index );
3601 }
3602 }
3603
3604 return( eigen_value );
3605 }
3606
3607
3622 template < class T1, class T2, class Allocator1, class Allocator2 >
3623 const matrix< T2, Allocator2 >& svd( const matrix< T1, Allocator1 > &a, matrix< T1, Allocator1 > &u, matrix< T2, Allocator2 > &s, matrix< T1, Allocator1 > &vt, matrix_style::style style = matrix_style::ge )
3624 {
3625 try
3626 {
3627 matrix< T1, Allocator1 > a_( a );
3628 return( __svd__::__svd__< __numeric__::is_complex< T1 >::value >::svd( a_, u, s, vt, style ) );
3629 }
3630 catch( const numerical_exception &e )
3631 {
3632 throw( e );
3633 }
3634 }
3635
3636
3652 template < class T1, class T2, class Allocator1, class Allocator2 >
3653 const matrix< T2, Allocator2 >& svd( const matrix< T1, Allocator1 > &a, matrix< T2, Allocator2 > &s, matrix< T1, Allocator1 > &vt, matrix_style::style style = matrix_style::ge )
3654 {
3655 try
3656 {
3657 matrix< T1, Allocator1 > a_( a );
3658 return( __svd__::__svd__< __numeric__::is_complex< T1 >::value >::svd( a_, s, vt, style ) );
3659 }
3660 catch( const numerical_exception &e )
3661 {
3662 throw( e );
3663 }
3664 }
3665
3666
3668 // 行列演算グループの終わり
3669
3670
3671 // mist名前空間の終わり
3672 _MIST_END
3673
3674
3675 #endif // __INCLUDE_NUMERIC__
3676