Is there any known inconsistent behavior in C ++ std :: vector multiplication?
I have been translating Python code to C ++.
I have a problem in the below snippet. These two codes should be the same, but they give me different results.
I'm lost. The error starts when k=1
. What's going on, where is the error?
If it matters: I am compiling C ++ code using Eclipse Parallel Mars IDE , Windows 10, MinGW.
Python 3.5:
import numpy as np
a = np.array(([2,-1,0,0],[-1,1.5,-0.5,0],[0,-0.5,0.75,-0.25],[0,0,-0.25,0.25]))
b = np.array(([0,0,0,1]))
for k in range(len(b)-1,-1,-1):
p = 0
print("")
print("k = ", k)
print("b[k] = ", b[k])
for m in range(k+1, len(b)):
print("m = ", m)
print("b[m] = ", b[m])
print("a[k,m] = ",a[k,m])
p += a[k,m] * b[m];
print("p = ", p)
b[k] = (b[k] - p)/a[k,k];
C ++ 11:
#include <iostream>
#include <vector>
int main(){
std::vector< std::vector<double> > a = { {2,-1,0,0}, {-1,1.5,-0.5,0},
{0,-0.5,0.75,-0.25}, {0,0,-0.25,0.25} };
std::vector<double> b = {0,0,0,1};
double p;
for (int k = b.size()-1; k >= 0; --k) {
p = 0;
std::cout << std::endl << "k = " << k << std::endl;
std::cout << "b[k] = " << b[k] << std::endl;
for (size_t m = k+1; m < b.size(); ++m) {
std::cout << "m = " << m << std::endl;
std::cout << "b[m] = " << b[m] << std::endl;
std::cout << "a[k][m] = " << a[k][m] << std::endl;
p += a[k][m] * b[m];
std::cout << "p = " << p << std::endl;
}
b[k] = (b[k] - p)/a[k][k];
}
return 0;
}
Output
Python 3.5:
k = 3 b[k] = 1 k = 2 b[k] = 0 m = 3 b[m] = 4 a[k,m] = -0.25 p = -1.0 k = 1 b[k] = 0 m = 2 b[m] = 1 a[k,m] = -0.5 p = -0.5 m = 3 b[m] = 4 a[k,m] = 0.0 p = -0.5 k = 0 b[k] = 0 m = 1 b[m] = 0 a[k,m] = -1.0 p = 0.0 m = 2 b[m] = 1 a[k,m] = 0.0 p = 0.0 m = 3 b[m] = 4 a[k,m] = 0.0 p = 0.0
C ++ 11:
k = 3
b[k] = 1
k = 2
b[k] = 0
m = 3
b[m] = 4
a[k][m] = -0.25
p = -1
k = 1
b[k] = 0
m = 2
b[m] = 1.33333
a[k][m] = -0.5
p = -0.666667
m = 3
b[m] = 4
a[k][m] = 0
p = -0.666667
k = 0
b[k] = 0
m = 1
b[m] = 0.444444
a[k][m] = -1
p = -0.444444
m = 2
b[m] = 1.33333
a[k][m] = 0
p = -0.444444
m = 3
b[m] = 4
a[k][m] = 0
p = -0.444444
source to share
Your Python code is faulty. This is a truncation of numbers, resulting in integer values where you would expect a float with a fractional component.
In particular, it np.array(([0,0,0,1]))
creates a numpy array with a built-in data type, which means that the b[k]
floating point value is truncated to an integer when assigned . From the docs for a numpy.array()
relatively optional argument dtype
(emphasis mine):
The required data type for the array. If not specified, the type will be determined as the minimum type required to store objects in sequence.
Since you only provided integer values in the input array, numpy indicates that you want to create an array of integers.
The C ++ code is correct.
When I fix your Python code to use floating point values everywhere, the result matches the C ++ version:
import numpy as np
a = np.array(([2.0,-1.0,0.0,0.0],[-1.0,1.5,-0.5,0.0],[0.0,-0.5,0.75,-0.25],[0.0,0.0,-0.25,0.25]))
b = np.array(([0.0,0.0,0.0,1.0]))
for k in range(len(b)-1,-1,-1):
p = 0
print("")
print("k = ", k)
print("b[k] = ", b[k])
for m in range(k+1, len(b)):
print("m = ", m)
print("b[m] = ", b[m])
print("a[k,m] = ",a[k,m])
p += a[k,m] * b[m];
print("p = ", p)
b[k] = (b[k] - p)/a[k,k];
source to share
If you specify dtype
as 'float64'
when creating the array b
:
b = np.array(([0,0,0,1]), dtype='float64')
The result is as follows, this matches your C ++ program:
k = 3 b[k] = 1.0 k = 2 b[k] = 0.0 m = 3 b[m] = 4.0 a[k,m] = -0.25 p = -1.0 k = 1 b[k] = 0.0 m = 2 b[m] = 1.3333333333333333 a[k,m] = -0.5 p = -0.66666666666666663 m = 3 b[m] = 4.0 a[k,m] = 0.0 p = -0.66666666666666663 k = 0 b[k] = 0.0 m = 1 b[m] = 0.44444444444444442 a[k,m] = -1.0 p = -0.44444444444444442 m = 2 b[m] = 1.3333333333333333 a[k,m] = 0.0 p = -0.44444444444444442 m = 3 b[m] = 4.0 a[k,m] = 0.0 p = -0.44444444444444442
source to share