/*
*   Java Numerical Toolkit (JNT) for Linear Algebra
*
*   BETA VERASION INCOMPLETE AND SUBJECT TO CHANGE
*
*   Mathematical and Computational Sciences Division
*   National Institute of Standards and Technology
*
*/

/**
* Optimized kernels of the Basic Linear Algebra Subprograms (BLAS)
* in Java.  Currently has only various unrolled versions of
* daxpy() and ddot() kernels.  Final implementation will 
* include Level 1, 2, and 3 kernels.
*
* @version 0.20  11 Sep 1997
* @author Roldan Pozo
*/


import java.math.*;


public class JBLAS
{
                           
    public static final void daxpy(int N, double a, double x[], int x_off, double y[], int y_off)
    {

        for (int i=0; i<N; i++)
            y[i + y_off] += a * x[i + x_off];
    }

    public static final double ddot(int N, double x[], int x_off, double y[], int y_off)
    {
        double sum = 0.0;
        for (int i=0; i<N; i++)
            sum += x[i + x_off] * y[i + y_off];
        return sum;
    }
    
    public static final void daxpy_unroll4(
            int N, double a, double x[], int x_off, double y[], int y_off)
    {
        if (N < 4)    // don't bother for small loops
        {
            for (int i=0; i<N; i++)
            y[i + y_off] += a * x[i + x_off];
        }
        else
        {
            int m = N & 3;      // m = N % 4
            if (m > 0)
            {
                for (int i=0; i<m; i++)
                    y[i + y_off] += a * x[i + x_off];
            }
            int Nm1 = N-1;
            for (int i=m; i<Nm1; i+=4)
            {
                y[i + y_off] += a * x[i + x_off];
                y[i+1 + y_off] += a * x[i+1 + x_off];
                y[i+2 + y_off] += a * x[i+2 + x_off];
                y[i+3 + y_off] += a * x[i+3 + x_off];
            }
        }
    }
    public static final double ddot_unroll4(
            int N, double x[], int x_off, double y[], int y_off)
    {
        double sum = 0.0;
        if (N < 4)    // don't bother for small loops
        {
            for (int i=0; i<N; i++)
                sum += x[i + x_off] * y[i + y_off];
        }
        else
        {
            double sum2 = 0.0;
            int m = N & 3;      // m = N % 4
            if (m > 0)
            {
                for (int i=0; i<m; i++)
                    sum += x[i] * y[i + y_off];
            }
            int Nm1 = N-1;
            for (int i=m; i<Nm1; i+=4)
            {
                sum += x[i + x_off] * y[i + y_off];
                sum2 += x[i+1 + x_off] * y[i+1 + y_off];
                sum += x[i+2 + x_off] * y[i+2 + y_off];
                sum2 += x[i+3 + x_off] * y[i+3 + y_off];
            }
            sum += sum2;
        }
        return sum;
    }

    public static final double ddot_unroll8_sum4(
            int N, double x[], int x_off, double y[], int y_off)
    {
        double sum = 0.0;
        
        if (N < 8)    // don't bother for small loops
        {
            for (int i=0; i<N; i++)
                sum += x[i + x_off]*y[i + y_off]; 
        }
        else
        {
            double sum2 = 0;
            double sum3 = 0;
            double sum4 = 0;
            
            int m = N & 7;      // m = N % 4
            if (m > 0)
            {
                for (int i=0; i<m; i++)
                    sum += x[i + x_off]*y[i + y_off]; 
            }
            int Nm1 = N-1;
            for (int i=m; i<Nm1; )
            {
                sum += x[i + x_off]*y[i + y_off];
                i++;
                sum2 += x[i + x_off]*y[i + y_off];
                i++;
                sum3 += x[i + x_off]*y[i + y_off];
                i++;
                sum4 += x[i + x_off]*y[i + y_off];
                i++;
                sum += x[i + x_off]*y[i + y_off];
                i++;
                sum2 += x[i + x_off]*y[i + y_off];
                i++;
                sum3 += x[i + x_off]*y[i + y_off];
                i++;
                sum4 += x[i + x_off]*y[i + y_off];
                i++;
            }   
        
            sum += sum2 + sum3 + sum4;      
        }
        
        return sum;
    }


    public static final void daxpy_unroll8(
            int N, double a, double x[], int x_off, double y[], int y_off)
    {
        if (N < 8)    // don't bother for small loops
        {
            for (int i=0; i<N; i++)
            y[i + y_off] += a * x[i + x_off];
        }
        else
        {
            int m = N & 7;      // m = N % 4
            if (m > 0)
            {
                for (int i=0; i<m; i++)
                    y[i + y_off] += a * x[i + x_off];
            }
            int Nm1 = N-1;
            for (int i=m; i<Nm1; i+=8)
            {
                y[i + y_off] += a * x[i + x_off];
                y[i+1 + y_off] += a * x[i+1 + x_off];
                y[i+2 + y_off] += a * x[i+2 + x_off];
                y[i+3 + y_off] += a * x[i+3 + x_off];
                y[i+4 + y_off] += a * x[i+4 + x_off];
                y[i+5 + y_off] += a * x[i+5 + x_off];
                y[i+6 + y_off] += a * x[i+6 + x_off];
                y[i+7 + y_off] += a * x[i+7 + x_off];
            }
        }
    }

    public static final void daxpy_unroll8_by1(
            int N, double a, double x[], int x_off, double y[], int y_off)
    {
        if (N < 8)    // don't bother for small loops
        {
            for (int i=0; i<N; i++)
            y[i + y_off] += a * x[i + x_off];
        }
        else
        {
            int m = N & 7;      // m = N % 4
            if (m > 0)
            {
                for (int i=0; i<m; i++)
                    y[i + y_off] += a * x[i + x_off];
            }
            int Nm1 = N-1;
            for (int i=m; i<Nm1; )
            {
                y[i + y_off] += a * x[i + x_off];
                i++;
                y[i + y_off] += a * x[i + x_off];
                i++;
                y[i + y_off] += a * x[i + x_off];
                i++;
                y[i + y_off] += a * x[i + x_off];
                i++;
                y[i + y_off] += a * x[i + x_off];
                i++;
                y[i + y_off] += a * x[i + x_off];
                i++;
                y[i + y_off] += a * x[i + x_off];
                i++;
                y[i + y_off] += a * x[i + x_off];
                i++;
            }
        }
    }

            
            
            
            
/*
    public static final void daxpy(int N, double a, double x[], int xi,
        double y[], int yi)
    {

        for (int i=0; i<N; i++)
            y[yi+i] += a * x[xi+i];
    }
*/
        
    public static final void daxpy(int N, double a, double x[], 
        int x_offset, int incx, double y[], int y_offset, int incy)
    {
        if (incx ==1 && incy==1)
        {
            if (x_offset == 0 && y_offset==0)
            {
                for (int i=0; i<N; i++)
                    y[i] += a * x[i];
            }
            else
            {
                for (int i=0; i<N; i++)
                    y[y_offset+i] += a * x[x_offset+i];
            }
        }
        else  /* the most general case */
        {
            for (int i=0; i<N; i++)
                y[y_offset + i*incy] += a * x[x_offset + i*incx];
        }
    }



}
/* end class BLAS */

.