C lang

  /********************
   *      prop_W      *
   ********************/

  for(int ip = 0; ip < 4; ++ip) { // propagator type

    // create source
    for(int i = 0; i < 2*3*4*XYZTnodeSites; i++) source[i] = 0.0;

    // only on the iT time-slice
    int iT = t_S;
    if (t0 <= iT && iT < t1) {
      int it = iT % TnodeSites;

      Float q_z, *prop;
      // notation: (pm|pi) where pm (pi) is the intermediate (initial) src momentum
      if      (ip == 0) { prop = prop_Q[0]; q_z = - 2.0 * M_PI / Zsites; } // (-p|+p)
      else if (ip == 1) { prop = prop_Q[1]; q_z = + 2.0 * M_PI / Zsites; } // (+p|-p)
      else if (ip == 2) { prop = prop_Q[1]; q_z = 0.0;                   } // ( 0|-p)
      else if (ip == 3) { prop = prop_Q[2]; q_z = - 2.0 * M_PI / Zsites; } // (-p| 0)

      for(    int iz = 0; iz < ZnodeSites; ++iz) {

        int iZ = iz + ZnodeCoor * ZnodeSites;
        COMPLEX Factor = cos(q_z * iZ) + COMPLEX_I * sin(q_z * iZ);

        for(  int iy = 0; iy < YnodeSites; ++iy) {
          for(int ix = 0; ix < XnodeSites; ++ix) {

            int ixyzt = ix + XnodeSites*(iy + YnodeSites*(iz + ZnodeSites*it));

            // indices of the intermediate (m for middle) source pion
            for(  int Dm = 0; Dm < 4; Dm++) {
              for(int Cm = 0; Cm < 3; Cm++) {

              Source(Cm,Dm, ixyzt) = Factor * ZGM(Dm,5) * Prop(Cm,IGM(Dm,5), ixyzt);

            } }
      } } }
    }

    // call the solver
    int iter = invoke_solver_dirichlet_bc(lat,cgarg,source,tmp,is_fix_gauge,N_Dirichlet_BC,ary_Dirichlet_BC,&true_res);