# include <stdlib.h>
# include <stdio.h>
# include <math.h>
# include <time.h>
# include <complex.h>
# include "blas0.h"
# include "blas1_d.h"
/******************************************************************************/
double dasum ( int n, double x[], int incx )
/******************************************************************************/
/*
Purpose:
DASUM takes the sum of the absolute values of a vector.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
30 March 2007
Author:
Original FORTRAN77 version by Charles Lawson, Richard Hanson,
David Kincaid, Fred Krogh.
C version by John Burkardt.
Reference:
Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
LINPACK User's Guide,
SIAM, 1979.
Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
Basic Linear Algebra Subprograms for Fortran Usage,
Algorithm 539,
ACM Transactions on Mathematical Software,
Volume 5, Number 3, September 1979, pages 308-323.
Parameters:
Input, int N, the number of entries in the vector.
Input, double X[*], the vector to be examined.
Input, int INCX, the increment between successive entries of X.
INCX must not be negative.
Output, double DASUM, the sum of the absolute values of X.
*/
{
int i;
int j;
double value;
value = 0.0;
j = 0;
for ( i = 0; i < n; i++ )
{
value = value + fabs ( x[j] );
j = j + incx;
}
return value;
}
/******************************************************************************/
void daxpy ( int n, double da, double dx[], int incx, double dy[], int incy )
/******************************************************************************/
/*
Purpose:
DAXPY computes constant times a vector plus a vector.
Discussion:
This routine uses unrolled loops for increments equal to one.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
30 March 2007
Author:
Original FORTRAN77 version by Charles Lawson, Richard Hanson,
David Kincaid, Fred Krogh.
C version by John Burkardt.
Reference:
Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
LINPACK User's Guide,
SIAM, 1979.
Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
Basic Linear Algebra Subprograms for Fortran Usage,
Algorithm 539,
ACM Transactions on Mathematical Software,
Volume 5, Number 3, September 1979, pages 308-323.
Parameters:
Input, int N, the number of elements in DX and DY.
Input, double DA, the multiplier of DX.
Input, double DX[*], the first vector.
Input, int INCX, the increment between successive entries of DX.
Input/output, double DY[*], the second vector.
On output, DY[*] has been replaced by DY[*] + DA * DX[*].
Input, int INCY, the increment between successive entries of DY.
*/
{
int i;
int ix;
int iy;
int m;
if ( n <= 0 )
{
return;
}
if ( da == 0.0 )
{
return;
}
/*
Code for unequal increments or equal increments
not equal to 1.
*/
if ( incx != 1 || incy != 1 )
{
if ( 0 <= incx )
{
ix = 0;
}
else
{
ix = ( - n + 1 ) * incx;
}
if ( 0 <= incy )
{
iy = 0;
}
else
{
iy = ( - n + 1 ) * incy;
}
for ( i = 0; i < n; i++ )
{
dy[iy] = dy[iy] + da * dx[ix];
ix = ix + incx;
iy = iy + incy;
}
}
/*
Code for both increments equal to 1.
*/
else
{
m = n % 4;
for ( i = 0; i < m; i++ )
{
dy[i] = dy[i] + da * dx[i];
}
for ( i = m; i < n; i = i + 4 )
{
dy[i ] = dy[i ] + da * dx[i ];
dy[i+1] = dy[i+1] + da * dx[i+1];
dy[i+2] = dy[i+2] + da * dx[i+2];
dy[i+3] = dy[i+3] + da * dx[i+3];
}
}
return;
}
/******************************************************************************/
void dcopy ( int n, double dx[], int incx, double dy[], int incy )
/******************************************************************************/
/*
Purpose:
DCOPY copies a vector X to a vector Y.
Discussion:
The routine uses unrolled loops for increments equal to one.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
30 March 2007
Author:
Original FORTRAN77 version by Charles Lawson, Richard Hanson,
David Kincaid, Fred Krogh.
C version by John Burkardt.
Reference:
Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
LINPACK User's Guide,
SIAM, 1979.
Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
Basic Linear Algebra Subprograms for Fortran Usage,
Algorithm 539,
ACM Transactions on Mathematical Software,
Volume 5, Number 3, September 1979, pages 308-323.
Parameters:
Input, int N, the number of elements in DX and DY.
Input, double DX[*], the first vector.
Input, int INCX, the increment between successive entries of DX.
Output, double DY[*], the second vector.
Input, int INCY, the increment between successive entries of DY.
*/
{
int i;
int ix;
int iy;
int m;
if ( n <= 0 )
{
return;
}
if ( incx == 1 && incy == 1 )
{
m = n % 7;
if ( m != 0 )
{
for ( i = 0; i < m; i++ )
{
dy[i] = dx[i];
}
}
for ( i = m; i < n; i = i + 7 )
{
dy[i] = dx[i];
dy[i + 1] = dx[i + 1];
dy[i + 2] = dx[i + 2];
dy[i + 3] = dx[i + 3];
dy[i + 4] = dx[i + 4];
dy[i + 5] = dx[i + 5];
dy[i + 6] = dx[i + 6];
}
}
else
{
if ( 0 <= incx )
{
ix = 0;
}
else
{
ix = ( -n + 1 ) * incx;
}
if ( 0 <= incy )
{
iy = 0;
}
else
{
iy = ( -n + 1 ) * incy;
}
for ( i = 0; i < n; i++ )
{
dy[iy] = dx[ix];
ix = ix + incx;
iy = iy + incy;
}
}
return;
}
/******************************************************************************/
double ddot ( int n, double dx[], int incx, double dy[], int incy )
/******************************************************************************/
/*
Purpose:
DDOT forms the dot product of two vectors.
Discussion:
This routine uses unrolled loops for increments equal to one.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
30 March 2007
Author:
Original FORTRAN77 version by Charles Lawson, Richard Hanson,
David Kincaid, Fred Krogh.
C version by John Burkardt.
Reference:
Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
LINPACK User's Guide,
SIAM, 1979.
Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
Basic Linear Algebra Subprograms for Fortran Usage,
Algorithm 539,
ACM Transactions on Mathematical Software,
Volume 5, Number 3, September 1979, pages 308-323.
Parameters:
Input, int N, the number of entries in the vectors.
Input, double DX[*], the first vector.
Input, int INCX, the increment between successive entries in DX.
Input, double DY[*], the second vector.
Input, int INCY, the increment between successive entries in DY.
Output, double DDOT, the sum of the product of the corresponding
entries of DX and DY.
*/
{
double dtemp;
int i;
int ix;
int iy;
int m;
dtemp = 0.0;
if ( n <= 0 )
{
return dtemp;
}
/*
Code for unequal increments or equal increments
not equal to 1.
*/
if ( incx != 1 || incy != 1 )
{
if ( 0 <= incx )
{
ix = 0;
}
else
{
ix = ( - n + 1 ) * incx;
}
if ( 0 <= incy )
{
iy = 0;
}
else
{
iy = ( - n + 1 ) * incy;
}
for ( i = 0; i < n; i++ )
{
dtemp = dtemp + dx[ix] * dy[iy];
ix = ix + incx;
iy = iy + incy;
}
}
/*
Code for both increments equal to 1.
*/
else
{
m = n % 5;
for ( i = 0; i < m; i++ )
{
dtemp = dtemp + dx[i]