Skip to content

Commit c8b236f

Browse files
bosilcabosilca
authored
Merge pull request #147 from therault/fix-sumsq-updates
Fix missing corner cases when computing sum-square updates --
2 parents 771f834 + d3a8916 commit c8b236f

8 files changed

Lines changed: 213 additions & 221 deletions

File tree

‎src/cores/core_zgessq.c‎

Lines changed: 6 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -15,18 +15,7 @@
1515
#include <math.h>
1616
#include <lapacke.h>
1717
#include "common.h"
18-
19-
#define COMPLEX
20-
21-
#define UPDATE( __nb, __value ) \
22-
if (__value != 0. ){ \
23-
if ( *scale < __value ) { \
24-
*sumsq = __nb + (*sumsq) * ( *scale / __value ) * ( *scale / __value ); \
25-
*scale = __value; \
26-
} else { \
27-
*sumsq = *sumsq + __nb * ( __value / *scale ) * ( __value / *scale ); \
28-
} \
29-
}
18+
#include "sumsq_update.h"
3019

3120
/*****************************************************************************
3221
*
@@ -91,19 +80,16 @@ int CORE_zgessq(int M, int N,
9180
double *scale, double *sumsq)
9281
{
9382
int i, j;
94-
double tmp;
9583
double *ptr;
9684

9785
for(j=0; j<N; j++) {
9886
ptr = (double*) ( A + j * LDA );
9987
for(i=0; i<M; i++, ptr++) {
100-
tmp = fabs(*ptr);
101-
UPDATE( 1., tmp );
102-
103-
#ifdef COMPLEX
104-
ptr++;
105-
tmp = fabs(*ptr);
106-
UPDATE( 1., tmp );
88+
sumsq_update( 1, scale, sumsq, ptr );
89+
90+
#if defined(PRECISION_z) || defined(PRECISION_c)
91+
ptr++;
92+
sumsq_update( 1, scale, sumsq, ptr );
10793
#endif
10894
}
10995
}

‎src/cores/core_zhessq.c‎

Lines changed: 20 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -15,18 +15,7 @@
1515
#include <math.h>
1616
#include <lapacke.h>
1717
#include "common.h"
18-
19-
#define COMPLEX
20-
21-
#define UPDATE( __nb, __value ) \
22-
if (__value != 0. ){ \
23-
if ( *scale < __value ) { \
24-
*sumsq = __nb + (*sumsq) * ( *scale / __value ) * ( *scale / __value ); \
25-
*scale = __value; \
26-
} else { \
27-
*sumsq = *sumsq + __nb * ( __value / *scale ) * ( __value / *scale ); \
28-
} \
29-
}
18+
#include "sumsq_update.h"
3019

3120
/*****************************************************************************
3221
*
@@ -97,31 +86,26 @@ int CORE_zhessq(PLASMA_enum uplo, int N,
9786
double *scale, double *sumsq)
9887
{
9988
int i, j;
100-
double tmp;
10189
double *ptr;
10290

10391
if ( uplo == PlasmaUpper ) {
10492
for(j=0; j<N; j++) {
10593
ptr = (double*) ( A + j * LDA );
10694

10795
for(i=0; i<j; i++, ptr++) {
108-
109-
tmp = fabs(*ptr);
110-
UPDATE( 2., tmp );
111-
112-
#ifdef COMPLEX
113-
ptr++;
114-
tmp = fabs(*ptr);
115-
UPDATE( 2., tmp );
96+
sumsq_update( 2, scale, sumsq, ptr );
97+
98+
#if defined(PRECISION_z) || defined(PRECISION_c)
99+
ptr++;
100+
sumsq_update( 2, scale, sumsq, ptr );
116101
#endif
117102
}
118103

119104
/* Diagonal */
120-
tmp = fabs(*ptr);
121-
UPDATE( 1., tmp );
122-
123-
#ifdef COMPLEX
124-
ptr++;
105+
sumsq_update( 1, scale, sumsq, ptr );
106+
107+
#if defined(PRECISION_z) || defined(PRECISION_c)
108+
ptr++;
125109
#endif
126110
}
127111
} else {
@@ -130,23 +114,19 @@ int CORE_zhessq(PLASMA_enum uplo, int N,
130114
ptr = (double*) ( A + j * LDA + j);
131115

132116
/* Diagonal */
133-
tmp = fabs(*ptr);
134-
UPDATE( 1., tmp );
135-
ptr++;
136-
137-
#ifdef COMPLEX
138-
ptr++;
117+
sumsq_update( 1, scale, sumsq, ptr );
118+
ptr++;
119+
120+
#if defined(PRECISION_z) || defined(PRECISION_c)
121+
ptr++;
139122
#endif
140123

141124
for(i=j+1; i<N; i++, ptr++) {
142-
143-
tmp = fabs(*ptr);
144-
UPDATE( 2., tmp );
145-
146-
#ifdef COMPLEX
147-
ptr++;
148-
tmp = fabs(*ptr);
149-
UPDATE( 2., tmp );
125+
sumsq_update( 2, scale, sumsq, ptr );
126+
127+
#if defined(PRECISION_z) || defined(PRECISION_c)
128+
ptr++;
129+
sumsq_update( 2, scale, sumsq, ptr );
150130
#endif
151131
}
152132
}

‎src/cores/core_zsyssq.c‎

Lines changed: 22 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -15,18 +15,7 @@
1515
#include <math.h>
1616
#include <lapacke.h>
1717
#include "common.h"
18-
19-
#define COMPLEX
20-
21-
#define UPDATE( __nb, __value ) \
22-
if (__value != 0. ){ \
23-
if ( *scale < __value ) { \
24-
*sumsq = __nb + (*sumsq) * ( *scale / __value ) * ( *scale / __value ); \
25-
*scale = __value; \
26-
} else { \
27-
*sumsq = *sumsq + __nb * ( __value / *scale ) * ( __value / *scale ); \
28-
} \
29-
}
18+
#include "sumsq_update.h"
3019

3120
/*****************************************************************************
3221
*
@@ -97,33 +86,27 @@ int CORE_zsyssq(PLASMA_enum uplo, int N,
9786
double *scale, double *sumsq)
9887
{
9988
int i, j;
100-
double tmp;
10189
double *ptr;
10290

10391
if ( uplo == PlasmaUpper ) {
10492
for(j=0; j<N; j++) {
10593
ptr = (double*) ( A + j * LDA );
10694

10795
for(i=0; i<j; i++, ptr++) {
108-
109-
tmp = fabs(*ptr);
110-
UPDATE( 2., tmp );
111-
112-
#ifdef COMPLEX
113-
ptr++;
114-
tmp = fabs(*ptr);
115-
UPDATE( 2., tmp );
96+
sumsq_update( 2, scale, sumsq, ptr );
97+
98+
#if defined(PRECISION_z) || defined(PRECISION_c)
99+
ptr++;
100+
sumsq_update( 2, scale, sumsq, ptr );
116101
#endif
117102
}
118103

119104
/* Diagonal */
120-
tmp = fabs(*ptr);
121-
UPDATE( 1., tmp );
122-
123-
#ifdef COMPLEX
124-
ptr++;
125-
tmp = fabs(*ptr);
126-
UPDATE( 1., tmp );
105+
sumsq_update( 1, scale, sumsq, ptr );
106+
107+
#if defined(PRECISION_z) || defined(PRECISION_c)
108+
ptr++;
109+
sumsq_update( 1, scale, sumsq, ptr );
127110
#endif
128111
}
129112
} else {
@@ -132,25 +115,20 @@ int CORE_zsyssq(PLASMA_enum uplo, int N,
132115
ptr = (double*) ( A + j * LDA + j);
133116

134117
/* Diagonal */
135-
tmp = fabs(*ptr);
136-
UPDATE( 1., tmp );
137-
ptr++;
138-
139-
#ifdef COMPLEX
140-
tmp = fabs(*ptr);
141-
UPDATE( 1., tmp );
142-
ptr++;
118+
sumsq_update( 1, scale, sumsq, ptr );
119+
ptr++;
120+
121+
#if defined(PRECISION_z) || defined(PRECISION_c)
122+
sumsq_update( 1, scale, sumsq, ptr );
123+
ptr++;
143124
#endif
144125

145126
for(i=j+1; i<N; i++, ptr++) {
146-
147-
tmp = fabs(*ptr);
148-
UPDATE( 2., tmp );
149-
150-
#ifdef COMPLEX
151-
ptr++;
152-
tmp = fabs(*ptr);
153-
UPDATE( 2., tmp );
127+
sumsq_update( 2, scale, sumsq, ptr );
128+
129+
#if defined(PRECISION_z) || defined(PRECISION_c)
130+
ptr++;
131+
sumsq_update( 2, scale, sumsq, ptr );
154132
#endif
155133
}
156134
}

‎src/cores/core_ztrssq.c‎

Lines changed: 10 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -15,18 +15,7 @@
1515
#include <math.h>
1616
#include <lapacke.h>
1717
#include "common.h"
18-
19-
#define COMPLEX
20-
21-
#define UPDATE( __nb, __value ) \
22-
if (__value != 0. ){ \
23-
if ( *scale < __value ) { \
24-
*sumsq = __nb + (*sumsq) * ( *scale / __value ) * ( *scale / __value ); \
25-
*scale = __value; \
26-
} else { \
27-
*sumsq = *sumsq + __nb * ( __value / *scale ) * ( __value / *scale ); \
28-
} \
29-
}
18+
#include "sumsq_update.h"
3019

3120
/*****************************************************************************
3221
*
@@ -97,7 +86,7 @@ int CORE_ztrssq(PLASMA_enum uplo, PLASMA_enum diag, int M, int N,
9786

9887
if ( diag == PlasmaUnit ){
9988
tmp = sqrt( min(M, N) );
100-
UPDATE( 1., tmp );
89+
sumsq_update( 1, scale, sumsq, &tmp );
10190
}
10291

10392
if (uplo == PlasmaUpper ) {
@@ -108,13 +97,11 @@ int CORE_ztrssq(PLASMA_enum uplo, PLASMA_enum diag, int M, int N,
10897
imax = min(j+1-idiag, M);
10998

11099
for(i=0; i<imax; i++, ptr++) {
111-
tmp = fabs(*ptr);
112-
UPDATE( 1., tmp );
113-
114-
#ifdef COMPLEX
100+
sumsq_update( 1, scale, sumsq, ptr );
101+
102+
#if defined(PRECISION_z) || defined(PRECISION_c)
115103
ptr++;
116-
tmp = fabs(*ptr);
117-
UPDATE( 1., tmp );
104+
sumsq_update( 1, scale, sumsq, ptr );
118105
#endif
119106
}
120107
}
@@ -126,13 +113,11 @@ int CORE_ztrssq(PLASMA_enum uplo, PLASMA_enum diag, int M, int N,
126113
ptr = (double*) ( A + j * (LDA+1) + idiag );
127114

128115
for(i=j+idiag; i<M; i++, ptr++) {
129-
tmp = fabs(*ptr);
130-
UPDATE( 1., tmp );
131-
132-
#ifdef COMPLEX
116+
sumsq_update( 1, scale, sumsq, ptr );
117+
118+
#if defined(PRECISION_z) || defined(PRECISION_c)
133119
ptr++;
134-
tmp = fabs(*ptr);
135-
UPDATE( 1., tmp );
120+
sumsq_update( 1, scale, sumsq, ptr );
136121
#endif
137122
}
138123
}

0 commit comments

Comments
 (0)