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_cbdsqr( char *uplo, lapack_int *n, lapack_int *ncvt,
00055 lapack_int *nru, lapack_int *ncc,
00056 lapack_int *ldvt, lapack_int *ldu,
00057 lapack_int *ldc );
00058 static void init_d( lapack_int size, float *d );
00059 static void init_e( lapack_int size, float *e );
00060 static void init_vt( lapack_int size, lapack_complex_float *vt );
00061 static void init_u( lapack_int size, lapack_complex_float *u );
00062 static void init_c( lapack_int size, lapack_complex_float *c );
00063 static void init_work( lapack_int size, float *work );
00064 static int compare_cbdsqr( float *d, float *d_i, float *e, float *e_i,
00065 lapack_complex_float *vt, lapack_complex_float *vt_i,
00066 lapack_complex_float *u, lapack_complex_float *u_i,
00067 lapack_complex_float *c, lapack_complex_float *c_i,
00068 lapack_int info, lapack_int info_i, lapack_int ldc,
00069 lapack_int ldu, lapack_int ldvt, lapack_int n,
00070 lapack_int ncc, lapack_int ncvt, lapack_int nru );
00071
00072 int main(void)
00073 {
00074
00075 char uplo, uplo_i;
00076 lapack_int n, n_i;
00077 lapack_int ncvt, ncvt_i;
00078 lapack_int nru, nru_i;
00079 lapack_int ncc, ncc_i;
00080 lapack_int ldvt, ldvt_i;
00081 lapack_int ldvt_r;
00082 lapack_int ldu, ldu_i;
00083 lapack_int ldu_r;
00084 lapack_int ldc, ldc_i;
00085 lapack_int ldc_r;
00086 lapack_int info, info_i;
00087 lapack_int i;
00088 int failed;
00089
00090
00091 float *d = NULL, *d_i = NULL;
00092 float *e = NULL, *e_i = NULL;
00093 lapack_complex_float *vt = NULL, *vt_i = NULL;
00094 lapack_complex_float *u = NULL, *u_i = NULL;
00095 lapack_complex_float *c = NULL, *c_i = NULL;
00096 float *work = NULL, *work_i = NULL;
00097 float *d_save = NULL;
00098 float *e_save = NULL;
00099 lapack_complex_float *vt_save = NULL;
00100 lapack_complex_float *u_save = NULL;
00101 lapack_complex_float *c_save = NULL;
00102 lapack_complex_float *vt_r = NULL;
00103 lapack_complex_float *u_r = NULL;
00104 lapack_complex_float *c_r = NULL;
00105
00106
00107 init_scalars_cbdsqr( &uplo, &n, &ncvt, &nru, &ncc, &ldvt, &ldu, &ldc );
00108 ldvt_r = ncvt+2;
00109 ldu_r = n+2;
00110 ldc_r = ncc+2;
00111 uplo_i = uplo;
00112 n_i = n;
00113 ncvt_i = ncvt;
00114 nru_i = nru;
00115 ncc_i = ncc;
00116 ldvt_i = ldvt;
00117 ldu_i = ldu;
00118 ldc_i = ldc;
00119
00120
00121 d = (float *)LAPACKE_malloc( n * sizeof(float) );
00122 e = (float *)LAPACKE_malloc( n * sizeof(float) );
00123 vt = (lapack_complex_float *)
00124 LAPACKE_malloc( ldvt*ncvt * sizeof(lapack_complex_float) );
00125 u = (lapack_complex_float *)
00126 LAPACKE_malloc( ldu*n * sizeof(lapack_complex_float) );
00127 c = (lapack_complex_float *)
00128 LAPACKE_malloc( ldc*ncc * sizeof(lapack_complex_float) );
00129 work = (float *)LAPACKE_malloc( 4*n * sizeof(float) );
00130
00131
00132 d_i = (float *)LAPACKE_malloc( n * sizeof(float) );
00133 e_i = (float *)LAPACKE_malloc( n * sizeof(float) );
00134 vt_i = (lapack_complex_float *)
00135 LAPACKE_malloc( ldvt*ncvt * sizeof(lapack_complex_float) );
00136 u_i = (lapack_complex_float *)
00137 LAPACKE_malloc( ldu*n * sizeof(lapack_complex_float) );
00138 c_i = (lapack_complex_float *)
00139 LAPACKE_malloc( ldc*ncc * sizeof(lapack_complex_float) );
00140 work_i = (float *)LAPACKE_malloc( 4*n * sizeof(float) );
00141
00142
00143 d_save = (float *)LAPACKE_malloc( n * sizeof(float) );
00144 e_save = (float *)LAPACKE_malloc( n * sizeof(float) );
00145 vt_save = (lapack_complex_float *)
00146 LAPACKE_malloc( ldvt*ncvt * sizeof(lapack_complex_float) );
00147 u_save = (lapack_complex_float *)
00148 LAPACKE_malloc( ldu*n * sizeof(lapack_complex_float) );
00149 c_save = (lapack_complex_float *)
00150 LAPACKE_malloc( ldc*ncc * sizeof(lapack_complex_float) );
00151
00152
00153 vt_r = (lapack_complex_float *)
00154 LAPACKE_malloc( n*(ncvt+2) * sizeof(lapack_complex_float) );
00155 u_r = (lapack_complex_float *)
00156 LAPACKE_malloc( nru*(n+2) * sizeof(lapack_complex_float) );
00157 c_r = (lapack_complex_float *)
00158 LAPACKE_malloc( n*(ncc+2) * sizeof(lapack_complex_float) );
00159
00160
00161 init_d( n, d );
00162 init_e( n, e );
00163 init_vt( ldvt*ncvt, vt );
00164 init_u( ldu*n, u );
00165 init_c( ldc*ncc, c );
00166 init_work( 4*n, work );
00167
00168
00169 for( i = 0; i < n; i++ ) {
00170 d_save[i] = d[i];
00171 }
00172 for( i = 0; i < n; i++ ) {
00173 e_save[i] = e[i];
00174 }
00175 for( i = 0; i < ldvt*ncvt; i++ ) {
00176 vt_save[i] = vt[i];
00177 }
00178 for( i = 0; i < ldu*n; i++ ) {
00179 u_save[i] = u[i];
00180 }
00181 for( i = 0; i < ldc*ncc; i++ ) {
00182 c_save[i] = c[i];
00183 }
00184
00185
00186 cbdsqr_( &uplo, &n, &ncvt, &nru, &ncc, d, e, vt, &ldvt, u, &ldu, c, &ldc,
00187 work, &info );
00188
00189
00190
00191 for( i = 0; i < n; i++ ) {
00192 d_i[i] = d_save[i];
00193 }
00194 for( i = 0; i < n; i++ ) {
00195 e_i[i] = e_save[i];
00196 }
00197 for( i = 0; i < ldvt*ncvt; i++ ) {
00198 vt_i[i] = vt_save[i];
00199 }
00200 for( i = 0; i < ldu*n; i++ ) {
00201 u_i[i] = u_save[i];
00202 }
00203 for( i = 0; i < ldc*ncc; i++ ) {
00204 c_i[i] = c_save[i];
00205 }
00206 for( i = 0; i < 4*n; i++ ) {
00207 work_i[i] = work[i];
00208 }
00209 info_i = LAPACKE_cbdsqr_work( LAPACK_COL_MAJOR, uplo_i, n_i, ncvt_i, nru_i,
00210 ncc_i, d_i, e_i, vt_i, ldvt_i, u_i, ldu_i,
00211 c_i, ldc_i, work_i );
00212
00213 failed = compare_cbdsqr( d, d_i, e, e_i, vt, vt_i, u, u_i, c, c_i, info,
00214 info_i, ldc, ldu, ldvt, n, ncc, ncvt, nru );
00215 if( failed == 0 ) {
00216 printf( "PASSED: column-major middle-level interface to cbdsqr\n" );
00217 } else {
00218 printf( "FAILED: column-major middle-level interface to cbdsqr\n" );
00219 }
00220
00221
00222
00223 for( i = 0; i < n; i++ ) {
00224 d_i[i] = d_save[i];
00225 }
00226 for( i = 0; i < n; i++ ) {
00227 e_i[i] = e_save[i];
00228 }
00229 for( i = 0; i < ldvt*ncvt; i++ ) {
00230 vt_i[i] = vt_save[i];
00231 }
00232 for( i = 0; i < ldu*n; i++ ) {
00233 u_i[i] = u_save[i];
00234 }
00235 for( i = 0; i < ldc*ncc; i++ ) {
00236 c_i[i] = c_save[i];
00237 }
00238 for( i = 0; i < 4*n; i++ ) {
00239 work_i[i] = work[i];
00240 }
00241 info_i = LAPACKE_cbdsqr( LAPACK_COL_MAJOR, uplo_i, n_i, ncvt_i, nru_i,
00242 ncc_i, d_i, e_i, vt_i, ldvt_i, u_i, ldu_i, c_i,
00243 ldc_i );
00244
00245 failed = compare_cbdsqr( d, d_i, e, e_i, vt, vt_i, u, u_i, c, c_i, info,
00246 info_i, ldc, ldu, ldvt, n, ncc, ncvt, nru );
00247 if( failed == 0 ) {
00248 printf( "PASSED: column-major high-level interface to cbdsqr\n" );
00249 } else {
00250 printf( "FAILED: column-major high-level interface to cbdsqr\n" );
00251 }
00252
00253
00254
00255 for( i = 0; i < n; i++ ) {
00256 d_i[i] = d_save[i];
00257 }
00258 for( i = 0; i < n; i++ ) {
00259 e_i[i] = e_save[i];
00260 }
00261 for( i = 0; i < ldvt*ncvt; i++ ) {
00262 vt_i[i] = vt_save[i];
00263 }
00264 for( i = 0; i < ldu*n; i++ ) {
00265 u_i[i] = u_save[i];
00266 }
00267 for( i = 0; i < ldc*ncc; i++ ) {
00268 c_i[i] = c_save[i];
00269 }
00270 for( i = 0; i < 4*n; i++ ) {
00271 work_i[i] = work[i];
00272 }
00273
00274 if( ncvt != 0 ) {
00275 LAPACKE_cge_trans( LAPACK_COL_MAJOR, n, ncvt, vt_i, ldvt, vt_r,
00276 ncvt+2 );
00277 }
00278 if( nru != 0 ) {
00279 LAPACKE_cge_trans( LAPACK_COL_MAJOR, nru, n, u_i, ldu, u_r, n+2 );
00280 }
00281 if( ncc != 0 ) {
00282 LAPACKE_cge_trans( LAPACK_COL_MAJOR, n, ncc, c_i, ldc, c_r, ncc+2 );
00283 }
00284 info_i = LAPACKE_cbdsqr_work( LAPACK_ROW_MAJOR, uplo_i, n_i, ncvt_i, nru_i,
00285 ncc_i, d_i, e_i, vt_r, ldvt_r, u_r, ldu_r,
00286 c_r, ldc_r, work_i );
00287
00288 if( ncvt != 0 ) {
00289 LAPACKE_cge_trans( LAPACK_ROW_MAJOR, n, ncvt, vt_r, ncvt+2, vt_i,
00290 ldvt );
00291 }
00292 if( nru != 0 ) {
00293 LAPACKE_cge_trans( LAPACK_ROW_MAJOR, nru, n, u_r, n+2, u_i, ldu );
00294 }
00295 if( ncc != 0 ) {
00296 LAPACKE_cge_trans( LAPACK_ROW_MAJOR, n, ncc, c_r, ncc+2, c_i, ldc );
00297 }
00298
00299 failed = compare_cbdsqr( d, d_i, e, e_i, vt, vt_i, u, u_i, c, c_i, info,
00300 info_i, ldc, ldu, ldvt, n, ncc, ncvt, nru );
00301 if( failed == 0 ) {
00302 printf( "PASSED: row-major middle-level interface to cbdsqr\n" );
00303 } else {
00304 printf( "FAILED: row-major middle-level interface to cbdsqr\n" );
00305 }
00306
00307
00308
00309 for( i = 0; i < n; i++ ) {
00310 d_i[i] = d_save[i];
00311 }
00312 for( i = 0; i < n; i++ ) {
00313 e_i[i] = e_save[i];
00314 }
00315 for( i = 0; i < ldvt*ncvt; i++ ) {
00316 vt_i[i] = vt_save[i];
00317 }
00318 for( i = 0; i < ldu*n; i++ ) {
00319 u_i[i] = u_save[i];
00320 }
00321 for( i = 0; i < ldc*ncc; i++ ) {
00322 c_i[i] = c_save[i];
00323 }
00324 for( i = 0; i < 4*n; i++ ) {
00325 work_i[i] = work[i];
00326 }
00327
00328
00329 if( ncvt != 0 ) {
00330 LAPACKE_cge_trans( LAPACK_COL_MAJOR, n, ncvt, vt_i, ldvt, vt_r,
00331 ncvt+2 );
00332 }
00333 if( nru != 0 ) {
00334 LAPACKE_cge_trans( LAPACK_COL_MAJOR, nru, n, u_i, ldu, u_r, n+2 );
00335 }
00336 if( ncc != 0 ) {
00337 LAPACKE_cge_trans( LAPACK_COL_MAJOR, n, ncc, c_i, ldc, c_r, ncc+2 );
00338 }
00339 info_i = LAPACKE_cbdsqr( LAPACK_ROW_MAJOR, uplo_i, n_i, ncvt_i, nru_i,
00340 ncc_i, d_i, e_i, vt_r, ldvt_r, u_r, ldu_r, c_r,
00341 ldc_r );
00342
00343 if( ncvt != 0 ) {
00344 LAPACKE_cge_trans( LAPACK_ROW_MAJOR, n, ncvt, vt_r, ncvt+2, vt_i,
00345 ldvt );
00346 }
00347 if( nru != 0 ) {
00348 LAPACKE_cge_trans( LAPACK_ROW_MAJOR, nru, n, u_r, n+2, u_i, ldu );
00349 }
00350 if( ncc != 0 ) {
00351 LAPACKE_cge_trans( LAPACK_ROW_MAJOR, n, ncc, c_r, ncc+2, c_i, ldc );
00352 }
00353
00354 failed = compare_cbdsqr( d, d_i, e, e_i, vt, vt_i, u, u_i, c, c_i, info,
00355 info_i, ldc, ldu, ldvt, n, ncc, ncvt, nru );
00356 if( failed == 0 ) {
00357 printf( "PASSED: row-major high-level interface to cbdsqr\n" );
00358 } else {
00359 printf( "FAILED: row-major high-level interface to cbdsqr\n" );
00360 }
00361
00362
00363 if( d != NULL ) {
00364 LAPACKE_free( d );
00365 }
00366 if( d_i != NULL ) {
00367 LAPACKE_free( d_i );
00368 }
00369 if( d_save != NULL ) {
00370 LAPACKE_free( d_save );
00371 }
00372 if( e != NULL ) {
00373 LAPACKE_free( e );
00374 }
00375 if( e_i != NULL ) {
00376 LAPACKE_free( e_i );
00377 }
00378 if( e_save != NULL ) {
00379 LAPACKE_free( e_save );
00380 }
00381 if( vt != NULL ) {
00382 LAPACKE_free( vt );
00383 }
00384 if( vt_i != NULL ) {
00385 LAPACKE_free( vt_i );
00386 }
00387 if( vt_r != NULL ) {
00388 LAPACKE_free( vt_r );
00389 }
00390 if( vt_save != NULL ) {
00391 LAPACKE_free( vt_save );
00392 }
00393 if( u != NULL ) {
00394 LAPACKE_free( u );
00395 }
00396 if( u_i != NULL ) {
00397 LAPACKE_free( u_i );
00398 }
00399 if( u_r != NULL ) {
00400 LAPACKE_free( u_r );
00401 }
00402 if( u_save != NULL ) {
00403 LAPACKE_free( u_save );
00404 }
00405 if( c != NULL ) {
00406 LAPACKE_free( c );
00407 }
00408 if( c_i != NULL ) {
00409 LAPACKE_free( c_i );
00410 }
00411 if( c_r != NULL ) {
00412 LAPACKE_free( c_r );
00413 }
00414 if( c_save != NULL ) {
00415 LAPACKE_free( c_save );
00416 }
00417 if( work != NULL ) {
00418 LAPACKE_free( work );
00419 }
00420 if( work_i != NULL ) {
00421 LAPACKE_free( work_i );
00422 }
00423
00424 return 0;
00425 }
00426
00427
00428 static void init_scalars_cbdsqr( char *uplo, lapack_int *n, lapack_int *ncvt,
00429 lapack_int *nru, lapack_int *ncc,
00430 lapack_int *ldvt, lapack_int *ldu,
00431 lapack_int *ldc )
00432 {
00433 *uplo = 'U';
00434 *n = 4;
00435 *ncvt = 4;
00436 *nru = 6;
00437 *ncc = 0;
00438 *ldvt = 8;
00439 *ldu = 8;
00440 *ldc = 8;
00441
00442 return;
00443 }
00444
00445
00446 static void init_d( lapack_int size, float *d ) {
00447 lapack_int i;
00448 for( i = 0; i < size; i++ ) {
00449 d[i] = 0;
00450 }
00451 d[0] = -3.087005138e+000;
00452 d[1] = 2.066039324e+000;
00453 d[2] = 1.873128653e+000;
00454 d[3] = 2.002182961e+000;
00455 }
00456 static void init_e( lapack_int size, float *e ) {
00457 lapack_int i;
00458 for( i = 0; i < size; i++ ) {
00459 e[i] = 0;
00460 }
00461 e[0] = 2.112571001e+000;
00462 e[1] = 1.262810111e+000;
00463 e[2] = -1.612633824e+000;
00464 e[3] = 0.000000000e+000;
00465 }
00466 static void init_vt( lapack_int size, lapack_complex_float *vt ) {
00467 lapack_int i;
00468 for( i = 0; i < size; i++ ) {
00469 vt[i] = lapack_make_complex_float( 0.0f, 0.0f );
00470 }
00471 vt[0] = lapack_make_complex_float( 1.000000000e+000, 0.000000000e+000 );
00472 vt[8] = lapack_make_complex_float( 0.000000000e+000, 0.000000000e+000 );
00473 vt[16] = lapack_make_complex_float( 0.000000000e+000, 0.000000000e+000 );
00474 vt[24] = lapack_make_complex_float( 0.000000000e+000, 0.000000000e+000 );
00475 vt[1] = lapack_make_complex_float( 0.000000000e+000, 0.000000000e+000 );
00476 vt[9] = lapack_make_complex_float( -2.312345505e-001, -5.404263735e-001 );
00477 vt[17] = lapack_make_complex_float( 1.786240637e-001, -5.887280703e-001 );
00478 vt[25] = lapack_make_complex_float( -4.047983885e-001, -3.348147273e-001 );
00479 vt[2] = lapack_make_complex_float( 0.000000000e+000, 0.000000000e+000 );
00480 vt[10] = lapack_make_complex_float( 4.962260723e-001, 2.648075521e-001 );
00481 vt[18] = lapack_make_complex_float( -3.556231856e-001, -6.888166666e-001 );
00482 vt[26] = lapack_make_complex_float( 2.756552696e-001, -8.194240928e-002 );
00483 vt[3] = lapack_make_complex_float( 0.000000000e+000, 0.000000000e+000 );
00484 vt[11] = lapack_make_complex_float( 2.917790413e-001, 5.029628277e-001 );
00485 vt[19] = lapack_make_complex_float( 4.850706458e-002, 1.349202693e-001 );
00486 vt[27] = lapack_make_complex_float( -7.155819535e-001, -3.595546782e-001 );
00487 }
00488 static void init_u( lapack_int size, lapack_complex_float *u ) {
00489 lapack_int i;
00490 for( i = 0; i < size; i++ ) {
00491 u[i] = lapack_make_complex_float( 0.0f, 0.0f );
00492 }
00493 u[0] = lapack_make_complex_float( -3.109810352e-001, 2.623902261e-001 );
00494 u[8] = lapack_make_complex_float( -6.521001458e-001, -5.532329679e-001 );
00495 u[16] = lapack_make_complex_float( -4.266516492e-002, -3.605512530e-002 );
00496 u[24] = lapack_make_complex_float( 2.634342611e-001, 7.411344349e-002 );
00497 u[1] = lapack_make_complex_float( 3.174597919e-001, -6.413983703e-001 );
00498 u[9] = lapack_make_complex_float( -3.487945795e-001, -7.205679268e-002 );
00499 u[17] = lapack_make_complex_float( -2.287386060e-001, -6.908839103e-003 );
00500 u[25] = lapack_make_complex_float( -1.101401895e-001, 3.261797875e-002 );
00501 u[2] = lapack_make_complex_float( -2.008419037e-001, 1.490117311e-001 );
00502 u[10] = lapack_make_complex_float( 3.102515340e-001, -2.303087153e-002 );
00503 u[18] = lapack_make_complex_float( -1.854601800e-001, 1.817280650e-001 );
00504 u[26] = lapack_make_complex_float( 2.956088483e-001, -5.647819638e-001 );
00505 u[3] = lapack_make_complex_float( 1.198572665e-001, -1.230966449e-001 );
00506 u[11] = lapack_make_complex_float( 4.591666162e-003, 4.640752450e-004 );
00507 u[19] = lapack_make_complex_float( 3.304663301e-001, -4.820711017e-001 );
00508 u[27] = lapack_make_complex_float( 6.749325991e-002, -3.463969827e-001 );
00509 u[4] = lapack_make_complex_float( -2.688689828e-001, -1.652086675e-001 );
00510 u[12] = lapack_make_complex_float( -1.793801785e-001, 5.860745907e-002 );
00511 u[20] = lapack_make_complex_float( 5.235373974e-001, 2.579777837e-001 );
00512 u[28] = lapack_make_complex_float( -3.926950991e-001, -1.449706554e-001 );
00513 u[5] = lapack_make_complex_float( -3.498536646e-001, 9.070278704e-002 );
00514 u[13] = lapack_make_complex_float( -8.288498223e-002, 5.058868229e-002 );
00515 u[21] = lapack_make_complex_float( -3.202392757e-001, -3.037969768e-001 );
00516 u[29] = lapack_make_complex_float( -3.173573017e-001, -3.241353333e-001 );
00517 }
00518 static void init_c( lapack_int size, lapack_complex_float *c ) {
00519 lapack_int i;
00520 for( i = 0; i < size; i++ ) {
00521 c[i] = lapack_make_complex_float( 0.0f, 0.0f );
00522 }
00523 }
00524 static void init_work( lapack_int size, float *work ) {
00525 lapack_int i;
00526 for( i = 0; i < size; i++ ) {
00527 work[i] = 0;
00528 }
00529 }
00530
00531
00532
00533 static int compare_cbdsqr( float *d, float *d_i, float *e, float *e_i,
00534 lapack_complex_float *vt, lapack_complex_float *vt_i,
00535 lapack_complex_float *u, lapack_complex_float *u_i,
00536 lapack_complex_float *c, lapack_complex_float *c_i,
00537 lapack_int info, lapack_int info_i, lapack_int ldc,
00538 lapack_int ldu, lapack_int ldvt, lapack_int n,
00539 lapack_int ncc, lapack_int ncvt, lapack_int nru )
00540 {
00541 lapack_int i;
00542 int failed = 0;
00543 for( i = 0; i < n; i++ ) {
00544 failed += compare_floats(d[i],d_i[i]);
00545 }
00546 for( i = 0; i < n; i++ ) {
00547 failed += compare_floats(e[i],e_i[i]);
00548 }
00549 if( ncvt != 0 ) {
00550 for( i = 0; i < ldvt*ncvt; i++ ) {
00551 failed += compare_complex_floats(vt[i],vt_i[i]);
00552 }
00553 }
00554 if( nru != 0 ) {
00555 for( i = 0; i < ldu*n; i++ ) {
00556 failed += compare_complex_floats(u[i],u_i[i]);
00557 }
00558 }
00559 if( ncc != 0 ) {
00560 for( i = 0; i < ldc*ncc; i++ ) {
00561 failed += compare_complex_floats(c[i],c_i[i]);
00562 }
00563 }
00564 failed += (info == info_i) ? 0 : 1;
00565 if( info != 0 || info_i != 0 ) {
00566 printf( "info=%d, info_i=%d\n",(int)info,(int)info_i );
00567 }
00568
00569 return failed;
00570 }