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_dpbrfs( char *uplo, lapack_int *n, lapack_int *kd,
00055 lapack_int *nrhs, lapack_int *ldab,
00056 lapack_int *ldafb, lapack_int *ldb,
00057 lapack_int *ldx );
00058 static void init_ab( lapack_int size, double *ab );
00059 static void init_afb( lapack_int size, double *afb );
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_dpbrfs( 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 kd, kd_i;
00076 lapack_int nrhs, nrhs_i;
00077 lapack_int ldab, ldab_i;
00078 lapack_int ldab_r;
00079 lapack_int ldafb, ldafb_i;
00080 lapack_int ldafb_r;
00081 lapack_int ldb, ldb_i;
00082 lapack_int ldb_r;
00083 lapack_int ldx, ldx_i;
00084 lapack_int ldx_r;
00085 lapack_int info, info_i;
00086 lapack_int i;
00087 int failed;
00088
00089
00090 double *ab = NULL, *ab_i = NULL;
00091 double *afb = NULL, *afb_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 *ab_r = NULL;
00102 double *afb_r = NULL;
00103 double *b_r = NULL;
00104 double *x_r = NULL;
00105
00106
00107 init_scalars_dpbrfs( &uplo, &n, &kd, &nrhs, &ldab, &ldafb, &ldb, &ldx );
00108 ldab_r = n+2;
00109 ldafb_r = n+2;
00110 ldb_r = nrhs+2;
00111 ldx_r = nrhs+2;
00112 uplo_i = uplo;
00113 n_i = n;
00114 kd_i = kd;
00115 nrhs_i = nrhs;
00116 ldab_i = ldab;
00117 ldafb_i = ldafb;
00118 ldb_i = ldb;
00119 ldx_i = ldx;
00120
00121
00122 ab = (double *)LAPACKE_malloc( ldab*n * sizeof(double) );
00123 afb = (double *)LAPACKE_malloc( ldafb*n * sizeof(double) );
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 ab_i = (double *)LAPACKE_malloc( ldab*n * sizeof(double) );
00133 afb_i = (double *)LAPACKE_malloc( ldafb*n * sizeof(double) );
00134 b_i = (double *)LAPACKE_malloc( ldb*nrhs * sizeof(double) );
00135 x_i = (double *)LAPACKE_malloc( ldx*nrhs * sizeof(double) );
00136 ferr_i = (double *)LAPACKE_malloc( nrhs * sizeof(double) );
00137 berr_i = (double *)LAPACKE_malloc( nrhs * sizeof(double) );
00138 work_i = (double *)LAPACKE_malloc( 3*n * sizeof(double) );
00139 iwork_i = (lapack_int *)LAPACKE_malloc( n * sizeof(lapack_int) );
00140
00141
00142 x_save = (double *)LAPACKE_malloc( ldx*nrhs * sizeof(double) );
00143 ferr_save = (double *)LAPACKE_malloc( nrhs * sizeof(double) );
00144 berr_save = (double *)LAPACKE_malloc( nrhs * sizeof(double) );
00145
00146
00147 ab_r = (double *)LAPACKE_malloc( (kd+1)*(n+2) * sizeof(double) );
00148 afb_r = (double *)LAPACKE_malloc( (kd+1)*(n+2) * sizeof(double) );
00149 b_r = (double *)LAPACKE_malloc( n*(nrhs+2) * sizeof(double) );
00150 x_r = (double *)LAPACKE_malloc( n*(nrhs+2) * sizeof(double) );
00151
00152
00153 init_ab( ldab*n, ab );
00154 init_afb( ldafb*n, afb );
00155 init_b( ldb*nrhs, b );
00156 init_x( ldx*nrhs, x );
00157 init_ferr( nrhs, ferr );
00158 init_berr( nrhs, berr );
00159 init_work( 3*n, work );
00160 init_iwork( n, iwork );
00161
00162
00163 for( i = 0; i < ldx*nrhs; i++ ) {
00164 x_save[i] = x[i];
00165 }
00166 for( i = 0; i < nrhs; i++ ) {
00167 ferr_save[i] = ferr[i];
00168 }
00169 for( i = 0; i < nrhs; i++ ) {
00170 berr_save[i] = berr[i];
00171 }
00172
00173
00174 dpbrfs_( &uplo, &n, &kd, &nrhs, ab, &ldab, afb, &ldafb, b, &ldb, x, &ldx,
00175 ferr, berr, work, iwork, &info );
00176
00177
00178
00179 for( i = 0; i < ldab*n; i++ ) {
00180 ab_i[i] = ab[i];
00181 }
00182 for( i = 0; i < ldafb*n; i++ ) {
00183 afb_i[i] = afb[i];
00184 }
00185 for( i = 0; i < ldb*nrhs; i++ ) {
00186 b_i[i] = b[i];
00187 }
00188 for( i = 0; i < ldx*nrhs; i++ ) {
00189 x_i[i] = x_save[i];
00190 }
00191 for( i = 0; i < nrhs; i++ ) {
00192 ferr_i[i] = ferr_save[i];
00193 }
00194 for( i = 0; i < nrhs; i++ ) {
00195 berr_i[i] = berr_save[i];
00196 }
00197 for( i = 0; i < 3*n; i++ ) {
00198 work_i[i] = work[i];
00199 }
00200 for( i = 0; i < n; i++ ) {
00201 iwork_i[i] = iwork[i];
00202 }
00203 info_i = LAPACKE_dpbrfs_work( LAPACK_COL_MAJOR, uplo_i, n_i, kd_i, nrhs_i,
00204 ab_i, ldab_i, afb_i, ldafb_i, b_i, ldb_i, x_i,
00205 ldx_i, ferr_i, berr_i, work_i, iwork_i );
00206
00207 failed = compare_dpbrfs( x, x_i, ferr, ferr_i, berr, berr_i, info, info_i,
00208 ldx, nrhs );
00209 if( failed == 0 ) {
00210 printf( "PASSED: column-major middle-level interface to dpbrfs\n" );
00211 } else {
00212 printf( "FAILED: column-major middle-level interface to dpbrfs\n" );
00213 }
00214
00215
00216
00217 for( i = 0; i < ldab*n; i++ ) {
00218 ab_i[i] = ab[i];
00219 }
00220 for( i = 0; i < ldafb*n; i++ ) {
00221 afb_i[i] = afb[i];
00222 }
00223 for( i = 0; i < ldb*nrhs; i++ ) {
00224 b_i[i] = b[i];
00225 }
00226 for( i = 0; i < ldx*nrhs; i++ ) {
00227 x_i[i] = x_save[i];
00228 }
00229 for( i = 0; i < nrhs; i++ ) {
00230 ferr_i[i] = ferr_save[i];
00231 }
00232 for( i = 0; i < nrhs; i++ ) {
00233 berr_i[i] = berr_save[i];
00234 }
00235 for( i = 0; i < 3*n; i++ ) {
00236 work_i[i] = work[i];
00237 }
00238 for( i = 0; i < n; i++ ) {
00239 iwork_i[i] = iwork[i];
00240 }
00241 info_i = LAPACKE_dpbrfs( LAPACK_COL_MAJOR, uplo_i, n_i, kd_i, nrhs_i, ab_i,
00242 ldab_i, afb_i, ldafb_i, b_i, ldb_i, x_i, ldx_i,
00243 ferr_i, berr_i );
00244
00245 failed = compare_dpbrfs( x, x_i, ferr, ferr_i, berr, berr_i, info, info_i,
00246 ldx, nrhs );
00247 if( failed == 0 ) {
00248 printf( "PASSED: column-major high-level interface to dpbrfs\n" );
00249 } else {
00250 printf( "FAILED: column-major high-level interface to dpbrfs\n" );
00251 }
00252
00253
00254
00255 for( i = 0; i < ldab*n; i++ ) {
00256 ab_i[i] = ab[i];
00257 }
00258 for( i = 0; i < ldafb*n; i++ ) {
00259 afb_i[i] = afb[i];
00260 }
00261 for( i = 0; i < ldb*nrhs; i++ ) {
00262 b_i[i] = b[i];
00263 }
00264 for( i = 0; i < ldx*nrhs; i++ ) {
00265 x_i[i] = x_save[i];
00266 }
00267 for( i = 0; i < nrhs; i++ ) {
00268 ferr_i[i] = ferr_save[i];
00269 }
00270 for( i = 0; i < nrhs; i++ ) {
00271 berr_i[i] = berr_save[i];
00272 }
00273 for( i = 0; i < 3*n; i++ ) {
00274 work_i[i] = work[i];
00275 }
00276 for( i = 0; i < n; i++ ) {
00277 iwork_i[i] = iwork[i];
00278 }
00279
00280 LAPACKE_dge_trans( LAPACK_COL_MAJOR, kd+1, n, ab_i, ldab, ab_r, n+2 );
00281 LAPACKE_dge_trans( LAPACK_COL_MAJOR, kd+1, n, afb_i, ldafb, afb_r, n+2 );
00282 LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, nrhs, b_i, ldb, b_r, nrhs+2 );
00283 LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, nrhs, x_i, ldx, x_r, nrhs+2 );
00284 info_i = LAPACKE_dpbrfs_work( LAPACK_ROW_MAJOR, uplo_i, n_i, kd_i, nrhs_i,
00285 ab_r, ldab_r, afb_r, ldafb_r, b_r, ldb_r, x_r,
00286 ldx_r, ferr_i, berr_i, work_i, iwork_i );
00287
00288 LAPACKE_dge_trans( LAPACK_ROW_MAJOR, n, nrhs, x_r, nrhs+2, x_i, ldx );
00289
00290 failed = compare_dpbrfs( x, x_i, ferr, ferr_i, berr, berr_i, info, info_i,
00291 ldx, nrhs );
00292 if( failed == 0 ) {
00293 printf( "PASSED: row-major middle-level interface to dpbrfs\n" );
00294 } else {
00295 printf( "FAILED: row-major middle-level interface to dpbrfs\n" );
00296 }
00297
00298
00299
00300 for( i = 0; i < ldab*n; i++ ) {
00301 ab_i[i] = ab[i];
00302 }
00303 for( i = 0; i < ldafb*n; i++ ) {
00304 afb_i[i] = afb[i];
00305 }
00306 for( i = 0; i < ldb*nrhs; i++ ) {
00307 b_i[i] = b[i];
00308 }
00309 for( i = 0; i < ldx*nrhs; i++ ) {
00310 x_i[i] = x_save[i];
00311 }
00312 for( i = 0; i < nrhs; i++ ) {
00313 ferr_i[i] = ferr_save[i];
00314 }
00315 for( i = 0; i < nrhs; i++ ) {
00316 berr_i[i] = berr_save[i];
00317 }
00318 for( i = 0; i < 3*n; i++ ) {
00319 work_i[i] = work[i];
00320 }
00321 for( i = 0; i < n; i++ ) {
00322 iwork_i[i] = iwork[i];
00323 }
00324
00325
00326 LAPACKE_dge_trans( LAPACK_COL_MAJOR, kd+1, n, ab_i, ldab, ab_r, n+2 );
00327 LAPACKE_dge_trans( LAPACK_COL_MAJOR, kd+1, n, afb_i, ldafb, afb_r, n+2 );
00328 LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, nrhs, b_i, ldb, b_r, nrhs+2 );
00329 LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, nrhs, x_i, ldx, x_r, nrhs+2 );
00330 info_i = LAPACKE_dpbrfs( LAPACK_ROW_MAJOR, uplo_i, n_i, kd_i, nrhs_i, ab_r,
00331 ldab_r, afb_r, ldafb_r, b_r, ldb_r, x_r, ldx_r,
00332 ferr_i, berr_i );
00333
00334 LAPACKE_dge_trans( LAPACK_ROW_MAJOR, n, nrhs, x_r, nrhs+2, x_i, ldx );
00335
00336 failed = compare_dpbrfs( x, x_i, ferr, ferr_i, berr, berr_i, info, info_i,
00337 ldx, nrhs );
00338 if( failed == 0 ) {
00339 printf( "PASSED: row-major high-level interface to dpbrfs\n" );
00340 } else {
00341 printf( "FAILED: row-major high-level interface to dpbrfs\n" );
00342 }
00343
00344
00345 if( ab != NULL ) {
00346 LAPACKE_free( ab );
00347 }
00348 if( ab_i != NULL ) {
00349 LAPACKE_free( ab_i );
00350 }
00351 if( ab_r != NULL ) {
00352 LAPACKE_free( ab_r );
00353 }
00354 if( afb != NULL ) {
00355 LAPACKE_free( afb );
00356 }
00357 if( afb_i != NULL ) {
00358 LAPACKE_free( afb_i );
00359 }
00360 if( afb_r != NULL ) {
00361 LAPACKE_free( afb_r );
00362 }
00363 if( b != NULL ) {
00364 LAPACKE_free( b );
00365 }
00366 if( b_i != NULL ) {
00367 LAPACKE_free( b_i );
00368 }
00369 if( b_r != NULL ) {
00370 LAPACKE_free( b_r );
00371 }
00372 if( x != NULL ) {
00373 LAPACKE_free( x );
00374 }
00375 if( x_i != NULL ) {
00376 LAPACKE_free( x_i );
00377 }
00378 if( x_r != NULL ) {
00379 LAPACKE_free( x_r );
00380 }
00381 if( x_save != NULL ) {
00382 LAPACKE_free( x_save );
00383 }
00384 if( ferr != NULL ) {
00385 LAPACKE_free( ferr );
00386 }
00387 if( ferr_i != NULL ) {
00388 LAPACKE_free( ferr_i );
00389 }
00390 if( ferr_save != NULL ) {
00391 LAPACKE_free( ferr_save );
00392 }
00393 if( berr != NULL ) {
00394 LAPACKE_free( berr );
00395 }
00396 if( berr_i != NULL ) {
00397 LAPACKE_free( berr_i );
00398 }
00399 if( berr_save != NULL ) {
00400 LAPACKE_free( berr_save );
00401 }
00402 if( work != NULL ) {
00403 LAPACKE_free( work );
00404 }
00405 if( work_i != NULL ) {
00406 LAPACKE_free( work_i );
00407 }
00408 if( iwork != NULL ) {
00409 LAPACKE_free( iwork );
00410 }
00411 if( iwork_i != NULL ) {
00412 LAPACKE_free( iwork_i );
00413 }
00414
00415 return 0;
00416 }
00417
00418
00419 static void init_scalars_dpbrfs( char *uplo, lapack_int *n, lapack_int *kd,
00420 lapack_int *nrhs, lapack_int *ldab,
00421 lapack_int *ldafb, lapack_int *ldb,
00422 lapack_int *ldx )
00423 {
00424 *uplo = 'L';
00425 *n = 4;
00426 *kd = 1;
00427 *nrhs = 2;
00428 *ldab = 9;
00429 *ldafb = 9;
00430 *ldb = 8;
00431 *ldx = 8;
00432
00433 return;
00434 }
00435
00436
00437 static void init_ab( lapack_int size, double *ab ) {
00438 lapack_int i;
00439 for( i = 0; i < size; i++ ) {
00440 ab[i] = 0;
00441 }
00442 ab[0] = 5.49000000000000020e+000;
00443 ab[9] = 5.62999999999999990e+000;
00444 ab[18] = 2.60000000000000010e+000;
00445 ab[27] = 5.16999999999999990e+000;
00446 ab[1] = 2.68000000000000020e+000;
00447 ab[10] = -2.39000000000000010e+000;
00448 ab[19] = -2.22000000000000020e+000;
00449 ab[28] = 0.00000000000000000e+000;
00450 }
00451 static void init_afb( lapack_int size, double *afb ) {
00452 lapack_int i;
00453 for( i = 0; i < size; i++ ) {
00454 afb[i] = 0;
00455 }
00456 afb[0] = 2.34307490277199640e+000;
00457 afb[9] = 2.07887720150650870e+000;
00458 afb[18] = 1.13061224833700420e+000;
00459 afb[27] = 1.14652471173422900e+000;
00460 afb[1] = 1.14379612740053750e+000;
00461 afb[10] = -1.14965905550747700e+000;
00462 afb[19] = -1.96353790016458380e+000;
00463 afb[28] = 0.00000000000000000e+000;
00464 }
00465 static void init_b( lapack_int size, double *b ) {
00466 lapack_int i;
00467 for( i = 0; i < size; i++ ) {
00468 b[i] = 0;
00469 }
00470 b[0] = 2.20900000000000000e+001;
00471 b[8] = 5.09999999999999960e+000;
00472 b[1] = 9.31000000000000050e+000;
00473 b[9] = 3.08099999999999990e+001;
00474 b[2] = -5.24000000000000020e+000;
00475 b[10] = -2.58200000000000000e+001;
00476 b[3] = 1.18300000000000000e+001;
00477 b[11] = 2.28999999999999990e+001;
00478 }
00479 static void init_x( lapack_int size, double *x ) {
00480 lapack_int i;
00481 for( i = 0; i < size; i++ ) {
00482 x[i] = 0;
00483 }
00484 x[0] = 5.00000000000000000e+000;
00485 x[8] = -2.00000000000000220e+000;
00486 x[1] = -2.00000000000000090e+000;
00487 x[9] = 6.00000000000000360e+000;
00488 x[2] = -3.00000000000000130e+000;
00489 x[10] = -9.99999999999996230e-001;
00490 x[3] = 9.99999999999999440e-001;
00491 x[11] = 4.00000000000000090e+000;
00492 }
00493 static void init_ferr( lapack_int size, double *ferr ) {
00494 lapack_int i;
00495 for( i = 0; i < size; i++ ) {
00496 ferr[i] = 0;
00497 }
00498 }
00499 static void init_berr( lapack_int size, double *berr ) {
00500 lapack_int i;
00501 for( i = 0; i < size; i++ ) {
00502 berr[i] = 0;
00503 }
00504 }
00505 static void init_work( lapack_int size, double *work ) {
00506 lapack_int i;
00507 for( i = 0; i < size; i++ ) {
00508 work[i] = 0;
00509 }
00510 }
00511 static void init_iwork( lapack_int size, lapack_int *iwork ) {
00512 lapack_int i;
00513 for( i = 0; i < size; i++ ) {
00514 iwork[i] = 0;
00515 }
00516 }
00517
00518
00519
00520 static int compare_dpbrfs( double *x, double *x_i, double *ferr, double *ferr_i,
00521 double *berr, double *berr_i, lapack_int info,
00522 lapack_int info_i, lapack_int ldx, lapack_int nrhs )
00523 {
00524 lapack_int i;
00525 int failed = 0;
00526 for( i = 0; i < ldx*nrhs; i++ ) {
00527 failed += compare_doubles(x[i],x_i[i]);
00528 }
00529 for( i = 0; i < nrhs; i++ ) {
00530 failed += compare_doubles(ferr[i],ferr_i[i]);
00531 }
00532 for( i = 0; i < nrhs; i++ ) {
00533 failed += compare_doubles(berr[i],berr_i[i]);
00534 }
00535 failed += (info == info_i) ? 0 : 1;
00536 if( info != 0 || info_i != 0 ) {
00537 printf( "info=%d, info_i=%d\n",(int)info,(int)info_i );
00538 }
00539
00540 return failed;
00541 }