Commit 16fbe0c0 authored by Daniel Brown's avatar Daniel Brown
Browse files

Due to both SuperLU and Finesse defining doublecomplex and complex_t, both of...

Due to both SuperLU and Finesse defining doublecomplex and complex_t, both of which are identical in memory layout terms just a difference in naming, I have changed the SuperLU definition of a double precision complex to use Finesse's struct. Obviously any future update to SuperLU will need to apply a similar refactoring.
parent 72a1bef5
......@@ -181,8 +181,8 @@ cgstrs(trans_t trans, SuperMatrix *L, SuperMatrix *U,
for (i = 0; i < nrow; i++) {
irow = L_SUB(iptr);
c_sub(&rhs_work[irow], &rhs_work[irow], &work_col[i]);
work_col[i].r = 0.0;
work_col[i].i = 0.0;
work_col[i].re = 0.0;
work_col[i].im = 0.0;
iptr++;
}
}
......@@ -197,8 +197,8 @@ cgstrs(trans_t trans, SuperMatrix *L, SuperMatrix *U,
for (i = 0; i < nrow; i++) {
irow = L_SUB(iptr);
c_sub(&rhs_work[irow], &rhs_work[irow], &work[i]);
work[i].r = 0.0;
work[i].i = 0.0;
work[i].re = 0.0;
work[i].im = 0.0;
iptr++;
}
}
......
......@@ -85,8 +85,8 @@ clacon_(int *n, complex *v, complex *x, float *est, int *kase)
safmin = slamch_("Safe minimum");
if ( *kase == 0 ) {
for (i = 0; i < *n; ++i) {
x[i].r = 1. / (float) (*n);
x[i].i = 0.;
x[i].re = 1. / (float) (*n);
x[i].im = 0.;
}
*kase = 1;
jump = 1;
......@@ -116,8 +116,8 @@ clacon_(int *n, complex *v, complex *x, float *est, int *kase)
d__1 = c_abs(&x[i]);
if (d__1 > safmin) {
d__1 = 1 / d__1;
x[i].r *= d__1;
x[i].i *= d__1;
x[i].re *= d__1;
x[i].im *= d__1;
} else {
x[i] = one;
}
......@@ -157,8 +157,8 @@ L90:
d__1 = c_abs(&x[i]);
if (d__1 > safmin) {
d__1 = 1 / d__1;
x[i].r *= d__1;
x[i].i *= d__1;
x[i].re *= d__1;
x[i].im *= d__1;
} else {
x[i] = one;
}
......@@ -173,7 +173,7 @@ L110:
jlast = j;
j = icmax1_(n, &x[0], &c__1);
--j;
if (x[jlast].r != (d__1 = x[j].r, fabs(d__1)) && iter < 5) {
if (x[jlast].re != (d__1 = x[j].re, fabs(d__1)) && iter < 5) {
++iter;
goto L50;
}
......@@ -182,8 +182,8 @@ L110:
L120:
altsgn = 1.;
for (i = 1; i <= *n; ++i) {
x[i-1].r = altsgn * ((float) (i - 1) / (float) (*n - 1) + 1.);
x[i-1].i = 0.;
x[i-1].re = altsgn * ((float) (i - 1) / (float) (*n - 1) + 1.);
x[i-1].im = 0.;
altsgn = -altsgn;
}
*kase = 1;
......
......@@ -41,7 +41,7 @@ cband(int n, int b, int nonz, complex **nzval, int **rowind, int **colptr)
ilow = SUPERLU_MAX(0, j - ub);
ihigh = SUPERLU_MIN(n-1, j + lb);
for (i = ilow; i <= ihigh; ++i) {
val[i-ilow].r = dlaran_(iseed);
val[i-ilow].re = dlaran_(iseed);
row[i-ilow] = i;
}
lasta += ihigh - ilow + 1;
......@@ -83,7 +83,7 @@ cblockdiag(int nb, /* number of blocks */
val = &a[lasta];
row = &asub[lasta];
for (i = 0; i < bs; ++i) {
val[i].r = dlaran_(iseed);
val[i].re = dlaran_(iseed);
row[i] = i + rstart;
}
lasta += bs;
......
......@@ -101,8 +101,8 @@ int cReadValues(FILE *fp, int n, complex *destination, int perline, int persize)
pair = 1;
} else {
/* The value is imaginary part */
destination[i].r = realpart;
destination[i++].i = atof(&buf[s]);
destination[i].re = realpart;
destination[i++].im = atof(&buf[s]);
pair = 0;
}
buf[(j+1)*persize] = tmp; /* recover the char at that place */
......
......@@ -162,8 +162,8 @@ static int cReadValues(FILE *fp, int n, complex *destination, int perline, int p
pair = 1;
} else {
/* The value is imaginary part */
destination[i].r = realpart;
destination[i++].i = atof(&buf[s]);
destination[i].re = realpart;
destination[i++].im = atof(&buf[s]);
pair = 0;
}
buf[(j+1)*persize] = tmp; /* recover the char at that place */
......
......@@ -21,27 +21,27 @@ void z_div(doublecomplex *c, doublecomplex *a, doublecomplex *b)
double ratio, den;
double abr, abi, cr, ci;
if( (abr = b->r) < 0.)
if( (abr = b->re) < 0.)
abr = - abr;
if( (abi = b->i) < 0.)
if( (abi = b->im) < 0.)
abi = - abi;
if( abr <= abi ) {
if (abi == 0) {
fprintf(stderr, "z_div: division by zero\n");
exit(-1);
}
ratio = b->r / b->i ;
den = b->i * (1 + ratio*ratio);
cr = (a->r*ratio + a->i) / den;
ci = (a->i*ratio - a->r) / den;
ratio = b->re / b->im ;
den = b->im * (1 + ratio*ratio);
cr = (a->re*ratio + a->im) / den;
ci = (a->im*ratio - a->re) / den;
} else {
ratio = b->i / b->r ;
den = b->r * (1 + ratio*ratio);
cr = (a->r + a->i*ratio) / den;
ci = (a->i - a->r*ratio) / den;
ratio = b->im / b->re ;
den = b->re * (1 + ratio*ratio);
cr = (a->re + a->im*ratio) / den;
ci = (a->im - a->re*ratio) / den;
}
c->r = cr;
c->i = ci;
c->re = cr;
c->im = ci;
}
......@@ -49,8 +49,8 @@ void z_div(doublecomplex *c, doublecomplex *a, doublecomplex *b)
double z_abs(doublecomplex *z)
{
double temp;
double real = z->r;
double imag = z->i;
double real = z->re;
double imag = z->im;
if (real < 0) real = -real;
if (imag < 0) imag = -imag;
......@@ -71,8 +71,8 @@ double z_abs(doublecomplex *z)
/* Returns abs(z.r) + abs(z.i) */
double z_abs1(doublecomplex *z)
{
double real = z->r;
double imag = z->i;
double real = z->re;
double imag = z->im;
if (real < 0) real = -real;
if (imag < 0) imag = -imag;
......@@ -85,22 +85,22 @@ void z_exp(doublecomplex *r, doublecomplex *z)
{
double expx;
expx = exp(z->r);
r->r = expx * cos(z->i);
r->i = expx * sin(z->i);
expx = exp(z->re);
r->re = expx * cos(z->im);
r->im = expx * sin(z->im);
}
/* Return the complex conjugate */
void d_cnjg(doublecomplex *r, doublecomplex *z)
{
r->r = z->r;
r->i = -z->i;
r->re = z->re;
r->im = -z->im;
}
/* Return the imaginary part */
double d_imag(doublecomplex *z)
{
return (z->i);
return (z->im);
}
......@@ -68,17 +68,17 @@ int icmax1_(int *n, complex *cx, int *incx)
/* CODE FOR INCREMENT NOT EQUAL TO 1 */
ix = 1;
smax = (r__1 = CX(1).r, fabs(r__1));
smax = (r__1 = CX(1).re, fabs(r__1));
ix += *incx;
i__1 = *n;
for (i = 2; i <= *n; ++i) {
i__2 = ix;
if ((r__1 = CX(ix).r, fabs(r__1)) <= smax) {
if ((r__1 = CX(ix).re, fabs(r__1)) <= smax) {
goto L10;
}
ret_val = i;
i__2 = ix;
smax = (r__1 = CX(ix).r, fabs(r__1));
smax = (r__1 = CX(ix).re, fabs(r__1));
L10:
ix += *incx;
/* L20: */
......@@ -88,16 +88,16 @@ L10:
/* CODE FOR INCREMENT EQUAL TO 1 */
L30:
smax = (r__1 = CX(1).r, fabs(r__1));
smax = (r__1 = CX(1).re, fabs(r__1));
i__1 = *n;
for (i = 2; i <= *n; ++i) {
i__2 = i;
if ((r__1 = CX(i).r, fabs(r__1)) <= smax) {
if ((r__1 = CX(i).re, fabs(r__1)) <= smax) {
goto L40;
}
ret_val = i;
i__2 = i;
smax = (r__1 = CX(i).r, fabs(r__1));
smax = (r__1 = CX(i).re, fabs(r__1));
L40:
;
}
......
......@@ -62,17 +62,17 @@ izmax1_(int *n, doublecomplex *cx, int *incx)
/* CODE FOR INCREMENT NOT EQUAL TO 1 */
ix = 1;
smax = (d__1 = CX(1).r, fabs(d__1));
smax = (d__1 = CX(1).re, fabs(d__1));
ix += *incx;
i__1 = *n;
for (i = 2; i <= *n; ++i) {
i__2 = ix;
if ((d__1 = CX(ix).r, fabs(d__1)) <= smax) {
if ((d__1 = CX(ix).re, fabs(d__1)) <= smax) {
goto L10;
}
ret_val = i;
i__2 = ix;
smax = (d__1 = CX(ix).r, fabs(d__1));
smax = (d__1 = CX(ix).re, fabs(d__1));
L10:
ix += *incx;
/* L20: */
......@@ -82,16 +82,16 @@ L10:
/* CODE FOR INCREMENT EQUAL TO 1 */
L30:
smax = (d__1 = CX(1).r, fabs(d__1));
smax = (d__1 = CX(1).re, fabs(d__1));
i__1 = *n;
for (i = 2; i <= *n; ++i) {
i__2 = i;
if ((d__1 = CX(i).r, fabs(d__1)) <= smax) {
if ((d__1 = CX(i).re, fabs(d__1)) <= smax) {
goto L40;
}
ret_val = i;
i__2 = i;
smax = (d__1 = CX(i).r, fabs(d__1));
smax = (d__1 = CX(i).re, fabs(d__1));
L40:
;
}
......
......@@ -319,7 +319,7 @@ if ( jcol==BADCOL )
#else
for (kcol = fsupc; kcol <= krep; ++kcol) {
if ( (dense_col[inv_perm_r[kcol]].r != 0.0) || (dense_col[inv_perm_r[kcol]].i != 0.0) ) {
if ( (dense_col[inv_perm_r[kcol]].re != 0.0) || (dense_col[inv_perm_r[kcol]].im != 0.0) ) {
repfnz_col[krep] = kcol;
break; /* Found the leading nonzero in the U-segment */
}
......
......@@ -298,7 +298,7 @@ ccheck_zero_vec(int pnum, char *msg, int n, complex *vec)
nonzero = FALSE;
for (i = 0; i < n; ++i) {
if ((vec[i].r != 0.0) || (vec[i].i != 0.0))
if ((vec[i].re != 0.0) || (vec[i].im != 0.0))
{
printf("(%d) vec[%d] = %.10e; should be zero!\n",
pnum, i, vec[i]);
......@@ -318,8 +318,8 @@ cGenXtrue(int n, int nrhs, complex *x, int ldx)
int i, j;
for (j = 0; j < nrhs; ++j) {
for (i = 0; i < n; ++i) {
x[i + j*ldx].r = 1.0;
x[i + j*ldx].i = 0.0;
x[i + j*ldx].re = 1.0;
x[i + j*ldx].im = 0.0;
}
}
}
......@@ -458,7 +458,7 @@ int print_complex_vec(char *what, int n, int *ind, complex *vec)
{
int i;
printf("%s: n %d\n", what, n);
for (i = 0; i < n; ++i) printf("%d\t%f%f\n", ind[i], vec[i].r, vec[i].i);
for (i = 0; i < n; ++i) printf("%d\t%f%f\n", ind[i], vec[i].re, vec[i].im);
return 0;
}
......@@ -319,7 +319,7 @@ if ( jcol==BADCOL )
#else
for (kcol = fsupc; kcol <= krep; ++kcol) {
if ( (dense_col[inv_perm_r[kcol]].r != 0.0) || (dense_col[inv_perm_r[kcol]].i != 0.0) ) {
if ( (dense_col[inv_perm_r[kcol]].re != 0.0) || (dense_col[inv_perm_r[kcol]].im != 0.0) ) {
repfnz_col[krep] = kcol;
break; /* Found the leading nonzero in the U-segment */
}
......
......@@ -298,7 +298,7 @@ zcheck_zero_vec(int pnum, char *msg, int n, doublecomplex *vec)
nonzero = FALSE;
for (i = 0; i < n; ++i) {
if ((vec[i].r != 0.0) || (vec[i].i != 0.0))
if ((vec[i].re != 0.0) || (vec[i].im != 0.0))
{
printf("(%d) vec[%d] = %.10e; should be zero!\n",
pnum, i, vec[i]);
......@@ -318,8 +318,8 @@ zGenXtrue(int n, int nrhs, doublecomplex *x, int ldx)
int i, j;
for (j = 0; j < nrhs; ++j) {
for (i = 0; i < n; ++i) {
x[i + j*ldx].r = 1.0;
x[i + j*ldx].i = 0.0;
x[i + j*ldx].re = 1.0;
x[i + j*ldx].im = 0.0;
}
}
}
......@@ -458,7 +458,7 @@ int print_doublecomplex_vec(char *what, int n, int *ind, doublecomplex *vec)
{
int i;
printf("%s: n %d\n", what, n);
for (i = 0; i < n; ++i) printf("%d\t%f%f\n", ind[i], vec[i].r, vec[i].i);
for (i = 0; i < n; ++i) printf("%d\t%f%f\n", ind[i], vec[i].re, vec[i].im);
return 0;
}
......@@ -21,27 +21,27 @@ void c_div(complex *c, complex *a, complex *b)
float ratio, den;
float abr, abi, cr, ci;
if( (abr = b->r) < 0.)
if( (abr = b->re) < 0.)
abr = - abr;
if( (abi = b->i) < 0.)
if( (abi = b->im) < 0.)
abi = - abi;
if( abr <= abi ) {
if (abi == 0) {
fprintf(stderr, "c_div: division by zero\n");
exit(-1);
}
ratio = b->r / b->i ;
den = b->i * (1 + ratio*ratio);
cr = (a->r*ratio + a->i) / den;
ci = (a->i*ratio - a->r) / den;
ratio = b->re / b->im ;
den = b->im * (1 + ratio*ratio);
cr = (a->re*ratio + a->im) / den;
ci = (a->im*ratio - a->re) / den;
} else {
ratio = b->i / b->r ;
den = b->r * (1 + ratio*ratio);
cr = (a->r + a->i*ratio) / den;
ci = (a->i - a->r*ratio) / den;
ratio = b->im / b->re ;
den = b->re * (1 + ratio*ratio);
cr = (a->re + a->im*ratio) / den;
ci = (a->im - a->re*ratio) / den;
}
c->r = cr;
c->i = ci;
c->re = cr;
c->im = ci;
}
......@@ -49,8 +49,8 @@ void c_div(complex *c, complex *a, complex *b)
double c_abs(complex *z)
{
float temp;
float real = z->r;
float imag = z->i;
float real = z->re;
float imag = z->im;
if (real < 0) real = -real;
if (imag < 0) imag = -imag;
......@@ -71,8 +71,8 @@ double c_abs(complex *z)
/* Returns abs(z.r) + abs(z.i) */
double c_abs1(complex *z)
{
float real = z->r;
float imag = z->i;
float real = z->re;
float imag = z->im;
if (real < 0) real = -real;
if (imag < 0) imag = -imag;
......@@ -85,22 +85,22 @@ void c_exp(complex *r, complex *z)
{
float expx;
expx = exp(z->r);
r->r = expx * cos(z->i);
r->i = expx * sin(z->i);
expx = exp(z->re);
r->re = expx * cos(z->im);
r->im = expx * sin(z->im);
}
/* Return the complex conjugate */
void r_cnjg(complex *r, complex *z)
{
r->r = z->r;
r->i = -z->i;
r->re = z->re;
r->im = -z->im;
}
/* Return the imaginary part */
double r_imag(complex *z)
{
return (z->i);
return (z->im);
}
......@@ -15,39 +15,40 @@
#ifndef DCOMPLEX_INCLUDE
#define DCOMPLEX_INCLUDE
typedef struct { double r, i; } doublecomplex;
typedef struct { double re, im; } complex_t;
typedef complex_t doublecomplex;
/* Macro definitions */
/* Complex Addition c = a + b */
#define z_add(c, a, b) { (c)->r = (a)->r + (b)->r; \
(c)->i = (a)->i + (b)->i; }
#define z_add(c, a, b) { (c)->re = (a)->re + (b)->re; \
(c)->im = (a)->im + (b)->im; }
/* Complex Subtraction c = a - b */
#define z_sub(c, a, b) { (c)->r = (a)->r - (b)->r; \
(c)->i = (a)->i - (b)->i; }
#define z_sub(c, a, b) { (c)->re = (a)->re - (b)->re; \
(c)->im = (a)->im - (b)->im; }
/* Complex-Double Multiplication */
#define zd_mult(c, a, b) { (c)->r = (a)->r * (b); \
(c)->i = (a)->i * (b); }
#define zd_mult(c, a, b) { (c)->re = (a)->re * (b); \
(c)->im = (a)->im * (b); }
/* Complex-Complex Multiplication */
#define zz_mult(c, a, b) { \
double cr, ci; \
cr = (a)->r * (b)->r - (a)->i * (b)->i; \
ci = (a)->i * (b)->r + (a)->r * (b)->i; \
(c)->r = cr; \
(c)->i = ci; \
cr = (a)->re * (b)->re - (a)->im * (b)->im; \
ci = (a)->im * (b)->re + (a)->re * (b)->im; \
(c)->re = cr; \
(c)->im = ci; \
}
#define zz_conj(a, b) { \
(a)->r = (b)->r; \
(a)->i = -((b)->i); \
(a)->re = (b)->re; \
(a)->im = -((b)->im); \
}
/* Complex equality testing */
#define z_eq(a, b) ( (a)->r == (b)->r && (a)->i == (b)->i )
#define z_eq(a, b) ( (a)->re == (b)->re && (a)->im == (b)->im )
#ifdef __cplusplus
......
......@@ -15,39 +15,39 @@
#ifndef SCOMPLEX_INCLUDE
#define SCOMPLEX_INCLUDE
typedef struct { float r, i; } complex;
typedef struct { float re, im; } complex;
/* Macro definitions */
/* Complex Addition c = a + b */
#define c_add(c, a, b) { (c)->r = (a)->r + (b)->r; \
(c)->i = (a)->i + (b)->i; }
#define c_add(c, a, b) { (c)->re = (a)->re + (b)->re; \
(c)->im = (a)->im + (b)->im; }
/* Complex Subtraction c = a - b */
#define c_sub(c, a, b) { (c)->r = (a)->r - (b)->r; \
(c)->i = (a)->i - (b)->i; }
#define c_sub(c, a, b) { (c)->re = (a)->re - (b)->re; \
(c)->im = (a)->im - (b)->im; }
/* Complex-Double Multiplication */
#define cs_mult(c, a, b) { (c)->r = (a)->r * (b); \
(c)->i = (a)->i * (b); }
#define cs_mult(c, a, b) { (c)->re = (a)->re * (b); \
(c)->im = (a)->im * (b); }
/* Complex-Complex Multiplication */
#define cc_mult(c, a, b) { \
float cr, ci; \
cr = (a)->r * (b)->r - (a)->i * (b)->i; \
ci = (a)->i * (b)->r + (a)->r * (b)->i; \
(c)->r = cr; \
(c)->i = ci; \
cr = (a)->re * (b)->re - (a)->im * (b)->im; \
ci = (a)->im * (b)->re + (a)->re * (b)->im; \
(c)->re = cr; \
(c)->im = ci; \
}
#define cc_conj(a, b) { \
(a)->r = (b)->r; \
(a)->i = -((b)->i); \
(a)->re = (b)->re; \
(a)->im = -((b)->im); \
}
/* Complex equality testing */
#define c_eq(a, b) ( (a)->r == (b)->r && (a)->i == (b)->i )
#define c_eq(a, b) ( (a)->re == (b)->re && (a)->im == (b)->im )
#ifdef __cplusplus
......
......@@ -181,8 +181,8 @@ zgstrs(trans_t trans, SuperMatrix *L, SuperMatrix *U,
for (i = 0; i < nrow; i++) {
irow = L_SUB(iptr);
z_sub(&rhs_work[irow], &rhs_work[irow], &work_col[i]);
work_col[i].r = 0.0;
work_col[i].i = 0.0;
work_col[i].re = 0.0;
work_col[i].im = 0.0;
iptr++;
}
}
......@@ -197,8 +197,8 @@ zgstrs(trans_t trans, SuperMatrix *L, SuperMatrix *U,
for (i = 0; i < nrow; i++) {
irow = L_SUB(iptr);
z_sub(&rhs_work[irow], &rhs_work[irow], &work[i]);
work[i].r = 0.0;
work[i].i = 0.0;
work[i].re = 0.0;
work[i].im = 0.0;
iptr++;
}
}
......
......@@ -85,8 +85,8 @@ zlacon_(int *n, doublecomplex *v, doublecomplex *x, double *est, int *kase)
safmin = dlamch_("Safe minimum");
if ( *kase == 0 ) {
for (i = 0; i < *n; ++i) {
x[i].r = 1. / (double) (*n);
x[i].i = 0.;
x[i].re = 1. / (double) (*n);
x[i].im =