use rand::{rngs::StdRng, thread_rng, Rng, SeedableRng};
use std::mem;
extern crate ndarray;
use ndarray::prelude::*;
pub mod error;
use error::SvdLibError;
pub trait SMat {
fn nrows(&self) -> usize;
fn ncols(&self) -> usize;
fn nnz(&self) -> usize;
fn svd_opa(&self, x: &[f64], y: &mut [f64], transposed: bool); }
#[derive(Debug, Clone, PartialEq)]
pub struct SvdRec {
pub d: usize,
pub ut: Array2<f64>,
pub s: Array1<f64>,
pub vt: Array2<f64>,
pub diagnostics: Diagnostics,
}
#[derive(Debug, Clone, PartialEq)]
pub struct Diagnostics {
pub non_zero: usize,
pub dimensions: usize,
pub iterations: usize,
pub transposed: bool,
pub lanczos_steps: usize,
pub ritz_values_stabilized: usize,
pub significant_values: usize,
pub singular_values: usize,
pub end_interval: [f64; 2],
pub kappa: f64,
pub random_seed: u32,
}
#[allow(non_snake_case)]
pub fn svd(A: &dyn SMat) -> Result<SvdRec, SvdLibError> {
svdLAS2(A, 0, 0, &[-1.0e-30, 1.0e-30], 1.0e-6, 0)
}
#[allow(non_snake_case)]
pub fn svd_dim(A: &dyn SMat, dimensions: usize) -> Result<SvdRec, SvdLibError> {
svdLAS2(A, dimensions, 0, &[-1.0e-30, 1.0e-30], 1.0e-6, 0)
}
#[allow(non_snake_case)]
pub fn svd_dim_seed(A: &dyn SMat, dimensions: usize, random_seed: u32) -> Result<SvdRec, SvdLibError> {
svdLAS2(A, dimensions, 0, &[-1.0e-30, 1.0e-30], 1.0e-6, random_seed)
}
#[allow(clippy::redundant_field_names)]
#[allow(non_snake_case)]
pub fn svdLAS2(
A: &dyn SMat,
dimensions: usize,
iterations: usize,
end_interval: &[f64; 2],
kappa: f64,
random_seed: u32,
) -> Result<SvdRec, SvdLibError> {
let random_seed = match random_seed > 0 {
true => random_seed,
false => thread_rng().gen::<_>(),
};
let min_nrows_ncols = A.nrows().min(A.ncols());
let dimensions = match dimensions {
n if n == 0 || n > min_nrows_ncols => min_nrows_ncols,
_ => dimensions,
};
let iterations = match iterations {
n if n == 0 || n > min_nrows_ncols => min_nrows_ncols,
n if n < dimensions => dimensions,
_ => iterations,
};
if dimensions < 2 {
return Err(SvdLibError::Las2Error(format!(
"svdLAS2: insufficient dimensions: {dimensions}"
)));
}
assert!(dimensions > 1 && dimensions <= min_nrows_ncols);
assert!(iterations >= dimensions && iterations <= min_nrows_ncols);
let transposed = A.ncols() as f64 >= (A.nrows() as f64 * 1.2);
let nrows = if transposed { A.ncols() } else { A.nrows() };
let ncols = if transposed { A.nrows() } else { A.ncols() };
let mut wrk = WorkSpace::new(nrows, ncols, transposed, iterations)?;
let mut store = Store::new(ncols)?;
let mut neig = 0;
let steps = lanso(
A,
dimensions,
iterations,
end_interval,
&mut wrk,
&mut neig,
&mut store,
random_seed,
)?;
let kappa = kappa.abs().max(eps34());
let mut R = ritvec(A, dimensions, kappa, &mut wrk, steps, neig, &mut store)?;
if transposed {
mem::swap(&mut R.Ut, &mut R.Vt);
}
Ok(SvdRec {
d: R.d,
ut: Array::from_shape_vec((R.d, R.Ut.cols), R.Ut.value)?,
s: Array::from_shape_vec(R.d, R.S)?,
vt: Array::from_shape_vec((R.d, R.Vt.cols), R.Vt.value)?,
diagnostics: Diagnostics {
non_zero: A.nnz(),
dimensions: dimensions,
iterations: iterations,
transposed: transposed,
lanczos_steps: steps + 1,
ritz_values_stabilized: neig,
significant_values: R.d,
singular_values: R.nsig,
end_interval: *end_interval,
kappa: kappa,
random_seed: random_seed,
},
})
}
const MAXLL: usize = 2;
fn eps34() -> f64 {
f64::EPSILON.powf(0.75) }
struct Store {
n: usize,
vecs: Vec<Vec<f64>>,
}
impl Store {
fn new(n: usize) -> Result<Self, SvdLibError> {
Ok(Self { n, vecs: vec![] })
}
fn storq(&mut self, idx: usize, v: &[f64]) {
while idx + MAXLL >= self.vecs.len() {
self.vecs.push(vec![0.0; self.n]);
}
self.vecs[idx + MAXLL].copy_from_slice(v);
}
fn storp(&mut self, idx: usize, v: &[f64]) {
while idx >= self.vecs.len() {
self.vecs.push(vec![0.0; self.n]);
}
self.vecs[idx].copy_from_slice(v);
}
fn retrq(&mut self, idx: usize) -> &[f64] {
&self.vecs[idx + MAXLL]
}
fn retrp(&mut self, idx: usize) -> &[f64] {
&self.vecs[idx]
}
}
struct WorkSpace {
nrows: usize,
ncols: usize,
transposed: bool,
w0: Vec<f64>, w1: Vec<f64>, w2: Vec<f64>, w3: Vec<f64>, w4: Vec<f64>, w5: Vec<f64>, alf: Vec<f64>, eta: Vec<f64>, oldeta: Vec<f64>, bet: Vec<f64>, bnd: Vec<f64>, ritz: Vec<f64>, temp: Vec<f64>, }
impl WorkSpace {
fn new(nrows: usize, ncols: usize, transposed: bool, iterations: usize) -> Result<Self, SvdLibError> {
Ok(Self {
nrows,
ncols,
transposed,
w0: vec![0.0; ncols],
w1: vec![0.0; ncols],
w2: vec![0.0; ncols],
w3: vec![0.0; ncols],
w4: vec![0.0; ncols],
w5: vec![0.0; ncols],
alf: vec![0.0; iterations],
eta: vec![0.0; iterations],
oldeta: vec![0.0; iterations],
bet: vec![0.0; 1 + iterations],
ritz: vec![0.0; 1 + iterations],
bnd: vec![f64::MAX; 1 + iterations],
temp: vec![0.0; nrows],
})
}
}
#[derive(Debug, Clone, PartialEq)]
struct DMat {
cols: usize,
value: Vec<f64>,
}
#[allow(non_snake_case)]
#[derive(Debug, Clone, PartialEq)]
struct SVDRawRec {
d: usize,
nsig: usize,
Ut: DMat,
S: Vec<f64>,
Vt: DMat,
}
fn compare(computed: f64, expected: f64) -> bool {
(expected - computed).abs() < f64::EPSILON
}
fn insert_sort<T: PartialOrd>(n: usize, array1: &mut [T], array2: &mut [T]) {
for i in 1..n {
for j in (1..i + 1).rev() {
if array1[j - 1] <= array1[j] {
break;
}
array1.swap(j - 1, j);
array2.swap(j - 1, j);
}
}
}
#[allow(non_snake_case)]
#[rustfmt::skip]
fn svd_opb(A: &dyn SMat, x: &[f64], y: &mut [f64], temp: &mut [f64], transposed: bool) {
let nrows = if transposed { A.ncols() } else { A.nrows() };
let ncols = if transposed { A.nrows() } else { A.ncols() };
assert_eq!(x.len(), ncols, "svd_opb: x must be A.ncols() in length, x = {}, A.ncols = {}", x.len(), ncols);
assert_eq!(y.len(), ncols, "svd_opb: y must be A.ncols() in length, y = {}, A.ncols = {}", y.len(), ncols);
assert_eq!(temp.len(), nrows, "svd_opa: temp must be A.nrows() in length, temp = {}, A.nrows = {}", temp.len(), nrows);
A.svd_opa(x, temp, transposed); A.svd_opa(temp, y, !transposed); }
fn svd_daxpy(da: f64, x: &[f64], y: &mut [f64]) {
for (xval, yval) in x.iter().zip(y.iter_mut()) {
*yval += da * xval
}
}
fn svd_idamax(n: usize, x: &[f64]) -> usize {
assert!(n > 0, "svd_idamax: unexpected inputs!");
match n {
1 => 0,
_ => {
let mut imax = 0;
for (i, xval) in x.iter().enumerate().take(n).skip(1) {
if xval.abs() > x[imax].abs() {
imax = i;
}
}
imax
}
}
}
fn svd_fsign(a: f64, b: f64) -> f64 {
match a >= 0.0 && b >= 0.0 || a < 0.0 && b < 0.0 {
true => a,
false => -a,
}
}
fn svd_pythag(a: f64, b: f64) -> f64 {
match a.abs().max(b.abs()) {
n if n > 0.0 => {
let mut p = n;
let mut r = (a.abs().min(b.abs()) / p).powi(2);
let mut t = 4.0 + r;
while !compare(t, 4.0) {
let s = r / t;
let u = 1.0 + 2.0 * s;
p *= u;
r *= (s / u).powi(2);
t = 4.0 + r;
}
p
}
_ => 0.0,
}
}
fn svd_ddot(x: &[f64], y: &[f64]) -> f64 {
x.iter().zip(y).map(|(a, b)| a * b).sum()
}
fn svd_norm(x: &[f64]) -> f64 {
svd_ddot(x, x).sqrt()
}
fn svd_datx(d: f64, x: &[f64], y: &mut [f64]) {
for (i, xval) in x.iter().enumerate() {
y[i] = d * xval;
}
}
fn svd_dscal(d: f64, x: &mut [f64]) {
for elem in x.iter_mut() {
*elem *= d;
}
}
fn svd_dcopy(n: usize, offset: usize, x: &[f64], y: &mut [f64]) {
if n > 0 {
let start = n - 1;
for i in 0..n {
y[offset + start - i] = x[offset + i];
}
}
}
fn imtqlb(n: usize, d: &mut [f64], e: &mut [f64], bnd: &mut [f64]) -> Result<(), SvdLibError> {
if n == 1 {
return Ok(());
}
bnd[0] = 1.0;
let last = n - 1;
for i in 1..=last {
bnd[i] = 0.0;
e[i - 1] = e[i];
}
e[last] = 0.0;
let mut i = 0;
for l in 0..=last {
let mut iteration = 0;
while iteration <= 30 {
let mut m = l;
while m < n {
if m == last {
break;
}
let test = d[m].abs() + d[m + 1].abs();
if compare(test, test + e[m].abs()) {
break; }
m += 1;
}
let mut p = d[l];
let mut f = bnd[l];
if m == l {
let mut exchange = true;
if l > 0 {
i = l;
while i >= 1 && exchange {
if p < d[i - 1] {
d[i] = d[i - 1];
bnd[i] = bnd[i - 1];
i -= 1;
} else {
exchange = false;
}
}
}
if exchange {
i = 0;
}
d[i] = p;
bnd[i] = f;
iteration = 31;
} else {
if iteration == 30 {
return Err(SvdLibError::ImtqlbError(
"imtqlb no convergence to an eigenvalue after 30 iterations".to_string(),
));
}
iteration += 1;
let mut g = (d[l + 1] - p) / (2.0 * e[l]);
let mut r = svd_pythag(g, 1.0);
g = d[m] - p + e[l] / (g + svd_fsign(r, g));
let mut s = 1.0;
let mut c = 1.0;
p = 0.0;
assert!(m > 0, "imtqlb: expected 'm' to be non-zero");
i = m - 1;
let mut underflow = false;
while !underflow && i >= l {
f = s * e[i];
let b = c * e[i];
r = svd_pythag(f, g);
e[i + 1] = r;
if compare(r, 0.0) {
underflow = true;
break;
}
s = f / r;
c = g / r;
g = d[i + 1] - p;
r = (d[i] - g) * s + 2.0 * c * b;
p = s * r;
d[i + 1] = g + p;
g = c * r - b;
f = bnd[i + 1];
bnd[i + 1] = s * bnd[i] + c * f;
bnd[i] = c * bnd[i] - s * f;
if i == 0 {
break;
}
i -= 1;
}
if underflow {
d[i + 1] -= p;
} else {
d[l] -= p;
e[l] = g;
}
e[m] = 0.0;
}
}
}
Ok(())
}
#[allow(non_snake_case)]
fn startv(
A: &dyn SMat,
wrk: &mut WorkSpace,
step: usize,
store: &mut Store,
random_seed: u32,
) -> Result<f64, SvdLibError> {
let mut rnm2 = svd_ddot(&wrk.w0, &wrk.w0);
for id in 0..3 {
if id > 0 || step > 0 || compare(rnm2, 0.0) {
let mut bytes = [0; 32];
for (i, b) in random_seed.to_le_bytes().iter().enumerate() {
bytes[i] = *b;
}
let mut seeded_rng = StdRng::from_seed(bytes);
wrk.w0.fill_with(|| seeded_rng.gen_range(-1.0..1.0));
}
wrk.w3.copy_from_slice(&wrk.w0);
svd_opb(A, &wrk.w3, &mut wrk.w0, &mut wrk.temp, wrk.transposed);
wrk.w3.copy_from_slice(&wrk.w0);
rnm2 = svd_ddot(&wrk.w3, &wrk.w3);
if rnm2 > 0.0 {
break;
}
}
if rnm2 <= 0.0 {
return Err(SvdLibError::StartvError(format!("rnm2 <= 0.0, rnm2 = {}", rnm2)));
}
if step > 0 {
for i in 0..step {
let v = store.retrq(i);
svd_daxpy(-svd_ddot(&wrk.w3, v), v, &mut wrk.w0);
}
svd_daxpy(-svd_ddot(&wrk.w4, &wrk.w0), &wrk.w2, &mut wrk.w0);
wrk.w3.copy_from_slice(&wrk.w0);
rnm2 = match svd_ddot(&wrk.w3, &wrk.w3) {
dot if dot <= f64::EPSILON * rnm2 => 0.0,
dot => dot,
}
}
Ok(rnm2.sqrt())
}
#[allow(non_snake_case)]
fn stpone(A: &dyn SMat, wrk: &mut WorkSpace, store: &mut Store, random_seed: u32) -> Result<(f64, f64), SvdLibError> {
let mut rnm = startv(A, wrk, 0, store, random_seed)?;
if compare(rnm, 0.0) {
return Err(SvdLibError::StponeError("rnm == 0.0".to_string()));
}
svd_datx(rnm.recip(), &wrk.w0, &mut wrk.w1);
svd_dscal(rnm.recip(), &mut wrk.w3);
svd_opb(A, &wrk.w3, &mut wrk.w0, &mut wrk.temp, wrk.transposed);
wrk.alf[0] = svd_ddot(&wrk.w0, &wrk.w3);
svd_daxpy(-wrk.alf[0], &wrk.w1, &mut wrk.w0);
let t = svd_ddot(&wrk.w0, &wrk.w3);
wrk.alf[0] += t;
svd_daxpy(-t, &wrk.w1, &mut wrk.w0);
wrk.w4.copy_from_slice(&wrk.w0);
rnm = svd_norm(&wrk.w4);
let anorm = rnm + wrk.alf[0].abs();
Ok((rnm, f64::EPSILON.sqrt() * anorm))
}
#[allow(non_snake_case)]
#[allow(clippy::too_many_arguments)]
fn lanczos_step(
A: &dyn SMat,
wrk: &mut WorkSpace,
first: usize,
last: usize,
ll: &mut usize,
enough: &mut bool,
rnm: &mut f64,
tol: &mut f64,
store: &mut Store,
) -> Result<usize, SvdLibError> {
let eps1 = f64::EPSILON * (wrk.ncols as f64).sqrt();
let mut j = first;
while j < last {
mem::swap(&mut wrk.w1, &mut wrk.w2);
mem::swap(&mut wrk.w3, &mut wrk.w4);
store.storq(j - 1, &wrk.w2);
if j - 1 < MAXLL {
store.storp(j - 1, &wrk.w4);
}
wrk.bet[j] = *rnm;
if compare(*rnm, 0.0) {
*rnm = startv(A, wrk, j, store, 0)?;
if compare(*rnm, 0.0) {
*enough = true;
}
}
if *enough {
mem::swap(&mut wrk.w1, &mut wrk.w2);
break;
}
svd_datx(rnm.recip(), &wrk.w0, &mut wrk.w1);
svd_dscal(rnm.recip(), &mut wrk.w3);
svd_opb(A, &wrk.w3, &mut wrk.w0, &mut wrk.temp, wrk.transposed);
svd_daxpy(-*rnm, &wrk.w2, &mut wrk.w0);
wrk.alf[j] = svd_ddot(&wrk.w0, &wrk.w3);
svd_daxpy(-wrk.alf[j], &wrk.w1, &mut wrk.w0);
if j <= MAXLL && wrk.alf[j - 1].abs() > 4.0 * wrk.alf[j].abs() {
*ll = j;
}
for i in 0..(j - 1).min(*ll) {
let v1 = store.retrp(i);
let t = svd_ddot(v1, &wrk.w0);
let v2 = store.retrq(i);
svd_daxpy(-t, v2, &mut wrk.w0);
wrk.eta[i] = eps1;
wrk.oldeta[i] = eps1;
}
let t = svd_ddot(&wrk.w0, &wrk.w4);
svd_daxpy(-t, &wrk.w2, &mut wrk.w0);
if wrk.bet[j] > 0.0 {
wrk.bet[j] += t;
}
let t = svd_ddot(&wrk.w0, &wrk.w3);
svd_daxpy(-t, &wrk.w1, &mut wrk.w0);
wrk.alf[j] += t;
wrk.w4.copy_from_slice(&wrk.w0);
*rnm = svd_norm(&wrk.w4);
let anorm = wrk.bet[j] + wrk.alf[j].abs() + *rnm;
*tol = f64::EPSILON.sqrt() * anorm;
ortbnd(wrk, j, *rnm, eps1);
purge(wrk.ncols, *ll, wrk, j, rnm, *tol, store);
if *rnm <= *tol {
*rnm = 0.0;
}
j += 1;
}
Ok(j)
}
fn purge(n: usize, ll: usize, wrk: &mut WorkSpace, step: usize, rnm: &mut f64, tol: f64, store: &mut Store) {
if step < ll + 2 {
return;
}
let reps = f64::EPSILON.sqrt();
let eps1 = f64::EPSILON * (n as f64).sqrt();
let k = svd_idamax(step - (ll + 1), &wrk.eta) + ll;
if wrk.eta[k].abs() > reps {
let reps1 = eps1 / reps;
let mut iteration = 0;
let mut flag = true;
while iteration < 2 && flag {
if *rnm > tol {
let mut tq = 0.0;
let mut tr = 0.0;
for i in ll..step {
let v = store.retrq(i);
let t = svd_ddot(v, &wrk.w3);
tq += t.abs();
svd_daxpy(-t, v, &mut wrk.w1);
let t = svd_ddot(v, &wrk.w4);
tr += t.abs();
svd_daxpy(-t, v, &mut wrk.w0);
}
wrk.w3.copy_from_slice(&wrk.w1);
let t = svd_ddot(&wrk.w0, &wrk.w3);
tr += t.abs();
svd_daxpy(-t, &wrk.w1, &mut wrk.w0);
wrk.w4.copy_from_slice(&wrk.w0);
*rnm = svd_norm(&wrk.w4);
if tq <= reps1 && tr <= *rnm * reps1 {
flag = false;
}
}
iteration += 1;
}
for i in ll..=step {
wrk.eta[i] = eps1;
wrk.oldeta[i] = eps1;
}
}
}
fn ortbnd(wrk: &mut WorkSpace, step: usize, rnm: f64, eps1: f64) {
if step < 1 {
return;
}
if !compare(rnm, 0.0) && step > 1 {
wrk.oldeta[0] =
(wrk.bet[1] * wrk.eta[1] + (wrk.alf[0] - wrk.alf[step]) * wrk.eta[0] - wrk.bet[step] * wrk.oldeta[0]) / rnm
+ eps1;
if step > 2 {
for i in 1..=step - 2 {
wrk.oldeta[i] = (wrk.bet[i + 1] * wrk.eta[i + 1]
+ (wrk.alf[i] - wrk.alf[step]) * wrk.eta[i]
+ wrk.bet[i] * wrk.eta[i - 1]
- wrk.bet[step] * wrk.oldeta[i])
/ rnm
+ eps1;
}
}
}
wrk.oldeta[step - 1] = eps1;
mem::swap(&mut wrk.oldeta, &mut wrk.eta);
wrk.eta[step] = eps1;
}
fn error_bound(
enough: &mut bool,
endl: f64,
endr: f64,
ritz: &mut [f64],
bnd: &mut [f64],
step: usize,
tol: f64,
) -> usize {
assert!(step > 0, "error_bound: expected 'step' to be non-zero");
let mid = svd_idamax(step + 1, bnd);
let mut i = ((step + 1) + (step - 1)) / 2;
while i > mid + 1 {
if (ritz[i - 1] - ritz[i]).abs() < eps34() * ritz[i].abs() && bnd[i] > tol && bnd[i - 1] > tol {
bnd[i - 1] = (bnd[i].powi(2) + bnd[i - 1].powi(2)).sqrt();
bnd[i] = 0.0;
}
i -= 1;
}
let mut i = ((step + 1) - (step - 1)) / 2;
while i + 1 < mid {
if (ritz[i + 1] - ritz[i]).abs() < eps34() * ritz[i].abs() && bnd[i] > tol && bnd[i + 1] > tol {
bnd[i + 1] = (bnd[i].powi(2) + bnd[i + 1].powi(2)).sqrt();
bnd[i] = 0.0;
}
i += 1;
}
let mut neig = 0;
let mut gapl = ritz[step] - ritz[0];
for i in 0..=step {
let mut gap = gapl;
if i < step {
gapl = ritz[i + 1] - ritz[i];
}
gap = gap.min(gapl);
if gap > bnd[i] {
bnd[i] *= bnd[i] / gap;
}
if bnd[i] <= 16.0 * f64::EPSILON * ritz[i].abs() {
neig += 1;
if !*enough {
*enough = endl < ritz[i] && ritz[i] < endr;
}
}
}
neig
}
fn imtql2(nm: usize, n: usize, d: &mut [f64], e: &mut [f64], z: &mut [f64]) -> Result<(), SvdLibError> {
if n == 1 {
return Ok(());
}
assert!(n > 1, "imtql2: expected 'n' to be > 1");
let last = n - 1;
for i in 1..n {
e[i - 1] = e[i];
}
e[last] = 0.0;
let nnm = n * nm;
for l in 0..n {
let mut iteration = 0;
while iteration <= 30 {
let mut m = l;
while m < n {
if m == last {
break;
}
let test = d[m].abs() + d[m + 1].abs();
if compare(test, test + e[m].abs()) {
break; }
m += 1;
}
if m == l {
break;
}
if iteration == 30 {
return Err(SvdLibError::Imtql2Error(
"imtql2 no convergence to an eigenvalue after 30 iterations".to_string(),
));
}
iteration += 1;
let mut g = (d[l + 1] - d[l]) / (2.0 * e[l]);
let mut r = svd_pythag(g, 1.0);
g = d[m] - d[l] + e[l] / (g + svd_fsign(r, g));
let mut s = 1.0;
let mut c = 1.0;
let mut p = 0.0;
assert!(m > 0, "imtql2: expected 'm' to be non-zero");
let mut i = m - 1;
let mut underflow = false;
while !underflow && i >= l {
let mut f = s * e[i];
let b = c * e[i];
r = svd_pythag(f, g);
e[i + 1] = r;
if compare(r, 0.0) {
underflow = true;
} else {
s = f / r;
c = g / r;
g = d[i + 1] - p;
r = (d[i] - g) * s + 2.0 * c * b;
p = s * r;
d[i + 1] = g + p;
g = c * r - b;
for k in (0..nnm).step_by(n) {
let index = k + i;
f = z[index + 1];
z[index + 1] = s * z[index] + c * f;
z[index] = c * z[index] - s * f;
}
if i == 0 {
break;
}
i -= 1;
}
}
if underflow {
d[i + 1] -= p;
} else {
d[l] -= p;
e[l] = g;
}
e[m] = 0.0;
}
}
for l in 1..n {
let i = l - 1;
let mut k = i;
let mut p = d[i];
for (j, item) in d.iter().enumerate().take(n).skip(l) {
if *item < p {
k = j;
p = *item;
}
}
if k != i {
d[k] = d[i];
d[i] = p;
for j in (0..nnm).step_by(n) {
z.swap(j + i, j + k);
}
}
}
Ok(())
}
fn rotate_array(a: &mut [f64], x: usize) {
let n = a.len();
let mut j = 0;
let mut start = 0;
let mut t1 = a[0];
for _ in 0..n {
j = match j >= x {
true => j - x,
false => j + n - x,
};
let t2 = a[j];
a[j] = t1;
if j == start {
j += 1;
start = j;
t1 = a[j];
} else {
t1 = t2;
}
}
}
#[allow(non_snake_case)]
fn ritvec(
A: &dyn SMat,
dimensions: usize,
kappa: f64,
wrk: &mut WorkSpace,
steps: usize,
neig: usize,
store: &mut Store,
) -> Result<SVDRawRec, SvdLibError> {
let js = steps + 1;
let jsq = js * js;
let mut s = vec![0.0; jsq];
for i in (0..jsq).step_by(js + 1) {
s[i] = 1.0;
}
let mut Vt = DMat {
cols: wrk.ncols,
value: vec![0.0; wrk.ncols * dimensions],
};
svd_dcopy(js, 0, &wrk.alf, &mut Vt.value);
svd_dcopy(steps, 1, &wrk.bet, &mut wrk.w5);
imtql2(js, js, &mut Vt.value, &mut wrk.w5, &mut s)?;
let mut nsig = 0;
let mut x = 0;
let mut id2 = jsq - js;
for k in 0..js {
if wrk.bnd[k] <= kappa * wrk.ritz[k].abs() && k + 1 > js - neig {
x = match x {
0 => dimensions - 1,
_ => x - 1,
};
let offset = x * Vt.cols;
Vt.value[offset..offset + Vt.cols].fill(0.0);
let mut idx = id2 + js;
for i in 0..js {
idx -= js;
if s[idx] != 0.0 {
for (j, item) in store.retrq(i).iter().enumerate().take(Vt.cols) {
Vt.value[j + offset] += s[idx] * item;
}
}
}
nsig += 1;
}
id2 += 1;
}
if x > 0 {
rotate_array(&mut Vt.value, x * Vt.cols);
}
let d = dimensions.min(nsig);
let mut S = vec![0.0; d];
let mut Ut = DMat {
cols: wrk.nrows,
value: vec![0.0; wrk.nrows * d],
};
Vt.value.resize(Vt.cols * d, 0.0);
let mut tmp_vec = vec![0.0; Vt.cols];
for (i, sval) in S.iter_mut().enumerate() {
let vt_offset = i * Vt.cols;
let ut_offset = i * Ut.cols;
let vt_vec = &Vt.value[vt_offset..vt_offset + Vt.cols];
let ut_vec = &mut Ut.value[ut_offset..ut_offset + Ut.cols];
svd_opb(A, vt_vec, &mut tmp_vec, &mut wrk.temp, wrk.transposed);
let t = svd_ddot(vt_vec, &tmp_vec);
*sval = t.sqrt();
svd_daxpy(-t, vt_vec, &mut tmp_vec);
wrk.bnd[js] = svd_norm(&tmp_vec) * sval.recip();
A.svd_opa(vt_vec, ut_vec, wrk.transposed);
svd_dscal(sval.recip(), ut_vec);
}
Ok(SVDRawRec {
d,
nsig,
Ut,
S,
Vt,
})
}
#[allow(non_snake_case)]
#[allow(clippy::too_many_arguments)]
fn lanso(
A: &dyn SMat,
dim: usize,
iterations: usize,
end_interval: &[f64; 2],
wrk: &mut WorkSpace,
neig: &mut usize,
store: &mut Store,
random_seed: u32,
) -> Result<usize, SvdLibError> {
let (endl, endr) = (end_interval[0], end_interval[1]);
let rnm_tol = stpone(A, wrk, store, random_seed)?;
let mut rnm = rnm_tol.0;
let mut tol = rnm_tol.1;
let eps1 = f64::EPSILON * (wrk.ncols as f64).sqrt();
wrk.eta[0] = eps1;
wrk.oldeta[0] = eps1;
let mut ll = 0;
let mut first = 1;
let mut last = iterations.min(dim.max(8) + dim);
let mut enough = false;
let mut j = 0;
let mut intro = 0;
while !enough {
if rnm <= tol {
rnm = 0.0;
}
let steps = lanczos_step(A, wrk, first, last, &mut ll, &mut enough, &mut rnm, &mut tol, store)?;
j = match enough {
true => steps - 1,
false => last - 1,
};
first = j + 1;
wrk.bet[first] = rnm;
let mut l = 0;
for _ in 0..j {
if l > j {
break;
}
let mut i = l;
while i <= j {
if compare(wrk.bet[i + 1], 0.0) {
break;
}
i += 1;
}
i = i.min(j);
let sz = i - l;
svd_dcopy(sz + 1, l, &wrk.alf, &mut wrk.ritz);
svd_dcopy(sz, l + 1, &wrk.bet, &mut wrk.w5);
imtqlb(sz + 1, &mut wrk.ritz[l..], &mut wrk.w5[l..], &mut wrk.bnd[l..])?;
for m in l..=i {
wrk.bnd[m] = rnm * wrk.bnd[m].abs();
}
l = i + 1;
}
insert_sort(j + 1, &mut wrk.ritz, &mut wrk.bnd);
*neig = error_bound(&mut enough, endl, endr, &mut wrk.ritz, &mut wrk.bnd, j, tol);
if *neig < dim {
if *neig == 0 {
last = first + 9;
intro = first;
} else {
last = first + 3.max(1 + ((j - intro) * (dim - *neig)) / *neig);
}
last = last.min(iterations);
} else {
enough = true
}
enough = enough || first >= iterations;
}
store.storq(j, &wrk.w1);
Ok(j)
}
impl SvdRec {
pub fn recompose(&self) -> Array2<f64> {
let sdiag = Array2::from_diag(&self.s);
self.ut.t().dot(&sdiag).dot(&self.vt)
}
}
#[rustfmt::skip]
impl SMat for nalgebra_sparse::csc::CscMatrix<f64> {
fn nrows(&self) -> usize { self.nrows() }
fn ncols(&self) -> usize { self.ncols() }
fn nnz(&self) -> usize { self.nnz() }
fn svd_opa(&self, x: &[f64], y: &mut [f64], transposed: bool) {
let nrows = if transposed { self.ncols() } else { self.nrows() };
let ncols = if transposed { self.nrows() } else { self.ncols() };
assert_eq!(x.len(), ncols, "svd_opa: x must be A.ncols() in length, x = {}, A.ncols = {}", x.len(), ncols);
assert_eq!(y.len(), nrows, "svd_opa: y must be A.nrows() in length, y = {}, A.nrows = {}", y.len(), nrows);
let (major_offsets, minor_indices, values) = self.csc_data();
y.fill(0.0);
if transposed {
for (i, yval) in y.iter_mut().enumerate() {
for j in major_offsets[i]..major_offsets[i + 1] {
*yval += values[j] * x[minor_indices[j]];
}
}
} else {
for (i, xval) in x.iter().enumerate() {
for j in major_offsets[i]..major_offsets[i + 1] {
y[minor_indices[j]] += values[j] * xval;
}
}
}
}
}
#[rustfmt::skip]
impl SMat for nalgebra_sparse::csr::CsrMatrix<f64> {
fn nrows(&self) -> usize { self.nrows() }
fn ncols(&self) -> usize { self.ncols() }
fn nnz(&self) -> usize { self.nnz() }
fn svd_opa(&self, x: &[f64], y: &mut [f64], transposed: bool) {
let nrows = if transposed { self.ncols() } else { self.nrows() };
let ncols = if transposed { self.nrows() } else { self.ncols() };
assert_eq!(x.len(), ncols, "svd_opa: x must be A.ncols() in length, x = {}, A.ncols = {}", x.len(), ncols);
assert_eq!(y.len(), nrows, "svd_opa: y must be A.nrows() in length, y = {}, A.nrows = {}", y.len(), nrows);
let (major_offsets, minor_indices, values) = self.csr_data();
y.fill(0.0);
if !transposed {
for (i, yval) in y.iter_mut().enumerate() {
for j in major_offsets[i]..major_offsets[i + 1] {
*yval += values[j] * x[minor_indices[j]];
}
}
} else {
for (i, xval) in x.iter().enumerate() {
for j in major_offsets[i]..major_offsets[i + 1] {
y[minor_indices[j]] += values[j] * xval;
}
}
}
}
}
#[rustfmt::skip]
impl SMat for nalgebra_sparse::coo::CooMatrix<f64> {
fn nrows(&self) -> usize { self.nrows() }
fn ncols(&self) -> usize { self.ncols() }
fn nnz(&self) -> usize { self.nnz() }
fn svd_opa(&self, x: &[f64], y: &mut [f64], transposed: bool) {
let nrows = if transposed { self.ncols() } else { self.nrows() };
let ncols = if transposed { self.nrows() } else { self.ncols() };
assert_eq!(x.len(), ncols, "svd_opa: x must be A.ncols() in length, x = {}, A.ncols = {}", x.len(), ncols);
assert_eq!(y.len(), nrows, "svd_opa: y must be A.nrows() in length, y = {}, A.nrows = {}", y.len(), nrows);
y.fill(0.0);
if transposed {
for (i, j, v) in self.triplet_iter() {
y[j] += v * x[i];
}
} else {
for (i, j, v) in self.triplet_iter() {
y[i] += v * x[j];
}
}
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
#[rustfmt::skip]
fn basic_2x2() {
let mut coo = nalgebra_sparse::coo::CooMatrix::<f64>::new(2, 2);
coo.push(0, 0, 4.0);
coo.push(1, 0, 3.0);
coo.push(1, 1, -5.0);
let csc = nalgebra_sparse::csc::CscMatrix::from(&coo);
let svd = svd(&csc).unwrap();
assert_eq!(svd.d, svd.ut.nrows());
assert_eq!(svd.d, svd.s.dim());
assert_eq!(svd.d, svd.vt.nrows());
let matrix_approximation = svd.ut.t().dot(&Array2::from_diag(&svd.s)).dot(&svd.vt);
assert_eq!(svd.recompose(), matrix_approximation);
let epsilon = 1.0e-12;
assert_eq!(svd.d, 2);
assert!((matrix_approximation[[0, 0]] - 4.0).abs() < epsilon);
assert!((matrix_approximation[[0, 1]] - 0.0).abs() < epsilon);
assert!((matrix_approximation[[1, 0]] - 3.0).abs() < epsilon);
assert!((matrix_approximation[[1, 1]] - -5.0).abs() < epsilon);
assert!((svd.s[0] - 6.3245553203368).abs() < epsilon);
assert!((svd.s[1] - 3.1622776601684).abs() < epsilon);
}
#[test]
#[rustfmt::skip]
fn identity_3x3() {
let mut coo = nalgebra_sparse::coo::CooMatrix::<f64>::new(3, 3);
coo.push(0, 0, 1.0);
coo.push(1, 1, 1.0);
coo.push(2, 2, 1.0);
let csc = nalgebra_sparse::csc::CscMatrix::from(&coo);
let svd = svd(&csc).unwrap();
assert_eq!(svd.d, svd.ut.nrows());
assert_eq!(svd.d, svd.s.dim());
assert_eq!(svd.d, svd.vt.nrows());
let epsilon = 1.0e-12;
assert_eq!(svd.d, 1);
assert!((svd.s[0] - 1.0).abs() < epsilon);
}
}