sdot.c

Go to the documentation of this file.

/* forms the dot product of two vectors
   uses unrolled loops for increments equal to one
   Ver.1.0, May,23,1988                            */
 
#include "arith.h"
 
REAL sdot(n,sx,isx,jsx,incx,sy,isy,jsy,incy)
int n,incx,incy,isx,jsx,isy,jsy;
MATRIX sx,sy;
{
    REAL stemp;
    int i,ix,iy,m;
 
    stemp=0.0;
 
    if(n <= 0)
        return;
 
    if((incx != 1) || (incy != 1)){
 
        /* code for unequal increments or equal increments
           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++){
            stemp= stemp+sx[ix][1]*sy[iy][1];
            ix= ix+incx;
            iy= iy+incy;
        }
 
        return(stemp);
*/
    }
 
    /* code for both increments equal to 1 */
 
    /* clean-up loop */
 
    m=n%5;
    if(m != 0){
        for(i=0; i<m; i++)
            stemp= stemp+sx[isx+i][jsx]*sy[isy+i][jsy];
    }
    if((m == 0) || (n >= 5)){
        for(i=m; i<n; i += 5)
            stemp= stemp+sx[isx+i][jsx]*sy[isy+i][jsy]
                        +sx[isx+i+1][jsx]*sy[isy+i+1][jsy]
                        +sx[isx+i+2][jsx]*sy[isy+i+2][jsy]
                        +sx[isx+i+3][jsx]*sy[isy+i+3][jsy]
                        +sx[isx+i+4][jsx]*sy[isy+i+4][jsy];
    }
 
    return(stemp);
}