public static void daxpy_f77 (int n, double da, double dx[], int incx, double
                      dy[], int incy) {

      int i,ix,iy,m;

      if (n <= 0) return;
      if (da == 0.0) return;

      if ((incx == 1) && (incy == 1)) {

//   both increments equal to 1

         m = n%4;

         for (i = 1; i <= m; i++) {

            dy[i] += da*dx[i];

         }

         for (i = m+1; i <= n; i += 4) {

            dy[i]   += da*dx[i];
            dy[i+1] += da*dx[i+1];
            dy[i+2] += da*dx[i+2];
            dy[i+3] += da*dx[i+3];

         }

         return;

      } else {

//   at least one increment not equal to 1

         ix = 1;
         iy = 1;

         if (incx < 0) ix = (-n+1)*incx + 1;
         if (incy < 0) iy = (-n+1)*incy + 1;

         for (i = 1; i <= n; i++) {

            dy[iy] += da*dx[ix];

            ix += incx;
            iy += incy;

         }

         return;

      }

   }