SparseMat

有限要素法で使う疎行列用ライブラリです。本ライブラリのstaticライブラリファイルは 「SparseMat.lib」 です。

EigenのEigen::SparseMatrixをベースとして作成されており、Eigenの疎行列ラッパーとして機能します。 double型、complex型の両方に対応し、疎行列の四則演算なども用意されています。ソルバとしては、オリジナルの加速係数つきICCGソルバが用意されています。その他、Eigenが用意しているICCGとBiCGstabのラッパも用意しています。

内部のintは、slv_intとして別の型として定義しています。

・using slv_int = int;　(SparseMatTMPL.hppの冒頭)

基本はこのままでいいと思いますが、超大規模にしたい場合はunsingedやlong intに書き換えてください

SparseMatクラス

実数用の疎行列クラスです。基本的には、まず初期化した後、add関数で行列要素を追加していきます。追加が完了したら、fixメソッドで疎行列を確定させてください。確定後にメソッドや四則演算が行えます。fix前に各種操作を行うとエラーになる恐れがあるので注意してください。「SparseMat/SparseMat.hpp」をインクルードして使ってください。

* 基本操作

SparseMat();

コンストラクタ：疎行列を初期化します。インスタンスが作成されるだけで内部は何も初期化されません。この状態でaddをしてもエラーになります。

SparseMat(slv_int x);

コンストラクタ：疎行列を初期化します。行数サイズをxで初期化します。（内部でtempInitializeを呼び、初期化します）

SparseMat(const SparseMat& mat);

コンストラクタ：別の疎行列を代入して初期化します。本方法で初期化すると最初からfix状態になります。

bool isFixed()

呼び出し元インスタンスが確定済み(fixを呼び出した後か)かどうかを返します。

bool isEmpty() const

行列の中身が空かどうかを返します。

void tempInitialize()

疎行列の状況を初期化し、fixされていない状態にします。

void fix()

疎行列を確定させます。

void resetMat()

fix済みの疎行列の０にします。（疎行列の位置はそのまま、値だけゼロにします）

void add(slv_int gyo, slv_int retu, double val)

疎行列の指定の位置に値を加えます。行の位置をgyoで、列の位置をretuで指定し、そこにValを加えます。 fix前の場合は任意の位置にaddできます。 fix後は、すでに値のある位置以外を指定するとエラーになります。

slv_int isInclude(slv_int gyo, slv_int target_r)const

疎行列のi行目に、target_r列があるかどうかを判定します。あったらそのindexを返します。

void getTargetRowVal(slv_int target, std::vector<slv_int>& row_pos, std::vector& row_val)const

指定した列の非ゼロの行位置と値をvector「row_val」にコピーします。

void getTargetColVal(slv_int target, std::vector<slv_int>& col_pos, std::vector& col_val)const

指定した行の非ゼロの列位置と値をvector「col_val」にコピーします。

void printMat(const std::string& str="Mat.csv")

疎行列の内容をstrで指定したファイルに書き出します。

void print()

疎行列をコンソールに表示します。

* オペレータ群

double* operator(const double vec) const;

dcomplex* operator(const dcomplex vec) const;

疎行列に列ベクトルvecを掛けて、結果ベクトルを返します。列ベクトルはC++のデフォルトの配列で渡します。（vecは実数、複素数どちらでも対応）

Eigen::VectorXd operator*(const Eigen::VectorXd& vec) const;

Eigen::VectorXcd operator*(const Eigen::VectorXcd& vec) const;

疎行列に列ベクトルvecを掛けて、結果ベクトルを返します。列ベクトルはEigenのVectorで渡します。（vecは実数、複素数どちらでも対応）

void operator*=(const double x)

自身の全要素にxを掛けます。

SparseMat operator*(const double x) const;

SparseMatC operator*(const dcomplex x) const;

自身の全要素にxを掛けた後の疎行列を返します。（実数、複素数どちらでも対応）

SparseMat operator*(const SparseMat& mat) const;

SparseMatC operator*(const SparseMatC& mat) const;

自身と行列matの行列積を行い、結果の疎行列を返します。

SparseMat operator+(const SparseMat& mat) const;

SparseMatC operator+(const SparseMatC& mat) const;

自身と行列matを足し、結果の疎行列を返します。

* その他の行列オペレータ群

SparseMat trans() const

自身の転置行列を返します。

SparseMat makeSubMat(slv_int range1a, slv_int range1b, slv_int range2a, slv_int range2b)

範囲[r1a, r1b]×[r2a, r2b]の部分を抜き出して部分行列を作り、得られた部分行列を返します。

void MatDiv(SparseMat& matK11, SparseMat& matK12, slv_int rangeA, slv_int rangeB)

列rangeBを境に、行列を2つに分け、行rangeAより下は削除します。

SparseMat getMatLower() const;

下三角行列を作成して返します。

SparseMat getMatUpper() const

上三角行列を作成して返します。

SparseMat inv() const

逆行列（密アルゴリズム・大行列でやるとメモリが吹き飛ぶ！）

void round()

行列の要素を四捨五入して丸めます。

SparseMatCクラス

複素数用の疎行列クラスです。基本的な操作は全て実数用のSparseMatと同じです。「SparseMat/SparseMatC.hpp」をインクルードして使ってください。

SparseMatOperatorsクラス

複数の疎行列の操作を行うstaticなメソッドを管理するクラスです。「SparseMat/SparseMatOperators.hpp」をインクルードして使ってください。

static SparseMat plusShift(const SparseMat& mat1, const SparseMat& mat2, const double a1, const double a2, const slv_int pos1, const slv_int pos2)

static SparseMatC plusShift(const SparseMatC& mat1, const SparseMatC& mat2, const double a1, const double a2, const slv_int pos1, const slv_int pos2)

(a1mat1 + a2mat2)を計算して、結果の疎行列を返します。 mat2に足す位置をpos1,pos2でずらすことが可能です。

static double* dotVecMat2(const SparseMat& mat1, const SparseMat& mat2, const double* vec1);

static dcomplex* dotVecMat2(const SparseMat& mat1, const SparseMatC& mat2, const double* vec1);

static dcomplex* dotVecMat2(const SparseMatC& mat1, const SparseMat& mat2, const double* vec1);

static dcomplex* dotVecMat2(const SparseMatC& mat1, const SparseMatC& mat2, const double* vec1);

static dcomplex* dotVecMat2(const SparseMat& mat1, const SparseMat& mat2, const dcomplex* vec1);

static dcomplex* dotVecMat2(const SparseMat& mat1, const SparseMatC& mat2, const dcomplex* vec1);

static dcomplex* dotVecMat2(const SparseMatC& mat1, const SparseMat& mat2, const dcomplex* vec1);

static dcomplex* dotVecMat2(const SparseMatC& mat1, const SparseMatC& mat2, const dcomplex* vec1);

(mat1 + mat2)*vecBを計算し、得られたベクトルを返します。

static SparseMat dotMats(const SparseMat& matA, const SparseMat& matB, const SparseMat& matC);

static SparseMatC dotMats(const SparseMat& matA, const SparseMat& matB, const SparseMatC& matC);

static SparseMatC dotMats(const SparseMat& matA, const SparseMatC& matB, const SparseMatC& matC);

static SparseMatC dotMats(const SparseMatC& matA, const SparseMatC& matB, const SparseMatC& matC);

疎行列AとBとCをかけて結果を返します。

static void plusFix(SparseMat& matAB, const SparseMat& matA, const SparseMat& matB, double a1=1.0, double a2=1.0, slv_int pos1=0, slv_int pos2=0);

static void plusFix(SparseMatC& matAB, const SparseMat& matA, const SparseMatC& matB, double a1=1.0, double a2=1.0, slv_int pos1=0, slv_int pos2=0);

static void plusFix(SparseMatC& matAB, const SparseMatC& matA, const SparseMatC& matB, double a1=1.0, double a2=1.0, slv_int pos1=0, slv_int pos2=0);

fix済みの行列matABCに、(a1matA+a2matB)の結果を代入します。 matBの開始位置はposだけずらすことが可能です。

static void plusFix(SparseMat& matABC, const SparseMat& matA, const SparseMat& matB, const SparseMat& matC, double a1=1.0, double a2=1.0, double a3=1.0);

static void plusFix(SparseMatC& matABC, const SparseMatC& matA, const SparseMatC& matB, const SparseMatC& matC, double a1=1.0, double a2=1.0, double a3=1.0);

fix済みの行列matABCに、(a1A+a2B+a3*C)の結果を代入します。

static void dotFix(SparseMat& matAB, const SparseMat& matA, const SparseMat& matB);

static void dotFix(SparseMatC& matAB, const SparseMat& matA, const SparseMatC& matB);

static void dotFix(SparseMatC& matAB, const SparseMatC& matA, const SparseMatC& matB);

fix済みの行列matABCに、AとBの行列積の結果を代入します。

static void dotFix(SparseMat& matAB, const SparseMat& matA, const SparseMat& matB, const SparseMat& matC);

static void dotFix(SparseMatC& matAB, const SparseMat& matA, const SparseMat& matB, const SparseMatC& matC);

static void dotFix(SparseMatC& matAB, const SparseMat& matA, const SparseMatC& matB, const SparseMatC& matC);

static void dotFix(SparseMatC& matAB, const SparseMatC& matA, const SparseMatC& matB, const SparseMatC& matC);

fix済みの行列matABCに、AとBとCの行列積の結果を代入します。

MatSolversクラス

疎行列用のソルバを提供するクラスです。全てstaticメソッドです。「SparseMat/MatSolvers.hpp」をインクルードして使ってください。

static SparseMat plusShift(const SparseMat& mat1, const SparseMat& mat2, const double a1, const double a2, const slv_int pos1, const slv_int pos2)

static SparseMatC plusShift(const SparseMatC& mat1, const SparseMatC& mat2, const double a1, const double a2, const slv_int pos1, const slv_int pos2)

(a1mat1 + a2mat2)を計算して、結果の疎行列を返します。 mat2に足す位置をpos1,pos2でずらすことが可能です。

static bool solveICCG(const slv_int size0, const double conv_cri, const int max_ite, const double accera, const SparseMat& matA, const double vecB, double results, bool init=false)

static bool solveICCG(const slv_int size0, const double conv_cri, const int max_ite, const double accera, const double normB, const SparseMat& matA, const double vecB, double results, bool init=false)

static bool solveICCG(const slv_int size0, const double conv_cri, const int max_ite, const double accera, const SparseMat& matA, const Eigen::VectorXd& vecB, double *results, bool init=false)

自作のICCG法です。右辺をC++配列、Eigen::VectorXdで与えるかで2パターンあります。また、引く数normBは、右辺bのノルムを事前に与えるか、内部で計算するかです。（normBは、収束判定で使います。"r/normB < conv_cri"となると収束とします。）

size0：行列の行数です。 conv_cri：収束判定値です。10^-6など。 max_ite：最大反復数です。 accera：加速係数です。 matA：Ax=bのAです。対称行列としてください。 vecB:Ax=bのbです。size0のベクトルです。 results：解ベクトルです（Ax=bのx）。メモリは事前にsize0だけ確保しておいてください。 init：解ベクトルxをこの処理の前にゼロ初期化するかです。過去の値がxに入っていて流用したい場合はfalseにしてください。それ以外はtrueにすれば内部でゼロ初期化してから始めます。

static bool solveICCG(slv_int size0, const double conv_cri, const int max_ite, const double accera, const SparseMatC& matA, const dcomplex vecB, dcomplex results, bool init=false);

static bool solveICCG(slv_int size0, const double conv_cri, const int max_ite, const double accera, const double normB, const SparseMatC& matA, const dcomplex vecB, dcomplex results, bool init=false);

static bool solveICCG(slv_int size0, const double conv_cri, const int max_ite, const double accera, const SparseMatC& matA, const Eigen::VectorXcd& vecB, dcomplex *results, bool init=false);

上記ICCG法の複素数版です。使い方は同じです。ただしmatAは対称行列にのみ対応しています（エルミートなどは未対応）

static bool solveEigenICCG(const slv_int size0, const double conv_cri, const int max_ite, const SparseMat& matA, const Eigen::VectorXd& vecB, Eigen::VectorXd& results, bool init=false);

static bool solveEigenICCG(const slv_int size0, const double conv_cri, const int max_ite, const SparseMat& matA, const double* vecB, double* results, bool init=true);

static bool solveEigenBiCGstab(const slv_int size0, const double conv_cri, const int max_ite, const SparseMat& matA, const Eigen::VectorXd& vecB, Eigen::VectorXd& results, bool init=false);

static bool solveEigenBiCGstab(const slv_int size0, const double conv_cri, const int max_ite, const SparseMat& matA, const double* vecB, double* results, bool init=true);

static bool solveEigenICCG(const slv_int size0, const double conv_cri, const int max_ite, const SparseMatC& matA, const Eigen::VectorXcd& vecB, Eigen::VectorXcd& results, bool init=false);

static bool solveEigenICCG(const slv_int size0, const double conv_cri, const int max_ite, const SparseMatC& matA, const dcomplex* vecB, dcomplex* results, bool init=true);

static bool solveEigenBiCGstab(const slv_int size0, const double conv_cri, const int max_ite, const SparseMatC& matA, const Eigen::VectorXcd& vecB, Eigen::VectorXcd& results, bool init=false);

static bool solveEigenBiCGstab(const slv_int size0, const double conv_cri, const int max_ite, const SparseMatC& matA, const dcomplex* vecB, dcomplex* results, bool init=false);

Eigenに付属している疎行列ソルバのラッパーです。使い方は上記ICCG法の引数と同じです。

SparseMat

SparseMatクラス

* 基本操作

SparseMat();

SparseMat(slv_int x);

SparseMat(const SparseMat& mat);

bool isFixed()

bool isEmpty() const

void tempInitialize()

void fix()

void resetMat()

void add(slv_int gyo, slv_int retu, double val)

slv_int isInclude(slv_int gyo, slv_int target_r)const

void getTargetRowVal(slv_int target, std::vector<slv_int>& row_pos, std::vector& row_val)const

void getTargetColVal(slv_int target, std::vector<slv_int>& col_pos, std::vector& col_val)const

void printMat(const std::string& str="Mat.csv")

void print()

* オペレータ群

double* operator*(const double* vec) const;

dcomplex* operator*(const dcomplex* vec) const;

Eigen::VectorXd operator*(const Eigen::VectorXd& vec) const;

Eigen::VectorXcd operator*(const Eigen::VectorXcd& vec) const;

void operator*=(const double x)

SparseMat operator*(const double x) const;

SparseMatC operator*(const dcomplex x) const;

SparseMat operator*(const SparseMat& mat) const;

SparseMatC operator*(const SparseMatC& mat) const;

SparseMat operator+(const SparseMat& mat) const;

SparseMatC operator+(const SparseMatC& mat) const;

* その他の行列オペレータ群

SparseMat trans() const

SparseMat makeSubMat(slv_int range1a, slv_int range1b, slv_int range2a, slv_int range2b)

void MatDiv(SparseMat& matK11, SparseMat& matK12, slv_int rangeA, slv_int rangeB)

SparseMat getMatLower() const;

SparseMat getMatUpper() const

SparseMat inv() const

void round()

SparseMatCクラス

SparseMatOperatorsクラス

static SparseMat plusShift(const SparseMat& mat1, const SparseMat& mat2, const double a1, const double a2, const slv_int pos1, const slv_int pos2)

static SparseMatC plusShift(const SparseMatC& mat1, const SparseMatC& mat2, const double a1, const double a2, const slv_int pos1, const slv_int pos2)

static double* dotVecMat2(const SparseMat& mat1, const SparseMat& mat2, const double* vec1);

static dcomplex* dotVecMat2(const SparseMat& mat1, const SparseMatC& mat2, const double* vec1);

static dcomplex* dotVecMat2(const SparseMatC& mat1, const SparseMat& mat2, const double* vec1);

static dcomplex* dotVecMat2(const SparseMatC& mat1, const SparseMatC& mat2, const double* vec1);

static dcomplex* dotVecMat2(const SparseMat& mat1, const SparseMat& mat2, const dcomplex* vec1);

static dcomplex* dotVecMat2(const SparseMat& mat1, const SparseMatC& mat2, const dcomplex* vec1);

static dcomplex* dotVecMat2(const SparseMatC& mat1, const SparseMat& mat2, const dcomplex* vec1);

static dcomplex* dotVecMat2(const SparseMatC& mat1, const SparseMatC& mat2, const dcomplex* vec1);

static SparseMat dotMats(const SparseMat& matA, const SparseMat& matB, const SparseMat& matC);

static SparseMatC dotMats(const SparseMat& matA, const SparseMat& matB, const SparseMatC& matC);

static SparseMatC dotMats(const SparseMat& matA, const SparseMatC& matB, const SparseMatC& matC);

static SparseMatC dotMats(const SparseMatC& matA, const SparseMatC& matB, const SparseMatC& matC);

static void plusFix(SparseMat& matAB, const SparseMat& matA, const SparseMat& matB, double a1=1.0, double a2=1.0, slv_int pos1=0, slv_int pos2=0);

static void plusFix(SparseMatC& matAB, const SparseMat& matA, const SparseMatC& matB, double a1=1.0, double a2=1.0, slv_int pos1=0, slv_int pos2=0);

static void plusFix(SparseMatC& matAB, const SparseMatC& matA, const SparseMatC& matB, double a1=1.0, double a2=1.0, slv_int pos1=0, slv_int pos2=0);

static void plusFix(SparseMat& matABC, const SparseMat& matA, const SparseMat& matB, const SparseMat& matC, double a1=1.0, double a2=1.0, double a3=1.0);

static void plusFix(SparseMatC& matABC, const SparseMatC& matA, const SparseMatC& matB, const SparseMatC& matC, double a1=1.0, double a2=1.0, double a3=1.0);

static void dotFix(SparseMat& matAB, const SparseMat& matA, const SparseMat& matB);

static void dotFix(SparseMatC& matAB, const SparseMat& matA, const SparseMatC& matB);

static void dotFix(SparseMatC& matAB, const SparseMatC& matA, const SparseMatC& matB);

static void dotFix(SparseMat& matAB, const SparseMat& matA, const SparseMat& matB, const SparseMat& matC);

static void dotFix(SparseMatC& matAB, const SparseMat& matA, const SparseMat& matB, const SparseMatC& matC);

static void dotFix(SparseMatC& matAB, const SparseMat& matA, const SparseMatC& matB, const SparseMatC& matC);

static void dotFix(SparseMatC& matAB, const SparseMatC& matA, const SparseMatC& matB, const SparseMatC& matC);

MatSolversクラス

static SparseMat plusShift(const SparseMat& mat1, const SparseMat& mat2, const double a1, const double a2, const slv_int pos1, const slv_int pos2)

static SparseMatC plusShift(const SparseMatC& mat1, const SparseMatC& mat2, const double a1, const double a2, const slv_int pos1, const slv_int pos2)

static bool solveICCG(const slv_int size0, const double conv_cri, const int max_ite, const double accera, const SparseMat& matA, const double *vecB, double *results, bool init=false)

static bool solveICCG(const slv_int size0, const double conv_cri, const int max_ite, const double accera, const double normB, const SparseMat& matA, const double *vecB, double *results, bool init=false)

static bool solveICCG(const slv_int size0, const double conv_cri, const int max_ite, const double accera, const SparseMat& matA, const Eigen::VectorXd& vecB, double *results, bool init=false)

static bool solveICCG(slv_int size0, const double conv_cri, const int max_ite, const double accera, const SparseMatC& matA, const dcomplex *vecB, dcomplex *results, bool init=false);

static bool solveICCG(slv_int size0, const double conv_cri, const int max_ite, const double accera, const double normB, const SparseMatC& matA, const dcomplex *vecB, dcomplex *results, bool init=false);

static bool solveICCG(slv_int size0, const double conv_cri, const int max_ite, const double accera, const SparseMatC& matA, const Eigen::VectorXcd& vecB, dcomplex *results, bool init=false);

static bool solveEigenICCG(const slv_int size0, const double conv_cri, const int max_ite, const SparseMat& matA, const Eigen::VectorXd& vecB, Eigen::VectorXd& results, bool init=false);

static bool solveEigenICCG(const slv_int size0, const double conv_cri, const int max_ite, const SparseMat& matA, const double* vecB, double* results, bool init=true);

static bool solveEigenBiCGstab(const slv_int size0, const double conv_cri, const int max_ite, const SparseMat& matA, const Eigen::VectorXd& vecB, Eigen::VectorXd& results, bool init=false);

static bool solveEigenBiCGstab(const slv_int size0, const double conv_cri, const int max_ite, const SparseMat& matA, const double* vecB, double* results, bool init=true);

static bool solveEigenICCG(const slv_int size0, const double conv_cri, const int max_ite, const SparseMatC& matA, const Eigen::VectorXcd& vecB, Eigen::VectorXcd& results, bool init=false);

static bool solveEigenICCG(const slv_int size0, const double conv_cri, const int max_ite, const SparseMatC& matA, const dcomplex* vecB, dcomplex* results, bool init=true);

double* operator(const double vec) const;

dcomplex* operator(const dcomplex vec) const;

static bool solveICCG(const slv_int size0, const double conv_cri, const int max_ite, const double accera, const SparseMat& matA, const double vecB, double results, bool init=false)

static bool solveICCG(const slv_int size0, const double conv_cri, const int max_ite, const double accera, const double normB, const SparseMat& matA, const double vecB, double results, bool init=false)

static bool solveICCG(slv_int size0, const double conv_cri, const int max_ite, const double accera, const SparseMatC& matA, const dcomplex vecB, dcomplex results, bool init=false);

static bool solveICCG(slv_int size0, const double conv_cri, const int max_ite, const double accera, const double normB, const SparseMatC& matA, const dcomplex vecB, dcomplex results, bool init=false);