EigenFuncs.cpp

// EigenFuncs.cpp : DLL アプリケーション用にエクスポートされる関数を定義します。
//

#include "stdafx.h"

/*----------------*/
/* EigenFuncs.cpp */
/*----------------*/
// EigenFuncs.cpp : DLL アプリケーション用にエクスポートされる関数を定義します。
//

#define EIGEN_NO_DEBUG // コード内のassertを無効化．
#define EIGEN_DONT_PARALLELIZE // 並列を無効化．
#define EIGEN_UNROLLING_LIMIT 200 // 260415Cl 追加 デフォルト110から増やしミニ行列演算のループ展開を促進
#define EIGENFUNCS_EXPORTS

//#include <complex>
#include <stdio.h>
#include <string.h>
#include <vector>
#include "math.h"

#include <../Eigen/Core>
#include <../Eigen/Dense>
#include <../Eigen/Geometry>
#include <../Eigen/LU>
#include <../Eigen/Eigenvalues>
#include <../unsupported/Eigen/MatrixFunctions>

#include "EigenFuncs.h"

using namespace Eigen;
using namespace std;

const std::complex<double> two_pi_i = complex<double>(0, 2 * 3.141592653589793238462643383279);
//const size_t sizeComplex = sizeof(complex<double>);

#define Mat MatrixXcd
//縦ベクトル
#define Vec VectorXcd
//横ベクトル
#define VecR RowVectorXcd


extern "C" {
	//行列c0~c3をr0~r3の割合でブレンドする
	EIGEN_FUNCS_API void _Blend(int dim, double c0[], double c1[], double c2[], double c3[], double r0, double r1, double r2, double r3, double result[])
	{
		auto _c0 = Map<VectorXd>(c0, dim);
		auto _c1 = Map<VectorXd>(c1, dim);
		auto _c2 = Map<VectorXd>(c2, dim);
		auto _c3 = Map<VectorXd>(c3, dim);
		Map<VectorXd>(result, dim).noalias() = r0 * _c0 + r1 * _c1 + r2 * _c2 + r3 * _c3;
	}

	//行列c0~c3をr0~r3の割合でブレンドし、自己共役を得る
	EIGEN_FUNCS_API void _BlendAndConjugate(int dim, double c0[], double c1[], double c2[], double c3[], double r0, double r1, double r2, double r3, double result[])
	{
		auto _c0 = Map<Vec>((dcomplex*)c0, dim);
		auto _c1 = Map<Vec>((dcomplex*)c1, dim);
		auto _c2 = Map<Vec>((dcomplex*)c2, dim);
		auto _c3 = Map<Vec>((dcomplex*)c3, dim);
		Map<Vec>((dcomplex*)result, dim).noalias() = r0 * _c0.conjugate() + r1 * _c1.conjugate() + r2 * _c2.conjugate() + r3 * _c3.conjugate();

	}

	//STEMの非弾性散乱電子強度の計算用の特殊関数
	EIGEN_FUNCS_API void _BlendAdJointMul_Mul_Mul(int dim, double c0[], double c1[], double c2[], double c3[], double r0, double r1, double r2, double r3, double mat2[], double mat3[], double result[])
	{
		auto _c0 = Map<Mat>((dcomplex*)c0, dim, dim);
		auto _c1 = Map<Mat>((dcomplex*)c1, dim, dim);
		auto _c2 = Map<Mat>((dcomplex*)c2, dim, dim);
		auto _c3 = Map<Mat>((dcomplex*)c3, dim, dim);
		auto _m2 = Map<Mat>((dcomplex*)mat2, dim, dim);
		auto _m3 = Map<Mat>((dcomplex*)mat3, dim, dim);
		Map<Mat>((dcomplex*)result, dim, dim).noalias() = (r0 * _c0 + r1 * _c1 + r2 * _c2 + r3 * _c3).adjoint() * _m2 * _m3;
	}

	//STEMの非弾性散乱電子強度の計算用の特殊関数
	EIGEN_FUNCS_API void _AdJointMul_Mul_Mul(int dim, double mat1[], double mat2[], double mat3[], double result[])
	{
		auto _mat1 = Map<Mat>((dcomplex*)mat1, dim, dim);
		auto _mat2 = Map<Mat>((dcomplex*)mat2, dim, dim);
		auto _mat3 = Map<Mat>((dcomplex*)mat3, dim, dim);
		Map<Mat>((dcomplex*)result, dim, dim).noalias() = _mat1.adjoint() * _mat2 * _mat3;
	}

	//STEM用の特殊関数. 透過係数を求める
	EIGEN_FUNCS_API void _GenerateTC1(int dim, double thickness, double _kg_z[], double _val[], double _vec[], double _tc_k[])
	{
		auto val = (dcomplex*)_val;

		Vec exp_kgz = Vec(dim);
		Vec exp_val = Vec(dim);

		for (int i = 0; i < dim; ++i)
		{
			exp_kgz[i] = exp(two_pi_i * thickness * _kg_z[i]);
			exp_val[i] = exp(two_pi_i * thickness * val[i]);
		}
		auto m = Map<Mat>((dcomplex*)_vec, dim, dim);
		Map<Vec>((dcomplex*)_tc_k, dim).noalias() = exp_kgz.asDiagonal() * m * exp_val;
	}

	//STEM用の特殊関数. 透過係数を求める (260610Cl 追加: _GenerateTC1 の改良版)
	// 旧 _GenerateTC1 は呼び出しごとに Eigen Vec を2本 heap 確保し、exp_kgz.asDiagonal()*m*exp_val の
	// 左結合評価で一時行列が生じ得た。thread_local workspace 再利用 + a⊙(M·b) の評価順に変更 (数学的に同一、丸めのみ僅差)
	EIGEN_FUNCS_API void _GenerateTC1B(int dim, double thickness, double _kg_z[], double _val[], double _vec[], double _tc_k[])
	{
		static thread_local std::vector<dcomplex> work;
		if ((int)work.size() < dim)
			work.resize(dim);

		auto val = (dcomplex*)_val;
		for (int i = 0; i < dim; ++i)
			work[i] = exp(two_pi_i * thickness * val[i]);

		auto tc = (dcomplex*)_tc_k;
		Map<Vec>(tc, dim).noalias() = Map<Mat>((dcomplex*)_vec, dim, dim) * Map<Vec>(work.data(), dim);
		for (int i = 0; i < dim; ++i)
			tc[i] *= exp(two_pi_i * thickness * _kg_z[i]);
	}

	//STEMの非弾性散乱を q 単位で一括計算 (260610Cl 追加)
	// W = Uq(dim×dim) × TC_sub(dim×entryCount) を列タイル GEMM で計算し、エントリ e ごとに
	//   tmp = Σ_m r4[4e+m]·dot(conj(tcAll列 n4[4e+m]), W列e)   (Eigen の dot は第1引数を共役する)
	// を取り、sumD[d] += tmp·lenzAll[e·lDim+d] まで native 内で完結する。
	// 旧経路 (_BlendAndConjugate + _RowVec_SqMat_ColVec をエントリごとに P/Invoke 2回) と
	// 総和順序の入れ替えのみで数学的に等価。GEMV (メモリバウンド) の列挙を GEMM (高計算強度) に置き換える。
	EIGEN_FUNCS_API void _STEM_InelasticQ(int dim, int lDim, int entryCount,
		double _Uq[], double _tcAll[], int kIndices[], int n4[], double r4[], double _lenzAll[], double _sumD[])
	{
		constexpr int TILE = 64;
		static thread_local std::vector<dcomplex> tcTile, wTile;
		if ((int)tcTile.size() < dim * TILE)
		{
			tcTile.resize((size_t)dim * TILE);
			wTile.resize((size_t)dim * TILE);
		}

		auto U = Map<Mat>((dcomplex*)_Uq, dim, dim);
		auto tcAll = (dcomplex*)_tcAll;
		auto lenzAll = (dcomplex*)_lenzAll;
		auto sumD = (dcomplex*)_sumD;

		for (int e0 = 0; e0 < entryCount; e0 += TILE)
		{
			const int cols = entryCount - e0 < TILE ? entryCount - e0 : TILE;
			for (int c = 0; c < cols; ++c)//対象 k 列を連続バッファへ pack (コスト O(dim·cols) ≪ GEMM O(dim²·cols))
				memcpy(tcTile.data() + (size_t)c * dim, tcAll + (size_t)kIndices[e0 + c] * dim, sizeof(dcomplex) * dim);

			Map<Mat> TC(tcTile.data(), dim, cols);
			Map<Mat> W(wTile.data(), dim, cols);
			W.noalias() = U * TC;

			for (int c = 0; c < cols; ++c)
			{
				const int e = e0 + c;
				auto w = W.col(c);
				dcomplex tmp(0, 0);
				for (int m = 0; m < 4; ++m)
				{
					auto tcn = Map<Vec>(tcAll + (size_t)n4[e * 4 + m] * dim, dim);
					tmp += r4[e * 4 + m] * tcn.dot(w);//dot = Σ conj(tc_n)·w
				}
				for (int d = 0; d < lDim; ++d)
					sumD[d] += tmp * lenzAll[(size_t)e * lDim + d];
			}
		}
	}

	//STEM用の特殊関数. 透過係数を求める
	EIGEN_FUNCS_API void _GenerateTC2(int dim, double thickness, double _kg_z[], double _val[], double _vec[], double _tc_k[], double _tc_kq[])
	{
		auto val = (dcomplex*)_val;

		Vec exp_kgz = Vec(dim);
		Vec exp_val = Vec(dim);

		for (int i = 0; i < dim; ++i)
		{
			exp_kgz[i] = exp(two_pi_i * thickness * _kg_z[i]);
			exp_val[i] = exp(two_pi_i * thickness * val[i]);
		}
		auto m = Map<Mat>((dcomplex*)_vec, dim, dim);
		auto res = exp_kgz.asDiagonal() * m * exp_val;

		Map<Vec>((dcomplex*)_tc_k, dim).noalias() = res;
		Map<Vec>((dcomplex*)_tc_kq, dim).noalias() = res.conjugate();
	}

	//横ベクトル×正方行列×縦ベクトルの掛算. STEMの非弾性散乱を求めるときに使用
	EIGEN_FUNCS_API void _RowVec_SqMat_ColVec(int dim, double rowVec[], double sqMat[], double colVec[], double _result[])
	{
		auto rV = Map<VecR>((dcomplex*)rowVec, dim);
		auto m = Map<Mat>((dcomplex*)sqMat, dim, dim);
		auto cV = Map<Vec>((dcomplex*)colVec, dim);
		Map<Mat>((dcomplex*)_result, 1, 1).noalias() = rV * m * cV;
	}

	//横ベクトル×正方行列×縦ベクトルの掛算. STEMの非弾性散乱を求めるときに使用
	EIGEN_FUNCS_API void _SqMat_ColVec(int dim, double sqMat[], double colVec[], double _result[])
	{
		auto m = Map<Mat>((dcomplex*)sqMat, dim, dim);
		auto cV = Map<Vec>((dcomplex*)colVec, dim);
		Map<Vec>((dcomplex*)_result, dim).noalias() = m * cV;
	}

	EIGEN_FUNCS_API void _STEM_INEL1(int dim, double rowVec[], int n[], double r[], double sqMat[], double colVec[], double _result[])
	{
		auto rV1 = Map<VecR>((dcomplex*)(rowVec + n[0] * dim * 2), dim);
		auto rV2 = Map<VecR>((dcomplex*)(rowVec + n[1] * dim * 2), dim);
		auto rV3 = Map<VecR>((dcomplex*)(rowVec + n[2] * dim * 2), dim);
		auto rV4 = Map<VecR>((dcomplex*)(rowVec + n[3] * dim * 2), dim);
		auto m = Map<Mat>((dcomplex*)sqMat, dim, dim);
		auto cV = Map<Vec>((dcomplex*)colVec, dim);
		// 260602Cl 修正: 従来は blend に共役が無く、managed の BlendAndConjugate (各 c を共役してブレンド) + RowVec_SqMat_ColVec
		// 経路と数値不一致だった。r は実数なので conj(Σ r_k rV_k) = Σ r_k conj(rV_k)。各行ベクトルを共役してブレンドし、
		// _BlendAndConjugate (= r0*c0.conjugate()+…) と同形にする。これで 2 つの P/Invoke を 1 つに融合できる。
		// Map<Mat>((dcomplex*)_result, 1, 1).noalias() = (r[0] * rV1 + r[1] * rV2 + r[2] * rV3 + r[3] * rV4) * m * cV; // 260602Cl 変更前 (共役なし)
		Map<Mat>((dcomplex*)_result, 1, 1).noalias() = (r[0] * rV1.conjugate() + r[1] * rV2.conjugate() + r[2] * rV3.conjugate() + r[3] * rV4.conjugate()) * m * cV;
	}

	//複素非対称行列のmat1とmat2の要素ごとの掛算(アダマール積)を取る
	EIGEN_FUNCS_API void _PointwiseMultiply(int dim, double mat1[], double mat2[], double result[])
	{
		auto m1 = Map<Mat>((dcomplex*)mat1, dim, dim);
		auto m2 = Map<Mat>((dcomplex*)mat2, dim, dim);
		Map<Mat>((dcomplex*)result, dim, dim).noalias() = m1.cwiseProduct(m2);
	}

	//複素行列のmat1を共役転値して、mat2に掛ける
	EIGEN_FUNCS_API void _AdjointAndMultiply(int dim, double mat1[], double mat2[], double result[])
	{
		auto m1 = Map<Mat>((dcomplex*)mat1, dim, dim);
		auto m2 = Map<Mat>((dcomplex*)mat2, dim, dim);
		Map<Mat>((dcomplex*)result, dim, dim).noalias() = m1.adjoint() * m2;
	}

	//複素行列同士の乗算を求める
	EIGEN_FUNCS_API void _MultiplyMM(int dim, double mat1[], double mat2[], double result[])
	{
		auto m1 = Map<Mat>((dcomplex*)mat1, dim, dim);
		auto m2 = Map<Mat>((dcomplex*)mat2, dim, dim);
		Map<Mat>((dcomplex*)result, dim, dim).noalias() = m1 * m2;
	}

	//複素行列同士の乗算を求める
	EIGEN_FUNCS_API void _MultiplyMMM(int dim, double mat1[], double mat2[], double mat3[], double result[])
	{
		auto m1 = Map<Mat>((dcomplex*)mat1, dim, dim);
		auto m2 = Map<Mat>((dcomplex*)mat2, dim, dim);
		auto m3 = Map<Mat>((dcomplex*)mat3, dim, dim);
		Map<Mat>((dcomplex*)result, dim, dim).noalias() = m1 * m2 * m3;
	}

	//実数行列同士の乗算を求める
	EIGEN_FUNCS_API void _MultiplyMM_Real(int dim, double mat1[], double mat2[], double result[])
	{
		auto m1 = Map<MatrixXd>((double*)mat1, dim, dim);
		auto m2 = Map<MatrixXd>((double*)mat2, dim, dim);
		Map<MatrixXd>((double*)result, dim, dim).noalias() = m1 * m2;
	}

	//複素行列と複素ベクトルの乗算を求める
	EIGEN_FUNCS_API void _MultiplyMV(int dim, double mat[], double vec[], double result[])
	{
		auto m = Map<Mat>((dcomplex*)mat, dim, dim);
		auto v = Map<Vec>((dcomplex*)vec, dim);
		Map<Vec>((dcomplex*)result, dim).noalias() = m * v;
	}

	//複素ベクトル同士の乗算を求める
	EIGEN_FUNCS_API void _MultiplyVV(int dim, double vec1[], double vec2[], double result[])
	{
		auto v1 = Map<VecR>((dcomplex*)vec1, dim);
		auto v2 = Map<Vec>((dcomplex*)vec2, dim);
		Map<Mat>((dcomplex*)result, 1, 1).noalias() = v1 * v2;
	}

	//複素数と複素ベクトルの乗算を求める
	EIGEN_FUNCS_API void _MultiplySV(int dim, double real, double imag, double vec[], double result[])
	{
		auto s = dcomplex(real, imag);
		auto v = Map<Vec>((dcomplex*)vec, dim);
		Map<Vec>((dcomplex*)result, dim).noalias() = s * v;
	}

	//ベクトル同士の加算を求める. 実数と複素数の両方いける
	EIGEN_FUNCS_API void _AddVV(int dim, double vec1[], double vec2[], double result[])
	{
		auto v1 = Map<VectorXd>(vec1, dim);
		auto v2 = Map<VectorXd>(vec2, dim);
		Map<VectorXd>(result, dim).noalias() = v1 + v2;
	}

	//複素ベクトル同士の要素ごとの除算を求める
	EIGEN_FUNCS_API void _DivideVV(int dim, double vec1[], double vec2[], double result[])
	{
		auto v1 = Map<Vec>((dcomplex*)vec1, dim);
		auto v2 = Map<Vec>((dcomplex*)vec2, dim);

		auto res = (dcomplex*)result;

		for (int i = 0; i < dim; i++)
			res[i] = v1[i] / v2[i];
	}

	//複素ベクトル同士の減算を求める. 実数と複素数の両方いける
	EIGEN_FUNCS_API void _SubtractVV(int dim, double vec1[], double vec2[], double result[])
	{
		auto v1 = Map<VectorXd>(vec1, dim);
		auto v2 = Map<VectorXd>(vec2, dim);
		Map<VectorXd>(result, dim).noalias() = v1 - v2;
	}

	//複素非対称行列の逆行列を求める
	EIGEN_FUNCS_API void _Inverse(int dim, double mat[], double inverse[])
	{
		Map<Mat>((dcomplex*)inverse, dim, dim).noalias() = Map<Mat>((dcomplex*)mat, dim, dim).inverse();
	}

	//PartialPivLuSolve
	EIGEN_FUNCS_API void _PartialPivLuSolve(int dim, double mat[], double vec[], double result[])
	{
		auto m = Map<Mat>((dcomplex*)mat, dim, dim);
		auto v = Map<Vec>((dcomplex*)vec, dim);
		Map<Vec>((dcomplex*)result, dim).noalias() = m.partialPivLu().solve(v);
	}

	//複素非対称行列の固有値、固有ベクトルを求める
	EIGEN_FUNCS_API void _EigenSolver(int dim, double mat[], double eigenValues[], double eigenVectors[])
	{
		ComplexEigenSolver<Mat> solver(Map<Mat>((dcomplex*)mat, dim, dim));
		Map<Vec>((dcomplex*)eigenValues, dim).noalias() = solver.eigenvalues();
		Map<Mat>((dcomplex*)eigenVectors, dim, dim).noalias() = solver.eigenvectors();
	}

	//複素非対称行列の行列指数を求める
	EIGEN_FUNCS_API void _MatrixExponential(int dim, double mat[], double results[])
	{
		Map<Mat>((dcomplex*)results, dim, dim).noalias() = Map<MatrixXcd>((dcomplex*)mat, dim, dim).exp();
		//memcpy(results, Map<MatrixXcd>((dcomplex*)mat, dim, dim).exp().eval().data(), sizeComplex * dim * dim);//上の書き方でいいのか、未検証
	}

	//CBEDソルバー
	EIGEN_FUNCS_API void _CBEDSolver_Eigen(int dim, double potential[], double psi0[], int tDim, double thickness[], double result[])
	{
		ComplexEigenSolver<Mat> solver(Map<Mat>((dcomplex*)potential, dim, dim));
		auto values = solver.eigenvalues();
		Mat vectors = solver.eigenvectors() * solver.eigenvectors().partialPivLu().solve(Map<Vec>((dcomplex*)psi0, dim)).asDiagonal();

		Mat gamma = Mat(dim, tDim);
		for (int t = 0; t < tDim; ++t)
		{
			const auto coeff2 = two_pi_i * thickness[t];
			for (int g = 0; g < dim; ++g)
				gamma(g, t) = exp(values[g] * coeff2);
		}

		Map<Mat>((dcomplex*)result, dim, tDim).noalias() = vectors * gamma;
	}

	//CBEDソルバー. 固有値、固有ベクトル、αも返す
	EIGEN_FUNCS_API void _CBEDSolver_Eigen2(int dim, double potential[], double psi0[], int tDim, double thickness[], double Values[], double Vectors[], double Alphas[], double Tg[])
	{
		auto vals = Map<Vec>((dcomplex*)Values, dim);
		auto vecs = Map<Mat>((dcomplex*)Vectors, dim, dim);
		auto alphas = Map<Vec>((dcomplex*)Alphas, dim);
		auto tg = Map<Mat>((dcomplex*)Tg, dim, tDim);

		ComplexEigenSolver<Mat> solver(Map<Mat>((dcomplex*)potential, dim, dim));
		vals.noalias() = solver.eigenvalues();
		vecs.noalias() = solver.eigenvectors();
		alphas.noalias() = vecs.partialPivLu().solve(Map<Vec>((dcomplex*)psi0, dim));
		Vec gamma_alpha = Vec(dim);
		for (int t = 0; t < tDim; ++t)
		{
			const auto coeff2 = two_pi_i * thickness[t];
			for (int g = 0; g < dim; ++g)
				gamma_alpha[g] = exp(vals[g] * coeff2) * alphas[g];
			tg.col(t).noalias() = vecs * gamma_alpha;
		}
	}

	// 260704Cl 追加: Hermitian化 + 吸収摂動 CBED ソルバー (ReciPro_BetheMethod高速化提案.md §8 Phase H0 検証用)
	// 260829Cl 変更: managed getEigenMatrix の P 除算が列割り M = B̃D_P⁻¹ から行割り A = D_P⁻¹B̃ へ修正されたのに追随し、
	//   √P 変換の向きを反転 (At の si/sj、ψ0 側の ×√P、出力側の ÷√P の 3 箇所)。対称化後の S = D_P^{-1/2}B̃D_P^{-1/2}
	//   自体は両規約で同一なので、H/H'/摂動機構は無変更。
	// 入力は現行 getEigenMatrix の出力そのまま A = D_P⁻¹·(M+D_Q) (column-major, 行割り, 非Hermitian) と対角 P (>0 必須, 呼び出し側で保証)。
	// 対称化 Ã(i,j) = A(i,j)·√(P_i/P_j) は H + iH' (H, H' とも厳密に Hermitian; U(-g)=conj(U(g)) による) に分解でき、
	// SelfAdjointEigenSolver(H) の実固有値 γ_j とユニタリ W に吸収 H' を摂動として載せる。
	// ψ(t) = D_P^{-1/2} W diag(exp(2πiλ_j t)) α,  α = W^H D_P^{+1/2} ψ0 (W ユニタリのため LU 不要)。
	// perturb: 0 = 摂動なし (吸収完全無視, 比較用)
	//          1 = 対角 1 次摂動 λ_j = γ_j + i·w_j^H H' w_j (準縮退で破綻し得る)
	//          2 = 準縮退クラスター摂動: E = W^H H' W を作り、ペア条件 clusterTol·|E_jk| > γ_k−γ_j (結合強度 vs gap)
	//              を満たす固有値区間をクラスター化し、クラスター内 m×m 行列 diag(γ_C) + i·E_C を厳密対角化
	//              (縮退摂動論の標準処方)。W_C ← W_C·V_C, α_C ← V_C⁻¹·α_C で固有対も更新 (V_C は非ユニタリになり得る)
	//          3 = 2 に加えクラスタ外ペアの固有ベクトル 1 次補正 + 固有値 2 次補正 (RS 摂動論):
	//              U = blockdiag(V_C) + C, C_kj = (iE)_kj/(γ_j−γ_k) (別クラスタペアのみ)、
	//              λ_j −= Σ_k |E_jk|²/(γ_j−γ_k) (単独列のみ, 実シフト)。V = W·U, α = U⁻¹ W^H D^{-1/2}ψ0 (U は N×N LU 1本)。
	//              残差誤差は (|E|/gap) の 1 次 → 2 次に落ちる
	// stats (perturb>=2 のとき書き込み, nullptr 可): [0] = クラスター化後の残差危険度 max |E_jk|/gap (クラスタ外ペア),
	//              [1] = m>1 クラスター数, [2] = 最大クラスターサイズ
	// 出力契約は _CBEDSolver_Eigen と同一 (result = dim×tDim column-major)。
	EIGEN_FUNCS_API void _CBEDSolverHermitian(int dim, double potential[], double p[], double psi0[], int tDim, double thickness[], double result[], int perturb, double clusterTol, double stats[])
	{
		auto A = Map<Mat>((dcomplex*)potential, dim, dim);
		VectorXd sq(dim);
		for (int i = 0; i < dim; ++i)
			sq[i] = sqrt(p[i]);

		Mat At(dim, dim);
		for (int j = 0; j < dim; ++j)
		{
			const double sj = sq[j];
			for (int i = 0; i < dim; ++i)
				At(i, j) = A(i, j) * (sq[i] / sj); // 260829Cl 変更: 行割り入力 A に対し Ã = D^{+1/2}AD^{-1/2}。変更前: A(i, j) * (sj / sq[i])
		}

		Mat H = (At + At.adjoint()) * 0.5;
		SelfAdjointEigenSolver<Mat> solver(H);
		const auto& gamma = solver.eigenvalues(); // 実固有値 (昇順)
		Mat W = solver.eigenvectors();            // ユニタリ (クラスター処理で列を更新するためコピー)

		// α = W^H D_P^{+1/2} ψ0 (クラスター処理前の W で計算し、後段で α_C ← V_C⁻¹ α_C を適用)
		Vec dpsi0(dim);
		for (int i = 0; i < dim; ++i)
			dpsi0[i] = ((dcomplex*)psi0)[i] * sq[i]; // 260829Cl 変更: 行割り規約では ×√P。変更前: / sq[i]
		Vec alpha;
		alpha.noalias() = W.adjoint() * dpsi0;

		Vec lambda(dim);
		if (perturb == 0)
			for (int j = 0; j < dim; ++j)
				lambda[j] = gamma[j];
		else
		{
			Mat Habs = (At - At.adjoint()) * dcomplex(0, -0.5); // H' (Hermitian, 吸収)
			Mat T;
			T.noalias() = Habs * W; // O(N³) GEMM
			if (perturb == 1)
			{
				for (int j = 0; j < dim; ++j)
					lambda[j] = dcomplex(gamma[j], W.col(j).dot(T.col(j)).real()); // E_jj (Hermitian 二次形式なので実)
			}
			else
			{
				Mat E;
				E.noalias() = W.adjoint() * T; // E = W^H H' W (Hermitian)。O(N³) GEMM もう1本 (EVD 比では小)
				for (int j = 0; j < dim; ++j)
					lambda[j] = dcomplex(gamma[j], E(j, j).real());

				// ペア条件による区間併合: clusterTol·|E_jk| > gap なら [j,k] を同一クラスタへ (γ はソート済み)
				std::vector<int> reach(dim);
				for (int j = 0; j < dim; ++j)
					reach[j] = j;
				for (int j = 0; j < dim; ++j)
					for (int k = j + 1; k < dim; ++k)
						if (clusterTol * abs(E(j, k)) > gamma[k] - gamma[j] && k > reach[j])
							reach[j] = k;

				std::vector<int> clusterOf(dim);
				Mat U = Mat::Identity(dim, dim); // 基底 W 上の合成変換 (クラスタ回転 + ベクトル1次補正)
				bool hasU = false;
				int nClusters = 0, maxCluster = 1;
				int s = 0;
				while (s < dim)
				{
					int t2 = reach[s];
					for (int q = s; q <= t2; ++q) // 推移閉包 (区間の連鎖拡張)
						t2 = max(t2, reach[q]);
					const int m = t2 - s + 1;
					for (int q = s; q <= t2; ++q)
						clusterOf[q] = s;
					if (m > 1)
					{
						++nClusters;
						maxCluster = max(maxCluster, m);
						hasU = true;
						// クラスター内 m×m: B = diag(γ_C) + i·E_C を厳密対角化
						Mat B = dcomplex(0, 1) * E.block(s, s, m, m);
						for (int k = 0; k < m; ++k)
							B(k, k) += gamma[s + k];
						ComplexEigenSolver<Mat> cs(B);
						lambda.segment(s, m) = cs.eigenvalues();
						U.block(s, s, m, m) = cs.eigenvectors();
					}
					s = t2 + 1;
				}

				if (perturb >= 3)
				{   // クラスタ外ペアの RS 摂動: 固有ベクトル 1 次 C_kj = (iE)_kj/(γ_j−γ_k)、単独列の固有値 2 次 (実シフト)
					for (int j = 0; j < dim; ++j)
					{
						double shift2 = 0;
						for (int k = 0; k < dim; ++k)
							if (clusterOf[k] != clusterOf[j])
							{
								const double dg = gamma[j] - gamma[k];
								U(k, j) += dcomplex(0, 1) * E(k, j) / dg;
								shift2 -= (E(j, k) * E(k, j)).real() / dg; // = −|E_jk|²/dg (E は Hermitian)
							}
						if (clusterOf[j] == j && (j + 1 >= dim || clusterOf[j + 1] != j)) // 単独列のみ (クラスタ列の外部 2 次は高次扱いで省略)
							lambda[j] += shift2;
					}
					hasU = true;
				}

				if (hasU)
				{
					alpha = U.partialPivLu().solve(alpha.eval()); // α = U⁻¹ (W^H D^{-1/2} ψ0)
					W = (W * U).eval();                           // V = W·U
				}

				if (stats != nullptr)
				{
					double resid = 0; // クラスター化されずに残ったペアの max |E_jk|/gap = 1次近似の危険度
					for (int j = 0; j < dim; ++j)
						for (int k = j + 1; k < dim; ++k)
							if (clusterOf[j] != clusterOf[k])
								resid = max(resid, abs(E(j, k)) / (gamma[k] - gamma[j]));
					stats[0] = resid;
					stats[1] = nClusters;
					stats[2] = maxCluster;
				}
			}
		}

		auto res = Map<Mat>((dcomplex*)result, dim, tDim);
		Vec ga(dim);
		for (int t = 0; t < tDim; ++t)
		{
			const auto coeff = two_pi_i * thickness[t];
			for (int j = 0; j < dim; ++j)
				ga[j] = exp(lambda[j] * coeff) * alpha[j];
			res.col(t).noalias() = W * ga;
			for (int g = 0; g < dim; ++g)
				res(g, t) /= sq[g]; // ψ = D_P^{-1/2}·(W·ga)。260829Cl 変更 (行割り規約)。変更前: *= sq[g]
		}
	}

	//CBEDソルバー
	EIGEN_FUNCS_API void _CBEDSolver_MtxExp(int dim, double potential[], double psi0[], int tDim, double tStart, double tStep, double result[])
	{
		auto res = Map<Mat>((dcomplex*)result, dim, tDim);
		Mat matExp = (two_pi_i * tStart * Map<Mat>((dcomplex*)potential, dim, dim)).exp().eval();
		Vec vec = matExp * Map<Vec>((dcomplex*)psi0, dim);
		res.col(0).noalias() = vec;
		if (tStep == 0)
			return;

		if (tStart != tStep)
			matExp = (two_pi_i * tStep * Map<Mat>((dcomplex*)potential, dim, dim)).exp().eval();
		for (int t = 1; t < tDim; ++t)
		{
			vec = matExp * vec;
			res.col(t).noalias() = vec;
		}
	}

	// EBSD強度計算ソルバー
	EIGEN_FUNCS_API void _EBSDSolver(
		int bLen, int nAtoms, int tLen,
		double eigenValues[],
		double eigenVectors[],
		double alpha[],
		double phaseNG[],
		double sigma[],
		double thicknesses[],
		double intensity[])
	{
		auto vals = Map<Vec>((dcomplex*)eigenValues, bLen);
		auto C = Map<Mat>((dcomplex*)eigenVectors, bLen, bLen);
		auto a = Map<Vec>((dcomplex*)alpha, bLen);
		auto P = Map<Mat>((dcomplex*)phaseNG, nAtoms, bLen);

		// β = P × C × diag(α)
		Mat B(nAtoms, bLen);
		B.noalias() = P * C * a.asDiagonal();

		Mat sigB(nAtoms, bLen);
		for (int n = 0; n < nAtoms; ++n)
			sigB.row(n).noalias() = sigma[n] * B.row(n);
		Mat S(bLen, bLen);
		S.noalias() = B.adjoint() * sigB;

		// λ_{jj'} = 2πi(γ_j - conj(γ_j'))
		Mat lam(bLen, bLen);
		for (int j = 0; j < bLen; ++j)
			for (int jp = 0; jp < bLen; ++jp)
				lam(j, jp) = two_pi_i * (vals[j] - conj(vals[jp]));

		// 260420Cl S, F, lam, expAccum は全て Hermitian → sum = diag の real 和 + 上三角の 2·Re 和 で内側ループを半減

		// 等間隔判定
		bool isUniform = (tLen >= 2) && (abs(thicknesses[1] - 2.0 * thicknesses[0]) < 1e-10);

		if (isUniform && tLen > 1)
		{
			double dt = thicknesses[0];
			Mat expStep(bLen, bLen);
			// 260420Cl expStep も Hermitian → 上三角計算 + 下三角は conj。後段 Hadamard 累積で両側要るため完全に埋める
			for (int j = 0; j < bLen; ++j)
			{
				expStep(j, j) = exp(lam(j, j) * dt);
				for (int jp = j + 1; jp < bLen; ++jp)
				{
					auto v = exp(lam(j, jp) * dt);
					expStep(j, jp) = v;
					expStep(jp, j) = conj(v);
				}
			}

			Mat expAccum = expStep;

			for (int t = 0; t < tLen; ++t)
			{
				double thick = thicknesses[t];
				double sumReal = 0;
				for (int j = 0; j < bLen; ++j)
				{
					auto l = lam(j, j);
					dcomplex F = (abs(l) < 1e-15) ? dcomplex(thick, 0) : (expAccum(j, j) - 1.0) / l;
					sumReal += (S(j, j) * F).real();
				}
				for (int j = 0; j < bLen; ++j)
					for (int jp = j + 1; jp < bLen; ++jp)
					{
						auto l = lam(j, jp);
						dcomplex F = (abs(l) < 1e-15) ? dcomplex(thick, 0) : (expAccum(j, jp) - 1.0) / l;
						sumReal += 2.0 * (S(j, jp) * F).real();
					}
				intensity[t] = max(0.0, sumReal);

				if (t < tLen - 1)
					expAccum.array() *= expStep.array();
			}
		}
		else
		{
			for (int t = 0; t < tLen; ++t)
			{
				double thick = thicknesses[t];
				double sumReal = 0;
				for (int j = 0; j < bLen; ++j)
				{
					auto l = lam(j, j);
					dcomplex F = (abs(l) < 1e-15) ? dcomplex(thick, 0) : (exp(l * thick) - 1.0) / l;
					sumReal += (S(j, j) * F).real();
				}
				for (int j = 0; j < bLen; ++j)
					for (int jp = j + 1; jp < bLen; ++jp)
					{
						auto l = lam(j, jp);
						dcomplex F = (abs(l) < 1e-15) ? dcomplex(thick, 0) : (exp(l * thick) - 1.0) / l;
						sumReal += 2.0 * (S(j, jp) * F).real();
					}
				intensity[t] = max(0.0, sumReal);
			}
		}
	}

	// 260316Cl 追加
	// EBSD強度 + TDS バックグラウンド一括計算ソルバー
	EIGEN_FUNCS_API void _EBSDSolverWithTDS(
		int bLen, int nAtoms, int tLen,
		double eigenValues[],
		double eigenVectors[],
		double alpha[],
		double phaseNG[],
		double sigma[],
		double muBack[],
		double tdsCoeff,
		double thicknesses[],
		double intensity[],
		double tdsIntensity[])
	{
		auto vals = Map<Vec>((dcomplex*)eigenValues, bLen);
		auto C = Map<Mat>((dcomplex*)eigenVectors, bLen, bLen);
		auto a = Map<Vec>((dcomplex*)alpha, bLen);
		auto P = Map<Mat>((dcomplex*)phaseNG, nAtoms, bLen);
		auto U = Map<Mat>((dcomplex*)muBack, bLen, bLen);

		// --- 弾性: S = B† diag(σ) B ---
		Mat B(nAtoms, bLen);
		B.noalias() = P * C * a.asDiagonal();

		Mat sigB(nAtoms, bLen);
		for (int n = 0; n < nAtoms; ++n)
			sigB.row(n).noalias() = sigma[n] * B.row(n);
		Mat S(bLen, bLen);
		S.noalias() = B.adjoint() * sigB;

		// --- TDS: M = tdsCoeff × diag(α†) × C† × U × C × diag(α) ---
		Mat CUC(bLen, bLen);
		CUC.noalias() = C.adjoint() * U * C;
		for (int j = 0; j < bLen; ++j)
			for (int jp = 0; jp < bLen; ++jp)
				CUC(j, jp) *= tdsCoeff * conj(a[j]) * a[jp];

		// --- λ_{jj'} 事前計算 ---
		Mat lam(bLen, bLen);
		for (int j = 0; j < bLen; ++j)
			for (int jp = 0; jp < bLen; ++jp)
				lam(j, jp) = two_pi_i * (vals[j] - conj(vals[jp]));

		// 260420Cl S, CUC, lam, expAccum, F はいずれも Hermitian → 対角 + 上三角 2·Re で sumS/sumM を半減。CUC の a[j]a[jp] 係数も Hermitian 性を保存

		// 等間隔判定
		bool isUniform = (tLen >= 2) && (abs(thicknesses[1] - 2.0 * thicknesses[0]) < 1e-10);

		if (isUniform && tLen > 1)
		{
			double dt = thicknesses[0];
			Mat expStep(bLen, bLen);
			// 260420Cl expStep も Hermitian → 上三角計算 + 下三角 conj で埋める
			for (int j = 0; j < bLen; ++j)
			{
				expStep(j, j) = exp(lam(j, j) * dt);
				for (int jp = j + 1; jp < bLen; ++jp)
				{
					auto v = exp(lam(j, jp) * dt);
					expStep(j, jp) = v;
					expStep(jp, j) = conj(v);
				}
			}

			Mat expAccum = expStep;

			for (int t = 0; t < tLen; ++t)
			{
				double thick = thicknesses[t];
				double sumSReal = 0, sumMReal = 0;
				for (int j = 0; j < bLen; ++j)
				{
					auto l = lam(j, j);
					dcomplex F = (abs(l) < 1e-15) ? dcomplex(thick, 0) : (expAccum(j, j) - 1.0) / l;
					sumSReal += (S(j, j) * F).real();
					sumMReal += (CUC(j, j) * F).real();
				}
				for (int j = 0; j < bLen; ++j)
					for (int jp = j + 1; jp < bLen; ++jp)
					{
						auto l = lam(j, jp);
						dcomplex F = (abs(l) < 1e-15) ? dcomplex(thick, 0) : (expAccum(j, jp) - 1.0) / l;
						sumSReal += 2.0 * (S(j, jp) * F).real();
						sumMReal += 2.0 * (CUC(j, jp) * F).real();
					}
				intensity[t] = max(0.0, sumSReal);
				tdsIntensity[t] = max(0.0, sumMReal);

				if (t < tLen - 1)
					expAccum.array() *= expStep.array();
			}
		}
		else
		{
			for (int t = 0; t < tLen; ++t)
			{
				double thick = thicknesses[t];
				double sumSReal = 0, sumMReal = 0;
				for (int j = 0; j < bLen; ++j)
				{
					auto l = lam(j, j);
					dcomplex F = (abs(l) < 1e-15) ? dcomplex(thick, 0) : (exp(l * thick) - 1.0) / l;
					sumSReal += (S(j, j) * F).real();
					sumMReal += (CUC(j, j) * F).real();
				}
				for (int j = 0; j < bLen; ++j)
					for (int jp = j + 1; jp < bLen; ++jp)
					{
						auto l = lam(j, jp);
						dcomplex F = (abs(l) < 1e-15) ? dcomplex(thick, 0) : (exp(l * thick) - 1.0) / l;
						sumSReal += 2.0 * (S(j, jp) * F).real();
						sumMReal += 2.0 * (CUC(j, jp) * F).real();
					}
				intensity[t] = max(0.0, sumSReal);
				tdsIntensity[t] = max(0.0, sumMReal);
			}
		}
	}

} // extern "C"