-
c++ / Fortran90 issue
Ok, I am making a c++ object which does the job of this old Fortran code I have, wanting to integrate the functionality into my fresh shiney new c++ skillset. Massive issue is that though the code I am using looks identical (now, anyway I wrote it all c'd up so it looks like the fortran code as much as possible) the numbers come out all funky.
Is this a common problem, I tried calling the cmath cosl function with degrees and radians to check if that's the problem (it's not) the problem shows itself when a trig function is called.
I supply both sections of code for expert perusal and opinion:
c++ :
Code:
Matrix3ld transformMatrix::makeTransformationMatrix(Vector3ld angles, Vector3ld lengths)
{
cout.precision(19) ;
long double alpha, beta, gamma, a, ai, aj, ak,
b, bi, bj, bk,
c, ci, cj, ck ;
Matrix3ld returnMatrix ;
long double pi = 3.1415926535897932384626433832795 ;
alpha = pi*angles(0)/180.0 ;
beta = pi*angles(1)/180.0 ;
gamma = pi*angles(2)/180.0 ;
//alpha = angles(0) ;
//beta = angles(1) ;
//gamma = angles(2) ;
a = lengths(0) ;
b = lengths(1) ;
c = lengths(2) ;
ai = a ;
aj = 0 ;
ak = 0 ;
bi = b*cosl(gamma) ;
bj = b*sinl(gamma) ;
bk = 0 ;
ci = c*cosl(beta) ;
cj = (b*c*cosl(alpha) - bi*ci)/bj ;
ck = sqrtl(c*c-ci*ci-cj*cj) ;
returnMatrix(0,0) = ai ;
returnMatrix(0,1) = aj ;
returnMatrix(0,2) = ak ;
returnMatrix(1,0) = bi ;
returnMatrix(1,1) = bj ;
returnMatrix(1,2) = bk ;
returnMatrix(2,0) = ci ;
returnMatrix(2,1) = cj ;
returnMatrix(2,2) = ck ;
cout << "alpha: " << alpha*pi/180 << endl ;
cout << "beta: " << beta*pi/180 << endl ;
cout << "gamma: " << gamma*pi/180 << endl ;
cout << endl ;
cout << "a: " << a << endl ;
cout << "b: " << b << endl ;
cout << "c: " << c << endl ;
cout << endl ;
cout << "ai: " << ai << endl ;
cout << "bi: " << bi << ", bj: " << bj << endl ;
cout << "ci: " << ci << ", cj: " << cj << ", ck: " << ck << endl ;
cout << endl ;
return returnMatrix ;
}
and the fortran :
Code:
alpha_r = alpha*degtorad
beta_r = beta*degtorad
gamma_r = gamma*degtorad
Write(*,*) 'alpha: ', alpha_r
Write(*,*) 'beta: ', beta_r
Write(*,*) 'gamma: ', gamma_r
Write(*,*)
Write(*,*) 'a: ', a
Write(*,*) 'b: ', b
Write(*,*) 'c: ', c
Write(*,*)
ai = a
aj = 0
ak = 0
bi = b*cos(gamma_r)
bj = b*sin(gamma_r)
bk = 0
ci = c*cos(beta_r)
cj = (b*c*cos(alpha_r)-bi*ci)/bj
ck = sqrt(c**2-ci**2-cj**2)
lij%comp = 0
lij%comp(1,1) = ai
lij%comp(1,2) = bi
lij%comp(1,3) = ci
lij%comp(2,1) = aj
lij%comp(2,2) = bj
lij%comp(2,3) = cj
lij%comp(3,1) = ak
lij%comp(3,2) = bk
lij%comp(3,3) = ck
Write(*,*) 'ai: ', ai
Write(*,*) 'bi: ', bi, ', bj: ', bj
Write(*,*) 'ci: ', ci, ', cj: ', cj, ', ck: ', ck
Write(*,*)
Write(*,*) 'matrix lij'
Write(*,*) '----------'
Write(*,*) lij
Write(*,*)
lijinv=matrix_inverse(lij)
Write(*,*) 'matrix lijinv'
Write(*,*) '-------------'
Write(*,*) lijinv
Write(*,*)
please, please come up with some thing for me, the numbers are out by sign and magnitude (several orders of ten out) If this is a common problem then I am sorry but I can't find anything on tinterwebs. :(
Cheers for lookin if you pass me by ;)
-
OMG I am coming across as such an idiot in this forum, leaving code up due to need for ridicule at idocy........ can you spot my glaring error?
/headpalm
HOWEVER even though I built the matrix incorrectly (fixed now) the numbers are still diffferent and the form of the matrices are different also, as I shall demonstrate with the output now
output discrepancies
--------------------------
c++ : bi = 6.605e-16 ; cj = 1.7469e-15 upper diagonal matrix
fortran : bi = 5.035e-7 ; cj = -6.62032e-7 lower diagonal matrix
the matrix's I am using are transformations between orthogonal space and a crystal space so are fairly important to the operation of the program. I mean it's obvious that fortran references its matrix elements in the wrong way, but the numbers are all wrong still....
thx again guys who thought it was worth a look ;)
-
Hey guys,
after poking at my keyboard ineffectually for the past few hours trying to resolve this I am no further, any help at all would be appreciated. I have figured out that my c++ code is simpy getting closer to zero due to the number fed into the cosl function is closer to (should be exactly in this case but its not going to be that for every instance) PI/2 than the fortran version. This still does not explain the wildly different results I am getting when I multiply a vector (math name) by the matrix (again used in a math context).
I have read through the header files (eigen library) which define the operation of vector by matrix multiplication and checked that against the function in my fortran program which performs the same act, both match up (by multiplying the correct elements together) in functionality but give different results.
Please, am I missing something huge and glaring? If this is a RTFM moment pls let me know, I don't think it is....
cheers, morb
-
What is your current code?
-
The code for constructing the matrices is as above, the code for multiplication of the matrix in Fortran is :
Code:
! TRANSFORMING EUCLIDEAN RANDOM POSN VECTOR TO CRYSTAL BASIS (XPOINT -> XPOINTNEW)
xpointnew = xpoint*lijinv%comp(1,1) + ypoint*lijinv%comp(1,2) &
+ zpoint*lijinv%comp(1,3)
ypointnew = xpoint*lijinv%comp(2,1) + ypoint*lijinv%comp(2,2) &
+ zpoint*lijinv%comp(2,3)
zpointnew = xpoint*lijinv%comp(3,1) + ypoint*lijinv%comp(3,2) &
+ zpoint*lijinv%comp(3,3)
where xpoint, ypoint and zpoint are doubles and lijinv is the inverse of the matrix made in the above fortran code, using the usual determinant & matrix of cofactors approach (this is correct, I checked by multiplying the orig and the inv together, giving the unit matrix). The matrix is a true space transform martix, preserving the length of vectors on multiplication.
As you can see, through cunning use of the indexes a correct form of matrix*vector multiplication is acheived, though the matrix is the wrong way round (transform matrices should be upper diag) due to the way fortran indexes it matrices.
here is the c++ code for multiplying the vector and the matrix :
Code:
probeCrysBasis = euc_to_crys*probeEucBasis ;
I have, though, attempted to produce 'identical' code to the fortran with the same results, the vector is scaled and I have no idea why. I also checked to see whether the eigen library does something like fortran and it doesn't, it uses normal convention on indexing a matrix.
If you need to look at the whole source I could post it though it's not short....
cheers for the trouble laserlight
-
OK, I am off to hit the hay now, so I am boosting this to the top of the forum (whaaahahahaha, /evil_laugh) but here is the relevant bits of code in both fortran and c++ in reverse order specified above:
Code:
if(k != i) for(int t = 0 ; t <= trials-1 ; ++t)
{
phi = 2*pi*rnd.Uniform() ; //produce a random unit vector
theta = acosl(1 - 2*rnd.Uniform()) ;
probeEucBasis(0) = sinl(theta)*cosl(phi) ;
probeEucBasis(1) = sinl(theta)*sinl(phi) ;
probeEucBasis(2) = cosl(theta) ;
probeEucBasis = probeEucBasis*probeIdiameter ; //scale unit vector up by corresponding probe length
cout << "probeEuc Length: " << sqrtl(probeEucBasis.dot(probeEucBasis)) << endl ;
probeCrysBasis = euc_to_crys*probeEucBasis ;
cout << "unscaled probeCrys length: " << sqrtl(probeCrysBasis.dot(probeCrysBasis)) << " should equal above." << endl ;
probeCrysBasis(0) = probeCrysBasis(0)*basisLengths(0) ; //scale up by crystal lattice
probeCrysBasis(1) = probeCrysBasis(1)*basisLengths(1) ; //as basis change scales vector
probeCrysBasis(2) = probeCrysBasis(2)*basisLengths(2) ; //by 1/lattice consts (determinant)
cout << "scaled probeCrys length: " << sqrtl(probeCrysBasis.dot(probeCrysBasis)) << endl ;
and Code:
phi = ran0(seed)*twopi
costheta = 1 - ran0(seed) * 2.0
theta = Acos(costheta)
xpoint = sin(theta)*cos(phi)
ypoint = sin(theta)*sin(phi)
zpoint = costheta
xpoint=xpoint*atomsigma(i)/2.0
ypoint=ypoint*atomsigma(i)/2.0
zpoint=zpoint*atomsigma(i)/2.0
Write(*,*) 'probeEuc Length: ', sqrt(xpoint*xpoint + &
ypoint*ypoint + zpoint*zpoint)
xpointnew = xpoint*lijinv%comp(1,1) + ypoint*lijinv%comp(1,2) &
+ zpoint*lijinv%comp(1,3)
ypointnew = xpoint*lijinv%comp(2,1) + ypoint*lijinv%comp(2,2) &
+ zpoint*lijinv%comp(2,3)
zpointnew = xpoint*lijinv%comp(3,1) + ypoint*lijinv%comp(3,2) &
+ zpoint*lijinv%comp(3,3)
! checking lengths are preserved on lattice transform
Write(*,*) 'unscaled probeCrys Length : ', sqrt(xpoint*xpoint + &
ypoint*ypoint + zpoint*zpoint)
xpointnew = a*xpointnew
ypointnew = b*ypointnew
zpointnew = c*zpointnew
right, so, the problem is not that I am getting more and more accurate estimations of the value of cos(pi/2), it is that the vector is scaling after tranform in the c++ version of the code. I have spent far too long concentrating on this, which should be a simple program to write but has proven that I am c++'s .......... for the past 4 days.
Cheers to all who come to scorn me at my ineptitude.... and those who are helpful too I guess.
-
Sorry for wasting your time people, the problem was that I had not put the suffix 'new' onto my bug checking in the fortran code, as the vector name changed throuh the calculations (no, I don't know why, the prior vector names are not reused except in a new loop iteration, and it doesn't make it any easier to read, anyway, it's not my code)
Still getting weird output but I reckon I got this. Cheers anyway peeps.
