Skip to content

Commit af12b23

Browse files
termoshtttermoshtt
committed
cargo fmt --all
1 parent 55b46e9 commit af12b23

File tree

8 files changed

+201
-75
lines changed

8 files changed

+201
-75
lines changed

‎examples/eig.rs‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,4 +9,4 @@ fn main() {
99
let a_c: Array2<c64> = a.map(|f| c64::new(*f, 0.0));
1010
let av = a_c.dot(&vecs);
1111
println!("AV = \n{:?}", av);
12-
}
12+
}

‎src/eig.rs‎

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,9 @@
11
//! Eigenvalue decomposition for non-symmetric square matrices
22
3-
use ndarray::*;
43
use crate::error::*;
54
use crate::layout::*;
65
use crate::types::*;
6+
use ndarray::*;
77

88
/// Eigenvalue decomposition of general matrix reference
99
pub trait Eig {
@@ -27,7 +27,10 @@ where
2727
let layout = a.square_layout()?;
2828
let (s, t) = unsafe { A::eig(true, layout, a.as_allocated_mut()?)? };
2929
let (n, _) = layout.size();
30-
Ok((ArrayBase::from(s), ArrayBase::from(t).into_shape((n as usize, n as usize)).unwrap()))
30+
Ok((
31+
ArrayBase::from(s),
32+
ArrayBase::from(t).into_shape((n as usize, n as usize)).unwrap(),
33+
))
3134
}
3235
}
3336

@@ -49,4 +52,4 @@ where
4952
let (s, _) = unsafe { A::eig(true, a.square_layout()?, a.as_allocated_mut()?)? };
5053
Ok(ArrayBase::from(s))
5154
}
52-
}
55+
}

‎src/lapack/eig.rs‎

Lines changed: 71 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -17,13 +17,29 @@ pub trait Eig_: Scalar {
1717
macro_rules! impl_eig_complex {
1818
($scalar:ty, $ev:path) => {
1919
impl Eig_ for $scalar {
20-
unsafe fn eig(calc_v: bool, l: MatrixLayout, mut a: &mut [Self]) -> Result<(Vec<Self::Complex>, Vec<Self::Complex>)> {
20+
unsafe fn eig(
21+
calc_v: bool,
22+
l: MatrixLayout,
23+
mut a: &mut [Self],
24+
) -> Result<(Vec<Self::Complex>, Vec<Self::Complex>)> {
2125
let (n, _) = l.size();
2226
let jobvr = if calc_v { b'V' } else { b'N' };
2327
let mut w = vec![Self::Complex::zero(); n as usize];
2428
let mut vl = Vec::new();
2529
let mut vr = vec![Self::Complex::zero(); (n * n) as usize];
26-
let info = $ev(l.lapacke_layout(), b'N', jobvr, n, &mut a, n, &mut w, &mut vl, n, &mut vr, n);
30+
let info = $ev(
31+
l.lapacke_layout(),
32+
b'N',
33+
jobvr,
34+
n,
35+
&mut a,
36+
n,
37+
&mut w,
38+
&mut vl,
39+
n,
40+
&mut vr,
41+
n,
42+
);
2743
into_result(info, (w, vr))
2844
}
2945
}
@@ -33,49 +49,75 @@ macro_rules! impl_eig_complex {
3349
macro_rules! impl_eig_real {
3450
($scalar:ty, $ev:path) => {
3551
impl Eig_ for $scalar {
36-
unsafe fn eig(calc_v: bool, l: MatrixLayout, mut a: &mut [Self]) -> Result<(Vec<Self::Complex>, Vec<Self::Complex>)> {
52+
unsafe fn eig(
53+
calc_v: bool,
54+
l: MatrixLayout,
55+
mut a: &mut [Self],
56+
) -> Result<(Vec<Self::Complex>, Vec<Self::Complex>)> {
3757
let (n, _) = l.size();
3858
let jobvr = if calc_v { b'V' } else { b'N' };
3959
let mut wr = vec![Self::Real::zero(); n as usize];
4060
let mut wi = vec![Self::Real::zero(); n as usize];
4161
let mut vl = Vec::new();
4262
let mut vr = vec![Self::Real::zero(); (n * n) as usize];
43-
let info = $ev(l.lapacke_layout(), b'N', jobvr, n, &mut a, n, &mut wr, &mut wi, &mut vl, n, &mut vr, n);
44-
let w: Vec<Self::Complex> = wr.iter().zip(wi.iter()).map(|(&r, &i)| Self::Complex::new(r, i)).collect();
63+
let info = $ev(
64+
l.lapacke_layout(),
65+
b'N',
66+
jobvr,
67+
n,
68+
&mut a,
69+
n,
70+
&mut wr,
71+
&mut wi,
72+
&mut vl,
73+
n,
74+
&mut vr,
75+
n,
76+
);
77+
let w: Vec<Self::Complex> = wr
78+
.iter()
79+
.zip(wi.iter())
80+
.map(|(&r, &i)| Self::Complex::new(r, i))
81+
.collect();
4582
// If the j-th eigenvalue is real, then
4683
// eigenvector = [ vr[j], vr[j+n], vr[j+2*n], ... ].
4784
//
48-
// If the j-th and (j+1)-st eigenvalues form a complex conjugate pair,
85+
// If the j-th and (j+1)-st eigenvalues form a complex conjugate pair,
4986
// eigenvector(j) = [ vr[j] + i*vr[j+1], vr[j+n] + i*vr[j+n+1], vr[j+2*n] + i*vr[j+2*n+1], ... ] and
5087
// eigenvector(j+1) = [ vr[j] - i*vr[j+1], vr[j+n] - i*vr[j+n+1], vr[j+2*n] - i*vr[j+2*n+1], ... ].
51-
//
88+
//
5289
// Therefore, if eigenvector(j) is written as [ v_{j0}, v_{j1}, v_{j2}, ... ],
53-
// you have to make
90+
// you have to make
5491
// v = vec![ v_{00}, v_{10}, v_{20}, ..., v_{jk}, v_{(j+1)k}, v_{(j+2)k}, ... ] (v.len() = n*n)
5592
// based on wi and vr.
5693
// After that, v is converted to Array2 (see ../eig.rs).
5794
let n = n as usize;
5895
let mut flg = false;
59-
let conj: Vec<i8> = wi.iter().map(|&i| {
60-
if flg {
61-
flg = false;
62-
-1
63-
} else if i != 0.0 {
64-
flg = true;
65-
1
66-
} else {
67-
0
68-
}
69-
}).collect();
70-
let v: Vec<Self::Complex> = (0..n*n).map(|i| {
71-
let j = i % n;
72-
match conj[j] {
73-
1 => Self::Complex::new(vr[i], vr[i+1]),
74-
-1 => Self::Complex::new(vr[i-1], -vr[i]),
75-
_ => Self::Complex::new(vr[i], 0.0),
76-
}
77-
}).collect();
78-
96+
let conj: Vec<i8> = wi
97+
.iter()
98+
.map(|&i| {
99+
if flg {
100+
flg = false;
101+
-1
102+
} else if i != 0.0 {
103+
flg = true;
104+
1
105+
} else {
106+
0
107+
}
108+
})
109+
.collect();
110+
let v: Vec<Self::Complex> = (0..n * n)
111+
.map(|i| {
112+
let j = i % n;
113+
match conj[j] {
114+
1 => Self::Complex::new(vr[i], vr[i + 1]),
115+
-1 => Self::Complex::new(vr[i - 1], -vr[i]),
116+
_ => Self::Complex::new(vr[i], 0.0),
117+
}
118+
})
119+
.collect();
120+
79121
into_result(info, (w, v))
80122
}
81123
}
@@ -85,4 +127,4 @@ macro_rules! impl_eig_real {
85127
impl_eig_real!(f64, lapacke::dgeev);
86128
impl_eig_real!(f32, lapacke::sgeev);
87129
impl_eig_complex!(c64, lapacke::zgeev);
88-
impl_eig_complex!(c32, lapacke::cgeev);
130+
impl_eig_complex!(c32, lapacke::cgeev);

‎src/lapack/mod.rs‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,10 @@ use super::types::*;
2828
pub type Pivot = Vec<i32>;
2929

3030
/// Trait for primitive types which implements LAPACK subroutines
31-
pub trait Lapack: OperatorNorm_ + QR_ + SVD_ + SVDDC_ + Solve_ + Solveh_ + Cholesky_ + Eig_ + Eigh_ + Triangular_ {}
31+
pub trait Lapack:
32+
OperatorNorm_ + QR_ + SVD_ + SVDDC_ + Solve_ + Solveh_ + Cholesky_ + Eig_ + Eigh_ + Triangular_
33+
{
34+
}
3235

3336
impl Lapack for f32 {}
3437
impl Lapack for f64 {}

‎src/lobpcg/eig.rs‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
use super::lobpcg::{lobpcg, LobpcgResult, Order};
2-
use crate::{Lapack, Scalar, generate};
2+
use crate::{generate, Lapack, Scalar};
33
///! Implements truncated eigenvalue decomposition
44
///
55
use ndarray::prelude::*;

‎src/lobpcg/lobpcg.rs‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@ use crate::error::{LinalgError, Result};
66
use crate::{cholesky::*, close_l2, eigh::*, norm::*, triangular::*};
77
use crate::{Lapack, Scalar};
88
use ndarray::prelude::*;
9-
use ndarray::{OwnedRepr, ScalarOperand, Data};
9+
use ndarray::{Data, OwnedRepr, ScalarOperand};
1010
use num_traits::{Float, NumCast};
1111

1212
/// Find largest or smallest eigenvalues
@@ -324,7 +324,8 @@ pub fn lobpcg<
324324
//
325325
// first try to compute the eigenvalue decomposition of the span{R, X, P},
326326
// if this fails (or the algorithm was restarted), then just use span{R, X}
327-
let result = p_ap.as_ref()
327+
let result = p_ap
328+
.as_ref()
328329
.ok_or(LinalgError::Lapack { return_code: 1 })
329330
.and_then(|(active_p, active_ap)| {
330331
let xap = x.t().dot(active_ap);
@@ -352,7 +353,7 @@ pub fn lobpcg<
352353
stack![Axis(1), xp.t(), rp.t(), pp]
353354
]),
354355
size_x,
355-
&order
356+
&order,
356357
)
357358
})
358359
.or_else(|_| {
@@ -362,19 +363,18 @@ pub fn lobpcg<
362363
stack![Axis(0), stack![Axis(1), xax, xar], stack![Axis(1), xar.t(), rar]],
363364
Some(stack![Axis(0), stack![Axis(1), xx, xr], stack![Axis(1), xr.t(), rr]]),
364365
size_x,
365-
&order
366+
&order,
366367
)
367368
});
368369

369-
370370
// update eigenvalues and eigenvectors (lambda is also used in the next iteration)
371371
let eig_vecs;
372372
match result {
373373
Ok((x, y)) => {
374374
lambda = x;
375375
eig_vecs = y;
376-
},
377-
Err(x) => break Err(x)
376+
}
377+
Err(x) => break Err(x),
378378
}
379379

380380
// approximate eigenvector X and conjugate vectors P with solution of eigenproblem
@@ -432,8 +432,8 @@ mod tests {
432432
use super::LobpcgResult;
433433
use super::Order;
434434
use crate::close_l2;
435-
use crate::qr::*;
436435
use crate::generate;
436+
use crate::qr::*;
437437
use ndarray::prelude::*;
438438

439439
/// Test the `sorted_eigen` function

‎src/lobpcg/svd.rs‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@
33
///! This module computes the k largest/smallest singular values/vectors for a dense matrix.
44
use super::lobpcg::{lobpcg, LobpcgResult, Order};
55
use crate::error::Result;
6-
use crate::{Lapack, Scalar, generate};
6+
use crate::{generate, Lapack, Scalar};
77
use ndarray::prelude::*;
88
use ndarray::ScalarOperand;
99
use num_traits::{Float, NumCast};

0 commit comments

Comments
 (0)