1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
|
/*
-- MAGMA (version 2.9.0) --
Univ. of Tennessee, Knoxville
Univ. of California, Berkeley
Univ. of Colorado, Denver
@date January 2025
@author Mark Gates
*/
// Checks vector and matrix norms, for -m32 and -m64, with return as float or as double.
// In MacOS, -m64 must return double for both single and double precision
// functions, e.g., {s,d}dot, {s,d}nrm2, {s,d}lange, {s,d}lansy.
// Oddly, with -m32 both return float and double for single precision functions works.
// This is essentially a bug from an old f2c version of lapack (clapack).
//
// We work around this bug by putting replacement routines in the magma/lib/libblas_fix.a library.
// These correctly return float for the single precision functions that had issues.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
// hack to NOT include magma_lapack.h, by pre-defining MAGMA_LAPACK_H.
// we re-define the lapack prototypes below, so can't include that header.
#define MAGMA_LAPACK_H
#include "testings.h"
#include "magma_mangling.h"
// ------------------------------------------------------------
//#define LAPACK_RETURN_DOUBLE
#ifdef LAPACK_RETURN_DOUBLE
typedef double RETURN_FLOAT;
#else
typedef float RETURN_FLOAT;
#endif
// ------------------------------------------------------------
#ifdef __cplusplus
extern "C" {
#endif
#define blasf77_sdot FORTRAN_NAME( sdot, SDOT )
#define blasf77_snrm2 FORTRAN_NAME( snrm2, SNRM2 )
#define lapackf77_slange FORTRAN_NAME( slange, SLANGE )
#define lapackf77_slansy FORTRAN_NAME( slansy, SLANSY )
#define blasf77_ddot FORTRAN_NAME( ddot, DDOT )
#define blasf77_dnrm2 FORTRAN_NAME( dnrm2, DNRM2 )
#define lapackf77_dlange FORTRAN_NAME( dlange, DLANGE )
#define lapackf77_dlansy FORTRAN_NAME( dlansy, DLANSY )
RETURN_FLOAT
blasf77_sdot( const magma_int_t *n,
const float *x, const magma_int_t *incx,
const float *y, const magma_int_t *incy );
RETURN_FLOAT
blasf77_snrm2( const magma_int_t *n,
const float *x, const magma_int_t *incx );
RETURN_FLOAT
lapackf77_slange( const char *norm,
const magma_int_t *m, const magma_int_t *n,
const float *A, const magma_int_t *lda,
float *work );
RETURN_FLOAT
lapackf77_slansy( const char *norm, const char* uplo,
const magma_int_t *n,
const float *A, const magma_int_t *lda,
float *work );
double
blasf77_ddot( const magma_int_t *n,
const double *x, const magma_int_t *incx,
const double *y, const magma_int_t *incy );
double
blasf77_dnrm2( const magma_int_t *n,
const double *x, const magma_int_t *incx );
double
lapackf77_dlange( const char *norm,
const magma_int_t *m, const magma_int_t *n,
const double *A, const magma_int_t *lda,
double *work );
double
lapackf77_dlansy( const char *norm, const char* uplo,
const magma_int_t *n,
const double *A, const magma_int_t *lda,
double *work );
#ifdef __cplusplus
}
#endif
// ------------------------------------------------------------
// call matrix norms {s,d}lan{ge,sy}.
// return value, to check that the call stack isn't messed up.
float test( magma_int_t m, magma_int_t n )
{
#define sA(i,j) (sA + (i) + (j)*lda)
#define dA(i,j) (dA + (i) + (j)*lda)
float *sA, *swork;
float snorm_one, snorm_inf, snorm_fro, snorm_max;
double *dA, *dwork;
double dnorm_one, dnorm_inf, dnorm_fro, dnorm_max;
const magma_int_t ione = 1;
magma_int_t lda = max(m,n);
TESTING_CHECK( magma_smalloc_cpu( &&sA, lda*n ));
TESTING_CHECK( magma_dmalloc_cpu( &&dA, lda*n ));
TESTING_CHECK( magma_smalloc_cpu( &&swork, m ));
TESTING_CHECK( magma_dmalloc_cpu( &&dwork, m ));
for( magma_int_t j = 0; j < n; ++j ) {
for( magma_int_t i = 0; i < lda; ++i ) {
double tmp = rand() / (double)(RAND_MAX);
*sA(i,j) = tmp;
*dA(i,j) = tmp;
}
}
double error;
int status = 0;
// can repeat multiple times, but shows same results every time
status = 0;
for( magma_int_t i=0; i < 1; ++i ) {
snorm_one = blasf77_sdot( &m, sA, &ione, sA, &ione );
dnorm_one = blasf77_ddot( &m, dA, &ione, dA, &ione );
snorm_fro = blasf77_snrm2( &m, sA, &ione );
dnorm_fro = blasf77_dnrm2( &m, dA, &ione );
printf( "m %lld, sdot %12.8f, snrm2 %12.8f\n", (long long) m, snorm_one, snorm_fro );
printf( "m %lld, ddot %12.8f, dnrm2 %12.8f\n", (long long) m, dnorm_one, dnorm_fro );
error = fabs(snorm_one - dnorm_one) / dnorm_one;
status |= ! (error < 1e-6);
}
if ( status ) {
printf( "**** failed ****\n" );
}
else {
printf( "ok\n" );
}
printf( "\n" );
status = 0;
for( magma_int_t i=0; i < 1; ++i ) {
snorm_one = lapackf77_slange( "one", &m, &n, sA, &lda, swork );
snorm_inf = lapackf77_slange( "inf", &m, &n, sA, &lda, swork );
snorm_max = lapackf77_slange( "max", &m, &n, sA, &lda, swork );
snorm_fro = lapackf77_slange( "fro", &m, &n, sA, &lda, swork );
dnorm_one = lapackf77_dlange( "one", &m, &n, dA, &lda, dwork );
dnorm_inf = lapackf77_dlange( "inf", &m, &n, dA, &lda, dwork );
dnorm_max = lapackf77_dlange( "max", &m, &n, dA, &lda, dwork );
dnorm_fro = lapackf77_dlange( "fro", &m, &n, dA, &lda, dwork );
printf( "m %lld, n %lld, slange norm one %12.8f, inf %12.8f, max %12.8f, fro %12.8f\n",
(long long) m, (long long) n, snorm_one, snorm_inf, snorm_max, snorm_fro );
printf( "m %lld, n %lld, dlange norm one %12.8f, inf %12.8f, max %12.8f, fro %12.8f\n",
(long long) m, (long long) n, dnorm_one, dnorm_inf, dnorm_max, dnorm_fro );
error = fabs(snorm_one - dnorm_one) / dnorm_one;
status |= ! (error < 1e-6);
}
if ( status ) {
printf( "**** failed ****\n" );
}
else {
printf( "ok\n" );
}
printf( "\n" );
status = 0;
for( magma_int_t i=0; i < 1; ++i ) {
snorm_one = lapackf77_slansy( "one", "up", &n, sA, &lda, swork );
snorm_inf = lapackf77_slansy( "inf", "up", &n, sA, &lda, swork );
snorm_max = lapackf77_slansy( "max", "up", &n, sA, &lda, swork );
snorm_fro = lapackf77_slansy( "fro", "up", &n, sA, &lda, swork );
dnorm_one = lapackf77_dlansy( "one", "up", &n, dA, &lda, dwork );
dnorm_inf = lapackf77_dlansy( "inf", "up", &n, dA, &lda, dwork );
dnorm_max = lapackf77_dlansy( "max", "up", &n, dA, &lda, dwork );
dnorm_fro = lapackf77_dlansy( "fro", "up", &n, dA, &lda, dwork );
printf( "m %lld, n %lld, slansy norm one %12.8f, inf %12.8f, max %12.8f, fro %12.8f\n",
(long long) m, (long long) n, snorm_one, snorm_inf, snorm_max, snorm_fro );
printf( "m %lld, n %lld, dlansy norm one %12.8f, inf %12.8f, max %12.8f, fro %12.8f\n",
(long long) m, (long long) n, dnorm_one, dnorm_inf, dnorm_max, dnorm_fro );
error = fabs(snorm_one - dnorm_one) / dnorm_one;
status |= ! (error < 1e-6);
}
if ( status ) {
printf( "**** failed ****\n" );
}
else {
printf( "ok\n" );
}
printf( "\n" );
magma_free_cpu( sA );
magma_free_cpu( dA );
magma_free_cpu( swork );
magma_free_cpu( dwork );
return 1.125;
}
// ------------------------------------------------------------
int main( int argc, char** argv )
{
magma_int_t m = 100;
magma_int_t n = m;
if ( argc > 1 ) {
n = atoi( argv[1] );
}
float value;
printf( "--------------------\n" );
printf( "sizeof(void*) %lu, sizeof(RETURN_FLOAT) %lu\n\n",
sizeof(void*), sizeof(RETURN_FLOAT) );
// can repeat multiple times, but shows same results every time
for( magma_int_t i=0; i < 1; ++i ) {
value = test( m, n );
printf( "value %.4f\n\n", value );
}
return 0;
}
|
