mnddot.cxx

Go to the documentation of this file.

// @(#)root/minuit2:$Id$
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei   2003-2005

/**********************************************************************
 *                                                                    *
 * Copyright (c) 2005 LCG ROOT Math team,  CERN/PH-SFT                *
 *                                                                    *
 **********************************************************************/

/* ddot.f -- translated by f2c (version 20010320).
   You must link the resulting object file with the libraries:
   -lf2c -lm   (in that order)
*/

namespace ROOT {

   namespace Minuit2 {


double mnddot(unsigned int n, const double* dx, int incx, const double* dy,
              int incy) {
   /* System generated locals */
   int i__1;
   double ret_val;

   /* Local variables */
   int i__, m;
   double dtemp;
   int ix, iy, mp1;


   /*     forms the dot product of two vectors. */
   /*     uses unrolled loops for increments equal to one. */
   /*     jack dongarra, linpack, 3/11/78. */
   /*     modified 12/3/93, array(1) declarations changed to array(*) */


   /* Parameter adjustments */
   --dy;
   --dx;

   /* Function Body */
   ret_val = 0.;
   dtemp = 0.;
   if (n <= 0) {
      return ret_val;
   }
   if (incx == 1 && incy == 1) {
      goto L20;
   }

   /*        code for unequal increments or equal increments */
   /*          not equal to 1 */

   ix = 1;
   iy = 1;
   if (incx < 0) {
      ix = (-static_cast<int>(n) + 1) * incx + 1;
   }
   if (incy < 0) {
      iy = (-static_cast<int>(n) + 1) * incy + 1;
   }
   i__1 = n;
   for (i__ = 1; i__ <= i__1; ++i__) {
      dtemp += dx[ix] * dy[iy];
      ix += incx;
      iy += incy;
      /* L10: */
   }
   ret_val = dtemp;
   return ret_val;

   /*        code for both increments equal to 1 */


   /*        clean-up loop */

L20:
      m = n % 5;
   if (m == 0) {
      goto L40;
   }
   i__1 = m;
   for (i__ = 1; i__ <= i__1; ++i__) {
      dtemp += dx[i__] * dy[i__];
      /* L30: */
   }
   if (n < 5) {
      goto L60;
   }
L40:
      mp1 = m + 1;
   i__1 = n;
   for (i__ = mp1; i__ <= i__1; i__ += 5) {
      dtemp = dtemp + dx[i__] * dy[i__] + dx[i__ + 1] * dy[i__ + 1] + dx[
         i__ + 2] * dy[i__ + 2] + dx[i__ + 3] * dy[i__ + 3] + dx[i__ +
            4] * dy[i__ + 4];
      /* L50: */
   }
L60:
      ret_val = dtemp;
   return ret_val;
} /* ddot_ */


   }  // namespace Minuit2

}  // namespace ROOT