0

I'm not explain what is different in this two code 1 written in c++ and the second in Fortran... I didn't get the point of the difference..

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) += y.at(i) + aa.at(k) * x.at(ja.at(k)) ;
            k++;
         }
         while(k < ja.at(i+1)-1) ;

      }

}


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 : x)
         std::cout << i << ' ' ;   
    std::cout << std::endl;  
}

Fortran code, I know that c++ start the index at 0 and Fortran at 1 but I think that I've consider this! by the way the right one is the fortran code as follow :

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

I've solved I made a rough mistake ! using ja_ as index .. I forgot to subtrac 1 !

so the working function is : 

# 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++ )
      {

         auto k = ja.at(i)-1;

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

      }
      return y;
}

the difference of two codes is inside the use of ja_(index) as index of other vector without taking into account that in the Fortran code ja_(1) give us the first value of ja_ in the c++ the second ! @Vladimir F it is now clear ?

Drudox lebowsky
  • 1,020
  • 7
  • 21
  • Which difference? You must describe the difference, it is possible that even if we run the code we would see something different than what you see. Please show us your results and comment them. Read [ask] and [mcve]. – Vladimir F Героям слава Nov 18 '17 at 20:00
  • 2
    Your numbers in the initialized arrays are single precision, not double precision (plenty of duplicates). But really, you must show us the difference you are talking about. – Vladimir F Героям слава Nov 18 '17 at 20:02
  • One of the maybe duplicates: https://stackoverflow.com/questions/6146005/is-there-a-better-double-precision-assignment-in-fortran-90 See also https://groups.google.com/forum/#!topic/comp.lang.fortran/IoI5Mes2Se4[151-175] – Vladimir F Героям слава Nov 18 '17 at 20:07
  • This doesn't address the question, but if you want to avoid accessing the arrays out of bounds, check their sizes once at the beginning of the function. Then replace all those calls to `at(i)` with index operations. Checking bounds on every access to every array when you're using simple `for` loops wastes a bunch of time. – Pete Becker Nov 18 '17 at 20:39
  • You say "I know that c++ start the index at 0 and Fortran at 1 but I think that I've consider this!". You haven't done that correctly. `ja` is an array which stores indexes into another array, and you haven't accounted for the "off-by-one" associated with that. – Peter Nov 18 '17 at 22:02
  • I was asking about the differences **in the results** not for the explanation. Leave it now, but please, next time when you ask a question and you have some strange results, **show those results**. It is really, really, really important. Believe me, it is. Please do really read [ask]. – Vladimir F Героям слава Nov 19 '17 at 12:14

1 Answers1

1

There are several issues with the C++ code that I see.

First, the dot function does not return anything. It needs a return y; statement at the end.

Another is with the computation of y.at(i) in the innermost loop of that function. It should be either y.at(i) = y.at(i) + aa... or y.at(i) += aa... but not both.

A third problem is that in your output loop, you're using the wrong array (x instead of y).

There may be other problems. Running the programs in the debugger and comparing their executions can help spot them.

1201ProgramAlarm
  • 32,384
  • 7
  • 42
  • 56