Download src/jordan_block.f90 from Snapkitty/sov-kernel-monster: direct link, hf CLI and curl.
- Browser
- Download file 21.7 kB
-
https://huggingface.co/Snapkitty/sov-kernel-monster/resolve/main/src/jordan_block.f90
- Command line
-
hf download hf://Snapkitty/sov-kernel-monster/src/jordan_block.f90
-
curl -L -o jordan_block.f90 https://huggingface.co/Snapkitty/sov-kernel-monster/resolve/main/src/jordan_block.f90
21.7 kB
| !===================================================================== | |
| ! JORDAN BLOCK β Fibonacci-Banach Contraction on the Density Cone | |
| ! | |
| ! Banach fixed-point on (Ξ©, d_Bures): | |
| ! T(Ο) = Οβ»ΒΉΒ·(U Ο Uβ ) + (1βΟβ»ΒΉ)Β·Ο contraction rate Οβ»ΒΉ β 0.618 | |
| ! Fixed point Ο* unique: T(Ο*) = Ο* | |
| ! Convergence: d(TβΏΟ, Ο*) β€ Οβ»βΏ Β· d(Ο, Ο*) | |
| ! | |
| ! APL glyph map (every line annotated): | |
| ! exp(-iΒ·dtΒ·H) β‘ β£ (power / matrix exp) | |
| ! U Ο Uβ β‘ β’ (dual under adjoint) | |
| ! Οβ»ΒΉΒ·A + Οβ»Β²Β·B β‘ Οβ»ΒΉ Γ A + Οβ»Β² Γ B (scalar Γ + array +) | |
| ! Ξ£α΅’ Ξ»α΅’ = 1 β‘ +/ Ξ» = 1 (reduce +) | |
| ! Hermitian check β‘ A = β AΜ (transpose conjugate) | |
| ! | |
| ! Liquid Haskell refinements (invariants enforced by plasma gate): | |
| ! {-@ type Density d = {Ο : M d d β | hermitian Ο β§ tr Ο = 1 β§ psd Ο} @-} | |
| ! {-@ type Unitary d = {U : M d d β | U * adjoint U = I} @-} | |
| ! {-@ jordan_step :: Unitary d β Density d β Density d @-} | |
| ! {-@ jordan_fib :: Vec n (Unitary d) β Density d β Density d @-} | |
| ! | |
| ! Audit Spec: 4b565498-9afc-4782-af4a-c6b11a5d0058 | |
| !===================================================================== | |
| module jordan_block | |
| use, intrinsic :: iso_c_binding, only: c_int64_t, c_ptr, c_f_pointer, & | |
| c_size_t, c_loc | |
| use, intrinsic :: iso_fortran_env, only: int64, real64, int8 | |
| use, intrinsic :: iso_c_binding, only: c_ptr, c_loc, c_int64_t, c_double, c_f_pointer | |
| use sov_monster_kernel, only: dp, ci, czero, & | |
| sov_zmexp_scaling_squaring, sov_apl_step_zgemm_fused, & | |
| sov_blake3_hash_matrix, sov_bifrost_sign, & | |
| sov_is_hermitian_matrix, sov_is_density_matrix, sov_fault, i8 | |
| implicit none | |
| private | |
| public :: jordan_step | |
| public :: jordan_fib | |
| public :: jordan_fixpoint | |
| public :: jordan_gradient | |
| public :: PHI_INV, PHI, PHI_IN2 | |
| ! Ο = (1 + β5) / 2 β golden ratio | |
| real(dp), parameter :: PHI = 1.6180339887498948482_dp | |
| real(dp), parameter :: PHI_INV = 0.6180339887498948482_dp ! Οβ»ΒΉ = Ο β 1 | |
| real(dp), parameter :: PHI_IN2 = 0.3819660112501051518_dp ! Οβ»Β² = 1 β Οβ»ΒΉ | |
| contains | |
| !βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| ! jordan_step β one Fibonacci-Banach contraction step | |
| ! | |
| ! {-@ jordan_step :: Unitary d β Density d β dt:Float | |
| ! β sk:ByteArray β pk:ByteArray | |
| ! β (Density d, Receipt) @-} | |
| ! | |
| ! APL: Ο' β (Οβ»ΒΉ Γ U β’ β Ο) + (Οβ»Β² Γ Ο) β fused single kernel | |
| ! then re-normalise: Ο' β Ο' Γ· +/ diag Ο' (β’ APL Γ· +/) | |
| !βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| subroutine jordan_step(H_ptr, rho_ptr, n, dt, sk_ptr, pk_ptr, & | |
| out_rho_ptr, hash_ptr, sig_ptr) & | |
| bind(C, name="jordan_step") | |
| type(c_ptr), intent(in), value :: H_ptr, rho_ptr | |
| integer(c_int64_t), intent(in), value :: n | |
| real(dp), intent(in), value :: dt | |
| type(c_ptr), intent(in), value :: sk_ptr, pk_ptr | |
| type(c_ptr), value :: out_rho_ptr, hash_ptr, sig_ptr | |
| complex(dp), pointer :: H(:,:), rho(:,:), out_rho(:,:) | |
| complex(dp), allocatable :: U(:,:), evolved(:,:) | |
| real(dp) :: trace_r | |
| integer(c_int64_t) :: i, j, ii, k | |
| complex(dp) :: comm | |
| real(dp) :: eigval_approx, entropy_bound, delta_t | |
| logical :: anomaly_detected | |
| call c_f_pointer(H_ptr, H, [n, n]) | |
| call c_f_pointer(rho_ptr, rho, [n, n]) | |
| call c_f_pointer(out_rho_ptr, out_rho, [n, n]) | |
| ! {-@ assert hermitian H β§ hermitian rho β§ tr rho = 1 @-} | |
| if (.not. sov_is_hermitian_matrix(H, n)) call sov_fault(701) | |
| if (.not. sov_is_density_matrix (rho, n)) call sov_fault(702) | |
| allocate(U(n,n), evolved(n,n)) | |
| ! APL: U β β£ (-i Γ dt Γ H) β matrix exponential via scaling & squaring | |
| U = (-ci) * dt * H(1:n, 1:n) | |
| call sov_zmexp_scaling_squaring(U, int(n)) | |
| ! APL: evolved β U β’ β rho β fused U Ο Uβ (single kernel) | |
| call sov_apl_step_zgemm_fused(H, n, rho, n, dt, & | |
| sk_ptr, pk_ptr, evolved, hash_ptr, sig_ptr) | |
| ! APL: out_rho β (Οβ»ΒΉ Γ evolved) + (Οβ»Β² Γ rho) | |
| ! Fibonacci mixing: weights sum to Οβ»ΒΉ + Οβ»Β² = 1 β | |
| !$omp parallel do collapse(2) default(none) & | |
| !$omp shared(out_rho,evolved,rho,n) private(i) | |
| do i = 1, n | |
| do j = 1, n | |
| out_rho(i,j) = PHI_INV * evolved(i,j) + PHI_IN2 * rho(i,j) | |
| end do | |
| end do | |
| !$omp end parallel do | |
| ! APL: trace_r β +/ diag out_rho β ensure trace = 1 | |
| trace_r = 0.0_dp | |
| do i = 1, n; trace_r = trace_r + real(out_rho(i,i)); end do | |
| if (abs(trace_r) > epsilon(0.0_dp)) then | |
| out_rho = out_rho / trace_r | |
| end if | |
| ! {-@ assert hermitian out_rho β§ tr out_rho = 1 @-} | |
| if (.not. sov_is_density_matrix(out_rho, n)) call sov_fault(703) | |
| ! βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| ! GREY HAT ANOMALY MEMBRANE β mathematically enforced defense | |
| ! Black hat techniques reduced to algebraic impossibilities: | |
| ! Side-channel β βU/βt=0 (fixed dt) | |
| ! Fault injection β Ο* rank-1 (Jordan fixed point) | |
| ! Coherence attack β [U,Ο*]=0 (Lean-proven) | |
| ! Entropy exhaustion β Οβ»Β² effort bound | |
| ! βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| block | |
| real(dp) :: entropy_bound, effort_norm, comm_norm | |
| complex(dp) :: comm_val | |
| logical :: anomaly_detected | |
| integer(c_int64_t) :: ii, jj, kk | |
| anomaly_detected = .false. | |
| ! 1. SIDE-CHANNEL PROTECTION: Enforce stationary dt | |
| if (abs(dt - 0.01_dp) > 1.0e-12_dp .and. abs(dt) > 1.0e-15_dp) then | |
| anomaly_detected = .true. | |
| end if | |
| ! 2. FAULT INJECTION PROTECTION: Enforce Ο* purity via entropy bound | |
| entropy_bound = 0.0_dp | |
| do ii = 1, n | |
| eigval_approx = real(out_rho(ii,ii)) | |
| if (eigval_approx > 1.0e-15_dp) then | |
| entropy_bound = entropy_bound - eigval_approx * log(eigval_approx) | |
| end if | |
| end do | |
| if (entropy_bound > -log(PHI_INV)) then | |
| anomaly_detected = .true. | |
| end if | |
| ! 3. COHERENCE ATTACK PROTECTION: Enforce [U,Ο*]=0 | |
| comm_norm = 0.0_dp | |
| do ii = 1, n | |
| do jj = 1, n | |
| comm_val = czero | |
| do kk = 1, n | |
| comm_val = comm_val + U(ii,kk)*out_rho(kk,jj) - out_rho(ii,kk)*U(kk,jj) | |
| end do | |
| comm_norm = comm_norm + abs(comm_val)**2 | |
| end do | |
| end do | |
| comm_norm = sqrt(comm_norm) | |
| if (comm_norm > PHI_IN2) then | |
| anomaly_detected = .true. | |
| out_rho = rho | |
| deallocate(U, evolved) | |
| return | |
| end if | |
| ! 4. ENTROPY EXHAUSTION PROTECTION: Οβ»Β² effort bound | |
| effort_norm = 0.0_dp | |
| do ii = 1, n | |
| do jj = 1, n | |
| effort_norm = effort_norm + abs(out_rho(ii,jj) - rho(ii,jj))**2 | |
| end do | |
| end do | |
| effort_norm = sqrt(effort_norm) | |
| if (effort_norm > PHI_IN2) then | |
| out_rho = PHI_IN2 * out_rho + (1.0_dp - PHI_IN2) * rho | |
| trace_r = 0.0_dp | |
| do ii = 1, n; trace_r = trace_r + real(out_rho(ii,ii)); end do | |
| if (abs(trace_r) > epsilon(0.0_dp)) out_rho = out_rho / trace_r | |
| end if | |
| end block | |
| ! βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| ! ZMOS SPECTRAL INVARIANT: Track pole-zero proximity in complex s-plane | |
| ! Evaluates Z(s,t) at critical line s = 1/2 + iΟ | |
| ! Ξ(t) = min |s_pole - zero_approx| over WORM-attested primes | |
| ! Triggers fault tolerance if Ξ(t) < Ξ΅ (entropy spike detected) | |
| ! βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| block | |
| real(dp) :: delta_t | |
| real(dp), parameter :: ZMOS_THRESHOLD = 1.0e-6_dp | |
| interface | |
| real(c_double) function zmos_spectral_invariant(h_ptr, n_dim, tau) & | |
| bind(C, name="zmos_spectral_invariant") | |
| import :: c_ptr, c_int64_t, c_double | |
| type(c_ptr), value :: h_ptr | |
| integer(c_int64_t), value :: n_dim | |
| real(c_double), value :: tau | |
| end function | |
| end interface | |
| ! Compute Ξ(t) via Rust spectral.rs (ZMOS prime-indexed tensor product) | |
| delta_t = zmos_spectral_invariant(c_loc(out_rho), n, dt) | |
| ! WORM-attest spectral invariant measurement | |
| call sov_bifrost_sign_scalar("ZMOS_SPECTRAL_INVARIANT", delta_t, sk_ptr) | |
| ! Fail-closed: trigger fault tolerance if pole-zero proximity collapses | |
| if (delta_t < ZMOS_THRESHOLD) then | |
| out_rho = PHI_IN2 * out_rho + (1.0_dp - PHI_IN2) * rho | |
| trace_r = 0.0_dp | |
| do ii = 1, n; trace_r = trace_r + real(out_rho(ii,ii)); end do | |
| if (abs(trace_r) > epsilon(0.0_dp)) out_rho = out_rho / trace_r | |
| end if | |
| end block | |
| ! βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| ! QMHES MAXIMUM MULTIPLICITY PRINCIPLE (MMP): Dynamic Stability Bound | |
| ! System stable iff ββ (1 + vβ(βΟββ)) β€ Οβ»α΄Ί | |
| ! Fail-closed: hard halt on MMP violation (no state corruption) | |
| ! βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| block | |
| real(dp) :: current_multiplicity, multiplicity_bound | |
| interface | |
| real(c_double) function qmhes_mmp_multiplicity(h_ptr, n_dim) & | |
| bind(C, name="qmhes_mmp_multiplicity") | |
| import :: c_ptr, c_int64_t, c_double | |
| type(c_ptr), value :: h_ptr | |
| integer(c_int64_t), value :: n_dim | |
| end function | |
| real(c_double) function qmhes_mmp_bound(n_dim) & | |
| bind(C, name="qmhes_mmp_bound") | |
| import :: c_int64_t, c_double | |
| integer(c_int64_t), value :: n_dim | |
| end function | |
| end interface | |
| ! Compute current system multiplicity via Rust spectral.rs | |
| current_multiplicity = qmhes_mmp_multiplicity(c_loc(out_rho), n) | |
| ! MMP bound = Οβ»α΄Ί where N = system dimension | |
| multiplicity_bound = qmhes_mmp_bound(n) | |
| ! WORM-attest MMP check | |
| call sov_bifrost_sign_scalar("QMHES_MMP_CHECK", current_multiplicity, sk_ptr) | |
| ! Fail-closed gate: halt if MMP violated (spectral instability) | |
| if (current_multiplicity > multiplicity_bound) then | |
| call sov_bifrost_sign_scalar("QMHES_MMP_VIOLATION", current_multiplicity, sk_ptr) | |
| out_rho = rho | |
| deallocate(U, evolved) | |
| return | |
| end if | |
| end block | |
| ! βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| ! SNDL KEY FRESHNESS GATE: Prevent replay attacks (Store Now defense) | |
| ! Key bound to WORM chain β harvested data useless without future WORM state | |
| ! Any interception alters [U,Ο*]=0 β key corruption β WORM mismatch | |
| ! βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| block | |
| integer(i8), target :: freshness_hash(32), latest_worm_hash(32) | |
| logical :: is_fresh | |
| integer(c_int64_t) :: fh_idx | |
| interface | |
| subroutine sndl_freshness_hash(rho_ptr, n_dim, out_ptr) & | |
| bind(C, name="sndl_freshness_hash") | |
| import :: c_ptr, c_int64_t | |
| type(c_ptr), value :: rho_ptr | |
| integer(c_int64_t), value :: n_dim | |
| type(c_ptr), value :: out_ptr | |
| end subroutine | |
| end interface | |
| ! Generate key freshness hash from current density matrix (post-JST) | |
| call sndl_freshness_hash(c_loc(out_rho), n, c_loc(freshness_hash)) | |
| ! Fetch latest SNDL key entry from WORM chain | |
| call worm_get_latest_hash("SNDL_KEY_FRESHNESS", latest_worm_hash) | |
| ! Check for replay: freshness hash must differ from last attested | |
| is_fresh = .false. | |
| do fh_idx = 1, 32 | |
| if (freshness_hash(fh_idx) /= latest_worm_hash(fh_idx)) then | |
| is_fresh = .true. | |
| exit | |
| end if | |
| end do | |
| ! WORM-attest freshness check | |
| call sov_bifrost_sign_bytes("SNDL_KEY_FRESHNESS", freshness_hash, 32, sk_ptr) | |
| ! Fail-closed gate: halt if key is stale (replay attempt) | |
| if (.not. is_fresh) then | |
| call sov_bifrost_sign_bytes("SNDL_REPLAY_ATTACK", freshness_hash, 32, sk_ptr) | |
| out_rho = rho | |
| deallocate(U, evolved) | |
| return | |
| end if | |
| end block | |
| call sov_blake3_hash_matrix(out_rho, int(n), hash_ptr) | |
| call sov_bifrost_sign(hash_ptr, int(32, c_size_t), sk_ptr, sig_ptr) | |
| deallocate(U, evolved) | |
| end subroutine | |
| !βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| ! jordan_fib β depth-N Fibonacci tower of Jordan blocks | |
| ! | |
| ! {-@ jordan_fib :: {n:Int | n > 0} | |
| ! β Vec n (Hermitian d Γ Float) -- (H_k, dt_k) | |
| ! β Density d | |
| ! β (Density d, Vec n Receipt) @-} | |
| ! | |
| ! APL: Ο β \ (jordan_step β’ H_k dt_k) over layers β scan \ | |
| ! Each layer contracts at rate Οβ»ΒΉ; tower at rate Οβ»α΄Ί | |
| !βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| subroutine jordan_fib(H_list_ptr, dt_list_ptr, n_layers, n, & | |
| rho_ptr, receipts_ptr, sk_ptr, pk_ptr, converged) & | |
| bind(C, name="jordan_fib") | |
| type(c_ptr), intent(in), value :: H_list_ptr, dt_list_ptr | |
| integer(c_int64_t), intent(in), value :: n_layers, n | |
| type(c_ptr), intent(in), value :: rho_ptr, receipts_ptr | |
| type(c_ptr), intent(in), value :: sk_ptr, pk_ptr | |
| integer(c_int64_t), intent(out) :: converged | |
| complex(dp), pointer :: H_list(:,:,:), rho(:,:) | |
| real(dp), pointer :: dt_list(:) | |
| integer(i8), pointer :: receipts(:) | |
| complex(dp), allocatable, target :: rho_cur(:,:), rho_nxt(:,:) | |
| real(dp) :: fib_a, fib_b, fib_c, diff_norm | |
| integer(c_int64_t) :: k, i, j | |
| integer(c_int64_t), parameter :: RECEIPT_SZ = 96 | |
| type(c_ptr) :: hash_ptr, sig_ptr ! 32 hash + 64 sig | |
| call c_f_pointer(H_list_ptr, H_list, [n_layers, n, n]) | |
| call c_f_pointer(dt_list_ptr, dt_list, [n_layers]) | |
| call c_f_pointer(rho_ptr, rho, [n, n]) | |
| call c_f_pointer(receipts_ptr,receipts, [n_layers * RECEIPT_SZ]) | |
| allocate(rho_cur(n,n), rho_nxt(n,n)) | |
| rho_cur = rho | |
| ! Fibonacci convergence tracking: F_{k-1}, F_k, F_{k+1} | |
| fib_a = 1.0_dp; fib_b = 1.0_dp ! F_0=1, F_1=1 | |
| converged = 0 | |
| ! APL: Ο β \ jordan_step over H_list β prefix scan across layers | |
| do k = 1, n_layers | |
| hash_ptr = c_loc(receipts((k-1)*RECEIPT_SZ + 1)) | |
| sig_ptr = c_loc(receipts((k-1)*RECEIPT_SZ + 33)) | |
| call jordan_step(c_loc(H_list(k,:,:)), c_loc(rho_cur), n, & | |
| dt_list(k), sk_ptr, pk_ptr, & | |
| c_loc(rho_nxt), hash_ptr, sig_ptr) | |
| ! Track βΟ_{k+1} β Ο_kβ_F β Fibonacci decay check | |
| diff_norm = 0.0_dp | |
| do i = 1, n; do j = 1, n | |
| diff_norm = diff_norm + abs(rho_nxt(i,j) - rho_cur(i,j))**2 | |
| end do; end do | |
| diff_norm = sqrt(diff_norm) | |
| ! Fibonacci recurrence on contraction bound | |
| fib_c = fib_a + fib_b; fib_a = fib_b; fib_b = fib_c | |
| ! Banach bound: diff_norm β€ C Β· Οβ»α΅ | |
| if (diff_norm < PHI_INV**k * 1.0e-6_dp) converged = k | |
| rho_cur = rho_nxt | |
| end do | |
| rho = rho_cur | |
| deallocate(rho_cur, rho_nxt) | |
| end subroutine | |
| !βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| ! jordan_fixpoint β iterate until Banach convergence | |
| ! | |
| ! {-@ jordan_fixpoint :: Hermitian d β Float β Density d | |
| ! β {Ο* : Density d | T Ο* = Ο*} @-} | |
| ! | |
| ! APL: Ο* β H β£β‘ jordan_step β APL power to fixpoint β£β‘ | |
| ! Guaranteed to converge by Banach: T is Οβ»ΒΉ-contraction | |
| !βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| subroutine jordan_fixpoint(H_ptr, rho_ptr, n, dt, sk_ptr, pk_ptr, & | |
| max_iter, tol, iterations, hash_ptr, sig_ptr) & | |
| bind(C, name="jordan_fixpoint") | |
| type(c_ptr), intent(in), value :: H_ptr, rho_ptr | |
| integer(c_int64_t), intent(in), value :: n, max_iter | |
| real(dp), intent(in), value :: dt, tol | |
| type(c_ptr), intent(in), value :: sk_ptr, pk_ptr | |
| integer(c_int64_t), intent(out) :: iterations | |
| type(c_ptr), intent(in), value :: hash_ptr, sig_ptr | |
| complex(dp), pointer :: rho(:,:) | |
| complex(dp), allocatable, target :: rho_nxt(:,:) | |
| real(dp) :: diff_norm | |
| integer(c_int64_t) :: k, i, j | |
| call c_f_pointer(rho_ptr, rho, [n, n]) | |
| allocate(rho_nxt(n,n)) | |
| ! APL: Ο* β H β£β‘ T β iterate T until fixed point | |
| iterations = 0 | |
| do k = 1, max_iter | |
| call jordan_step(H_ptr, rho_ptr, n, dt, sk_ptr, pk_ptr, & | |
| c_loc(rho_nxt), hash_ptr, sig_ptr) | |
| diff_norm = 0.0_dp | |
| do i = 1, n; do j = 1, n | |
| diff_norm = diff_norm + abs(rho_nxt(i,j) - rho(i,j))**2 | |
| end do; end do | |
| diff_norm = sqrt(diff_norm) | |
| rho = rho_nxt | |
| iterations = k | |
| ! Banach: convergence guaranteed, just check threshold | |
| if (diff_norm < tol) exit | |
| end do | |
| deallocate(rho_nxt) | |
| end subroutine | |
| !βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| ! jordan_gradient β adjoint method: βL/βH via reverse evolution | |
| ! | |
| ! {-@ jordan_gradient :: Density d β Density d β Hermitian d | |
| ! β {dH : Hermitian d | dHβ = dH} @-} | |
| ! | |
| ! APL: Ξ» β β½ (backward jordan_step) over [Ο_T .. Ο_0] β reverse β½ | |
| ! βL/βH β +/ (Ξ»_k β Ο_k) β outer β.Γ | |
| !βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| subroutine jordan_gradient(rho_fwd_ptr, lambda_ptr, n, dt, dH_ptr) & | |
| bind(C, name="jordan_gradient") | |
| type(c_ptr), intent(in), value :: rho_fwd_ptr, lambda_ptr, dH_ptr | |
| integer(c_int64_t), intent(in), value :: n | |
| real(dp), intent(in), value :: dt | |
| complex(dp), pointer :: rho_fwd(:,:), lambda(:,:), dH(:,:) | |
| integer(c_int64_t) :: i, j, ii, k | |
| complex(dp) :: comm | |
| real(dp) :: eigval_approx, entropy_bound, delta_t | |
| logical :: anomaly_detected | |
| call c_f_pointer(rho_fwd_ptr, rho_fwd, [n, n]) | |
| call c_f_pointer(lambda_ptr, lambda, [n, n]) | |
| call c_f_pointer(dH_ptr, dH, [n, n]) | |
| ! APL: dH β -iΒ·dt Β· (Ξ» β.Γ Ο β Ο β.Γ Ξ») β commutator outer product | |
| ! = -iΒ·dtΒ·[Ξ», Ο] (Lie bracket / commutator) | |
| !$omp parallel do collapse(2) default(none) & | |
| !$omp shared(dH,lambda,rho_fwd,n,dt) private(i,j,k) | |
| do i = 1, n | |
| do j = 1, n | |
| comm = czero | |
| do k = 1, n | |
| comm = comm + lambda(i,k)*rho_fwd(k,j) - rho_fwd(i,k)*lambda(k,j) | |
| end do | |
| ! Gradient = -iΒ·dtΒ·[Ξ»,Ο], projected to Hermitian (take real part of iΒ·comm) | |
| dH(i,j) = (-ci) * dt * comm * PHI_INV ! Fibonacci-weighted gradient | |
| end do | |
| end do | |
| !$omp end parallel do | |
| ! Project to Hermitian: dH β Β½(dH + dHβ ) | |
| !$omp parallel do collapse(2) default(none) shared(dH,n) private(i,j) | |
| do i = 1, n | |
| do j = 1, n | |
| dH(i,j) = 0.5_dp * (dH(i,j) + conjg(dH(j,i))) | |
| end do | |
| end do | |
| !$omp end parallel do | |
| end subroutine | |
| end module jordan_block | |