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_dsyrfs( char *uplo, lapack_int *n, lapack_int *nrhs,
00055 lapack_int *lda, lapack_int *ldaf,
00056 lapack_int *ldb, lapack_int *ldx );
00057 static void init_a( lapack_int size, double *a );
00058 static void init_af( lapack_int size, double *af );
00059 static void init_ipiv( lapack_int size, lapack_int *ipiv );
00060 static void init_b( lapack_int size, double *b );
00061 static void init_x( lapack_int size, double *x );
00062 static void init_ferr( lapack_int size, double *ferr );
00063 static void init_berr( lapack_int size, double *berr );
00064 static void init_work( lapack_int size, double *work );
00065 static void init_iwork( lapack_int size, lapack_int *iwork );
00066 static int compare_dsyrfs( double *x, double *x_i, double *ferr, double *ferr_i,
00067 double *berr, double *berr_i, lapack_int info,
00068 lapack_int info_i, lapack_int ldx, lapack_int nrhs );
00069
00070 int main(void)
00071 {
00072
00073 char uplo, uplo_i;
00074 lapack_int n, n_i;
00075 lapack_int nrhs, nrhs_i;
00076 lapack_int lda, lda_i;
00077 lapack_int lda_r;
00078 lapack_int ldaf, ldaf_i;
00079 lapack_int ldaf_r;
00080 lapack_int ldb, ldb_i;
00081 lapack_int ldb_r;
00082 lapack_int ldx, ldx_i;
00083 lapack_int ldx_r;
00084 lapack_int info, info_i;
00085 lapack_int i;
00086 int failed;
00087
00088
00089 double *a = NULL, *a_i = NULL;
00090 double *af = NULL, *af_i = NULL;
00091 lapack_int *ipiv = NULL, *ipiv_i = NULL;
00092 double *b = NULL, *b_i = NULL;
00093 double *x = NULL, *x_i = NULL;
00094 double *ferr = NULL, *ferr_i = NULL;
00095 double *berr = NULL, *berr_i = NULL;
00096 double *work = NULL, *work_i = NULL;
00097 lapack_int *iwork = NULL, *iwork_i = NULL;
00098 double *x_save = NULL;
00099 double *ferr_save = NULL;
00100 double *berr_save = NULL;
00101 double *a_r = NULL;
00102 double *af_r = NULL;
00103 double *b_r = NULL;
00104 double *x_r = NULL;
00105
00106
00107 init_scalars_dsyrfs( &uplo, &n, &nrhs, &lda, &ldaf, &ldb, &ldx );
00108 lda_r = n+2;
00109 ldaf_r = n+2;
00110 ldb_r = nrhs+2;
00111 ldx_r = nrhs+2;
00112 uplo_i = uplo;
00113 n_i = n;
00114 nrhs_i = nrhs;
00115 lda_i = lda;
00116 ldaf_i = ldaf;
00117 ldb_i = ldb;
00118 ldx_i = ldx;
00119
00120
00121 a = (double *)LAPACKE_malloc( lda*n * sizeof(double) );
00122 af = (double *)LAPACKE_malloc( ldaf*n * sizeof(double) );
00123 ipiv = (lapack_int *)LAPACKE_malloc( n * sizeof(lapack_int) );
00124 b = (double *)LAPACKE_malloc( ldb*nrhs * sizeof(double) );
00125 x = (double *)LAPACKE_malloc( ldx*nrhs * sizeof(double) );
00126 ferr = (double *)LAPACKE_malloc( nrhs * sizeof(double) );
00127 berr = (double *)LAPACKE_malloc( nrhs * sizeof(double) );
00128 work = (double *)LAPACKE_malloc( 3*n * sizeof(double) );
00129 iwork = (lapack_int *)LAPACKE_malloc( n * sizeof(lapack_int) );
00130
00131
00132 a_i = (double *)LAPACKE_malloc( lda*n * sizeof(double) );
00133 af_i = (double *)LAPACKE_malloc( ldaf*n * sizeof(double) );
00134 ipiv_i = (lapack_int *)LAPACKE_malloc( n * sizeof(lapack_int) );
00135 b_i = (double *)LAPACKE_malloc( ldb*nrhs * sizeof(double) );
00136 x_i = (double *)LAPACKE_malloc( ldx*nrhs * sizeof(double) );
00137 ferr_i = (double *)LAPACKE_malloc( nrhs * sizeof(double) );
00138 berr_i = (double *)LAPACKE_malloc( nrhs * sizeof(double) );
00139 work_i = (double *)LAPACKE_malloc( 3*n * sizeof(double) );
00140 iwork_i = (lapack_int *)LAPACKE_malloc( n * sizeof(lapack_int) );
00141
00142
00143 x_save = (double *)LAPACKE_malloc( ldx*nrhs * sizeof(double) );
00144 ferr_save = (double *)LAPACKE_malloc( nrhs * sizeof(double) );
00145 berr_save = (double *)LAPACKE_malloc( nrhs * sizeof(double) );
00146
00147
00148 a_r = (double *)LAPACKE_malloc( n*(n+2) * sizeof(double) );
00149 af_r = (double *)LAPACKE_malloc( n*(n+2) * sizeof(double) );
00150 b_r = (double *)LAPACKE_malloc( n*(nrhs+2) * sizeof(double) );
00151 x_r = (double *)LAPACKE_malloc( n*(nrhs+2) * sizeof(double) );
00152
00153
00154 init_a( lda*n, a );
00155 init_af( ldaf*n, af );
00156 init_ipiv( n, ipiv );
00157 init_b( ldb*nrhs, b );
00158 init_x( ldx*nrhs, x );
00159 init_ferr( nrhs, ferr );
00160 init_berr( nrhs, berr );
00161 init_work( 3*n, work );
00162 init_iwork( n, iwork );
00163
00164
00165 for( i = 0; i < ldx*nrhs; i++ ) {
00166 x_save[i] = x[i];
00167 }
00168 for( i = 0; i < nrhs; i++ ) {
00169 ferr_save[i] = ferr[i];
00170 }
00171 for( i = 0; i < nrhs; i++ ) {
00172 berr_save[i] = berr[i];
00173 }
00174
00175
00176 dsyrfs_( &uplo, &n, &nrhs, a, &lda, af, &ldaf, ipiv, b, &ldb, x, &ldx, ferr,
00177 berr, work, iwork, &info );
00178
00179
00180
00181 for( i = 0; i < lda*n; i++ ) {
00182 a_i[i] = a[i];
00183 }
00184 for( i = 0; i < ldaf*n; i++ ) {
00185 af_i[i] = af[i];
00186 }
00187 for( i = 0; i < n; i++ ) {
00188 ipiv_i[i] = ipiv[i];
00189 }
00190 for( i = 0; i < ldb*nrhs; i++ ) {
00191 b_i[i] = b[i];
00192 }
00193 for( i = 0; i < ldx*nrhs; i++ ) {
00194 x_i[i] = x_save[i];
00195 }
00196 for( i = 0; i < nrhs; i++ ) {
00197 ferr_i[i] = ferr_save[i];
00198 }
00199 for( i = 0; i < nrhs; i++ ) {
00200 berr_i[i] = berr_save[i];
00201 }
00202 for( i = 0; i < 3*n; i++ ) {
00203 work_i[i] = work[i];
00204 }
00205 for( i = 0; i < n; i++ ) {
00206 iwork_i[i] = iwork[i];
00207 }
00208 info_i = LAPACKE_dsyrfs_work( LAPACK_COL_MAJOR, uplo_i, n_i, nrhs_i, a_i,
00209 lda_i, af_i, ldaf_i, ipiv_i, b_i, ldb_i, x_i,
00210 ldx_i, ferr_i, berr_i, work_i, iwork_i );
00211
00212 failed = compare_dsyrfs( x, x_i, ferr, ferr_i, berr, berr_i, info, info_i,
00213 ldx, nrhs );
00214 if( failed == 0 ) {
00215 printf( "PASSED: column-major middle-level interface to dsyrfs\n" );
00216 } else {
00217 printf( "FAILED: column-major middle-level interface to dsyrfs\n" );
00218 }
00219
00220
00221
00222 for( i = 0; i < lda*n; i++ ) {
00223 a_i[i] = a[i];
00224 }
00225 for( i = 0; i < ldaf*n; i++ ) {
00226 af_i[i] = af[i];
00227 }
00228 for( i = 0; i < n; i++ ) {
00229 ipiv_i[i] = ipiv[i];
00230 }
00231 for( i = 0; i < ldb*nrhs; i++ ) {
00232 b_i[i] = b[i];
00233 }
00234 for( i = 0; i < ldx*nrhs; i++ ) {
00235 x_i[i] = x_save[i];
00236 }
00237 for( i = 0; i < nrhs; i++ ) {
00238 ferr_i[i] = ferr_save[i];
00239 }
00240 for( i = 0; i < nrhs; i++ ) {
00241 berr_i[i] = berr_save[i];
00242 }
00243 for( i = 0; i < 3*n; i++ ) {
00244 work_i[i] = work[i];
00245 }
00246 for( i = 0; i < n; i++ ) {
00247 iwork_i[i] = iwork[i];
00248 }
00249 info_i = LAPACKE_dsyrfs( LAPACK_COL_MAJOR, uplo_i, n_i, nrhs_i, a_i, lda_i,
00250 af_i, ldaf_i, ipiv_i, b_i, ldb_i, x_i, ldx_i,
00251 ferr_i, berr_i );
00252
00253 failed = compare_dsyrfs( x, x_i, ferr, ferr_i, berr, berr_i, info, info_i,
00254 ldx, nrhs );
00255 if( failed == 0 ) {
00256 printf( "PASSED: column-major high-level interface to dsyrfs\n" );
00257 } else {
00258 printf( "FAILED: column-major high-level interface to dsyrfs\n" );
00259 }
00260
00261
00262
00263 for( i = 0; i < lda*n; i++ ) {
00264 a_i[i] = a[i];
00265 }
00266 for( i = 0; i < ldaf*n; i++ ) {
00267 af_i[i] = af[i];
00268 }
00269 for( i = 0; i < n; i++ ) {
00270 ipiv_i[i] = ipiv[i];
00271 }
00272 for( i = 0; i < ldb*nrhs; i++ ) {
00273 b_i[i] = b[i];
00274 }
00275 for( i = 0; i < ldx*nrhs; i++ ) {
00276 x_i[i] = x_save[i];
00277 }
00278 for( i = 0; i < nrhs; i++ ) {
00279 ferr_i[i] = ferr_save[i];
00280 }
00281 for( i = 0; i < nrhs; i++ ) {
00282 berr_i[i] = berr_save[i];
00283 }
00284 for( i = 0; i < 3*n; i++ ) {
00285 work_i[i] = work[i];
00286 }
00287 for( i = 0; i < n; i++ ) {
00288 iwork_i[i] = iwork[i];
00289 }
00290
00291 LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, n, a_i, lda, a_r, n+2 );
00292 LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, n, af_i, ldaf, af_r, n+2 );
00293 LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, nrhs, b_i, ldb, b_r, nrhs+2 );
00294 LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, nrhs, x_i, ldx, x_r, nrhs+2 );
00295 info_i = LAPACKE_dsyrfs_work( LAPACK_ROW_MAJOR, uplo_i, n_i, nrhs_i, a_r,
00296 lda_r, af_r, ldaf_r, ipiv_i, b_r, ldb_r, x_r,
00297 ldx_r, ferr_i, berr_i, work_i, iwork_i );
00298
00299 LAPACKE_dge_trans( LAPACK_ROW_MAJOR, n, nrhs, x_r, nrhs+2, x_i, ldx );
00300
00301 failed = compare_dsyrfs( x, x_i, ferr, ferr_i, berr, berr_i, info, info_i,
00302 ldx, nrhs );
00303 if( failed == 0 ) {
00304 printf( "PASSED: row-major middle-level interface to dsyrfs\n" );
00305 } else {
00306 printf( "FAILED: row-major middle-level interface to dsyrfs\n" );
00307 }
00308
00309
00310
00311 for( i = 0; i < lda*n; i++ ) {
00312 a_i[i] = a[i];
00313 }
00314 for( i = 0; i < ldaf*n; i++ ) {
00315 af_i[i] = af[i];
00316 }
00317 for( i = 0; i < n; i++ ) {
00318 ipiv_i[i] = ipiv[i];
00319 }
00320 for( i = 0; i < ldb*nrhs; i++ ) {
00321 b_i[i] = b[i];
00322 }
00323 for( i = 0; i < ldx*nrhs; i++ ) {
00324 x_i[i] = x_save[i];
00325 }
00326 for( i = 0; i < nrhs; i++ ) {
00327 ferr_i[i] = ferr_save[i];
00328 }
00329 for( i = 0; i < nrhs; i++ ) {
00330 berr_i[i] = berr_save[i];
00331 }
00332 for( i = 0; i < 3*n; i++ ) {
00333 work_i[i] = work[i];
00334 }
00335 for( i = 0; i < n; i++ ) {
00336 iwork_i[i] = iwork[i];
00337 }
00338
00339
00340 LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, n, a_i, lda, a_r, n+2 );
00341 LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, n, af_i, ldaf, af_r, n+2 );
00342 LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, nrhs, b_i, ldb, b_r, nrhs+2 );
00343 LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, nrhs, x_i, ldx, x_r, nrhs+2 );
00344 info_i = LAPACKE_dsyrfs( LAPACK_ROW_MAJOR, uplo_i, n_i, nrhs_i, a_r, lda_r,
00345 af_r, ldaf_r, ipiv_i, b_r, ldb_r, x_r, ldx_r,
00346 ferr_i, berr_i );
00347
00348 LAPACKE_dge_trans( LAPACK_ROW_MAJOR, n, nrhs, x_r, nrhs+2, x_i, ldx );
00349
00350 failed = compare_dsyrfs( x, x_i, ferr, ferr_i, berr, berr_i, info, info_i,
00351 ldx, nrhs );
00352 if( failed == 0 ) {
00353 printf( "PASSED: row-major high-level interface to dsyrfs\n" );
00354 } else {
00355 printf( "FAILED: row-major high-level interface to dsyrfs\n" );
00356 }
00357
00358
00359 if( a != NULL ) {
00360 LAPACKE_free( a );
00361 }
00362 if( a_i != NULL ) {
00363 LAPACKE_free( a_i );
00364 }
00365 if( a_r != NULL ) {
00366 LAPACKE_free( a_r );
00367 }
00368 if( af != NULL ) {
00369 LAPACKE_free( af );
00370 }
00371 if( af_i != NULL ) {
00372 LAPACKE_free( af_i );
00373 }
00374 if( af_r != NULL ) {
00375 LAPACKE_free( af_r );
00376 }
00377 if( ipiv != NULL ) {
00378 LAPACKE_free( ipiv );
00379 }
00380 if( ipiv_i != NULL ) {
00381 LAPACKE_free( ipiv_i );
00382 }
00383 if( b != NULL ) {
00384 LAPACKE_free( b );
00385 }
00386 if( b_i != NULL ) {
00387 LAPACKE_free( b_i );
00388 }
00389 if( b_r != NULL ) {
00390 LAPACKE_free( b_r );
00391 }
00392 if( x != NULL ) {
00393 LAPACKE_free( x );
00394 }
00395 if( x_i != NULL ) {
00396 LAPACKE_free( x_i );
00397 }
00398 if( x_r != NULL ) {
00399 LAPACKE_free( x_r );
00400 }
00401 if( x_save != NULL ) {
00402 LAPACKE_free( x_save );
00403 }
00404 if( ferr != NULL ) {
00405 LAPACKE_free( ferr );
00406 }
00407 if( ferr_i != NULL ) {
00408 LAPACKE_free( ferr_i );
00409 }
00410 if( ferr_save != NULL ) {
00411 LAPACKE_free( ferr_save );
00412 }
00413 if( berr != NULL ) {
00414 LAPACKE_free( berr );
00415 }
00416 if( berr_i != NULL ) {
00417 LAPACKE_free( berr_i );
00418 }
00419 if( berr_save != NULL ) {
00420 LAPACKE_free( berr_save );
00421 }
00422 if( work != NULL ) {
00423 LAPACKE_free( work );
00424 }
00425 if( work_i != NULL ) {
00426 LAPACKE_free( work_i );
00427 }
00428 if( iwork != NULL ) {
00429 LAPACKE_free( iwork );
00430 }
00431 if( iwork_i != NULL ) {
00432 LAPACKE_free( iwork_i );
00433 }
00434
00435 return 0;
00436 }
00437
00438
00439 static void init_scalars_dsyrfs( char *uplo, lapack_int *n, lapack_int *nrhs,
00440 lapack_int *lda, lapack_int *ldaf,
00441 lapack_int *ldb, lapack_int *ldx )
00442 {
00443 *uplo = 'L';
00444 *n = 4;
00445 *nrhs = 2;
00446 *lda = 8;
00447 *ldaf = 8;
00448 *ldb = 8;
00449 *ldx = 8;
00450
00451 return;
00452 }
00453
00454
00455 static void init_a( lapack_int size, double *a ) {
00456 lapack_int i;
00457 for( i = 0; i < size; i++ ) {
00458 a[i] = 0;
00459 }
00460 a[0] = 2.06999999999999980e+000;
00461 a[8] = 0.00000000000000000e+000;
00462 a[16] = 0.00000000000000000e+000;
00463 a[24] = 0.00000000000000000e+000;
00464 a[1] = 3.87000000000000010e+000;
00465 a[9] = -2.09999999999999990e-001;
00466 a[17] = 0.00000000000000000e+000;
00467 a[25] = 0.00000000000000000e+000;
00468 a[2] = 4.20000000000000020e+000;
00469 a[10] = 1.87000000000000010e+000;
00470 a[18] = 1.14999999999999990e+000;
00471 a[26] = 0.00000000000000000e+000;
00472 a[3] = -1.14999999999999990e+000;
00473 a[11] = 6.30000000000000000e-001;
00474 a[19] = 2.06000000000000010e+000;
00475 a[27] = -1.81000000000000010e+000;
00476 }
00477 static void init_af( lapack_int size, double *af ) {
00478 lapack_int i;
00479 for( i = 0; i < size; i++ ) {
00480 af[i] = 0;
00481 }
00482 af[0] = 2.06999999999999980e+000;
00483 af[8] = 0.00000000000000000e+000;
00484 af[16] = 0.00000000000000000e+000;
00485 af[24] = 0.00000000000000000e+000;
00486 af[1] = 4.20000000000000020e+000;
00487 af[9] = 1.14999999999999990e+000;
00488 af[17] = 0.00000000000000000e+000;
00489 af[25] = 0.00000000000000000e+000;
00490 af[2] = 2.23041384055834070e-001;
00491 af[10] = 8.11501032143910230e-001;
00492 af[18] = -2.59067708640519000e+000;
00493 af[26] = 0.00000000000000000e+000;
00494 af[3] = 6.53658376748910360e-001;
00495 af[11] = -5.95969723778629450e-001;
00496 af[19] = 3.03084679550618070e-001;
00497 af[27] = 4.07385198134887610e-001;
00498 }
00499 static void init_ipiv( lapack_int size, lapack_int *ipiv ) {
00500 lapack_int i;
00501 for( i = 0; i < size; i++ ) {
00502 ipiv[i] = 0;
00503 }
00504 ipiv[0] = -3;
00505 ipiv[1] = -3;
00506 ipiv[2] = 3;
00507 ipiv[3] = 4;
00508 }
00509 static void init_b( lapack_int size, double *b ) {
00510 lapack_int i;
00511 for( i = 0; i < size; i++ ) {
00512 b[i] = 0;
00513 }
00514 b[0] = -9.50000000000000000e+000;
00515 b[8] = 2.78500000000000010e+001;
00516 b[1] = -8.38000000000000080e+000;
00517 b[9] = 9.90000000000000040e+000;
00518 b[2] = -6.07000000000000030e+000;
00519 b[10] = 1.92500000000000000e+001;
00520 b[3] = -9.59999999999999960e-001;
00521 b[11] = 3.93000000000000020e+000;
00522 }
00523 static void init_x( lapack_int size, double *x ) {
00524 lapack_int i;
00525 for( i = 0; i < size; i++ ) {
00526 x[i] = 0;
00527 }
00528 x[0] = -4.00000000000000710e+000;
00529 x[8] = 1.00000000000000440e+000;
00530 x[1] = -1.00000000000000270e+000;
00531 x[9] = 4.00000000000000180e+000;
00532 x[2] = 2.00000000000000800e+000;
00533 x[10] = 2.99999999999999470e+000;
00534 x[3] = 5.00000000000001240e+000;
00535 x[11] = 1.99999999999999200e+000;
00536 }
00537 static void init_ferr( lapack_int size, double *ferr ) {
00538 lapack_int i;
00539 for( i = 0; i < size; i++ ) {
00540 ferr[i] = 0;
00541 }
00542 }
00543 static void init_berr( lapack_int size, double *berr ) {
00544 lapack_int i;
00545 for( i = 0; i < size; i++ ) {
00546 berr[i] = 0;
00547 }
00548 }
00549 static void init_work( lapack_int size, double *work ) {
00550 lapack_int i;
00551 for( i = 0; i < size; i++ ) {
00552 work[i] = 0;
00553 }
00554 }
00555 static void init_iwork( lapack_int size, lapack_int *iwork ) {
00556 lapack_int i;
00557 for( i = 0; i < size; i++ ) {
00558 iwork[i] = 0;
00559 }
00560 }
00561
00562
00563
00564 static int compare_dsyrfs( double *x, double *x_i, double *ferr, double *ferr_i,
00565 double *berr, double *berr_i, lapack_int info,
00566 lapack_int info_i, lapack_int ldx, lapack_int nrhs )
00567 {
00568 lapack_int i;
00569 int failed = 0;
00570 for( i = 0; i < ldx*nrhs; i++ ) {
00571 failed += compare_doubles(x[i],x_i[i]);
00572 }
00573 for( i = 0; i < nrhs; i++ ) {
00574 failed += compare_doubles(ferr[i],ferr_i[i]);
00575 }
00576 for( i = 0; i < nrhs; i++ ) {
00577 failed += compare_doubles(berr[i],berr_i[i]);
00578 }
00579 failed += (info == info_i) ? 0 : 1;
00580 if( info != 0 || info_i != 0 ) {
00581 printf( "info=%d, info_i=%d\n",(int)info,(int)info_i );
00582 }
00583
00584 return failed;
00585 }