dhseqr_1.c
Go to the documentation of this file.
00001 /*****************************************************************************
00002   Copyright (c) 2010, Intel Corp.
00003   All rights reserved.
00004 
00005   Redistribution and use in source and binary forms, with or without
00006   modification, are permitted provided that the following conditions are met:
00007 
00008     * Redistributions of source code must retain the above copyright notice,
00009       this list of conditions and the following disclaimer.
00010     * Redistributions in binary form must reproduce the above copyright
00011       notice, this list of conditions and the following disclaimer in the
00012       documentation and/or other materials provided with the distribution.
00013     * Neither the name of Intel Corporation nor the names of its contributors
00014       may be used to endorse or promote products derived from this software
00015       without specific prior written permission.
00016 
00017   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
00018   AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00019   IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
00020   ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
00021   LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
00022   CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
00023   SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
00024   INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
00025   CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
00026   ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
00027   THE POSSIBILITY OF SUCH DAMAGE.
00028 *****************************************************************************/
00029 /*  Contents: test routine for C interface to LAPACK
00030 *   Author: Intel Corporation
00031 *   Created in March, 2010
00032 *
00033 * Purpose
00034 *
00035 * dhseqr_1 is the test program for the C interface to LAPACK
00036 * routine dhseqr
00037 * The program doesn't require an input, the input data is hardcoded in the
00038 * test program.
00039 * The program tests the C interface in the four combinations:
00040 *   1) column-major layout, middle-level interface
00041 *   2) column-major layout, high-level interface
00042 *   3) row-major layout, middle-level interface
00043 *   4) row-major layout, high-level interface
00044 * The output of the C interface function is compared to those obtained from
00045 * the corresponiding LAPACK routine with the same input data, and the
00046 * comparison diagnostics is then printed on the standard output having PASSED
00047 * keyword if the test is passed, and FAILED keyword if the test isn't passed.
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_dhseqr( char *job, char *compz, lapack_int *n,
00055                                  lapack_int *ilo, lapack_int *ihi,
00056                                  lapack_int *ldh, lapack_int *ldz,
00057                                  lapack_int *lwork );
00058 static void init_h( lapack_int size, double *h );
00059 static void init_wr( lapack_int size, double *wr );
00060 static void init_wi( lapack_int size, double *wi );
00061 static void init_z( lapack_int size, double *z );
00062 static void init_work( lapack_int size, double *work );
00063 static int compare_dhseqr( double *h, double *h_i, double *wr, double *wr_i,
00064                            double *wi, double *wi_i, double *z, double *z_i,
00065                            lapack_int info, lapack_int info_i, char compz,
00066                            lapack_int ldh, lapack_int ldz, lapack_int n );
00067 
00068 int main(void)
00069 {
00070     /* Local scalars */
00071     char job, job_i;
00072     char compz, compz_i;
00073     lapack_int n, n_i;
00074     lapack_int ilo, ilo_i;
00075     lapack_int ihi, ihi_i;
00076     lapack_int ldh, ldh_i;
00077     lapack_int ldh_r;
00078     lapack_int ldz, ldz_i;
00079     lapack_int ldz_r;
00080     lapack_int lwork, lwork_i;
00081     lapack_int info, info_i;
00082     lapack_int i;
00083     int failed;
00084 
00085     /* Local arrays */
00086     double *h = NULL, *h_i = NULL;
00087     double *wr = NULL, *wr_i = NULL;
00088     double *wi = NULL, *wi_i = NULL;
00089     double *z = NULL, *z_i = NULL;
00090     double *work = NULL, *work_i = NULL;
00091     double *h_save = NULL;
00092     double *wr_save = NULL;
00093     double *wi_save = NULL;
00094     double *z_save = NULL;
00095     double *h_r = NULL;
00096     double *z_r = NULL;
00097 
00098     /* Iniitialize the scalar parameters */
00099     init_scalars_dhseqr( &job, &compz, &n, &ilo, &ihi, &ldh, &ldz, &lwork );
00100     ldh_r = n+2;
00101     ldz_r = n+2;
00102     job_i = job;
00103     compz_i = compz;
00104     n_i = n;
00105     ilo_i = ilo;
00106     ihi_i = ihi;
00107     ldh_i = ldh;
00108     ldz_i = ldz;
00109     lwork_i = lwork;
00110 
00111     /* Allocate memory for the LAPACK routine arrays */
00112     h = (double *)LAPACKE_malloc( ldh*n * sizeof(double) );
00113     wr = (double *)LAPACKE_malloc( n * sizeof(double) );
00114     wi = (double *)LAPACKE_malloc( n * sizeof(double) );
00115     z = (double *)LAPACKE_malloc( ldz*n * sizeof(double) );
00116     work = (double *)LAPACKE_malloc( lwork * sizeof(double) );
00117 
00118     /* Allocate memory for the C interface function arrays */
00119     h_i = (double *)LAPACKE_malloc( ldh*n * sizeof(double) );
00120     wr_i = (double *)LAPACKE_malloc( n * sizeof(double) );
00121     wi_i = (double *)LAPACKE_malloc( n * sizeof(double) );
00122     z_i = (double *)LAPACKE_malloc( ldz*n * sizeof(double) );
00123     work_i = (double *)LAPACKE_malloc( lwork * sizeof(double) );
00124 
00125     /* Allocate memory for the backup arrays */
00126     h_save = (double *)LAPACKE_malloc( ldh*n * sizeof(double) );
00127     wr_save = (double *)LAPACKE_malloc( n * sizeof(double) );
00128     wi_save = (double *)LAPACKE_malloc( n * sizeof(double) );
00129     z_save = (double *)LAPACKE_malloc( ldz*n * sizeof(double) );
00130 
00131     /* Allocate memory for the row-major arrays */
00132     h_r = (double *)LAPACKE_malloc( n*(n+2) * sizeof(double) );
00133     z_r = (double *)LAPACKE_malloc( n*(n+2) * sizeof(double) );
00134 
00135     /* Initialize input arrays */
00136     init_h( ldh*n, h );
00137     init_wr( n, wr );
00138     init_wi( n, wi );
00139     init_z( ldz*n, z );
00140     init_work( lwork, work );
00141 
00142     /* Backup the ouptut arrays */
00143     for( i = 0; i < ldh*n; i++ ) {
00144         h_save[i] = h[i];
00145     }
00146     for( i = 0; i < n; i++ ) {
00147         wr_save[i] = wr[i];
00148     }
00149     for( i = 0; i < n; i++ ) {
00150         wi_save[i] = wi[i];
00151     }
00152     for( i = 0; i < ldz*n; i++ ) {
00153         z_save[i] = z[i];
00154     }
00155 
00156     /* Call the LAPACK routine */
00157     dhseqr_( &job, &compz, &n, &ilo, &ihi, h, &ldh, wr, wi, z, &ldz, work,
00158              &lwork, &info );
00159 
00160     /* Initialize input data, call the column-major middle-level
00161      * interface to LAPACK routine and check the results */
00162     for( i = 0; i < ldh*n; i++ ) {
00163         h_i[i] = h_save[i];
00164     }
00165     for( i = 0; i < n; i++ ) {
00166         wr_i[i] = wr_save[i];
00167     }
00168     for( i = 0; i < n; i++ ) {
00169         wi_i[i] = wi_save[i];
00170     }
00171     for( i = 0; i < ldz*n; i++ ) {
00172         z_i[i] = z_save[i];
00173     }
00174     for( i = 0; i < lwork; i++ ) {
00175         work_i[i] = work[i];
00176     }
00177     info_i = LAPACKE_dhseqr_work( LAPACK_COL_MAJOR, job_i, compz_i, n_i, ilo_i,
00178                                   ihi_i, h_i, ldh_i, wr_i, wi_i, z_i, ldz_i,
00179                                   work_i, lwork_i );
00180 
00181     failed = compare_dhseqr( h, h_i, wr, wr_i, wi, wi_i, z, z_i, info, info_i,
00182                              compz, ldh, ldz, n );
00183     if( failed == 0 ) {
00184         printf( "PASSED: column-major middle-level interface to dhseqr\n" );
00185     } else {
00186         printf( "FAILED: column-major middle-level interface to dhseqr\n" );
00187     }
00188 
00189     /* Initialize input data, call the column-major high-level
00190      * interface to LAPACK routine and check the results */
00191     for( i = 0; i < ldh*n; i++ ) {
00192         h_i[i] = h_save[i];
00193     }
00194     for( i = 0; i < n; i++ ) {
00195         wr_i[i] = wr_save[i];
00196     }
00197     for( i = 0; i < n; i++ ) {
00198         wi_i[i] = wi_save[i];
00199     }
00200     for( i = 0; i < ldz*n; i++ ) {
00201         z_i[i] = z_save[i];
00202     }
00203     for( i = 0; i < lwork; i++ ) {
00204         work_i[i] = work[i];
00205     }
00206     info_i = LAPACKE_dhseqr( LAPACK_COL_MAJOR, job_i, compz_i, n_i, ilo_i,
00207                              ihi_i, h_i, ldh_i, wr_i, wi_i, z_i, ldz_i );
00208 
00209     failed = compare_dhseqr( h, h_i, wr, wr_i, wi, wi_i, z, z_i, info, info_i,
00210                              compz, ldh, ldz, n );
00211     if( failed == 0 ) {
00212         printf( "PASSED: column-major high-level interface to dhseqr\n" );
00213     } else {
00214         printf( "FAILED: column-major high-level interface to dhseqr\n" );
00215     }
00216 
00217     /* Initialize input data, call the row-major middle-level
00218      * interface to LAPACK routine and check the results */
00219     for( i = 0; i < ldh*n; i++ ) {
00220         h_i[i] = h_save[i];
00221     }
00222     for( i = 0; i < n; i++ ) {
00223         wr_i[i] = wr_save[i];
00224     }
00225     for( i = 0; i < n; i++ ) {
00226         wi_i[i] = wi_save[i];
00227     }
00228     for( i = 0; i < ldz*n; i++ ) {
00229         z_i[i] = z_save[i];
00230     }
00231     for( i = 0; i < lwork; i++ ) {
00232         work_i[i] = work[i];
00233     }
00234 
00235     LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, n, h_i, ldh, h_r, n+2 );
00236     if( LAPACKE_lsame( compz, 'i' ) || LAPACKE_lsame( compz, 'v' ) ) {
00237         LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, n, z_i, ldz, z_r, n+2 );
00238     }
00239     info_i = LAPACKE_dhseqr_work( LAPACK_ROW_MAJOR, job_i, compz_i, n_i, ilo_i,
00240                                   ihi_i, h_r, ldh_r, wr_i, wi_i, z_r, ldz_r,
00241                                   work_i, lwork_i );
00242 
00243     LAPACKE_dge_trans( LAPACK_ROW_MAJOR, n, n, h_r, n+2, h_i, ldh );
00244     if( LAPACKE_lsame( compz, 'i' ) || LAPACKE_lsame( compz, 'v' ) ) {
00245         LAPACKE_dge_trans( LAPACK_ROW_MAJOR, n, n, z_r, n+2, z_i, ldz );
00246     }
00247 
00248     failed = compare_dhseqr( h, h_i, wr, wr_i, wi, wi_i, z, z_i, info, info_i,
00249                              compz, ldh, ldz, n );
00250     if( failed == 0 ) {
00251         printf( "PASSED: row-major middle-level interface to dhseqr\n" );
00252     } else {
00253         printf( "FAILED: row-major middle-level interface to dhseqr\n" );
00254     }
00255 
00256     /* Initialize input data, call the row-major high-level
00257      * interface to LAPACK routine and check the results */
00258     for( i = 0; i < ldh*n; i++ ) {
00259         h_i[i] = h_save[i];
00260     }
00261     for( i = 0; i < n; i++ ) {
00262         wr_i[i] = wr_save[i];
00263     }
00264     for( i = 0; i < n; i++ ) {
00265         wi_i[i] = wi_save[i];
00266     }
00267     for( i = 0; i < ldz*n; i++ ) {
00268         z_i[i] = z_save[i];
00269     }
00270     for( i = 0; i < lwork; i++ ) {
00271         work_i[i] = work[i];
00272     }
00273 
00274     /* Init row_major arrays */
00275     LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, n, h_i, ldh, h_r, n+2 );
00276     if( LAPACKE_lsame( compz, 'i' ) || LAPACKE_lsame( compz, 'v' ) ) {
00277         LAPACKE_dge_trans( LAPACK_COL_MAJOR, n, n, z_i, ldz, z_r, n+2 );
00278     }
00279     info_i = LAPACKE_dhseqr( LAPACK_ROW_MAJOR, job_i, compz_i, n_i, ilo_i,
00280                              ihi_i, h_r, ldh_r, wr_i, wi_i, z_r, ldz_r );
00281 
00282     LAPACKE_dge_trans( LAPACK_ROW_MAJOR, n, n, h_r, n+2, h_i, ldh );
00283     if( LAPACKE_lsame( compz, 'i' ) || LAPACKE_lsame( compz, 'v' ) ) {
00284         LAPACKE_dge_trans( LAPACK_ROW_MAJOR, n, n, z_r, n+2, z_i, ldz );
00285     }
00286 
00287     failed = compare_dhseqr( h, h_i, wr, wr_i, wi, wi_i, z, z_i, info, info_i,
00288                              compz, ldh, ldz, n );
00289     if( failed == 0 ) {
00290         printf( "PASSED: row-major high-level interface to dhseqr\n" );
00291     } else {
00292         printf( "FAILED: row-major high-level interface to dhseqr\n" );
00293     }
00294 
00295     /* Release memory */
00296     if( h != NULL ) {
00297         LAPACKE_free( h );
00298     }
00299     if( h_i != NULL ) {
00300         LAPACKE_free( h_i );
00301     }
00302     if( h_r != NULL ) {
00303         LAPACKE_free( h_r );
00304     }
00305     if( h_save != NULL ) {
00306         LAPACKE_free( h_save );
00307     }
00308     if( wr != NULL ) {
00309         LAPACKE_free( wr );
00310     }
00311     if( wr_i != NULL ) {
00312         LAPACKE_free( wr_i );
00313     }
00314     if( wr_save != NULL ) {
00315         LAPACKE_free( wr_save );
00316     }
00317     if( wi != NULL ) {
00318         LAPACKE_free( wi );
00319     }
00320     if( wi_i != NULL ) {
00321         LAPACKE_free( wi_i );
00322     }
00323     if( wi_save != NULL ) {
00324         LAPACKE_free( wi_save );
00325     }
00326     if( z != NULL ) {
00327         LAPACKE_free( z );
00328     }
00329     if( z_i != NULL ) {
00330         LAPACKE_free( z_i );
00331     }
00332     if( z_r != NULL ) {
00333         LAPACKE_free( z_r );
00334     }
00335     if( z_save != NULL ) {
00336         LAPACKE_free( z_save );
00337     }
00338     if( work != NULL ) {
00339         LAPACKE_free( work );
00340     }
00341     if( work_i != NULL ) {
00342         LAPACKE_free( work_i );
00343     }
00344 
00345     return 0;
00346 }
00347 
00348 /* Auxiliary function: dhseqr scalar parameters initialization */
00349 static void init_scalars_dhseqr( char *job, char *compz, lapack_int *n,
00350                                  lapack_int *ilo, lapack_int *ihi,
00351                                  lapack_int *ldh, lapack_int *ldz,
00352                                  lapack_int *lwork )
00353 {
00354     *job = 'S';
00355     *compz = 'V';
00356     *n = 4;
00357     *ilo = 2;
00358     *ihi = 4;
00359     *ldh = 8;
00360     *ldz = 8;
00361     *lwork = 512;
00362 
00363     return;
00364 }
00365 
00366 /* Auxiliary functions: dhseqr array parameters initialization */
00367 static void init_h( lapack_int size, double *h ) {
00368     lapack_int i;
00369     for( i = 0; i < size; i++ ) {
00370         h[i] = 0;
00371     }
00372     h[0] = -4.00000000000000020e-001;  /* h[0,0] */
00373     h[8] = 3.20000000000000020e+000;  /* h[0,1] */
00374     h[16] = -9.22503046669019610e-001;  /* h[0,2] */
00375     h[24] = -7.69148803086154140e+000;  /* h[0,3] */
00376     h[1] = 0.00000000000000000e+000;  /* h[1,0] */
00377     h[9] = 2.00000000000000010e-001;  /* h[1,1] */
00378     h[17] = -4.38260274722939160e+000;  /* h[1,2] */
00379     h[25] = 4.67492684410557580e-001;  /* h[1,3] */
00380     h[2] = 0.00000000000000000e+000;  /* h[2,0] */
00381     h[10] = -2.94416371827383250e+000;  /* h[2,1] */
00382     h[18] = -8.93239118145844560e-001;  /* h[2,2] */
00383     h[26] = -6.79197286602599700e-001;  /* h[2,3] */
00384     h[3] = 0.00000000000000000e+000;  /* h[3,0] */
00385     h[11] = 7.26487042240654410e-001;  /* h[3,1] */
00386     h[19] = -2.13919728660259920e+000;  /* h[3,2] */
00387     h[27] = 6.69323911814584530e+000;  /* h[3,3] */
00388 }
00389 static void init_wr( lapack_int size, double *wr ) {
00390     lapack_int i;
00391     for( i = 0; i < size; i++ ) {
00392         wr[i] = 0;
00393     }
00394 }
00395 static void init_wi( lapack_int size, double *wi ) {
00396     lapack_int i;
00397     for( i = 0; i < size; i++ ) {
00398         wi[i] = 0;
00399     }
00400 }
00401 static void init_z( lapack_int size, double *z ) {
00402     lapack_int i;
00403     for( i = 0; i < size; i++ ) {
00404         z[i] = 0;
00405     }
00406     z[0] = 1.00000000000000000e+000;  /* z[0,0] */
00407     z[8] = 0.00000000000000000e+000;  /* z[0,1] */
00408     z[16] = 0.00000000000000000e+000;  /* z[0,2] */
00409     z[24] = 0.00000000000000000e+000;  /* z[0,3] */
00410     z[1] = 0.00000000000000000e+000;  /* z[1,0] */
00411     z[9] = 1.00000000000000000e+000;  /* z[1,1] */
00412     z[17] = 0.00000000000000000e+000;  /* z[1,2] */
00413     z[25] = 0.00000000000000000e+000;  /* z[1,3] */
00414     z[2] = 0.00000000000000000e+000;  /* z[2,0] */
00415     z[10] = 0.00000000000000000e+000;  /* z[2,1] */
00416     z[18] = -3.09086072337558140e-001;  /* z[2,2] */
00417     z[26] = -9.51034068730947980e-001;  /* z[2,3] */
00418     z[3] = 0.00000000000000000e+000;  /* z[3,0] */
00419     z[11] = 0.00000000000000000e+000;  /* z[3,1] */
00420     z[19] = -9.51034068730947980e-001;  /* z[3,2] */
00421     z[27] = 3.09086072337558360e-001;  /* z[3,3] */
00422 }
00423 static void init_work( lapack_int size, double *work ) {
00424     lapack_int i;
00425     for( i = 0; i < size; i++ ) {
00426         work[i] = 0;
00427     }
00428 }
00429 
00430 /* Auxiliary function: C interface to dhseqr results check */
00431 /* Return value: 0 - test is passed, non-zero - test is failed */
00432 static int compare_dhseqr( double *h, double *h_i, double *wr, double *wr_i,
00433                            double *wi, double *wi_i, double *z, double *z_i,
00434                            lapack_int info, lapack_int info_i, char compz,
00435                            lapack_int ldh, lapack_int ldz, lapack_int n )
00436 {
00437     lapack_int i;
00438     int failed = 0;
00439     for( i = 0; i < ldh*n; i++ ) {
00440         failed += compare_doubles(h[i],h_i[i]);
00441     }
00442     for( i = 0; i < n; i++ ) {
00443         failed += compare_doubles(wr[i],wr_i[i]);
00444     }
00445     for( i = 0; i < n; i++ ) {
00446         failed += compare_doubles(wi[i],wi_i[i]);
00447     }
00448     if( LAPACKE_lsame( compz, 'i' ) || LAPACKE_lsame( compz, 'v' ) ) {
00449         for( i = 0; i < ldz*n; i++ ) {
00450             failed += compare_doubles(z[i],z_i[i]);
00451         }
00452     }
00453     failed += (info == info_i) ? 0 : 1;
00454     if( info != 0 || info_i != 0 ) {
00455         printf( "info=%d, info_i=%d\n",(int)info,(int)info_i );
00456     }
00457 
00458     return failed;
00459 }


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:45