c++ / Fortran90 issue

This is a discussion on c++ / Fortran90 issue within the C++ Programming forums, part of the General Programming Boards category; Ok, I am making a c++ object which does the job of this old Fortran code I have, wanting to ...

  1. #1
    Registered User
    Join Date
    Nov 2010
    Posts
    15

    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

  2. #2
    Registered User
    Join Date
    Nov 2010
    Posts
    15
    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
    Last edited by morbidslug; 11-15-2010 at 10:11 AM. Reason: not fixed after all...

  3. #3
    Registered User
    Join Date
    Nov 2010
    Posts
    15
    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

  4. #4
    C++ Witch laserlight's Avatar
    Join Date
    Oct 2003
    Location
    Singapore
    Posts
    21,711
    What is your current code?
    C + C++ Compiler: MinGW port of GCC
    Version Control System: Bazaar

    Look up a C++ Reference and learn How To Ask Questions The Smart Way

  5. #5
    Registered User
    Join Date
    Nov 2010
    Posts
    15
    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

  6. #6
    Registered User
    Join Date
    Nov 2010
    Posts
    15
    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.
    Last edited by morbidslug; 11-15-2010 at 01:41 PM. Reason: dodgy first edit

  7. #7
    Registered User
    Join Date
    Nov 2010
    Posts
    15

    Lightbulb

    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.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. K&R Learning Issue
    By TimHarrington in forum C Programming
    Replies: 48
    Last Post: 09-06-2010, 04:33 AM
  2. bandwidth issue / network issue with wireless device communication
    By vlrk in forum Networking/Device Communication
    Replies: 0
    Last Post: 07-05-2010, 11:52 PM
  3. float calculation issue
    By George2 in forum C# Programming
    Replies: 1
    Last Post: 05-26-2008, 04:56 AM
  4. type safe issue
    By George2 in forum C++ Programming
    Replies: 4
    Last Post: 02-12-2008, 08:32 PM
  5. my first issue of GDM
    By DavidP in forum A Brief History of Cprogramming.com
    Replies: 0
    Last Post: 09-12-2002, 04:02 PM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21