556 lines
19 KiB
JavaScript
556 lines
19 KiB
JavaScript
const FLOAT64_BYTES = Float64Array.BYTES_PER_ELEMENT;
|
|
|
|
export const DEFAULT_HERMITIAN_JACOBI_MAX_SWEEPS = 32;
|
|
export const DEFAULT_MAX_HERMITIAN_EIGEN_WORK_BYTES = 16 * 1024 * 1024;
|
|
export const DEFAULT_MAX_HERMITIAN_JACOBI_PAIR_UPDATES = 1_000_000_000;
|
|
|
|
function assertPositiveInteger(value, label) {
|
|
if (!Number.isSafeInteger(value) || value < 1) {
|
|
throw new TypeError(`${label} must be a positive safe integer.`);
|
|
}
|
|
}
|
|
|
|
function assertMaximumSweeps(value) {
|
|
if (
|
|
!Number.isSafeInteger(value) ||
|
|
value < 0 ||
|
|
value > DEFAULT_HERMITIAN_JACOBI_MAX_SWEEPS
|
|
) {
|
|
throw new RangeError(
|
|
`maximumSweeps must be an integer from 0 through ${DEFAULT_HERMITIAN_JACOBI_MAX_SWEEPS}.`,
|
|
);
|
|
}
|
|
}
|
|
|
|
function checkedProduct(factors, label) {
|
|
let product = 1;
|
|
for (const factor of factors) {
|
|
if (factor !== 0 && product > Number.MAX_SAFE_INTEGER / factor) {
|
|
throw new RangeError(`${label} exceeds the safe integer range.`);
|
|
}
|
|
product *= factor;
|
|
}
|
|
return product;
|
|
}
|
|
|
|
function roundoffGamma(operationCount) {
|
|
const product = operationCount * (Number.EPSILON / 2);
|
|
if (!(product < 1)) {
|
|
throw new RangeError("roundoff bound is undefined for this workload.");
|
|
}
|
|
return product / (1 - product);
|
|
}
|
|
|
|
/**
|
|
* Estimate the bounded typed-array allocation and cyclic-Jacobi work before
|
|
* any decomposition workspace is allocated.
|
|
*/
|
|
export function estimateHermitianEigenWorkload({
|
|
order,
|
|
maximumSweeps = DEFAULT_HERMITIAN_JACOBI_MAX_SWEEPS,
|
|
}) {
|
|
assertPositiveInteger(order, "order");
|
|
assertMaximumSweeps(maximumSweeps);
|
|
|
|
const matrixElementCount = checkedProduct([order, order], "matrix element count");
|
|
const pairCount = checkedProduct([order, order - 1], "Jacobi pair count") / 2;
|
|
const worstCaseRotations = checkedProduct(
|
|
[pairCount, maximumSweeps],
|
|
"Jacobi rotation count",
|
|
);
|
|
// A rotation updates n-2 matrix row pairs and n eigenvector row pairs.
|
|
const pairUpdatesPerRotation = Math.max(0, 2 * order - 2);
|
|
const worstCasePairUpdates = checkedProduct(
|
|
[worstCaseRotations, pairUpdatesPerRotation],
|
|
"Jacobi pair-update count",
|
|
);
|
|
// Two split-complex matrices plus normalized and physical eigenvalue arrays.
|
|
const typedArrayElementCount = checkedProduct(
|
|
[4 * matrixElementCount + 2 * order, FLOAT64_BYTES],
|
|
"Hermitian eigensolver allocation",
|
|
) / FLOAT64_BYTES;
|
|
const allocationBytes = checkedProduct(
|
|
[typedArrayElementCount, FLOAT64_BYTES],
|
|
"Hermitian eigensolver allocation",
|
|
);
|
|
|
|
return Object.freeze({
|
|
order,
|
|
maximumSweeps,
|
|
matrixElementCount,
|
|
pairCount,
|
|
pairUpdatesPerRotation,
|
|
worstCaseRotations,
|
|
worstCasePairUpdates,
|
|
allocationBytes,
|
|
});
|
|
}
|
|
|
|
function inferOrder(real, imaginary, explicitOrder) {
|
|
if (real == null || imaginary == null) {
|
|
throw new TypeError("real and imaginary matrix components are required.");
|
|
}
|
|
const order = explicitOrder ?? Math.sqrt(real.length);
|
|
assertPositiveInteger(order, "order");
|
|
const expectedLength = checkedProduct([order, order], "matrix element count");
|
|
if (real.length !== expectedLength || imaginary.length !== expectedLength) {
|
|
throw new RangeError(
|
|
`Hermitian matrix components must each contain ${expectedLength} values.`,
|
|
);
|
|
}
|
|
return order;
|
|
}
|
|
|
|
function matrixFrobeniusNorm(real, imaginary, order) {
|
|
let sum = 0;
|
|
for (let index = 0; index < order * order; index += 1) {
|
|
sum += real[index] * real[index] + imaginary[index] * imaginary[index];
|
|
}
|
|
return Math.sqrt(sum);
|
|
}
|
|
|
|
function offDiagonalFrobeniusNorm(real, imaginary, order) {
|
|
let upperTriangleSum = 0;
|
|
for (let row = 0; row < order; row += 1) {
|
|
for (let column = row + 1; column < order; column += 1) {
|
|
const index = row * order + column;
|
|
upperTriangleSum +=
|
|
real[index] * real[index] + imaginary[index] * imaginary[index];
|
|
}
|
|
}
|
|
return Math.sqrt(2 * upperTriangleSum);
|
|
}
|
|
|
|
function rotatePair(matrixReal, matrixImaginary, vectorsReal, vectorsImaginary, order, p, q) {
|
|
const pq = p * order + q;
|
|
const offReal = matrixReal[pq];
|
|
const offImaginary = matrixImaginary[pq];
|
|
const magnitude = Math.hypot(offReal, offImaginary);
|
|
if (magnitude === 0) return false;
|
|
|
|
const diagonalP = matrixReal[p * order + p];
|
|
const diagonalQ = matrixReal[q * order + q];
|
|
const halfDifference = (diagonalQ - diagonalP) / 2;
|
|
const tangent = halfDifference === 0
|
|
? 1
|
|
: Math.sign(halfDifference) * magnitude /
|
|
(Math.abs(halfDifference) + Math.hypot(halfDifference, magnitude));
|
|
const cosine = 1 / Math.hypot(1, tangent);
|
|
const sine = tangent * cosine;
|
|
const phaseReal = offReal / magnitude;
|
|
const phaseImaginary = offImaginary / magnitude;
|
|
|
|
for (let row = 0; row < order; row += 1) {
|
|
if (row === p || row === q) continue;
|
|
const rowP = row * order + p;
|
|
const rowQ = row * order + q;
|
|
const pReal = matrixReal[rowP];
|
|
const pImaginary = matrixImaginary[rowP];
|
|
const qReal = matrixReal[rowQ];
|
|
const qImaginary = matrixImaginary[rowQ];
|
|
|
|
const conjugatePhaseTimesQReal =
|
|
phaseReal * qReal + phaseImaginary * qImaginary;
|
|
const conjugatePhaseTimesQImaginary =
|
|
phaseReal * qImaginary - phaseImaginary * qReal;
|
|
const phaseTimesPReal = phaseReal * pReal - phaseImaginary * pImaginary;
|
|
const phaseTimesPImaginary = phaseReal * pImaginary + phaseImaginary * pReal;
|
|
|
|
const nextPReal = cosine * pReal - sine * conjugatePhaseTimesQReal;
|
|
const nextPImaginary = cosine * pImaginary - sine * conjugatePhaseTimesQImaginary;
|
|
const nextQReal = sine * phaseTimesPReal + cosine * qReal;
|
|
const nextQImaginary = sine * phaseTimesPImaginary + cosine * qImaginary;
|
|
|
|
matrixReal[rowP] = nextPReal;
|
|
matrixImaginary[rowP] = nextPImaginary;
|
|
matrixReal[p * order + row] = nextPReal;
|
|
matrixImaginary[p * order + row] = -nextPImaginary;
|
|
matrixReal[rowQ] = nextQReal;
|
|
matrixImaginary[rowQ] = nextQImaginary;
|
|
matrixReal[q * order + row] = nextQReal;
|
|
matrixImaginary[q * order + row] = -nextQImaginary;
|
|
}
|
|
|
|
matrixReal[p * order + p] = diagonalP - tangent * magnitude;
|
|
matrixImaginary[p * order + p] = 0;
|
|
matrixReal[q * order + q] = diagonalQ + tangent * magnitude;
|
|
matrixImaginary[q * order + q] = 0;
|
|
matrixReal[pq] = 0;
|
|
matrixImaginary[pq] = 0;
|
|
matrixReal[q * order + p] = 0;
|
|
matrixImaginary[q * order + p] = 0;
|
|
|
|
for (let row = 0; row < order; row += 1) {
|
|
const rowP = row * order + p;
|
|
const rowQ = row * order + q;
|
|
const pReal = vectorsReal[rowP];
|
|
const pImaginary = vectorsImaginary[rowP];
|
|
const qReal = vectorsReal[rowQ];
|
|
const qImaginary = vectorsImaginary[rowQ];
|
|
const conjugatePhaseTimesQReal =
|
|
phaseReal * qReal + phaseImaginary * qImaginary;
|
|
const conjugatePhaseTimesQImaginary =
|
|
phaseReal * qImaginary - phaseImaginary * qReal;
|
|
const phaseTimesPReal = phaseReal * pReal - phaseImaginary * pImaginary;
|
|
const phaseTimesPImaginary = phaseReal * pImaginary + phaseImaginary * pReal;
|
|
vectorsReal[rowP] = cosine * pReal - sine * conjugatePhaseTimesQReal;
|
|
vectorsImaginary[rowP] = cosine * pImaginary - sine * conjugatePhaseTimesQImaginary;
|
|
vectorsReal[rowQ] = sine * phaseTimesPReal + cosine * qReal;
|
|
vectorsImaginary[rowQ] = sine * phaseTimesPImaginary + cosine * qImaginary;
|
|
}
|
|
return true;
|
|
}
|
|
|
|
function canonicalizeColumnPhases(real, imaginary, order) {
|
|
for (let column = 0; column < order; column += 1) {
|
|
let pivotRow = 0;
|
|
let pivotMagnitude = -1;
|
|
for (let row = 0; row < order; row += 1) {
|
|
const index = row * order + column;
|
|
const magnitude = Math.hypot(real[index], imaginary[index]);
|
|
if (magnitude > pivotMagnitude) {
|
|
pivotMagnitude = magnitude;
|
|
pivotRow = row;
|
|
}
|
|
}
|
|
const pivotIndex = pivotRow * order + column;
|
|
const pivotReal = real[pivotIndex];
|
|
const pivotImaginary = imaginary[pivotIndex];
|
|
const normalizedPivotMagnitude = Math.hypot(pivotReal, pivotImaginary);
|
|
if (!(normalizedPivotMagnitude > 0) || !Number.isFinite(normalizedPivotMagnitude)) {
|
|
throw new RangeError("Jacobi eigenvector phase canonicalization failed.");
|
|
}
|
|
const phaseReal = pivotReal / normalizedPivotMagnitude;
|
|
const phaseImaginary = -pivotImaginary / normalizedPivotMagnitude;
|
|
for (let row = 0; row < order; row += 1) {
|
|
const index = row * order + column;
|
|
const valueReal = real[index];
|
|
const valueImaginary = imaginary[index];
|
|
real[index] = valueReal * phaseReal - valueImaginary * phaseImaginary;
|
|
imaginary[index] = valueReal * phaseImaginary + valueImaginary * phaseReal;
|
|
}
|
|
real[pivotIndex] = normalizedPivotMagnitude;
|
|
imaginary[pivotIndex] = 0;
|
|
}
|
|
}
|
|
|
|
function projectedInputEntry(real, imaginary, order, row, column, scale) {
|
|
if (row === column) return [real[row * order + row] / scale, 0];
|
|
const forward = row * order + column;
|
|
const reverse = column * order + row;
|
|
return [
|
|
(real[forward] + real[reverse]) / (2 * scale),
|
|
(imaginary[forward] - imaginary[reverse]) / (2 * scale),
|
|
];
|
|
}
|
|
|
|
function residualFrobeniusNorm(
|
|
inputReal,
|
|
inputImaginary,
|
|
vectorsReal,
|
|
vectorsImaginary,
|
|
normalizedEigenvalues,
|
|
order,
|
|
scale,
|
|
) {
|
|
if (scale === 0) return 0;
|
|
let residualSquared = 0;
|
|
for (let column = 0; column < order; column += 1) {
|
|
const eigenvalue = normalizedEigenvalues[column];
|
|
for (let row = 0; row < order; row += 1) {
|
|
let productReal = 0;
|
|
let productImaginary = 0;
|
|
for (let inner = 0; inner < order; inner += 1) {
|
|
const [matrixReal, matrixImaginary] = projectedInputEntry(
|
|
inputReal,
|
|
inputImaginary,
|
|
order,
|
|
row,
|
|
inner,
|
|
scale,
|
|
);
|
|
const vectorIndex = inner * order + column;
|
|
const vectorReal = vectorsReal[vectorIndex];
|
|
const vectorImaginary = vectorsImaginary[vectorIndex];
|
|
productReal += matrixReal * vectorReal - matrixImaginary * vectorImaginary;
|
|
productImaginary += matrixReal * vectorImaginary + matrixImaginary * vectorReal;
|
|
}
|
|
const vectorIndex = row * order + column;
|
|
const differenceReal = productReal - eigenvalue * vectorsReal[vectorIndex];
|
|
const differenceImaginary =
|
|
productImaginary - eigenvalue * vectorsImaginary[vectorIndex];
|
|
residualSquared +=
|
|
differenceReal * differenceReal + differenceImaginary * differenceImaginary;
|
|
}
|
|
}
|
|
return Math.sqrt(residualSquared);
|
|
}
|
|
|
|
function orthogonalityError(real, imaginary, order) {
|
|
let errorSquared = 0;
|
|
for (let left = 0; left < order; left += 1) {
|
|
for (let right = 0; right < order; right += 1) {
|
|
let productReal = 0;
|
|
let productImaginary = 0;
|
|
for (let row = 0; row < order; row += 1) {
|
|
const leftIndex = row * order + left;
|
|
const rightIndex = row * order + right;
|
|
productReal +=
|
|
real[leftIndex] * real[rightIndex] +
|
|
imaginary[leftIndex] * imaginary[rightIndex];
|
|
productImaginary +=
|
|
real[leftIndex] * imaginary[rightIndex] -
|
|
imaginary[leftIndex] * real[rightIndex];
|
|
}
|
|
if (left === right) productReal -= 1;
|
|
errorSquared += productReal * productReal + productImaginary * productImaginary;
|
|
}
|
|
}
|
|
return Math.sqrt(errorSquared);
|
|
}
|
|
|
|
/**
|
|
* Deterministic full cyclic Jacobi decomposition of a complex Hermitian
|
|
* matrix. Split-complex matrices are row-major; eigenvectors are returned as
|
|
* row-major columns, so component `row` of eigenvector `column` is at
|
|
* `row * order + column`.
|
|
*
|
|
* The input is projected to its nearest Hermitian matrix only after its
|
|
* anti-Hermitian discrepancy has passed the reported roundoff-scale audit.
|
|
* No diagonal loading or eigenvalue flooring is performed.
|
|
*/
|
|
export function decomposeHermitian({
|
|
real,
|
|
imaginary,
|
|
order: explicitOrder,
|
|
maximumSweeps = DEFAULT_HERMITIAN_JACOBI_MAX_SWEEPS,
|
|
relativeTolerance,
|
|
maximumWorkBytes = DEFAULT_MAX_HERMITIAN_EIGEN_WORK_BYTES,
|
|
maximumPairUpdates = DEFAULT_MAX_HERMITIAN_JACOBI_PAIR_UPDATES,
|
|
}) {
|
|
const order = inferOrder(real, imaginary, explicitOrder);
|
|
const workload = estimateHermitianEigenWorkload({ order, maximumSweeps });
|
|
if (!(maximumWorkBytes >= 0) || !Number.isFinite(maximumWorkBytes)) {
|
|
throw new TypeError("maximumWorkBytes must be finite and nonnegative.");
|
|
}
|
|
if (!(maximumPairUpdates >= 0) || !Number.isFinite(maximumPairUpdates)) {
|
|
throw new TypeError("maximumPairUpdates must be finite and nonnegative.");
|
|
}
|
|
if (workload.allocationBytes > maximumWorkBytes) {
|
|
throw new RangeError(
|
|
`Hermitian eigensolver requires ${workload.allocationBytes} typed-array bytes; limit is ${maximumWorkBytes}.`,
|
|
);
|
|
}
|
|
if (workload.worstCasePairUpdates > maximumPairUpdates) {
|
|
throw new RangeError(
|
|
`Hermitian eigensolver worst-case work is ${workload.worstCasePairUpdates} pair updates; limit is ${maximumPairUpdates}.`,
|
|
);
|
|
}
|
|
|
|
let scale = 0;
|
|
for (let index = 0; index < workload.matrixElementCount; index += 1) {
|
|
const valueReal = real[index];
|
|
const valueImaginary = imaginary[index];
|
|
if (!Number.isFinite(valueReal) || !Number.isFinite(valueImaginary)) {
|
|
throw new TypeError("Hermitian matrix components must be finite.");
|
|
}
|
|
// Component scaling avoids overflow in hypot() for otherwise finite
|
|
// split-complex entries near Number.MAX_VALUE.
|
|
scale = Math.max(scale, Math.abs(valueReal), Math.abs(valueImaginary));
|
|
}
|
|
|
|
const inputHermitianTolerance =
|
|
scale === 0 ? 0 : 32 * Number.EPSILON * Math.max(1, order) * scale;
|
|
let inputHermitianError = 0;
|
|
for (let row = 0; row < order; row += 1) {
|
|
inputHermitianError = Math.max(
|
|
inputHermitianError,
|
|
Math.abs(imaginary[row * order + row]),
|
|
);
|
|
for (let column = row + 1; column < order; column += 1) {
|
|
const forward = row * order + column;
|
|
const reverse = column * order + row;
|
|
inputHermitianError = Math.max(
|
|
inputHermitianError,
|
|
Math.hypot(
|
|
real[forward] - real[reverse],
|
|
imaginary[forward] + imaginary[reverse],
|
|
),
|
|
);
|
|
}
|
|
}
|
|
if (inputHermitianError > inputHermitianTolerance) {
|
|
throw new RangeError(
|
|
`matrix is not Hermitian within roundoff: error ${inputHermitianError}, tolerance ${inputHermitianTolerance}.`,
|
|
);
|
|
}
|
|
|
|
const matrixReal = new Float64Array(workload.matrixElementCount);
|
|
const matrixImaginary = new Float64Array(workload.matrixElementCount);
|
|
const vectorsReal = new Float64Array(workload.matrixElementCount);
|
|
const vectorsImaginary = new Float64Array(workload.matrixElementCount);
|
|
const normalizedEigenvalues = new Float64Array(order);
|
|
const eigenvalues = new Float64Array(order);
|
|
|
|
if (scale > 0) {
|
|
for (let row = 0; row < order; row += 1) {
|
|
const diagonal = row * order + row;
|
|
matrixReal[diagonal] = real[diagonal] / scale;
|
|
for (let column = row + 1; column < order; column += 1) {
|
|
const forward = row * order + column;
|
|
const reverse = column * order + row;
|
|
const projectedReal = (real[forward] + real[reverse]) / (2 * scale);
|
|
const projectedImaginary =
|
|
(imaginary[forward] - imaginary[reverse]) / (2 * scale);
|
|
matrixReal[forward] = projectedReal;
|
|
matrixImaginary[forward] = projectedImaginary;
|
|
matrixReal[reverse] = projectedReal;
|
|
matrixImaginary[reverse] = -projectedImaginary;
|
|
}
|
|
}
|
|
}
|
|
for (let index = 0; index < order; index += 1) {
|
|
vectorsReal[index * order + index] = 1;
|
|
}
|
|
|
|
const matrixNormNormalized = matrixFrobeniusNorm(
|
|
matrixReal,
|
|
matrixImaginary,
|
|
order,
|
|
);
|
|
const effectiveRelativeTolerance = relativeTolerance ?? roundoffGamma(32 * order);
|
|
if (!(effectiveRelativeTolerance > 0) || !Number.isFinite(effectiveRelativeTolerance)) {
|
|
throw new TypeError("relativeTolerance must be finite and positive.");
|
|
}
|
|
const absoluteTolerance = effectiveRelativeTolerance * matrixNormNormalized;
|
|
let offDiagonalNormalized = offDiagonalFrobeniusNorm(
|
|
matrixReal,
|
|
matrixImaginary,
|
|
order,
|
|
);
|
|
let converged = offDiagonalNormalized <= absoluteTolerance;
|
|
let sweeps = 0;
|
|
let rotations = 0;
|
|
const pairTolerance = workload.pairCount === 0
|
|
? 0
|
|
: absoluteTolerance / Math.sqrt(2 * workload.pairCount);
|
|
|
|
while (!converged && sweeps < maximumSweeps) {
|
|
let sweepRotations = 0;
|
|
for (let p = 0; p < order; p += 1) {
|
|
for (let q = p + 1; q < order; q += 1) {
|
|
const pq = p * order + q;
|
|
if (Math.hypot(matrixReal[pq], matrixImaginary[pq]) <= pairTolerance) {
|
|
continue;
|
|
}
|
|
if (
|
|
rotatePair(
|
|
matrixReal,
|
|
matrixImaginary,
|
|
vectorsReal,
|
|
vectorsImaginary,
|
|
order,
|
|
p,
|
|
q,
|
|
)
|
|
) {
|
|
sweepRotations += 1;
|
|
}
|
|
}
|
|
}
|
|
rotations += sweepRotations;
|
|
sweeps += 1;
|
|
offDiagonalNormalized = offDiagonalFrobeniusNorm(
|
|
matrixReal,
|
|
matrixImaginary,
|
|
order,
|
|
);
|
|
converged = offDiagonalNormalized <= absoluteTolerance;
|
|
if (sweepRotations === 0) break;
|
|
}
|
|
|
|
const stableOrder = Array.from({ length: order }, (_, index) => ({
|
|
index,
|
|
value: matrixReal[index * order + index],
|
|
}));
|
|
stableOrder.sort((left, right) => {
|
|
if (left.value > right.value) return -1;
|
|
if (left.value < right.value) return 1;
|
|
return left.index - right.index;
|
|
});
|
|
for (let column = 0; column < order; column += 1) {
|
|
const sourceColumn = stableOrder[column].index;
|
|
const normalizedEigenvalue = stableOrder[column].value;
|
|
normalizedEigenvalues[column] = Object.is(normalizedEigenvalue, -0)
|
|
? 0
|
|
: normalizedEigenvalue;
|
|
const physicalEigenvalue = normalizedEigenvalues[column] * scale;
|
|
if (!Number.isFinite(physicalEigenvalue)) {
|
|
throw new RangeError("Hermitian eigenvalue exceeds the finite Float64 range.");
|
|
}
|
|
eigenvalues[column] = Object.is(physicalEigenvalue, -0) ? 0 : physicalEigenvalue;
|
|
for (let row = 0; row < order; row += 1) {
|
|
const destination = row * order + column;
|
|
const source = row * order + sourceColumn;
|
|
matrixReal[destination] = vectorsReal[source];
|
|
matrixImaginary[destination] = vectorsImaginary[source];
|
|
}
|
|
}
|
|
// A product of unitary Jacobi rotations already has unit-norm columns.
|
|
// Canonicalize phase only; do not hide orthogonality drift by renormalizing.
|
|
canonicalizeColumnPhases(matrixReal, matrixImaginary, order);
|
|
|
|
const residualNormalized = residualFrobeniusNorm(
|
|
real,
|
|
imaginary,
|
|
matrixReal,
|
|
matrixImaginary,
|
|
normalizedEigenvalues,
|
|
order,
|
|
scale,
|
|
);
|
|
const relativeResidual = matrixNormNormalized === 0
|
|
? 0
|
|
: residualNormalized / matrixNormNormalized;
|
|
const relativeOffDiagonal = matrixNormNormalized === 0
|
|
? 0
|
|
: offDiagonalNormalized / matrixNormNormalized;
|
|
const residualTolerance = roundoffGamma(128 * order);
|
|
const orthogonalityTolerance = residualTolerance * Math.sqrt(order);
|
|
const measuredOrthogonalityError = orthogonalityError(
|
|
matrixReal,
|
|
matrixImaginary,
|
|
order,
|
|
);
|
|
const numericallyVerified =
|
|
converged &&
|
|
relativeResidual <= residualTolerance &&
|
|
measuredOrthogonalityError <= orthogonalityTolerance;
|
|
|
|
return {
|
|
order,
|
|
eigenvalues,
|
|
eigenvectors: {
|
|
real: matrixReal,
|
|
imaginary: matrixImaginary,
|
|
layout: "row-major-columns",
|
|
},
|
|
converged,
|
|
sweeps,
|
|
rotations,
|
|
scale,
|
|
inputHermitianError,
|
|
inputHermitianTolerance,
|
|
relativeTolerance: effectiveRelativeTolerance,
|
|
offDiagonalFrobeniusNorm: offDiagonalNormalized * scale,
|
|
relativeOffDiagonal,
|
|
residualFrobeniusNorm: residualNormalized * scale,
|
|
relativeResidual,
|
|
residualTolerance,
|
|
orthogonalityError: measuredOrthogonalityError,
|
|
orthogonalityTolerance,
|
|
numericallyVerified,
|
|
workload,
|
|
};
|
|
}
|