Verified Computation for Partial differential equations

行列クラスの基礎

行列クラスの基礎では

--近似計算--
・double型
・kv::dd型
・kv::mpfr<N>型
・double型のBLASとLapackによる高速化

--精度保証付き計算--
・kv::interval< double >型
・kv::interval< kv::dd >型
・kv::interval< kv::mpfr< N > >型
・kv::interval< double >型のBLASとLapackによる近似計算部分の高速化(BLASは丸めモードに対応していなくて良い)
・kv::interval< double >型のBLASとLapackによる高速化(BLASが丸めモードに対応している必要あり)

以上の型を行列クラスのテンプレートに渡したときに使用できる関数を紹介します。 見方が少し特殊なので以下の手順で見て下さい。
1. まずはdouble型を見て下さい。使用できる関数がそのままソースコードにもなっているので、御試して下さい。
2. double型の説明には(1)などの数字が書かれています。
3. 他の型については場合、double型との差分が書かれています。
4. 他の型の数字の部分をdouble型の数字の部分と置き換えると使用できる関数一覧とソースコードにもなります。

double型

//(1)インクルードファイル

#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void){

//(2)型

//vcp::matrix< double >
//あるいは
//vcp::matrix< double, vcp::mats< double > >
//同じ意味になります。

//(3)変数の宣言

vcp::matrix< double > A;
vcp::matrix< double > B, C, D;
vcp::matrix< double > x, b;

//(4)行列の初期化

int n = 10;
int m = 5;
A.zeros(n); // n×nのゼロ行列を作成
A.zeros(n, m); // n×mのゼロ行列を作成
A.ones(n); // n×nの全要素1の行列を作成
A.ones(n, m); // n×mの全要素1の行列を作成
A.rand(n); // n×nのランダム行列を作成
A.rand(n, m); // n×mのランダム行列を作成
A.eye(n); // n×nの単位行列を作成
B.eye(n);

//(5)行列サイズの取得

int row = A.rowsize(); // 行列Aの行数を戻り値intで返す
int column = A.columnsize(); // 行列Aの列数を戻り値intで返す

//(6)行列の要素へのアクセス

A(0,0) = 10; // 0<= x < A.rowsize(), 0 <= y < A.columnsize, A(x,y)にアクセス可
A(5,3) = A(0,0);

//(7)行列の四則計算

C = 1 + A; // 全要素に1を加えてCに代入
C = A + 2; // 全要素に2を加えてCに代入
C = A + B; // 行列同士の和(サイズが一致していないとエラー)
C += A; // C = C + A, +=のほうがメモリ的に優位

C = 1 - A; // 全要素に1を引いてCに代入
C = A - 2; // 全要素に2を引いてCに代入
C = A - B; // 行列同士の差(サイズが一致していないとエラー)
C -= A; // C = C - A, -=のほうがメモリ的に優位

C = 1 * A; // 全要素に1を掛けてCに代入
C = A * 2; // 全要素に2を掛けてCに代入
C = A * B; // 内積・行列積(A:n×m, B:m×l, n,m,l >= 1 ならOK)
C *= A; // C = C * A, *=のメモリ的な優位はない...
C = ltransmul(A); // C = transpose(A)*A

A = A + 1;
C = 1 / A; // 全要素に1を割ってCに代入
C = A / 2; // 全要素を2で割ってCに代入

//(8)行列の数学関数

C = abs(A); // 全要素の絶対値を取ってCに代入
C = sqrt(A); // 全要素の平方根を取ってCに代入
C = sin(A); // 全要素の正弦を取ってCに代入
C = cos(A); // 全要素の余弦を取ってCに代入
C = exp(A); // 全要素の指数関数を取ってCに代入
C = log(A); // 全要素の対数関数を取ってCに代入

//(9)Matlab like関数

C = sum(A); // Aが行列n×m -> 各列の和からなる行ベクトル1×m, Aがベクトル -> 全成分の和をCに代入
C = diag(A); // Aが行列 -> 対角成分を並べた列ベクトル, Aがベクトル-> ベクトルを対角成分に並べた行列
C = transpose(A); // Aの転置行列をCに代入
C = max(A); // Aが行列 -> 各列の最大要素を含む行ベクトル1×m, Aがベクトル -> 全成分の最大値をCに代入
C = min(A); // Aが行列 -> 各列の最小要素を含む行ベクトル1×m, Aがベクトル -> 全成分の最小値をCに代入
C = normone(A); // ベクトル・行列のどちらにも対応した1ノルムを計算
C = norminf(A); // ベクトル・行列のどちらにも対応した最大値ノルムを計算

//(10)連立一次方程式 Ax = b

A.rand(n);
b.ones(n,m);
b = A*b;
x = lss(A,b); // Ax = bを満たすxを求める, Aはn×n, bはn×m, n, m >= 1

//(11)逆行列

C = inv(A); // Aの逆行列をCに代入

//(12)対称固有値問題 Ax = lambda x (固有値のみ)

A = transpose(A) + A;
eigsym(A, C); // Ax = lambda xを満たすlambdaを対角成分に持つ行列Cを返す

//(13)対称固有値問題 Ax = lambda x (固有値と固有ベクトル)

eigsym(A, C, D); // Ax = lambda xに対し行列Cの対角成分に固有値，行列Dに固有ベクトル

//(14)一般化対称固有値問題 Ax = lambda B x (B:正定値対称行列) (固有値のみ)

B.rand(n);
A = transpose(A) + A;
B = ltransmul(B);
eigsymge(A, B, C); // Ax = lambda B xを満たすlambdaを対角成分に持つ行列Cを返す

//(15)一般化対称固有値問題 Ax = lambda B x (B:正定値対称行列) (固有値と固有ベクトル)

eigsymge(A, B, C, D); // Ax = lambda B xに対し行列Cの対角成分に固有値，行列Dに固有ベクトル

//(16)行列のディスプレイ出力

std::cout << "Matrix C = \n" << C << std::endl; //cout, endlでそのまま使えます

//(17)行列のバイナリファイルへの入出力

save(C,"filename"); // 行列Cをfilenameという名前でバイナリファイルへ保存
load(B,"filename"); // filenameという名前のバイナリファイルから行列Bへ読み込み
std::cout << "Matrix B = \n" << B << std::endl;

//(18)行列の結合

D = vercat(A,B,C); // Matlab関数 D = [A;B;C] : 引数の数は任意
D = horzcat(A,B,C); // Matlab関数 D = [A,B,C] : 引数の数は任意
C.rand(n);
A = vercat(D, horzcat(A,B,C)); // Matlab関数 A = [D; [A,B,C]]; : 引数の数は任意

//(19)行列の解放

A.clear(); // 行列Aを解放

return 0;
}

kv::dd型 (条件 : kvライブラリ)

//(1)インクルードファイル

#include <kv/dd.hpp>
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void){

//(2)型

//vcp::matrix< kv::dd >
//あるいは
//vcp::matrix< kv::dd, vcp::mats< kv::dd > >
//同じ意味になります。

//(3)変数の宣言

vcp::matrix< kv::dd > A;
vcp::matrix< kv::dd > B, C, D;
vcp::matrix< kv::dd > x, b;

kv::mpfr<N>型 (条件 : kvライブラリ)

//(1)インクルードファイル

#include <kv/mpfr.hpp>
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void){

//(2)型

//vcp::matrix< kv::mpfr< 300 > >
//あるいは
//vcp::matrix< kv::mpfr< 300 >, vcp::mats< kv::mpfr< 300 > > >
//同じ意味になります。

//(3)変数の宣言

vcp::matrix< kv::mpfr< 300 > > A;
vcp::matrix< kv::mpfr< 300 > > B, C, D;
vcp::matrix< kv::mpfr< 300 > > x, b;

//(17)行列のバイナリファイルへの入出力は使えません...

//save(C,"filename"); // 行列Cをfilenameという名前でバイナリファイルへ保存
//load(B,"filename"); // filenameという名前のバイナリファイルから行列Bへ読み込み
//std::cout << "Matrix B = \n" << B << std::endl;

double型 + BLAS (条件 : BLASとLapack (各スレッドの丸めモード変更の有無は問わない))

//(1)インクルードファイル

#include <vcp/pdblas.hpp>
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void){

//(2)型

//vcp::matrix< double, vcp::pdblas >

//(3)変数の宣言

vcp::matrix< double, vcp::pdblas > A;
vcp::matrix< double, vcp::pdblas > B, C, D;
vcp::matrix< double, vcp::pdblas > x, b;

kv::interval< double >型 (条件 : kvライブラリ)

//(1)インクルードファイル

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
#include <vcp/imats.hpp>
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void){

//(2)型

//vcp::matrix< kv::interval< double >, vcp::imats< double > >

//(3)変数の宣言

vcp::matrix< kv::interval< double >, vcp::imats< double > > A;
vcp::matrix< kv::interval< double >, vcp::imats< double > > B, C, D;
vcp::matrix< kv::interval< double >, vcp::imats< double > > x, b;

//(11)逆行列は区間Gaussが発動するので使わないで下さい

//C = inv(A); // Aの逆行列をCに代入

//(13)対称固有値問題 Ax = lambda x (固有値と固有ベクトル)は使えません...

//eigsym(A, C, D); // Ax = lambda xに対し行列Cの対角成分に固有値，行列Dに固有ベクトル

//(15)一般化対称固有値問題 Ax = lambda B x (固有値と固有ベクトル)は使えません...

//eigsymge(A, B, C, D); // Ax = lambda B xに対し行列Cの対角成分に固有値，行列Dに固有ベクトル

kv::interval< kv::dd >型 (条件 : kvライブラリ)

//(1)インクルードファイル

#include <kv/interval.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>
#include <vcp/imats.hpp>
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void){

//(2)型

//vcp::matrix< kv::interval< kv::dd >, vcp::imats< kv::dd > >

//(3)変数の宣言

vcp::matrix< kv::interval< kv::dd >, vcp::imats< kv::dd > > A;
vcp::matrix< kv::interval< kv::dd >, vcp::imats< kv::dd > > B, C, D;
vcp::matrix< kv::interval< kv::dd >, vcp::imats< kv::dd > > x, b;

//(11)逆行列は区間Gaussが発動するので使わないで下さい

//C = inv(A); // Aの逆行列をCに代入

//(13)対称固有値問題 Ax = lambda x (固有値と固有ベクトル)は使えません...

//eigsym(A, C, D); // Ax = lambda xに対し行列Cの対角成分に固有値，行列Dに固有ベクトル

//(15)一般化対称固有値問題 Ax = lambda B x (固有値と固有ベクトル)は使えません...

//eigsymge(A, B, C, D); // Ax = lambda B xに対し行列Cの対角成分に固有値，行列Dに固有ベクトル

kv::interval< kv::mpfr< N > >型 (条件 : kvライブラリ)

//(1)インクルードファイル

#include <kv/interval.hpp>
#include <kv/mpfr.hpp>
#include <kv/rmpfr.hpp>
#include <vcp/imats.hpp>
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void){

//(2)型

//vcp::matrix< kv::interval< kv::mpfr< N > >, vcp::imats< kv::mpfr< N > > >

//(3)変数の宣言

vcp::matrix< kv::interval< kv::mpfr< 300 > >, vcp::imats< kv::mpfr< 300 > > > A;
vcp::matrix< kv::interval< kv::mpfr< 300 > >, vcp::imats< kv::mpfr< 300 > > > B, C, D;
vcp::matrix< kv::interval< kv::mpfr< 300 > >, vcp::imats< kv::mpfr< 300 > > > x, b;

//(11)逆行列は区間Gaussが発動するので使わないで下さい

//C = inv(A); // Aの逆行列をCに代入

//(13)対称固有値問題 Ax = lambda x (固有値と固有ベクトル)は使えません...

//eigsym(A, C, D); // Ax = lambda xに対し行列Cの対角成分に固有値，行列Dに固有ベクトル

//(15)一般化対称固有値問題 Ax = lambda B x (固有値と固有ベクトル)は使えません...

//eigsymge(A, B, C, D); // Ax = lambda B xに対し行列Cの対角成分に固有値，行列Dに固有ベクトル

//(17)行列のバイナリファイルへの入出力は使えません...

//save(C,"filename"); // 行列Cをfilenameという名前でバイナリファイルへ保存
//load(B,"filename"); // filenameという名前のバイナリファイルから行列Bへ読み込み
//std::cout << "Matrix B = \n" << B << std::endl;

kv::interval< double >型 + BLAS 1 (条件 : kvライブラリ、Lapack、(丸めモードを要求しない)BLAS)

//(1)インクルードファイル

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
#include <vcp/pdblas.hpp>
#include <vcp/imats.hpp>
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void){

//(2)型

//vcp::matrix< kv::interval< double >, vcp::imats< double, vcp::pdblas > >

//(3)変数の宣言

vcp::matrix< kv::interval< double >, vcp::imats< double, vcp::pdblas > > A;
vcp::matrix< kv::interval< double >, vcp::imats< double, vcp::pdblas > > B, C, D;
vcp::matrix< kv::interval< double >, vcp::imats< double, vcp::pdblas > > x, b;

//(11)逆行列は区間Gaussが発動するので使わないで下さい

//C = inv(A); // Aの逆行列をCに代入

//(13)対称固有値問題 Ax = lambda x (固有値と固有ベクトル)は使えません...

//eigsym(A, C, D); // Ax = lambda xに対し行列Cの対角成分に固有値，行列Dに固有ベクトル

//(15)一般化対称固有値問題 Ax = lambda B x (固有値と固有ベクトル)は使えません...

//eigsymge(A, B, C, D); // Ax = lambda B xに対し行列Cの対角成分に固有値，行列Dに固有ベクトル

kv::interval< double >型 + BLAS 2 (条件 : kvライブラリ、Lapack、(丸めモードを要求する)BLAS)

//(1)インクルードファイル

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
#include <vcp/pidblas.hpp>
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void){

//(2)型

//vcp::matrix< kv::interval< double >, vcp::pidblas >

//(3)変数の宣言

vcp::matrix< kv::interval< double >, vcp::pidblas > A;
vcp::matrix< kv::interval< double >, vcp::pidblas > B, C, D;
vcp::matrix< kv::interval< double >, vcp::pidblas > x, b;

//(11)逆行列は区間Gaussが発動するので使わないで下さい

//C = inv(A); // Aの逆行列をCに代入

//(13)対称固有値問題 Ax = lambda x (固有値と固有ベクトル)は使えません...

//eigsym(A, C, D); // Ax = lambda xに対し行列Cの対角成分に固有値，行列Dに固有ベクトル

//(15)一般化対称固有値問題 Ax = lambda B x (固有値と固有ベクトル)は使えません...

//eigsymge(A, B, C, D); // Ax = lambda B xに対し行列Cの対角成分に固有値，行列Dに固有ベクトル

行列クラスの変換

vcp::convert関数

vcp::convert関数はvcp::matrix<T1, P1>型をvcp::matrix<T2, P2>に変換する関数である。 精度保証付き数値計算を行う際はやや危険な関数になる可能性があります。 特にvcp::matrix<kv::interval<T1>,P1>からvcp::matrix<T2,P2>へ変換する場合、kvのmid関数を使って近似計算になることに注意して下さい。 また、kvに入っているスカラー型(kv::ddなど)もvcp::convertによって同様に変換されます。 仕様を理解してから使うようにして下さい。

//使用例1 vcp::matrix<double> と vcp::matrix<kv::mpfr<N> >の変換

#include <kv/interval.hpp> //近似でもinterval.hppが必要
#include <kv/rdouble.hpp> //近似でもrdouble.hppが必要
#include <kv/mpfr.hpp>
#include <kv/rmpfr.hpp> //近似でもrmpfr.hppが必要
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void){

int n = 5;
vcp::matrix< double > A;
vcp::matrix< kv::mpfr< 500 > > B;
A.rand(n);
std::cout << "A = \n" << A << std::endl;
vcp::convert(A, B); // vcp::matrix<double> から vcp::matrix<kv::mpfr<500> >に変換
std::cout << "B = \n" << B << std::endl;
A.clear();
vcp::convert(B, A); // vcp::matrix<kv::mpfr<500> > から vcp::matrix<double>に変換
std::cout << "A = \n" << A << std::endl;

return 0;
}

//使用例2 vcp::matrix<kv::dd> と vcp::matrix<kv::mpfr<N> >の変換

#include <kv/interval.hpp> //近似でもinterval.hppが必要
#include <kv/rdouble.hpp> //doubleを使わなくてもrdouble.hppが必要
#include <kv/dd.hpp>
#include <kv/rdd.hpp> //近似でもrdd.hppが必要
#include <kv/mpfr.hpp>
#include <kv/rmpfr.hpp> //近似でもrmpfr.hppが必要
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void){

int n = 5;
vcp::matrix< kv::dd > A;
vcp::matrix< kv::mpfr< 500 > > B;
A.rand(n);
std::cout << "A = \n" << A << std::endl;
vcp::convert(A, B); // vcp::matrix<kv::dd> から vcp::matrix<kv::mpfr<500> >に変換
std::cout << "B = \n" << B << std::endl;
A.clear();
vcp::convert(B, A); // vcp::matrix<kv::mpfr<500> > から vcp::matrix<kv::dd>に変換
std::cout << "A = \n" << A << std::endl;

return 0;
}

//使用例3 vcp::matrix<kv::dd> と vcp::matrix<double, vcp::pdblas >の変換

#include <kv/interval.hpp> //近似でもinterval.hppが必要
#include <kv/rdouble.hpp> //近似でもrdouble.hppが必要
#include <kv/dd.hpp>
#include <kv/rdd.hpp> //近似でもrdd.hppが必要
#include <vcp/pdblas.hpp>
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void){

int n = 5;
vcp::matrix< double, vcp::pdblas > A;
vcp::matrix< kv::dd > B;
A.rand(n);
std::cout << "A = \n" << A << std::endl;
vcp::convert(A, B); // vcp::matrix<double, vcp::pdblas> から vcp::matrix<kv::dd>に変換
std::cout << "B = \n" << B << std::endl;
A.clear();
vcp::convert(B, A); // vcp::matrix<kv::dd> から vcp::matrix<double, vcp::pdblas>に変換
std::cout << "A = \n" << A << std::endl;

return 0;
}

//使用例4 vcp::matrix< kv::interval< kv::dd >, vcp::imats< kv::dd > > と vcp::matrix< double, vcp::pdblas >の変換

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>
#include <vcp/pdblas.hpp>
#include <vcp/imats.hpp>
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void){

int n = 5;
vcp::matrix< double, vcp::pdblas > A;
vcp::matrix< kv::interval< kv::dd >, vcp::imats< kv::dd > > B;
A.rand(n);
std::cout << "A = \n" << A << std::endl;
vcp::convert(A, B); // double -> interval<double> -> interval<kv::dd>の順で変換
std::cout << "B = \n" << B << std::endl;
A.clear();
vcp::convert(B, A); // (危険要素) interval<kv::dd> -> kv::dd(mid関数) -> doubleの順で変換
std::cout << "A = \n" << A << std::endl;

return 0;
}

//使用例5 vcp::matrix< kv::interval< double >, vcp::pidblas > と vcp::matrix< kv::interval< kv::dd >, vcp::imats< kv::dd > > の変換

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>
#include <vcp/imats.hpp>
#include <vcp/pidblas.hpp>
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void){

int n = 5;
vcp::matrix< kv::interval< double >, vcp::pidblas > A;
vcp::matrix< kv::interval< kv::dd >, vcp::imats< kv::dd > > B;
A.rand(n);
std::cout << "A = \n" << A << std::endl;
vcp::convert(A, B); // interval<double> -> interval<kv::dd> (精度保証付き)
std::cout << "B = \n" << B << std::endl;
A.clear();
vcp::convert(B, A); // interval<kv::dd> -> interval<double> (精度保証付き)
std::cout << "A = \n" << A << std::endl;

return 0;
}

//使用例6 (危険要素)スカラーでも同様に使える kv::interval< double > と kv::interval< kv::dd > の変換

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void){

kv::interval< double > A;
kv::interval< kv::dd > B;
A = 3;
std::cout << "A = \n" << A << std::endl;
vcp::convert(A, B); // interval<double> -> interval<kv::dd> (精度保証付き)
std::cout << "B = \n" << B << std::endl;
vcp::convert(B, A); // interval<kv::dd> -> interval<double> (精度保証付き)
std::cout << "A = \n" << A << std::endl;

return 0;
}

//使用不可例7 vcp::matrix< kv::interval< double >, vcp::pidblas > と vcp::matrix< double, vcp::pdblas > の変換

//#include <kv/interval.hpp>
//#include <kv/rdouble.hpp>
//#include <kv/dd.hpp>
//#include <kv/rdd.hpp>
//#include <vcp/imats.hpp>
//#include <vcp/pidblas.hpp>
//#include <vcp/matrix.hpp>
//#include <vcp/matrix_assist.hpp>

//int main(void){

//int n = 5;
//vcp::matrix< kv::interval< double >, vcp::pidblas > A;
//vcp::matrix< double, vcp::pdblas > B;
//A.rand(n);
//std::cout << "A = \n" << A << std::endl;
//vcp::convert(A, B); // 次に説明するvcp::mid(A,B)を使用
//std::cout << "B = \n" << B << std::endl;
//A.clear();
//vcp::convert(B, A); // 次に説明するvcp::interval(A,B)を使用
//std::cout << "A = \n" << A << std::endl;

//return 0;
//}

vcp::mid関数、vcp::rad関数

vcp:matrix< kv::interval< T >, P1 >の中心を計算する関数midと半径を計算する関数radの紹介をします。 どちらも vcp:matrix< kv::interval< T >, P1 > から vcp:matrix< T, P2 > に変換します。 ポリシーは変化しても構いませんが、型Tは変更してはいけません。 また、mid関数とrad関数はkvのmid関数とrad関数を呼び出しているので、この関数を使って中心半径型の区間に変形することはできません。 中心半径型に変形したい場合は次に説明するmidrad関数を使って下さい。

//使用例1 vcp::matrix< kv::interval< double >, vcp::pidblas > と vcp::matrix< double, vcp::pdblas > の変換

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
#include <vcp/pidblas.hpp>
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void){

int n = 5;
vcp::matrix< kv::interval< double >, vcp::pidblas > A;
vcp::matrix< double, vcp::pdblas > B;
A.rand(n);
A=ltransmul(A);
std::cout << "A = \n" << A << std::endl;
vcp::mid(A, B); // Aの全成分のそれぞれの中心値をBに代入
std::cout << "B = \n" << B << std::endl;
vcp::rad(A, B); // Aの全成分のそれぞれの半径値をBに代入
std::cout << "B = \n" << B << std::endl;

return 0;
}

vcp::midrad関数

vcp:matrix< kv::interval< T >, P1 >を中心半径型に変換する関数vcp::midrad関数の紹介をします。 vcp:matrix< kv::interval< T >, P1 > から 中心vcp:matrix< T, P2 > 半径vcp:matrix< T, P2 >に変換します。 ポリシーは変化しても構いませんが、型Tは変更してはいけません。 また、midrad関数はkvのmidrad関数を呼び出しているので、この関数を使って中心半径型の区間に変形できます。

//使用例1 vcp::matrix< kv::interval< double >, vcp::pidblas > の中心半径型への変換

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
#include <vcp/pidblas.hpp>
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void){

int n = 5;
vcp::matrix< kv::interval< double >, vcp::pidblas > A;
vcp::matrix< double, vcp::pdblas > B, C;
A.rand(n);
A=ltransmul(A);
std::cout << "A = \n" << A << std::endl;
vcp::midrad(A, B, C); // Aの全成分のそれぞれの中心値をBに、半径値をCに代入
std::cout << "B = \n" << B << std::endl;
std::cout << "C = \n" << C << std::endl;

return 0;
}

vcp::interval関数

点行列vcp:matrix< T, P1 >を区間行列vcp:matrix< kv::interval< T >, P2 >に変換する関数vcp::intervalの紹介をします。 ポリシーは変化しても構いませんが、型Tは変更してはいけません。

//使用例1 vcp::matrix< double, vcp::pidblas > の区間行列への変換

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
#include <vcp/pidblas.hpp>
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void){

int n = 5;
vcp::matrix< double, vcp::pdblas > A;
vcp::matrix< kv::interval< double >, vcp::pidblas > B;
A.rand(n);
std::cout << "A = \n" << A << std::endl;
vcp::interval(A, B); // Aの型をinterval<double>に変換してBに代入
std::cout << "B = \n" << B << std::endl;

return 0;
}

行列クラスのbool型

行列同士の比較

行列同士の比較をする場合，同じ型，同じポリシー同士でなければ比較できません。 また、比較結果はvcp::matrix<bool>型として扱うことができます。

//使用例

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
#include <vcp/pidblas.hpp>
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void) {
vcp::matrix<kv::interval<double >,vcp::pidblas> A, B;
vcp::matrix<bool> AA;

A.rand(5);
B.rand(5);
std::cout << (A < B) << std::endl; // 成分毎に比較
std::cout << (A <= B) << std::endl; // 成分毎に比較
std::cout << (A > B) << std::endl; // 成分毎に比較
std::cout << (A >= B) << std::endl; // 成分毎に比較

AA = A < B // 成分毎に比較した結果はmatrix<bool<型
return 0;
}

matrix<bool>型の演算子

vcp::matrix<bool>型の演算子として&&と||と!が演算子多重定義されています。 これらの演算子は成分毎の論理和(&&)と論理積(||)と否定(!)を計算します。

//使用例

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
#include <vcp/pidblas.hpp>
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void) {
vcp::matrix<kv::interval<double >,vcp::pidblas> A, B, C;

A.rand(5);
B.rand(5);
C.rand(5);

std::cout << ((A < B) && (A >= C)) << std::endl;
std::cout << ((A < B) || (A >= C)) << std::endl;
std::cout << !(A < B) << std::endl;

return 0;
}

matrix<bool>型の関数

vcp::matrix<bool>型の関数としてall関数、any関数 、none関数が定義されています。 引数はvcp::matrix<bool>型で戻り値はbool型です。
all関数:全ての要素がTureならばbool型でTrue。
any関数:少なくとも一つTureならばbool型でTrue。
all関数:全ての要素がFalseならばbool型でTrue。

//使用例

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
#include <vcp/pidblas.hpp>
#include <vcp/matrix.hpp>
#include <vcp/matrix_assist.hpp>

int main(void) {
vcp::matrix<kv::interval<double >,vcp::pidblas> A, B, C;

A.rand(5);
B.rand(5);
C.rand(5);

std::cout << all(A < B) << std::endl;
std::cout << any(A < B) << std::endl;
std::cout << none(A < B) << std::endl;

return 0;
}

vcp::timeの使い方

VCP Libraryのvcp_timer.hppには実行時間計測用のクラスtimerとtimerクラスのインスタンスtimeがある． vcp_timer.hppを以下のようにincludeすると使用可能になる

#include <vcp/vcp_timer.hpp>

インスタンスtimeを使うことによって簡単に実行時間を計測することが可能である． timerクラスのメソッドは
・tic : 時間計測開始(startに現在時刻を代入)
・toc : 時間計測終了(endに現在時刻を代入，diffに開始時刻と終了時刻の差を代入，unsigned int型のoにdiffを変換し出力。)
・output : unsigned int型のoを戻り値として返す関数
である．

//使用例

#include <vcp/vcp_timer.hpp>
#include <vcp/matrix.hpp>

int main(void) {
vcp::matrix< double > A, B, C;
int n = 100;
unsigned int t;

A.rand(n);
B.rand(n);

vcp::time.tic();
C = A*B;
vcp::time.toc();
t = vcp::time.output();
std::cout << t/1000.0 << " [sec]" << std::endl;
std::cout << t/1000.0/60.0 << " [min]" << std::endl;

return 0;
}

ルジャンドル基底クラス

ldbase.hppとは?

VCP Libraryのldbase.hppには例えば中尾先生，木下先生の論文"On very accurate verification of solutions for boundary value problems by using spectral methods"で取り上げられているLegendre多項式による\( H^1_0(\Omega)\)の基底を生成するクラス Legendre_Bases_Generator< DataType1, DataType2, MatrixPolicy >があります. 使い方については「test_Emden.cpp」や「test_Gray-Scott.cpp」を見た方が早いと思います． クラスLegendre_Bases_Generatorは3つの設定をすることで，偏微分方程式の近似解および解の存在証明を行う上で必要となる積分値，ベクトル，行列を返すメソッド関数を持っています.

想定している方程式:
\(n\)次元領域\( \Omega = (0, 1)^n \)とし，楕円型偏微分方程式系のDirichlet境界値問題 \begin{align} \left\{ \begin{array}{ll} -\Delta u_1 = f_1(u_1, u_2, \cdots, u_m) ~ &\mbox{in} ~ \Omega \\ -\Delta u_2 = f_2(u_1, u_2, \cdots, u_m) ~ &\mbox{in} ~ \Omega \\ &\cdots \\ -\Delta u_m = f_m(u_1, u_2, \cdots, u_m) ~ &\mbox{in} ~ \Omega \\ u_1 = u_2 = \cdots = u_m = 0 ~ &\mbox{on} ~ \partial \Omega \end{array} \right. \end{align} の解\( u_1, ~ u_2, ~ \cdots, ~ u_m \in H^1_0(\Omega) \)を求める際に必要な素材を作成します. 但し，\(f_1, ~ \cdots, ~ f_m \)は多項式に限ります.

Legendre_Bases_Generator型の変数の宣言とテンプレート引数

クラスLegendre_Bases_Generatorには，DataType1, DataType2, MatrixPolicyの3つのテンプレート引数が必要です.

型名: Legendre_Bases_Generator< DataType1, DataType2, MatrixPolicy >
DataType1: 内部で利用する精度保証付き高精度型(例 kv::interval< kv::mpfr< 500 > >)
DataType2: Newton法や固有値問題を計算する型.精度保証の有無は問わない.(例 kv::dd )
MatrixPolicy: DataType2に対応するvcp::matrixのPolicy(例 mats< kv::dd > )

以下，詳細は随時アップロードします．