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_dsptrs( char *uplo, lapack_int *n, lapack_int *nrhs,
00055 lapack_int *ldb );
00056 static void init_ap( lapack_int size, double *ap );
00057 static void init_ipiv( lapack_int size, lapack_int *ipiv );
00058 static void init_b( lapack_int size, double *b );
00059 static int compare_dsptrs( double *b, double *b_i, lapack_int info,
00060 lapack_int info_i, lapack_int ldb, lapack_int nrhs );
00061
00062 int main(void)
00063 {
00064
00065 char uplo, uplo_i;
00066 lapack_int n, n_i;
00067 lapack_int nrhs, nrhs_i;
00068 lapack_int ldb, ldb_i;
00069 lapack_int ldb_r;
00070 lapack_int info, info_i;
00071 lapack_int i;
00072 int failed;
00073
00074
00075 double *ap = NULL, *ap_i = NULL;
00076 lapack_int *ipiv = NULL, *ipiv_i = NULL;
00077 double *b = NULL, *b_i = NULL;
00078 double *b_save = NULL;
00079 double *ap_r = NULL;
00080 double *b_r = NULL;
00081
00082
00083 init_scalars_dsptrs( &uplo, &n, &nrhs, &ldb );
00084 ldb_r = nrhs+2;
00085 uplo_i = uplo;
00086 n_i = n;
00087 nrhs_i = nrhs;
00088 ldb_i = ldb;
00089
00090
00091 ap = (double *)LAPACKE_malloc( ((n*(n+1)/2)) * sizeof(double) );
00092 ipiv = (lapack_int *)LAPACKE_malloc( n * sizeof(lapack_int) );
00093 b = (double *)LAPACKE_malloc( ldb*nrhs * sizeof(double) );
00094
00095
00096 ap_i = (double *)LAPACKE_malloc( ((n*(n+1)/2)) * sizeof(double) );
00097 ipiv_i = (lapack_int *)LAPACKE_malloc( n * sizeof(lapack_int) );
00098 b_i = (double *)LAPACKE_malloc( ldb*nrhs * sizeof(double) );
00099
00100
00101 b_save = (double *)LAPACKE_malloc( ldb*nrhs * sizeof(double) );
00102
00103
00104 ap_r = (double *)LAPACKE_malloc( n*(n+1)/2 * sizeof(double) );
00105 b_r = (double *)LAPACKE_malloc( n*(nrhs+2) * sizeof(double) );
00106
00107
00108 init_ap( (n*(n+1)/2), ap );
00109 init_ipiv( n, ipiv );
00110 init_b( ldb*nrhs, b );
00111
00112
00113 for( i = 0; i < ldb*nrhs; i++ ) {
00114 b_save[i] = b[i];
00115 }
00116
00117
00118 dsptrs_( &uplo, &n, &nrhs, ap, ipiv, b, &ldb, &info );
00119
00120
00121
00122 for( i = 0; i < (n*(n+1)/2); i++ ) {
00123 ap_i[i] = ap[i];
00124 }
00125 for( i = 0; i < n; i++ ) {
00126 ipiv_i[i] = ipiv[i];
00127 }
00128 for( i = 0; i < ldb*nrhs; i++ ) {
00129 b_i[i] = b_save[i];
00130 }
00131 info_i = LAPACKE_dsptrs_work( LAPACK_COL_MAJOR, uplo_i, n_i, nrhs_i, ap_i,
00132 ipiv_i, b_i, ldb_i );
00133
00134 failed = compare_dsptrs( b, b_i, info, info_i, ldb, nrhs );
00135 if( failed == 0 ) {
00136 printf( "PASSED: column-major middle-level interface to dsptrs\n" );
00137 } else {
00138 printf( "FAILED: column-major middle-level interface to dsptrs\n" );
00139 }
00140
00141
00142
00143 for( i = 0; i < (n*(n+1)/2); i++ ) {
00144 ap_i[i] = ap[i];
00145 }
00146 for( i = 0; i < n; i++ ) {
00147 ipiv_i[i] = ipiv[i];
00148 }
00149 for( i = 0; i < ldb*nrhs; i++ ) {
00150 b_i[i] = b_save[i];
00151 }
00152 info_i = LAPACKE_dsptrs( LAPACK_COL_MAJOR, uplo_i, n_i, nrhs_i, ap_i,
00153 ipiv_i, b_i, ldb_i );
00154
00155 failed = compare_dsptrs( b, b_i, info, info_i, ldb, nrhs );
00156 if( failed == 0 ) {
00157 printf( "PASSED: column-major high-level interface to dsptrs\n" );
00158 } else {
00159 printf( "FAILED: column-major high-level interface to dsptrs\n" );
00160 }
00161
00162
00163
00164 for( i = 0; i < (n*(n+1)/2); i++ ) {
00165 ap_i[i] = ap[i];
00166 }
00167 for( i = 0; i < n; i++ ) {
00168 ipiv_i[i] = ipiv[i];
00169 }
00170 for( i = 0; i < ldb*nrhs; i++ ) {
00171 b_i[i] = b_save[i];
00172 }
00173
00174 LAPACKE_dpp_trans( LAPACK_COL_MAJOR, uplo, n, ap_i, ap_r );
00175 LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, nrhs, b_i, ldb, b_r, nrhs+2 );
00176 info_i = LAPACKE_dsptrs_work( LAPACK_ROW_MAJOR, uplo_i, n_i, nrhs_i, ap_r,
00177 ipiv_i, b_r, ldb_r );
00178
00179 LAPACKE_dge_trans( LAPACK_ROW_MAJOR, n, nrhs, b_r, nrhs+2, b_i, ldb );
00180
00181 failed = compare_dsptrs( b, b_i, info, info_i, ldb, nrhs );
00182 if( failed == 0 ) {
00183 printf( "PASSED: row-major middle-level interface to dsptrs\n" );
00184 } else {
00185 printf( "FAILED: row-major middle-level interface to dsptrs\n" );
00186 }
00187
00188
00189
00190 for( i = 0; i < (n*(n+1)/2); i++ ) {
00191 ap_i[i] = ap[i];
00192 }
00193 for( i = 0; i < n; i++ ) {
00194 ipiv_i[i] = ipiv[i];
00195 }
00196 for( i = 0; i < ldb*nrhs; i++ ) {
00197 b_i[i] = b_save[i];
00198 }
00199
00200
00201 LAPACKE_dpp_trans( LAPACK_COL_MAJOR, uplo, n, ap_i, ap_r );
00202 LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, nrhs, b_i, ldb, b_r, nrhs+2 );
00203 info_i = LAPACKE_dsptrs( LAPACK_ROW_MAJOR, uplo_i, n_i, nrhs_i, ap_r,
00204 ipiv_i, b_r, ldb_r );
00205
00206 LAPACKE_dge_trans( LAPACK_ROW_MAJOR, n, nrhs, b_r, nrhs+2, b_i, ldb );
00207
00208 failed = compare_dsptrs( b, b_i, info, info_i, ldb, nrhs );
00209 if( failed == 0 ) {
00210 printf( "PASSED: row-major high-level interface to dsptrs\n" );
00211 } else {
00212 printf( "FAILED: row-major high-level interface to dsptrs\n" );
00213 }
00214
00215
00216 if( ap != NULL ) {
00217 LAPACKE_free( ap );
00218 }
00219 if( ap_i != NULL ) {
00220 LAPACKE_free( ap_i );
00221 }
00222 if( ap_r != NULL ) {
00223 LAPACKE_free( ap_r );
00224 }
00225 if( ipiv != NULL ) {
00226 LAPACKE_free( ipiv );
00227 }
00228 if( ipiv_i != NULL ) {
00229 LAPACKE_free( ipiv_i );
00230 }
00231 if( b != NULL ) {
00232 LAPACKE_free( b );
00233 }
00234 if( b_i != NULL ) {
00235 LAPACKE_free( b_i );
00236 }
00237 if( b_r != NULL ) {
00238 LAPACKE_free( b_r );
00239 }
00240 if( b_save != NULL ) {
00241 LAPACKE_free( b_save );
00242 }
00243
00244 return 0;
00245 }
00246
00247
00248 static void init_scalars_dsptrs( char *uplo, lapack_int *n, lapack_int *nrhs,
00249 lapack_int *ldb )
00250 {
00251 *uplo = 'L';
00252 *n = 4;
00253 *nrhs = 2;
00254 *ldb = 8;
00255
00256 return;
00257 }
00258
00259
00260 static void init_ap( lapack_int size, double *ap ) {
00261 lapack_int i;
00262 for( i = 0; i < size; i++ ) {
00263 ap[i] = 0;
00264 }
00265 ap[0] = 2.06999999999999980e+000;
00266 ap[1] = 4.20000000000000020e+000;
00267 ap[2] = 2.23041384055834070e-001;
00268 ap[3] = 6.53658376748910360e-001;
00269 ap[4] = 1.14999999999999990e+000;
00270 ap[5] = 8.11501032143910230e-001;
00271 ap[6] = -5.95969723778629450e-001;
00272 ap[7] = -2.59067708640519000e+000;
00273 ap[8] = 3.03084679550618070e-001;
00274 ap[9] = 4.07385198134887610e-001;
00275 }
00276 static void init_ipiv( lapack_int size, lapack_int *ipiv ) {
00277 lapack_int i;
00278 for( i = 0; i < size; i++ ) {
00279 ipiv[i] = 0;
00280 }
00281 ipiv[0] = -3;
00282 ipiv[1] = -3;
00283 ipiv[2] = 3;
00284 ipiv[3] = 4;
00285 }
00286 static void init_b( lapack_int size, double *b ) {
00287 lapack_int i;
00288 for( i = 0; i < size; i++ ) {
00289 b[i] = 0;
00290 }
00291 b[0] = -9.50000000000000000e+000;
00292 b[8] = 2.78500000000000010e+001;
00293 b[1] = -8.38000000000000080e+000;
00294 b[9] = 9.90000000000000040e+000;
00295 b[2] = -6.07000000000000030e+000;
00296 b[10] = 1.92500000000000000e+001;
00297 b[3] = -9.59999999999999960e-001;
00298 b[11] = 3.93000000000000020e+000;
00299 }
00300
00301
00302
00303 static int compare_dsptrs( double *b, double *b_i, lapack_int info,
00304 lapack_int info_i, lapack_int ldb, lapack_int nrhs )
00305 {
00306 lapack_int i;
00307 int failed = 0;
00308 for( i = 0; i < ldb*nrhs; i++ ) {
00309 failed += compare_doubles(b[i],b_i[i]);
00310 }
00311 failed += (info == info_i) ? 0 : 1;
00312 if( info != 0 || info_i != 0 ) {
00313 printf( "info=%d, info_i=%d\n",(int)info,(int)info_i );
00314 }
00315
00316 return failed;
00317 }