Skip to content

Commit 64d9ff2

Browse files
seto77claude
seto77
and
claude
committed
Harden OpenGL init and optimize Bethe/HRTEM hot paths
- GLControlAlpha: guard static constructor WMI / OpenGL checks so Win10 with a broken WMI repository no longer crashes at launch (issues #48 #55). Also null-guard MakeCurrent/SetClip/GenerateBitmap against the DisablingOpenGL fallback. - BetheMethod: hoist t-invariant factors, swap r/g loop order, replace LINQ with plain for-loops, and use per-qIndex locks. - NativeWrapper.HRTEM_Helper: drop AsParallel + SelectMany and fill the output arrays with a single pass to cut GC pressure. - EigenFuncs.cpp: exploit Hermitian symmetry in expStep / F / CUC accumulations to halve the work. - Version.cs: add ver4.925 entry. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
1 parent 4d5ddfa commit 64d9ff2

5 files changed

Lines changed: 216 additions & 51 deletions

File tree

‎Crystallography.Native/EigenFuncs.cpp‎

Lines changed: 83 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -363,31 +363,53 @@ extern "C" {
363363
for (int jp = 0; jp < bLen; ++jp)
364364
lam(j, jp) = two_pi_i * (vals[j] - conj(vals[jp]));
365365

366+
// 260420Cl 追加: S と F (および lam) はいずれも Hermitian なので、sum = Σ S(j,jp)·F(j,jp) は
367+
// 対角項の real 和 + 上三角 (j<jp) の real 和の 2 倍に分解でき、内側ループが約 1/2 で済む。
368+
// S = B†·diag(σ)·B (σ 実) → Hermitian、lam(j,jp) = 2πi(γ_j - conj(γ_jp)) → Hermitian、
369+
// expAccum(j,jp) = exp(lam(j,jp)·dt) も Hermitian、F = (expAccum-1)/lam も Hermitian。
370+
366371
// 等間隔判定
367372
bool isUniform = (tLen >= 2) && (abs(thicknesses[1] - 2.0 * thicknesses[0]) < 1e-10);
368373

369374
if (isUniform && tLen > 1)
370375
{
371376
double dt = thicknesses[0];
372377
Mat expStep(bLen, bLen);
378+
// 260420Cl 変更: expStep も Hermitian (lam が Hermitian で Hadamard exp) なので上三角だけ計算し、
379+
// 下三角は conjugate で埋める。expStep は後段の Hadamard product でも使うため完全に埋める必要あり。
373380
for (int j = 0; j < bLen; ++j)
374-
for (int jp = 0; jp < bLen; ++jp)
375-
expStep(j, jp) = exp(lam(j, jp) * dt);
381+
{
382+
expStep(j, j) = exp(lam(j, j) * dt);
383+
for (int jp = j + 1; jp < bLen; ++jp)
384+
{
385+
auto v = exp(lam(j, jp) * dt);
386+
expStep(j, jp) = v;
387+
expStep(jp, j) = conj(v);
388+
}
389+
}
376390

377391
Mat expAccum = expStep;
378392

379393
for (int t = 0; t < tLen; ++t)
380394
{
381395
double thick = thicknesses[t];
382-
dcomplex sum = 0;
396+
double sumReal = 0; // 260420Cl 変更: Hermitian 対称性から結果は実数 → 実数で直接累積
397+
// 対角項: Hermitian のため S(j,j), F(j,j) ともに実数値。積も実数。
383398
for (int j = 0; j < bLen; ++j)
384-
for (int jp = 0; jp < bLen; ++jp)
399+
{
400+
auto l = lam(j, j);
401+
dcomplex F = (abs(l) < 1e-15) ? dcomplex(thick, 0) : (expAccum(j, j) - 1.0) / l;
402+
sumReal += (S(j, j) * F).real();
403+
}
404+
// 非対角 (j<jp): S(jp,j)·F(jp,j) = conj(S(j,jp)·F(j,jp)) なので 2·Re を取って上三角だけ。
405+
for (int j = 0; j < bLen; ++j)
406+
for (int jp = j + 1; jp < bLen; ++jp)
385407
{
386408
auto l = lam(j, jp);
387409
dcomplex F = (abs(l) < 1e-15) ? dcomplex(thick, 0) : (expAccum(j, jp) - 1.0) / l;
388-
sum += S(j, jp) * F;
410+
sumReal += 2.0 * (S(j, jp) * F).real();
389411
}
390-
intensity[t] = max(0.0, sum.real());
412+
intensity[t] = max(0.0, sumReal);
391413

392414
if (t < tLen - 1)
393415
expAccum.array() *= expStep.array();
@@ -398,15 +420,21 @@ extern "C" {
398420
for (int t = 0; t < tLen; ++t)
399421
{
400422
double thick = thicknesses[t];
401-
dcomplex sum = 0;
423+
double sumReal = 0; // 260420Cl 変更: Hermitian 対称性から実数で直接累積
424+
for (int j = 0; j < bLen; ++j)
425+
{
426+
auto l = lam(j, j);
427+
dcomplex F = (abs(l) < 1e-15) ? dcomplex(thick, 0) : (exp(l * thick) - 1.0) / l;
428+
sumReal += (S(j, j) * F).real();
429+
}
402430
for (int j = 0; j < bLen; ++j)
403-
for (int jp = 0; jp < bLen; ++jp)
431+
for (int jp = j + 1; jp < bLen; ++jp)
404432
{
405433
auto l = lam(j, jp);
406434
dcomplex F = (abs(l) < 1e-15) ? dcomplex(thick, 0) : (exp(l * thick) - 1.0) / l;
407-
sum += S(j, jp) * F;
435+
sumReal += 2.0 * (S(j, jp) * F).real();
408436
}
409-
intensity[t] = max(0.0, sum.real());
437+
intensity[t] = max(0.0, sumReal);
410438
}
411439
}
412440
}
@@ -455,33 +483,57 @@ extern "C" {
455483
for (int jp = 0; jp < bLen; ++jp)
456484
lam(j, jp) = two_pi_i * (vals[j] - conj(vals[jp]));
457485

486+
// 260420Cl 追加: S, CUC, lam, expAccum, F はいずれも Hermitian のため sumS, sumM とも
487+
// 対角 (実数) + 上三角 (2·Re) に分解して内側ループを半減。
488+
// CUC の Hermitian 性: CUC_pre = C†·U·C (U Hermitian) は Hermitian。
489+
// CUC(j,jp) *= tdsCoeff·conj(a[j])·a[jp] を適用しても Hermitian 性を保存する
490+
// (CUC(jp,j) = conj(CUC_pre(j,jp))·tdsCoeff·conj(a[jp])·a[j] = conj(CUC(j,jp)))。
491+
458492
// 等間隔判定
459493
bool isUniform = (tLen >= 2) && (abs(thicknesses[1] - 2.0 * thicknesses[0]) < 1e-10);
460494

461495
if (isUniform && tLen > 1)
462496
{
463497
double dt = thicknesses[0];
464498
Mat expStep(bLen, bLen);
499+
// 260420Cl 変更: expStep は Hermitian。上三角を計算し下三角は conjugate で埋める。
500+
// 後段の Hadamard 累積で両側が必要なため完全に埋める。
465501
for (int j = 0; j < bLen; ++j)
466-
for (int jp = 0; jp < bLen; ++jp)
467-
expStep(j, jp) = exp(lam(j, jp) * dt);
502+
{
503+
expStep(j, j) = exp(lam(j, j) * dt);
504+
for (int jp = j + 1; jp < bLen; ++jp)
505+
{
506+
auto v = exp(lam(j, jp) * dt);
507+
expStep(j, jp) = v;
508+
expStep(jp, j) = conj(v);
509+
}
510+
}
468511

469512
Mat expAccum = expStep;
470513

471514
for (int t = 0; t < tLen; ++t)
472515
{
473516
double thick = thicknesses[t];
474-
dcomplex sumS = 0, sumM = 0;
517+
double sumSReal = 0, sumMReal = 0; // 260420Cl 変更: 対称性により結果は実数
518+
// 対角項
475519
for (int j = 0; j < bLen; ++j)
476-
for (int jp = 0; jp < bLen; ++jp)
520+
{
521+
auto l = lam(j, j);
522+
dcomplex F = (abs(l) < 1e-15) ? dcomplex(thick, 0) : (expAccum(j, j) - 1.0) / l;
523+
sumSReal += (S(j, j) * F).real();
524+
sumMReal += (CUC(j, j) * F).real();
525+
}
526+
// 非対角 (上三角のみ、2 倍で計上)
527+
for (int j = 0; j < bLen; ++j)
528+
for (int jp = j + 1; jp < bLen; ++jp)
477529
{
478530
auto l = lam(j, jp);
479531
dcomplex F = (abs(l) < 1e-15) ? dcomplex(thick, 0) : (expAccum(j, jp) - 1.0) / l;
480-
sumS += S(j, jp) * F;
481-
sumM += CUC(j, jp) * F;
532+
sumSReal += 2.0 * (S(j, jp) * F).real();
533+
sumMReal += 2.0 * (CUC(j, jp) * F).real();
482534
}
483-
intensity[t] = max(0.0, sumS.real());
484-
tdsIntensity[t] = max(0.0, sumM.real());
535+
intensity[t] = max(0.0, sumSReal);
536+
tdsIntensity[t] = max(0.0, sumMReal);
485537

486538
if (t < tLen - 1)
487539
expAccum.array() *= expStep.array();
@@ -492,17 +544,24 @@ extern "C" {
492544
for (int t = 0; t < tLen; ++t)
493545
{
494546
double thick = thicknesses[t];
495-
dcomplex sumS = 0, sumM = 0;
547+
double sumSReal = 0, sumMReal = 0; // 260420Cl 変更
548+
for (int j = 0; j < bLen; ++j)
549+
{
550+
auto l = lam(j, j);
551+
dcomplex F = (abs(l) < 1e-15) ? dcomplex(thick, 0) : (exp(l * thick) - 1.0) / l;
552+
sumSReal += (S(j, j) * F).real();
553+
sumMReal += (CUC(j, j) * F).real();
554+
}
496555
for (int j = 0; j < bLen; ++j)
497-
for (int jp = 0; jp < bLen; ++jp)
556+
for (int jp = j + 1; jp < bLen; ++jp)
498557
{
499558
auto l = lam(j, jp);
500559
dcomplex F = (abs(l) < 1e-15) ? dcomplex(thick, 0) : (exp(l * thick) - 1.0) / l;
501-
sumS += S(j, jp) * F;
502-
sumM += CUC(j, jp) * F;
560+
sumSReal += 2.0 * (S(j, jp) * F).real();
561+
sumMReal += 2.0 * (CUC(j, jp) * F).real();
503562
}
504-
intensity[t] = max(0.0, sumS.real());
505-
tdsIntensity[t] = max(0.0, sumM.real());
563+
intensity[t] = max(0.0, sumSReal);
564+
tdsIntensity[t] = max(0.0, sumMReal);
506565
}
507566
}
508567
}

‎Crystallography.OpenGL/GLControlAlpha.cs‎

Lines changed: 31 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
using System;
55
using System.Collections.Generic;
66
using System.ComponentModel;
7+
using System.Diagnostics; // 260420Cl 追加 静的コンストラクタでの例外ログ出力用
78
using System.Drawing;
89
using System.Linq;
910
using System.Management;
@@ -514,19 +515,41 @@ private static Graphics getGraphics(GLControl control)
514515
/// <summary>静的コンストラクタ. バージョン情報などはここでチェック</summary>
515516
static GLControlAlpha()
516517
{
518+
// 260420Cl WMI/OpenGL 初期化を try/catch で囲む。
519+
// issues #48 #55 で Windows 10 環境において WMI リポジトリ破損等で
520+
// ManagementException("Classe non valide") が発生し、静的コンストラクタが
521+
// TypeInitializationException を投げてアプリ全体が起動不能になる事例が報告されている。
522+
// ここで例外を握り潰し、OpenGL 無効化フォールバック(DisablingOpenGL=true)へ倒す。
523+
517524
//ビデオカード検索
518-
using var searcher = new ManagementObjectSearcher(new SelectQuery("Win32_VideoController")); // (260320Ch) WMI リソースを確実に解放する
519-
using var videoControllers = searcher.Get();
520-
foreach (ManagementObject envVar in videoControllers)
521-
GraphicsInfo.Add((Convert.ToString(envVar["name"]) ?? string.Empty, Convert.ToString(envVar["DriverVersion"]) ?? string.Empty)); // (260320Ch) null 安全に GPU 情報を記録する
525+
try // 260420Cl 追加 WMI クエリを try/catch で保護
526+
{
527+
using var searcher = new ManagementObjectSearcher(new SelectQuery("Win32_VideoController")); // (260320Ch) WMI リソースを確実に解放する
528+
using var videoControllers = searcher.Get();
529+
foreach (ManagementObject envVar in videoControllers)
530+
GraphicsInfo.Add((Convert.ToString(envVar["name"]) ?? string.Empty, Convert.ToString(envVar["DriverVersion"]) ?? string.Empty)); // (260320Ch) null 安全に GPU 情報を記録する
531+
}
532+
catch (Exception e) // 260420Cl 追加 WMI リポジトリ破損や Win32_VideoController 未提供環境で落ちないようにする
533+
{
534+
Debug.WriteLine($"GLControlAlpha: WMI Win32_VideoController query failed: {e.Message}");
535+
}
522536

523537
//if(GraphicsInfo.Count>1 && GraphicsInfo[0].Product.Contains("Radeon"))
524538
// DisableTextRendering = true;
525539

526540
//using var glcontrol = new GLControl();
527541
//glcontrol.MakeCurrent();
528542
//var ver = CheckVersion();
529-
Version = checkSupportedVersion(FragShaders.PPLL, FragShaders.ZSORT); // (260319Ch) まず 4.3、次に 4.1 を試して実際に使える上限を判定
543+
try // 260420Cl 追加 OpenGL コンテキスト生成/シェーダ検査中の例外もフォールバック処理へ倒す
544+
{
545+
Version = checkSupportedVersion(FragShaders.PPLL, FragShaders.ZSORT); // (260319Ch) まず 4.3、次に 4.1 を試して実際に使える上限を判定
546+
}
547+
catch (Exception e) // 260420Cl 追加 ドライバ不良でもアプリ自体は起動を継続させる
548+
{
549+
Debug.WriteLine($"GLControlAlpha: OpenGL version check failed: {e.Message}");
550+
Version = 0;
551+
}
552+
530553
if (Version == 0)
531554
{
532555
DisablingOpenGL = true;
@@ -1107,7 +1130,7 @@ private void renderFullscreenQuad(int program)
11071130
#endregion
11081131

11091132
#region GlObjectの追加/削除
1110-
public void MakeCurrent() => glControl.MakeCurrent();
1133+
public void MakeCurrent() { if (glControl != null) glControl.MakeCurrent(); } // 260420Cl DisablingOpenGL 時は glControl が null のため null ガード
11111134

11121135
public void AddObjects(GLObject obj)
11131136
{
@@ -1191,6 +1214,7 @@ public void DeleteAllObjects()
11911214
/// <param name="clip"></param>
11921215
public void SetClip(Clip clip)
11931216
{
1217+
if (Program < 1 || glControl == null) return; // 260420Cl 追加 DisablingOpenGL 時の二次的 NullReferenceException を回避
11941218
glControl.MakeCurrent();
11951219
Clip = clip;
11961220
clip?.Generate(Program);
@@ -1619,6 +1643,7 @@ private void glControl_Resize(object sender, EventArgs e)
16191643

16201644
public Bitmap GenerateBitmap()
16211645
{
1646+
if (Program < 1 || glControl == null) return null; // 260420Cl 追加 DisablingOpenGL 時は null を返す (呼び出し側は null チェックするか操作失敗として扱う)
16221647
glControl.MakeCurrent();
16231648
var bmp = new Bitmap(glControl.ClientSize.Width, glControl.Height);
16241649
var data = bmp.LockBits(Rectangle.FromLTRB(0, 0, bmp.Width, bmp.Height), System.Drawing.Imaging.ImageLockMode.WriteOnly, System.Drawing.Imaging.PixelFormat.Format24bppRgb);

0 commit comments

Comments
 (0)