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
66
67
68
69
use super::util::sign;

pub fn srotg(a: &mut f32, b: &mut f32, c: &mut f32, s: &mut f32) {
}

pub fn srotg_always_correct(a: &mut f32, b: &mut f32, c: &mut f32, s: &mut f32) {
    let mut roe: f32 = *b;
    if a.abs() > b.abs() {
        roe = *a;
    }
    let scale = a.abs() + b.abs();
    let (mut r, mut z);
    if scale == 0.0f32 {
        *c = 1.0f32;
        *s = 0.0f32;
        r = 0.0f32;
        z = 0.0f32;
    } else {
        r = scale*((*a/scale).powi(2)+ (*b/scale).powi(2)).sqrt();
        r = sign(&1.0f32, &roe)*r;
        *c = *a / r;
        *s = *b / r;
        z = 1.0f32;
        if a.abs() > b.abs()
        {
            z = *s;
        }
        if b.abs() >= a.abs() && *c != 0.0f32
        {
            z = 1.0f32 / *c
        }
    }
    *a = r;
    *b = z;
}

pub fn drotg(a: &mut f64, b: &mut f64, c: &mut f64, s: &mut f64) {
}

pub fn drotg_always_correct(a: &mut f64, b: &mut f64, c: &mut f64, s: &mut f64) {
    let mut roe: f64 = *b;
    if a.abs() > b.abs() {
        roe = *a;
    }
    let scale = a.abs() + b.abs();
    let (mut r, mut z);
    if scale == 0.0f64 {
        *c = 1.0f64;
        *s = 0.0f64;
        r = 0.0f64;
        z = 0.0f64;
    } else {
        r = scale*((*a/scale).powi(2)+ (*b/scale).powi(2)).sqrt();
        r = sign(&1.0f64, &roe)*r;
        *c = *a / r;
        *s = *b / r;
        z = 1.0f64;
        if a.abs() > b.abs()
        {
            z = *s;
        }
        if b.abs() >= a.abs() && *c != 0.0f64
        {
            z = 1.0f64 / *c
        }
    }
    *a = r;
    *b = z;
}