Difference between Fortran code and is C++ version

[drudox84]

Hi all , I don't understand why this two doce give different result ... I know that Fortrn start index from 1 and C++ from 0 .. but i think that during the converion from fortran I've take into account this fact .. so the codes are :

C++

# include <iosfwd>
# include <vector>
# include <cmath>
# include <iomanip>
# include <iostream>

std::vector<double> dot( int size, std::vector<double> x, 
                         std::vector<double> aa, std::vector<int> ja)
{     

      std::vector<double> y(x.size());

      for(auto i = 0; i < size ; i++ )
            y.at(i) = aa.at(i) * x.at(i);

      for(auto i=0 ; i < size ; i++ )
      {
         //for(auto k=ja.at(i) ; k< ja.at(i+1)-1 ; k++ )
         auto k = ja.at(i);

         do
         {
            y.at(i) += aa.at(k) * x.at(ja.at(k)) ;
            k++;
         }
         while(k < ja.at(i+1)-1) ;

      }
      return y;
}


int main()
{
    std::vector<double> x = {0.,1.3,4.2,0.8} ;

    std::vector<double> aa = {1.01,4.07,6.08,9.9,0.,2.34,3.12,1.06,2.2};
    std::vector<int>    ja = {6,7,7,8,10,3,1,1,3};

    std::vector<double> y = dot(x.size(), x , aa , ja);

    for(auto& i : y)
         std::cout << i << ' ' ;   
    std::cout << std::endl;  
}

and Fortran that give me the right result :

MODULE MSR
 IMPLICIT NONE

CONTAINS
     subroutine amuxms (n, x, y, a,ja)
      real*8  x(*), y(*), a(*)
      integer n, ja(*)
      integer i, k
      do 10 i=1, n
        y(i) = a(i)*x(i)
 10     continue
      do 100 i = 1,n

         do 99 k=ja(i), ja(i+1)-1
            y(i) = y(i) + a(k) *x(ja(k))
 99      continue
 100  continue

      return

      end

END MODULE

PROGRAM MSRtest
USE MSR
IMPLICIT NONE
INTEGER :: i
REAL(KIND(0.D0)), DIMENSION(4) :: y, x= (/0.,1.3,4.2,0.8/)

REAL(KIND(0.D0)), DIMENSION(9) :: AA = (/ 1.01, 4.07, 6.08, 9.9, 0., 2.34, 3.12, 1.06, 2.2/) 
INTEGER , DIMENSION(9)         :: JA = (/6, 7, 7, 8, 10, 3, 1, 1, 3/)
CALL amuxms(4,x,y,aa,ja)

WRITE(6,FMT='(4F8.3)') (y(I), I=1,4)    

END PROGRAM

Fortran code return the right result : 9.828 5.291 25.536 17.160
c++ code return 4.056 6.669 26.914 9.68 that is wrong !!

could you help me ?

[jlb]

c++ code return 4.056 6.669 26.914 9.68 that is wrong !!

Really? What I get when I run the code is:

0 1.3 4.2 0.8

If I print the values of the y vector I get the values you reported however.

Did you follow your logic with your debugger to see where you are not following the Fortran logic?

[Peter87]

The values in JA are used as indices so when you move to C++ you need to subtract one from each of them. If you do that you will get the correct answer except for the last value. If you also remove the -1 from line 26 it seems to give the output that you want but I can't say I understand why so I'm not sure if it's a proper solution or not.

[lastchance]

As Peter87 says, they are used as indices. A quick fix is to add a -1 at the end of line 19, and before the last ) on line 23.

Alternatively, subtract 1 from every ja vector value and remove the -1 in line 26.

Neither C++ nor Fortran code is very readable - you aren't using array operations in Fortran, and a valarray might be a better vehicle here in c++.

[jonnin]

fortran is column based for 2-d and higher memory mappings. This is huge. C++ is row major.


also remember that the variable's name (starting letter) determines its type in F.

I don't think either of those will affect you here, but if there is more code we don't see, keep it in mind.

[lastchance]

@jonnin, the variable's name does not determine the type in fortran if implicit none is specified; there are also no 2d arrays here to confuse with row-major/column-major. If you want the c++ and fortran arrays to match indices then fortran does allow you to start the indexing from 0 (or anything else) if you want to: just change to
dimension(0:3)
for example. Then you could get a more direct comparison between the two languages.

[jonnin]

Did that change that? The fortran I know variables that start with I, j, or k are integer, everything else is 'real'. You can over-ride this by typing items, but that can trip up the unwary.

I know there was no 2d here, I said to watch for it if there was more code than what was shown.

Those and the already mentioned index from 0 vs index from 1 are the only 3 huge items I know of (granted the syntax is very different, but syntax can be fixed, these the 3 things that if you didn't know would drive you nuts trying to find the problem). There are also the fortan column oddities (what column the statements are in affect behaviors), but that would only affect really, really old code compiled on old compilers.

[TheIdeasMan]

there is also std::inner_product :+)

http://en.cppreference.com/w/cpp/algorithm/inner_product

[lastchance]

C++

#include <iostream>
#include <iomanip>
#include <vector>
using namespace std;


//======================================


void amuxms( int n, const vector<double> &x, vector<double> &y, const vector<double> &a, const vector<int> &ja )
{     
   for ( int i = 0; i < n; i++ )
   {
      y[i] = a[i] * x[i];
      for ( int k = ja[i]; k < ja[i+1]; k++ )
      {
         y[i] += a[k] * x[ ja[k] ];
      }
   }
}


//======================================


int main()
{
   vector<double> x  = { 0.0, 1.3, 4.2, 0.8 };
   int n = x.size();
   vector<double> y(n);
   vector<double> aa = { 1.01, 4.07, 6.08, 9.9, 0.0, 2.34, 3.12, 1.06, 2.2 };
   vector<int>    ja = { 6, 7, 7, 8, 10, 3, 1, 1, 3 };
   for ( int &j : ja ) j--;            // Convert ja to start indices at 0

   amuxms( n, x, y, aa, ja );

   for( double e : y ) cout << setw( 8 ) << setprecision( 3 ) << fixed << e;
   cout << endl;  
}


//====================================== 

   9.828   5.291  25.536  17.160

Fortran

subroutine amuxms( n, x, y, a, ja )
   implicit none
   integer, parameter :: rkind = kind( 0.0d0 )

   integer          , intent(in ) :: n
   real (kind=rkind), intent(in ) :: x(n), a(*)
   integer          , intent(in ) :: ja(*)
   real (kind=rkind), intent(out) :: y(n)
   integer i, k
   
   do i = 1, n
      y(i) = a(i) * x(i)
      do k = ja(i), ja(i+1) - 1
         y(i) = y(i) + a(k) * x( ja(k) )
      end do
   end do
   
end subroutine amuxms


!=======================================


program msrtest
   implicit none
   integer, parameter :: rkind = kind( 0.0d0 )

   integer, parameter :: n = 4
   real (kind=rkind) :: x(n)  = (/ 0.0, 1.3, 4.2, 0.8 /)
   real (kind=rkind) :: y(n)
   real (kind=rkind) :: aa(9) = (/ 1.01, 4.07, 6.08, 9.9, 0.0, 2.34, 3.12, 1.06, 2.2 /)
   integer           :: ja(9) = (/ 6, 7, 7, 8, 10, 3, 1, 1, 3 /)

   call amuxms( n, x, y, aa, ja )

   write( *, fmt = "( 4f8.3 )" ) y

end program msrtest


!=======================================

   9.828   5.291  25.536  17.160