I am trying to run a code that contains this function and I am getting the "*** stack smashing detected ***: terminated".

void Lapenta_Markidis( long double v[3], long double E[3], long double B[3], long double c2,  long double upart[3] ){
   
   long double   upartk[3], vbar[3];
   long double  tmp[3], Fk[3], C1[3], C2[3] ;
   long double  dupartk[3];
   long double gL = gamma_v(v, c2);
   

    for(int i=0;i<3;i  ){       
        upart[i] = v[i]*gL;       // momentum at start of time step
        upartk[i] = upart[i];
    }


    /* Start of the nonlinear cycle */

    long double     abserr = 1.0;
    long double     tol=1e-14;
    int             nkmax=10;
    int             nk = 0; //


    do{
        long double  J11, J12, J13,J21, J22,  J23, J31, J32, J33, Det;
        long double  iJ11, iJ12, iJ13,iJ21, iJ22,  iJ23, iJ31, iJ32, iJ33;
        long double  gL_new;

        nk     =  nk 1;
        gL_new = gamma_u(upartk, c2);

        for(int i=0;i<3;i  ){  
            vbar[i] = (upart[i]   upartk[i])/(gL_new   gL);
        }
        
        crossP(vbar,B,tmp);

        // Compute residual vector
        for(int i=0;i<3;i  ){ 
            Fk[i] = upartk[i] - upart[i] - q*dt/mp * (E[i]   tmp[i]);
        }

         // Compute auxiliary coefficients
        for(int i=0;i<3;i  ){
            C1[i] = (gL_new   gL - (upartk[i]*(upartk[i]   upart[i])) / (gL_new*c2) )/ pow((gL   gL_new),2) ;
            C2[i] = -( upartk[i] / (gL_new*c2) )/ ((gL_new   gL),2) ;
        }

    
        // Compute Jacobian
          J11 = 1. - (q*dt/mp) * (C2[1] * (upartk[2]   upart[2]) * B[3] - C2[1] * (upartk[3]   upart[3]) * B[2]) ;
          J12 = - (q*dt/mp)*(C1[2] * B[3] - C2[2] * (upartk[3]   upart[3]) * B[2]) ;
          J13 = - (q*dt/mp) * (C2[3] * (upartk[2]   upart[2]) * B[3] - C1[3] * B[2]) ;
          J21 = - q*dt/mp * (- C1[1] * B[3]   C2[1] * (upartk[3]   upart[3]) * B[1]) ;
          J22 = 1. - q*dt/mp * (- C2[2] * (upartk[1]   upart[1]) * B[3]   C2[2] * (upartk[3]   upart[3]) * B[1]) ;
          J23 = - q*dt/mp * (- C2[3] * (upartk[1]   upart[1]) * B[3]   C1[3] * B[1]) ;
          J31 = - q*dt/mp * (C1[1] * B[2] - C2[1] * (upartk[2]   upart[2]) * B[1]) ;
          J32 = - q*dt/mp * (C2[2] * (upartk[1]   upart[1]) * B[2] - C1[2] * B[1]) ;
          J33 = 1. - q*dt/mp * (C2[3] * (upartk[1]   upart[1]) * B[2] - C2[3] * (upartk[2]   upart[2]) * B[1]) ;

          // Compute inverse Jacobian
           Det = J11*J22*J33   J21*J32*J13   J31*J12*J23 - J11*J32*J23 - J31*J22*J13 - J21*J12*J33;
           iJ11 = (J22*J33 - J23*J32) / Det ;
           iJ12 = (J13*J32 - J12*J33) / Det ;
           iJ13 = (J12*J23 - J13*J22) / Det ;
           iJ21 = (J23*J31 - J21*J33) / Det ;
           iJ22 = (J11*J33 - J13*J31) / Det ;
           iJ23 = (J13*J21 - J11*J23) / Det ;
           iJ31 = (J21*J32 - J22*J31) / Det ;
           iJ32 = (J12*J31 - J11*J32) / Det ;
           iJ33 = (J11*J22 - J12*J21) / Det ;

        // Compute new upartk = upartk - J^(-1) * F(upartk)
          dupartk[1] = - (iJ11 * Fk[1]   iJ12 * Fk[2]   iJ13 * Fk[3]);
          dupartk[2] = - (iJ21 * Fk[1]   iJ22 * Fk[2]   iJ23 * Fk[3]);
          dupartk[3] = - (iJ31 * Fk[1]   iJ32 * Fk[2]   iJ33 * Fk[3]);

        // Check convergence
        for(int i=0;i<3;i  ){
          upartk[i]  =  dupartk[i] ;
        }
          abserr = sqrt(dotP(dupartk, dupartk));
    } while(abserr > tol && nk < nkmax); // End of while -> end of cycle

    // Update velocity 
    for(int i=0;i<3;i  ){
        upart[i] = upartk[i];
    }
   
    
}

I am trying to run a code that contains this function and I am getting the "*** stack smashing detected ***: terminated".

Any suggestions of what I am doing wrong? I am not that familiar with the C syntax, have I declared a variable, matrix in the wrong way?

CodePudding user response：

You mix 1-based indexing and 0-based indexing. But C arrays use 0-based indexing.

At several positions you use variable[3], where only variable[2] is allowed:

      dupartk[3] = - (iJ31 * Fk[1]   iJ32 * Fk[2]   iJ33 * Fk[3]);
//            ^                                               ^

Move all those indices by one and those accesses there be fine:

      dupartk[0] = - (iJ11 * Fk[0]   iJ12 * Fk[1]   iJ13 * Fk[2]);
      dupartk[1] = - (iJ21 * Fk[0]   iJ22 * Fk[1]   iJ23 * Fk[2]);
      dupartk[2] = - (iJ31 * Fk[0]   iJ32 * Fk[1]   iJ33 * Fk[2]);

But keep in mind that there are several other wrong indices, e.g. B[3] and C[3]. Check each index.

CodePudding user response：

The array indexes are from zero not one.

You write and read outside the bounds all of your arrays in your code.

You simply need to decrease all the indexes by 1.