Thread: 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

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

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

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.