# include <stdlib.h>
# include <stdio.h>
# include <math.h>
# include <time.h>
# include <complex.h>
# include "blas0.h"
# include "blas1_z.h"
/******************************************************************************/
double dzasum ( int n, double complex x[], int incx )
/******************************************************************************/
/*
Purpose:
DZASUM takes the sum of the absolute values of a vector.
Discussion:
This routine uses double precision complex arithmetic.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
31 March 2007
Author:
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,
ACM Transactions on Mathematical Software,
Volume 5, Number 3, pages 308-323, 1979.
Parameters:
Input, int N, the number of entries in the vector.
Input, double complex X[], the vector.
Input, int INCX, the increment between successive entries of X.
Output, double DZASUM, the sum of the absolute values.
*/
{
int i;
int ix;
double value;
value = 0.0;
if ( n <= 0 || incx <= 0 )
{
return value;
}
if ( incx == 1 )
{
for ( i = 0; i < n; i++ )
{
value = value + fabs ( creal ( x[i] ) )
+ fabs ( cimag ( x[i] ) );
}
}
else
{
ix = 0;
for ( i = 0; i < n; i++ )
{
value = value + fabs ( creal ( x[ix] ) )
+ fabs ( cimag ( x[ix] ) );
ix = ix + incx;
}
}
return value;
}
/******************************************************************************/
double dznrm2 ( int n, double complex x[], int incx )
/******************************************************************************/
/*
Purpose:
DZNRM2 returns the euclidean norm of a vector.
Discussion:
This routine uses double precision complex arithmetic.
DZNRM2 := sqrt ( sum ( conjg ( x(1:n) ) * x(1:n) ) )
= sqrt ( dot_product ( x(1:n), x(1:n) ) )
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
31 March 2007
Author:
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,
ACM Transactions on Mathematical Software,
Volume 5, Number 3, pages 308-323, 1979.
Parameters:
Input, int N, the number of entries in the vector.
Input, double complex X[], the vector.
Input, int INCX, the increment between successive entries of X.
Output, double DZNRM2, the norm of the vector.
*/
{
int i;
int ix;
double scale;
double ssq;
double temp;
double value;
if ( n < 1 || incx < 1 )
{
value = 0.0;
}
else
{
scale = 0.0;
ssq = 1.0;
ix = 0;
for ( i = 0; i < n; i++ )
{
if ( creal ( x[ix] ) != 0.0 )
{
temp = fabs ( creal ( x[ix] ) );
if ( scale < temp )
{
ssq = 1.0 + ssq * pow ( scale / temp, 2 );
scale = temp;
}
else
{
ssq = ssq + pow ( temp / scale, 2 );
}
}
if ( cimag ( x[ix] ) != 0.0 )
{
temp = fabs ( cimag ( x[ix] ) );
if ( scale < temp )
{
ssq = 1.0 + ssq * pow ( scale / temp, 2 );
scale = temp;
}
else
{
ssq = ssq + pow ( temp / scale, 2 );
}
}
ix = ix + incx;
}
value = scale * sqrt ( ssq );
}
return value;
}
/******************************************************************************/
int izamax ( int n, double complex x[], int incx )
/******************************************************************************/
/*
Purpose:
IZAMAX indexes the vector element of maximum absolute value.
Discussion:
This routine uses double precision complex arithmetic.
WARNING: This index is a 1-based index, not a 0-based index!
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
31 March 2007
Author:
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,
ACM Transactions on Mathematical Software,
Volume 5, Number 3, pages 308-323, 1979.
Parameters:
Input, int N, the number of entries in the vector.
Input, double complex X[], the vector.
Input, int INCX, the increment between successive entries of X.
Output, int IZAMAX, the index of the element of maximum
absolute value.
*/
{
int i;
int ix;
double smax;
int value;
value = 0;
if ( n < 1 || incx <= 0 )
{
return value;
}
value = 1;
if ( n == 1 )
{
return value;
}
if ( incx != 1 )
{
ix = 0;
smax = zabs1 ( x[0] );
ix = ix + incx;
for ( i = 1; i < n; i++ )
{
if ( smax < zabs1 ( x[ix] ) )
{
value = i + 1;
smax = zabs1 ( x[ix] );
}
ix = ix + incx;
}
}
else
{
smax = zabs1 ( x[0] );
for ( i = 1; i < n; i++ )
{
if ( smax < zabs1 ( x[i] ) )
{
value = i + 1;
smax = zabs1 ( x[i] );
}
}
}
return value;
}
/******************************************************************************/
void zaxpy ( int n, double complex ca, double complex cx[],
int incx, double complex cy[], int incy )
/******************************************************************************/
/*
Purpose:
ZAXPY computes a constant times a vector plus a vector.
Discussion:
This routine uses double precision complex arithmetic.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
31 March 2007
Author:
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 CX and CY.
Input, double complex CA, the multiplier of CX.
Input, double complex CX[], the first vector.
Input, int INCX, the increment between successive entries of CX.
Input/output, double complex CY[], the second vector.
On output, CY(*) has been replaced by CY(*) + CA * CX(*).
Input, int INCY, the increment between successive entries of CY.
*/
{
int i;
int ix;
int iy;
if ( n <= 0 )
{
return;
}
if ( zabs1 ( ca ) == 0.0 )
{
return;
}
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++ )
{
cy[iy] = cy[iy] + ca * cx[ix];
ix = ix + incx;
iy = iy + incy;
}
}
else
{
for ( i = 0; i < n; i++ )
{
cy[i] = cy[i] + ca * cx[i];
}
}
return;
}
/******************************************************************************/
void zcopy ( int n, double complex cx[], int incx, double complex cy[],
int incy )
/******************************************************************************/
/*
Purpose:
ZCOPY copies a vector X to a vector Y.
Discussion:
This routine uses double precision complex arithmetic.
Licensing:
This code is distribute