00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049 #include <stdio.h>
00050 #include "lapacke.h"
00051 #include "lapacke_utils.h"
00052 #include "test_utils.h"
00053
00054 static void init_scalars_dgehrd( lapack_int *n, lapack_int *ilo,
00055 lapack_int *ihi, lapack_int *lda,
00056 lapack_int *lwork );
00057 static void init_a( lapack_int size, double *a );
00058 static void init_tau( lapack_int size, double *tau );
00059 static void init_work( lapack_int size, double *work );
00060 static int compare_dgehrd( double *a, double *a_i, double *tau, double *tau_i,
00061 lapack_int info, lapack_int info_i, lapack_int lda,
00062 lapack_int n );
00063
00064 int main(void)
00065 {
00066
00067 lapack_int n, n_i;
00068 lapack_int ilo, ilo_i;
00069 lapack_int ihi, ihi_i;
00070 lapack_int lda, lda_i;
00071 lapack_int lda_r;
00072 lapack_int lwork, lwork_i;
00073 lapack_int info, info_i;
00074 lapack_int i;
00075 int failed;
00076
00077
00078 double *a = NULL, *a_i = NULL;
00079 double *tau = NULL, *tau_i = NULL;
00080 double *work = NULL, *work_i = NULL;
00081 double *a_save = NULL;
00082 double *tau_save = NULL;
00083 double *a_r = NULL;
00084
00085
00086 init_scalars_dgehrd( &n, &ilo, &ihi, &lda, &lwork );
00087 lda_r = n+2;
00088 n_i = n;
00089 ilo_i = ilo;
00090 ihi_i = ihi;
00091 lda_i = lda;
00092 lwork_i = lwork;
00093
00094
00095 a = (double *)LAPACKE_malloc( lda*n * sizeof(double) );
00096 tau = (double *)LAPACKE_malloc( (n-1) * sizeof(double) );
00097 work = (double *)LAPACKE_malloc( lwork * sizeof(double) );
00098
00099
00100 a_i = (double *)LAPACKE_malloc( lda*n * sizeof(double) );
00101 tau_i = (double *)LAPACKE_malloc( (n-1) * sizeof(double) );
00102 work_i = (double *)LAPACKE_malloc( lwork * sizeof(double) );
00103
00104
00105 a_save = (double *)LAPACKE_malloc( lda*n * sizeof(double) );
00106 tau_save = (double *)LAPACKE_malloc( (n-1) * sizeof(double) );
00107
00108
00109 a_r = (double *)LAPACKE_malloc( n*(n+2) * sizeof(double) );
00110
00111
00112 init_a( lda*n, a );
00113 init_tau( (n-1), tau );
00114 init_work( lwork, work );
00115
00116
00117 for( i = 0; i < lda*n; i++ ) {
00118 a_save[i] = a[i];
00119 }
00120 for( i = 0; i < (n-1); i++ ) {
00121 tau_save[i] = tau[i];
00122 }
00123
00124
00125 dgehrd_( &n, &ilo, &ihi, a, &lda, tau, work, &lwork, &info );
00126
00127
00128
00129 for( i = 0; i < lda*n; i++ ) {
00130 a_i[i] = a_save[i];
00131 }
00132 for( i = 0; i < (n-1); i++ ) {
00133 tau_i[i] = tau_save[i];
00134 }
00135 for( i = 0; i < lwork; i++ ) {
00136 work_i[i] = work[i];
00137 }
00138 info_i = LAPACKE_dgehrd_work( LAPACK_COL_MAJOR, n_i, ilo_i, ihi_i, a_i,
00139 lda_i, tau_i, work_i, lwork_i );
00140
00141 failed = compare_dgehrd( a, a_i, tau, tau_i, info, info_i, lda, n );
00142 if( failed == 0 ) {
00143 printf( "PASSED: column-major middle-level interface to dgehrd\n" );
00144 } else {
00145 printf( "FAILED: column-major middle-level interface to dgehrd\n" );
00146 }
00147
00148
00149
00150 for( i = 0; i < lda*n; i++ ) {
00151 a_i[i] = a_save[i];
00152 }
00153 for( i = 0; i < (n-1); i++ ) {
00154 tau_i[i] = tau_save[i];
00155 }
00156 for( i = 0; i < lwork; i++ ) {
00157 work_i[i] = work[i];
00158 }
00159 info_i = LAPACKE_dgehrd( LAPACK_COL_MAJOR, n_i, ilo_i, ihi_i, a_i, lda_i,
00160 tau_i );
00161
00162 failed = compare_dgehrd( a, a_i, tau, tau_i, info, info_i, lda, n );
00163 if( failed == 0 ) {
00164 printf( "PASSED: column-major high-level interface to dgehrd\n" );
00165 } else {
00166 printf( "FAILED: column-major high-level interface to dgehrd\n" );
00167 }
00168
00169
00170
00171 for( i = 0; i < lda*n; i++ ) {
00172 a_i[i] = a_save[i];
00173 }
00174 for( i = 0; i < (n-1); i++ ) {
00175 tau_i[i] = tau_save[i];
00176 }
00177 for( i = 0; i < lwork; i++ ) {
00178 work_i[i] = work[i];
00179 }
00180
00181 LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, n, a_i, lda, a_r, n+2 );
00182 info_i = LAPACKE_dgehrd_work( LAPACK_ROW_MAJOR, n_i, ilo_i, ihi_i, a_r,
00183 lda_r, tau_i, work_i, lwork_i );
00184
00185 LAPACKE_dge_trans( LAPACK_ROW_MAJOR, n, n, a_r, n+2, a_i, lda );
00186
00187 failed = compare_dgehrd( a, a_i, tau, tau_i, info, info_i, lda, n );
00188 if( failed == 0 ) {
00189 printf( "PASSED: row-major middle-level interface to dgehrd\n" );
00190 } else {
00191 printf( "FAILED: row-major middle-level interface to dgehrd\n" );
00192 }
00193
00194
00195
00196 for( i = 0; i < lda*n; i++ ) {
00197 a_i[i] = a_save[i];
00198 }
00199 for( i = 0; i < (n-1); i++ ) {
00200 tau_i[i] = tau_save[i];
00201 }
00202 for( i = 0; i < lwork; i++ ) {
00203 work_i[i] = work[i];
00204 }
00205
00206
00207 LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, n, a_i, lda, a_r, n+2 );
00208 info_i = LAPACKE_dgehrd( LAPACK_ROW_MAJOR, n_i, ilo_i, ihi_i, a_r, lda_r,
00209 tau_i );
00210
00211 LAPACKE_dge_trans( LAPACK_ROW_MAJOR, n, n, a_r, n+2, a_i, lda );
00212
00213 failed = compare_dgehrd( a, a_i, tau, tau_i, info, info_i, lda, n );
00214 if( failed == 0 ) {
00215 printf( "PASSED: row-major high-level interface to dgehrd\n" );
00216 } else {
00217 printf( "FAILED: row-major high-level interface to dgehrd\n" );
00218 }
00219
00220
00221 if( a != NULL ) {
00222 LAPACKE_free( a );
00223 }
00224 if( a_i != NULL ) {
00225 LAPACKE_free( a_i );
00226 }
00227 if( a_r != NULL ) {
00228 LAPACKE_free( a_r );
00229 }
00230 if( a_save != NULL ) {
00231 LAPACKE_free( a_save );
00232 }
00233 if( tau != NULL ) {
00234 LAPACKE_free( tau );
00235 }
00236 if( tau_i != NULL ) {
00237 LAPACKE_free( tau_i );
00238 }
00239 if( tau_save != NULL ) {
00240 LAPACKE_free( tau_save );
00241 }
00242 if( work != NULL ) {
00243 LAPACKE_free( work );
00244 }
00245 if( work_i != NULL ) {
00246 LAPACKE_free( work_i );
00247 }
00248
00249 return 0;
00250 }
00251
00252
00253 static void init_scalars_dgehrd( lapack_int *n, lapack_int *ilo,
00254 lapack_int *ihi, lapack_int *lda,
00255 lapack_int *lwork )
00256 {
00257 *n = 4;
00258 *ilo = 2;
00259 *ihi = 4;
00260 *lda = 8;
00261 *lwork = 512;
00262
00263 return;
00264 }
00265
00266
00267 static void init_a( lapack_int size, double *a ) {
00268 lapack_int i;
00269 for( i = 0; i < size; i++ ) {
00270 a[i] = 0;
00271 }
00272 a[0] = -4.00000000000000020e-001;
00273 a[8] = 3.20000000000000020e+000;
00274 a[16] = 7.59999999999999960e+000;
00275 a[24] = -1.50000000000000000e+000;
00276 a[1] = 0.00000000000000000e+000;
00277 a[9] = 2.00000000000000010e-001;
00278 a[17] = 9.10000000000000030e-001;
00279 a[25] = 4.31250000000000000e+000;
00280 a[2] = 0.00000000000000000e+000;
00281 a[10] = 9.10000000000000030e-001;
00282 a[18] = 5.13999999999999970e+000;
00283 a[26] = -4.09999999999999960e+000;
00284 a[3] = 0.00000000000000000e+000;
00285 a[11] = 2.79999999999999980e+000;
00286 a[19] = -2.64000000000000010e+000;
00287 a[27] = 6.60000000000000030e-001;
00288 }
00289 static void init_tau( lapack_int size, double *tau ) {
00290 lapack_int i;
00291 for( i = 0; i < size; i++ ) {
00292 tau[i] = 0;
00293 }
00294 }
00295 static void init_work( lapack_int size, double *work ) {
00296 lapack_int i;
00297 for( i = 0; i < size; i++ ) {
00298 work[i] = 0;
00299 }
00300 }
00301
00302
00303
00304 static int compare_dgehrd( double *a, double *a_i, double *tau, double *tau_i,
00305 lapack_int info, lapack_int info_i, lapack_int lda,
00306 lapack_int n )
00307 {
00308 lapack_int i;
00309 int failed = 0;
00310 for( i = 0; i < lda*n; i++ ) {
00311 failed += compare_doubles(a[i],a_i[i]);
00312 }
00313 for( i = 0; i < (n-1); i++ ) {
00314 failed += compare_doubles(tau[i],tau_i[i]);
00315 }
00316 failed += (info == info_i) ? 0 : 1;
00317 if( info != 0 || info_i != 0 ) {
00318 printf( "info=%d, info_i=%d\n",(int)info,(int)info_i );
00319 }
00320
00321 return failed;
00322 }