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 
Generated on Wed Nov 12 2014 19:44:21 for MIST by doxygen 1.8.1.2