1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
extern crate nalgebra as na;
use na::{DMatrix, DVector, DVectorSlice};
pub struct MGS {
pub basis: DMatrix<f64>,
}
impl MGS {
pub fn orthonormalize(basis: &mut DMatrix<f64>, start: usize) {
let end = basis.ncols();
for i in start..end {
for j in 0..i {
let proj = MGS::project(&basis.column(j), &basis.column(i));
basis.set_column(i, &(basis.column(i) - proj));
}
basis.set_column(i, &basis.column(i).normalize());
}
}
fn project(v1: &DVectorSlice<f64>, v2: &DVectorSlice<f64>) -> DVector<f64> {
let magnitud = v1.dot(&v2) / v1.dot(&v1);
return v1 * magnitud;
}
}
#[cfg(test)]
mod test {
extern crate nalgebra as na;
use na::DMatrix;
#[test]
fn test_gram_schmidt() {
let dim = 20;
let vectors = DMatrix::<f64>::new_random(dim, dim);
fun_test(vectors, 0);
}
#[test]
fn test_start_gram_schmidt() {
let arr = DMatrix::<f64>::from_vec(3, 3, vec![1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 2.0]);
fun_test(arr, 1);
}
fn fun_test(vectors: DMatrix<f64>, start: usize) {
let mut basis = vectors.clone();
super::MGS::orthonormalize(&mut basis, start);
let result = basis.transpose() * &basis;
assert!(result.is_identity(1e-8));
}
}