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_cpprfs( char *uplo, lapack_int *n, lapack_int *nrhs,
00055 lapack_int *ldb, lapack_int *ldx );
00056 static void init_ap( lapack_int size, lapack_complex_float *ap );
00057 static void init_afp( lapack_int size, lapack_complex_float *afp );
00058 static void init_b( lapack_int size, lapack_complex_float *b );
00059 static void init_x( lapack_int size, lapack_complex_float *x );
00060 static void init_ferr( lapack_int size, float *ferr );
00061 static void init_berr( lapack_int size, float *berr );
00062 static void init_work( lapack_int size, lapack_complex_float *work );
00063 static void init_rwork( lapack_int size, float *rwork );
00064 static int compare_cpprfs( lapack_complex_float *x, lapack_complex_float *x_i,
00065 float *ferr, float *ferr_i, float *berr,
00066 float *berr_i, lapack_int info, lapack_int info_i,
00067 lapack_int ldx, lapack_int nrhs );
00068
00069 int main(void)
00070 {
00071
00072 char uplo, uplo_i;
00073 lapack_int n, n_i;
00074 lapack_int nrhs, nrhs_i;
00075 lapack_int ldb, ldb_i;
00076 lapack_int ldb_r;
00077 lapack_int ldx, ldx_i;
00078 lapack_int ldx_r;
00079 lapack_int info, info_i;
00080 lapack_int i;
00081 int failed;
00082
00083
00084 lapack_complex_float *ap = NULL, *ap_i = NULL;
00085 lapack_complex_float *afp = NULL, *afp_i = NULL;
00086 lapack_complex_float *b = NULL, *b_i = NULL;
00087 lapack_complex_float *x = NULL, *x_i = NULL;
00088 float *ferr = NULL, *ferr_i = NULL;
00089 float *berr = NULL, *berr_i = NULL;
00090 lapack_complex_float *work = NULL, *work_i = NULL;
00091 float *rwork = NULL, *rwork_i = NULL;
00092 lapack_complex_float *x_save = NULL;
00093 float *ferr_save = NULL;
00094 float *berr_save = NULL;
00095 lapack_complex_float *ap_r = NULL;
00096 lapack_complex_float *afp_r = NULL;
00097 lapack_complex_float *b_r = NULL;
00098 lapack_complex_float *x_r = NULL;
00099
00100
00101 init_scalars_cpprfs( &uplo, &n, &nrhs, &ldb, &ldx );
00102 ldb_r = nrhs+2;
00103 ldx_r = nrhs+2;
00104 uplo_i = uplo;
00105 n_i = n;
00106 nrhs_i = nrhs;
00107 ldb_i = ldb;
00108 ldx_i = ldx;
00109
00110
00111 ap = (lapack_complex_float *)
00112 LAPACKE_malloc( ((n*(n+1)/2)) * sizeof(lapack_complex_float) );
00113 afp = (lapack_complex_float *)
00114 LAPACKE_malloc( ((n*(n+1)/2)) * sizeof(lapack_complex_float) );
00115 b = (lapack_complex_float *)
00116 LAPACKE_malloc( ldb*nrhs * sizeof(lapack_complex_float) );
00117 x = (lapack_complex_float *)
00118 LAPACKE_malloc( ldx*nrhs * sizeof(lapack_complex_float) );
00119 ferr = (float *)LAPACKE_malloc( nrhs * sizeof(float) );
00120 berr = (float *)LAPACKE_malloc( nrhs * sizeof(float) );
00121 work = (lapack_complex_float *)
00122 LAPACKE_malloc( 2*n * sizeof(lapack_complex_float) );
00123 rwork = (float *)LAPACKE_malloc( n * sizeof(float) );
00124
00125
00126 ap_i = (lapack_complex_float *)
00127 LAPACKE_malloc( ((n*(n+1)/2)) * sizeof(lapack_complex_float) );
00128 afp_i = (lapack_complex_float *)
00129 LAPACKE_malloc( ((n*(n+1)/2)) * sizeof(lapack_complex_float) );
00130 b_i = (lapack_complex_float *)
00131 LAPACKE_malloc( ldb*nrhs * sizeof(lapack_complex_float) );
00132 x_i = (lapack_complex_float *)
00133 LAPACKE_malloc( ldx*nrhs * sizeof(lapack_complex_float) );
00134 ferr_i = (float *)LAPACKE_malloc( nrhs * sizeof(float) );
00135 berr_i = (float *)LAPACKE_malloc( nrhs * sizeof(float) );
00136 work_i = (lapack_complex_float *)
00137 LAPACKE_malloc( 2*n * sizeof(lapack_complex_float) );
00138 rwork_i = (float *)LAPACKE_malloc( n * sizeof(float) );
00139
00140
00141 x_save = (lapack_complex_float *)
00142 LAPACKE_malloc( ldx*nrhs * sizeof(lapack_complex_float) );
00143 ferr_save = (float *)LAPACKE_malloc( nrhs * sizeof(float) );
00144 berr_save = (float *)LAPACKE_malloc( nrhs * sizeof(float) );
00145
00146
00147 ap_r = (lapack_complex_float *)
00148 LAPACKE_malloc( n*(n+1)/2 * sizeof(lapack_complex_float) );
00149 afp_r = (lapack_complex_float *)
00150 LAPACKE_malloc( n*(n+1)/2 * sizeof(lapack_complex_float) );
00151 b_r = (lapack_complex_float *)
00152 LAPACKE_malloc( n*(nrhs+2) * sizeof(lapack_complex_float) );
00153 x_r = (lapack_complex_float *)
00154 LAPACKE_malloc( n*(nrhs+2) * sizeof(lapack_complex_float) );
00155
00156
00157 init_ap( (n*(n+1)/2), ap );
00158 init_afp( (n*(n+1)/2), afp );
00159 init_b( ldb*nrhs, b );
00160 init_x( ldx*nrhs, x );
00161 init_ferr( nrhs, ferr );
00162 init_berr( nrhs, berr );
00163 init_work( 2*n, work );
00164 init_rwork( n, rwork );
00165
00166
00167 for( i = 0; i < ldx*nrhs; i++ ) {
00168 x_save[i] = x[i];
00169 }
00170 for( i = 0; i < nrhs; i++ ) {
00171 ferr_save[i] = ferr[i];
00172 }
00173 for( i = 0; i < nrhs; i++ ) {
00174 berr_save[i] = berr[i];
00175 }
00176
00177
00178 cpprfs_( &uplo, &n, &nrhs, ap, afp, b, &ldb, x, &ldx, ferr, berr, work,
00179 rwork, &info );
00180
00181
00182
00183 for( i = 0; i < (n*(n+1)/2); i++ ) {
00184 ap_i[i] = ap[i];
00185 }
00186 for( i = 0; i < (n*(n+1)/2); i++ ) {
00187 afp_i[i] = afp[i];
00188 }
00189 for( i = 0; i < ldb*nrhs; i++ ) {
00190 b_i[i] = b[i];
00191 }
00192 for( i = 0; i < ldx*nrhs; i++ ) {
00193 x_i[i] = x_save[i];
00194 }
00195 for( i = 0; i < nrhs; i++ ) {
00196 ferr_i[i] = ferr_save[i];
00197 }
00198 for( i = 0; i < nrhs; i++ ) {
00199 berr_i[i] = berr_save[i];
00200 }
00201 for( i = 0; i < 2*n; i++ ) {
00202 work_i[i] = work[i];
00203 }
00204 for( i = 0; i < n; i++ ) {
00205 rwork_i[i] = rwork[i];
00206 }
00207 info_i = LAPACKE_cpprfs_work( LAPACK_COL_MAJOR, uplo_i, n_i, nrhs_i, ap_i,
00208 afp_i, b_i, ldb_i, x_i, ldx_i, ferr_i, berr_i,
00209 work_i, rwork_i );
00210
00211 failed = compare_cpprfs( x, x_i, ferr, ferr_i, berr, berr_i, info, info_i,
00212 ldx, nrhs );
00213 if( failed == 0 ) {
00214 printf( "PASSED: column-major middle-level interface to cpprfs\n" );
00215 } else {
00216 printf( "FAILED: column-major middle-level interface to cpprfs\n" );
00217 }
00218
00219
00220
00221 for( i = 0; i < (n*(n+1)/2); i++ ) {
00222 ap_i[i] = ap[i];
00223 }
00224 for( i = 0; i < (n*(n+1)/2); i++ ) {
00225 afp_i[i] = afp[i];
00226 }
00227 for( i = 0; i < ldb*nrhs; i++ ) {
00228 b_i[i] = b[i];
00229 }
00230 for( i = 0; i < ldx*nrhs; i++ ) {
00231 x_i[i] = x_save[i];
00232 }
00233 for( i = 0; i < nrhs; i++ ) {
00234 ferr_i[i] = ferr_save[i];
00235 }
00236 for( i = 0; i < nrhs; i++ ) {
00237 berr_i[i] = berr_save[i];
00238 }
00239 for( i = 0; i < 2*n; i++ ) {
00240 work_i[i] = work[i];
00241 }
00242 for( i = 0; i < n; i++ ) {
00243 rwork_i[i] = rwork[i];
00244 }
00245 info_i = LAPACKE_cpprfs( LAPACK_COL_MAJOR, uplo_i, n_i, nrhs_i, ap_i, afp_i,
00246 b_i, ldb_i, x_i, ldx_i, ferr_i, berr_i );
00247
00248 failed = compare_cpprfs( x, x_i, ferr, ferr_i, berr, berr_i, info, info_i,
00249 ldx, nrhs );
00250 if( failed == 0 ) {
00251 printf( "PASSED: column-major high-level interface to cpprfs\n" );
00252 } else {
00253 printf( "FAILED: column-major high-level interface to cpprfs\n" );
00254 }
00255
00256
00257
00258 for( i = 0; i < (n*(n+1)/2); i++ ) {
00259 ap_i[i] = ap[i];
00260 }
00261 for( i = 0; i < (n*(n+1)/2); i++ ) {
00262 afp_i[i] = afp[i];
00263 }
00264 for( i = 0; i < ldb*nrhs; i++ ) {
00265 b_i[i] = b[i];
00266 }
00267 for( i = 0; i < ldx*nrhs; i++ ) {
00268 x_i[i] = x_save[i];
00269 }
00270 for( i = 0; i < nrhs; i++ ) {
00271 ferr_i[i] = ferr_save[i];
00272 }
00273 for( i = 0; i < nrhs; i++ ) {
00274 berr_i[i] = berr_save[i];
00275 }
00276 for( i = 0; i < 2*n; i++ ) {
00277 work_i[i] = work[i];
00278 }
00279 for( i = 0; i < n; i++ ) {
00280 rwork_i[i] = rwork[i];
00281 }
00282
00283 LAPACKE_cpp_trans( LAPACK_COL_MAJOR, uplo, n, ap_i, ap_r );
00284 LAPACKE_cpp_trans( LAPACK_COL_MAJOR, uplo, n, afp_i, afp_r );
00285 LAPACKE_cge_trans( LAPACK_COL_MAJOR, n, nrhs, b_i, ldb, b_r, nrhs+2 );
00286 LAPACKE_cge_trans( LAPACK_COL_MAJOR, n, nrhs, x_i, ldx, x_r, nrhs+2 );
00287 info_i = LAPACKE_cpprfs_work( LAPACK_ROW_MAJOR, uplo_i, n_i, nrhs_i, ap_r,
00288 afp_r, b_r, ldb_r, x_r, ldx_r, ferr_i, berr_i,
00289 work_i, rwork_i );
00290
00291 LAPACKE_cge_trans( LAPACK_ROW_MAJOR, n, nrhs, x_r, nrhs+2, x_i, ldx );
00292
00293 failed = compare_cpprfs( x, x_i, ferr, ferr_i, berr, berr_i, info, info_i,
00294 ldx, nrhs );
00295 if( failed == 0 ) {
00296 printf( "PASSED: row-major middle-level interface to cpprfs\n" );
00297 } else {
00298 printf( "FAILED: row-major middle-level interface to cpprfs\n" );
00299 }
00300
00301
00302
00303 for( i = 0; i < (n*(n+1)/2); i++ ) {
00304 ap_i[i] = ap[i];
00305 }
00306 for( i = 0; i < (n*(n+1)/2); i++ ) {
00307 afp_i[i] = afp[i];
00308 }
00309 for( i = 0; i < ldb*nrhs; i++ ) {
00310 b_i[i] = b[i];
00311 }
00312 for( i = 0; i < ldx*nrhs; i++ ) {
00313 x_i[i] = x_save[i];
00314 }
00315 for( i = 0; i < nrhs; i++ ) {
00316 ferr_i[i] = ferr_save[i];
00317 }
00318 for( i = 0; i < nrhs; i++ ) {
00319 berr_i[i] = berr_save[i];
00320 }
00321 for( i = 0; i < 2*n; i++ ) {
00322 work_i[i] = work[i];
00323 }
00324 for( i = 0; i < n; i++ ) {
00325 rwork_i[i] = rwork[i];
00326 }
00327
00328
00329 LAPACKE_cpp_trans( LAPACK_COL_MAJOR, uplo, n, ap_i, ap_r );
00330 LAPACKE_cpp_trans( LAPACK_COL_MAJOR, uplo, n, afp_i, afp_r );
00331 LAPACKE_cge_trans( LAPACK_COL_MAJOR, n, nrhs, b_i, ldb, b_r, nrhs+2 );
00332 LAPACKE_cge_trans( LAPACK_COL_MAJOR, n, nrhs, x_i, ldx, x_r, nrhs+2 );
00333 info_i = LAPACKE_cpprfs( LAPACK_ROW_MAJOR, uplo_i, n_i, nrhs_i, ap_r, afp_r,
00334 b_r, ldb_r, x_r, ldx_r, ferr_i, berr_i );
00335
00336 LAPACKE_cge_trans( LAPACK_ROW_MAJOR, n, nrhs, x_r, nrhs+2, x_i, ldx );
00337
00338 failed = compare_cpprfs( x, x_i, ferr, ferr_i, berr, berr_i, info, info_i,
00339 ldx, nrhs );
00340 if( failed == 0 ) {
00341 printf( "PASSED: row-major high-level interface to cpprfs\n" );
00342 } else {
00343 printf( "FAILED: row-major high-level interface to cpprfs\n" );
00344 }
00345
00346
00347 if( ap != NULL ) {
00348 LAPACKE_free( ap );
00349 }
00350 if( ap_i != NULL ) {
00351 LAPACKE_free( ap_i );
00352 }
00353 if( ap_r != NULL ) {
00354 LAPACKE_free( ap_r );
00355 }
00356 if( afp != NULL ) {
00357 LAPACKE_free( afp );
00358 }
00359 if( afp_i != NULL ) {
00360 LAPACKE_free( afp_i );
00361 }
00362 if( afp_r != NULL ) {
00363 LAPACKE_free( afp_r );
00364 }
00365 if( b != NULL ) {
00366 LAPACKE_free( b );
00367 }
00368 if( b_i != NULL ) {
00369 LAPACKE_free( b_i );
00370 }
00371 if( b_r != NULL ) {
00372 LAPACKE_free( b_r );
00373 }
00374 if( x != NULL ) {
00375 LAPACKE_free( x );
00376 }
00377 if( x_i != NULL ) {
00378 LAPACKE_free( x_i );
00379 }
00380 if( x_r != NULL ) {
00381 LAPACKE_free( x_r );
00382 }
00383 if( x_save != NULL ) {
00384 LAPACKE_free( x_save );
00385 }
00386 if( ferr != NULL ) {
00387 LAPACKE_free( ferr );
00388 }
00389 if( ferr_i != NULL ) {
00390 LAPACKE_free( ferr_i );
00391 }
00392 if( ferr_save != NULL ) {
00393 LAPACKE_free( ferr_save );
00394 }
00395 if( berr != NULL ) {
00396 LAPACKE_free( berr );
00397 }
00398 if( berr_i != NULL ) {
00399 LAPACKE_free( berr_i );
00400 }
00401 if( berr_save != NULL ) {
00402 LAPACKE_free( berr_save );
00403 }
00404 if( work != NULL ) {
00405 LAPACKE_free( work );
00406 }
00407 if( work_i != NULL ) {
00408 LAPACKE_free( work_i );
00409 }
00410 if( rwork != NULL ) {
00411 LAPACKE_free( rwork );
00412 }
00413 if( rwork_i != NULL ) {
00414 LAPACKE_free( rwork_i );
00415 }
00416
00417 return 0;
00418 }
00419
00420
00421 static void init_scalars_cpprfs( char *uplo, lapack_int *n, lapack_int *nrhs,
00422 lapack_int *ldb, lapack_int *ldx )
00423 {
00424 *uplo = 'L';
00425 *n = 4;
00426 *nrhs = 2;
00427 *ldb = 8;
00428 *ldx = 8;
00429
00430 return;
00431 }
00432
00433
00434 static void init_ap( lapack_int size, lapack_complex_float *ap ) {
00435 lapack_int i;
00436 for( i = 0; i < size; i++ ) {
00437 ap[i] = lapack_make_complex_float( 0.0f, 0.0f );
00438 }
00439 ap[0] = lapack_make_complex_float( 3.230000019e+000, 0.000000000e+000 );
00440 ap[1] = lapack_make_complex_float( 1.509999990e+000, 1.919999957e+000 );
00441 ap[2] = lapack_make_complex_float( 1.899999976e+000, -8.399999738e-001 );
00442 ap[3] = lapack_make_complex_float( 4.199999869e-001, -2.500000000e+000 );
00443 ap[4] = lapack_make_complex_float( 3.579999924e+000, 0.000000000e+000 );
00444 ap[5] = lapack_make_complex_float( -2.300000042e-001, -1.110000014e+000 );
00445 ap[6] = lapack_make_complex_float( -1.179999948e+000, -1.370000005e+000 );
00446 ap[7] = lapack_make_complex_float( 4.090000153e+000, 0.000000000e+000 );
00447 ap[8] = lapack_make_complex_float( 2.329999924e+000, 1.400000006e-001 );
00448 ap[9] = lapack_make_complex_float( 4.289999962e+000, 0.000000000e+000 );
00449 }
00450 static void init_afp( lapack_int size, lapack_complex_float *afp ) {
00451 lapack_int i;
00452 for( i = 0; i < size; i++ ) {
00453 afp[i] = lapack_make_complex_float( 0.0f, 0.0f );
00454 }
00455 afp[0] = lapack_make_complex_float( 1.797220111e+000, 0.000000000e+000 );
00456 afp[1] = lapack_make_complex_float( 8.401864767e-001, 1.068316579e+000 );
00457 afp[2] = lapack_make_complex_float( 1.057188272e+000, -4.673885107e-001 );
00458 afp[3] = lapack_make_complex_float( 2.336942554e-001, -1.391037226e+000 );
00459 afp[4] = lapack_make_complex_float( 1.316353440e+000, 0.000000000e+000 );
00460 afp[5] = lapack_make_complex_float( -4.701749682e-001, 3.130658865e-001 );
00461 afp[6] = lapack_make_complex_float( 8.335255831e-002, 3.676066920e-002 );
00462 afp[7] = lapack_make_complex_float( 1.560392976e+000, 0.000000000e+000 );
00463 afp[8] = lapack_make_complex_float( 9.359616637e-001, 9.899691939e-001 );
00464 afp[9] = lapack_make_complex_float( 6.603332758e-001, 0.000000000e+000 );
00465 }
00466 static void init_b( lapack_int size, lapack_complex_float *b ) {
00467 lapack_int i;
00468 for( i = 0; i < size; i++ ) {
00469 b[i] = lapack_make_complex_float( 0.0f, 0.0f );
00470 }
00471 b[0] = lapack_make_complex_float( 3.930000067e+000, -6.139999866e+000 );
00472 b[8] = lapack_make_complex_float( 1.480000019e+000, 6.579999924e+000 );
00473 b[1] = lapack_make_complex_float( 6.170000076e+000, 9.420000076e+000 );
00474 b[9] = lapack_make_complex_float( 4.650000095e+000, -4.750000000e+000 );
00475 b[2] = lapack_make_complex_float( -7.170000076e+000, -2.182999992e+001 );
00476 b[10] = lapack_make_complex_float( -4.909999847e+000, 2.289999962e+000 );
00477 b[3] = lapack_make_complex_float( 1.990000010e+000, -1.438000011e+001 );
00478 b[11] = lapack_make_complex_float( 7.639999866e+000, -1.078999996e+001 );
00479 }
00480 static void init_x( lapack_int size, lapack_complex_float *x ) {
00481 lapack_int i;
00482 for( i = 0; i < size; i++ ) {
00483 x[i] = lapack_make_complex_float( 0.0f, 0.0f );
00484 }
00485 x[0] = lapack_make_complex_float( 9.999954104e-001, -9.999989867e-001 );
00486 x[8] = lapack_make_complex_float( -9.999988675e-001, 2.000002623e+000 );
00487 x[1] = lapack_make_complex_float( 1.403683768e-006, 3.000000954e+000 );
00488 x[9] = lapack_make_complex_float( 3.000000238e+000, -4.000000954e+000 );
00489 x[2] = lapack_make_complex_float( -3.999997139e+000, -4.999999523e+000 );
00490 x[10] = lapack_make_complex_float( -2.000000000e+000, 2.999997854e+000 );
00491 x[3] = lapack_make_complex_float( 1.999998331e+000, 9.999975562e-001 );
00492 x[11] = lapack_make_complex_float( 3.999998569e+000, -4.999998569e+000 );
00493 }
00494 static void init_ferr( lapack_int size, float *ferr ) {
00495 lapack_int i;
00496 for( i = 0; i < size; i++ ) {
00497 ferr[i] = 0;
00498 }
00499 }
00500 static void init_berr( lapack_int size, float *berr ) {
00501 lapack_int i;
00502 for( i = 0; i < size; i++ ) {
00503 berr[i] = 0;
00504 }
00505 }
00506 static void init_work( lapack_int size, lapack_complex_float *work ) {
00507 lapack_int i;
00508 for( i = 0; i < size; i++ ) {
00509 work[i] = lapack_make_complex_float( 0.0f, 0.0f );
00510 }
00511 }
00512 static void init_rwork( lapack_int size, float *rwork ) {
00513 lapack_int i;
00514 for( i = 0; i < size; i++ ) {
00515 rwork[i] = 0;
00516 }
00517 }
00518
00519
00520
00521 static int compare_cpprfs( lapack_complex_float *x, lapack_complex_float *x_i,
00522 float *ferr, float *ferr_i, float *berr,
00523 float *berr_i, lapack_int info, lapack_int info_i,
00524 lapack_int ldx, lapack_int nrhs )
00525 {
00526 lapack_int i;
00527 int failed = 0;
00528 for( i = 0; i < ldx*nrhs; i++ ) {
00529 failed += compare_complex_floats(x[i],x_i[i]);
00530 }
00531 for( i = 0; i < nrhs; i++ ) {
00532 failed += compare_floats(ferr[i],ferr_i[i]);
00533 }
00534 for( i = 0; i < nrhs; i++ ) {
00535 failed += compare_floats(berr[i],berr_i[i]);
00536 }
00537 failed += (info == info_i) ? 0 : 1;
00538 if( info != 0 || info_i != 0 ) {
00539 printf( "info=%d, info_i=%d\n",(int)info,(int)info_i );
00540 }
00541
00542 return failed;
00543 }