-
Notifications
You must be signed in to change notification settings - Fork 6
/
xnrm2.cpp
96 lines (88 loc) · 1.74 KB
/
xnrm2.cpp
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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
#include <cmath>
#include "rt_nonfinite.h"
#include "processedout.h"
#include "xnrm2.h"
double b_xnrm2(const double x[1080], int ix0)
{
double y;
double scale;
int k;
double absxk;
double t;
y = 0.0;
scale = 3.3121686421112381E-170;
for (k = ix0; k <= ix0 + 359; k++) {
absxk = std::abs(x[k - 1]);
if (absxk > scale) {
t = scale / absxk;
y = 1.0 + y * t * t;
scale = absxk;
} else {
t = absxk / scale;
y += t * t;
}
}
return scale * std::sqrt(y);
}
double c_xnrm2(int n, const double x[1080], int ix0)
{
double y;
double scale;
int kend;
int k;
double absxk;
double t;
y = 0.0;
if (!(n < 1)) {
if (n == 1) {
y = std::abs(x[ix0 - 1]);
} else {
scale = 3.3121686421112381E-170;
kend = (ix0 + n) - 1;
for (k = ix0; k <= kend; k++) {
absxk = std::abs(x[k - 1]);
if (absxk > scale) {
t = scale / absxk;
y = 1.0 + y * t * t;
scale = absxk;
} else {
t = absxk / scale;
y += t * t;
}
}
y = scale * std::sqrt(y);
}
}
return y;
}
double xnrm2(int n, const emxArray_real_T *x, int ix0)
{
double y;
double scale;
int kend;
int k;
double absxk;
double t;
y = 0.0;
if (!(n < 1)) {
if (n == 1) {
y = std::abs(x->data[ix0 - 1]);
} else {
scale = 3.3121686421112381E-170;
kend = (ix0 + n) - 1;
for (k = ix0; k <= kend; k++) {
absxk = std::abs(x->data[k - 1]);
if (absxk > scale) {
t = scale / absxk;
y = 1.0 + y * t * t;
scale = absxk;
} else {
t = absxk / scale;
y += t * t;
}
}
y = scale * std::sqrt(y);
}
}
return y;
}