vslsbasic.c

使用INTEL矢量统计类库的程序,包括以下功能: &#61623 Raw and central moments up to 4th order &#61623 Kurtosis and

/*******************************************************************************
!                             INTEL CONFIDENTIAL
!  Copyright(C) 2007-2008 Intel Corporation. All Rights Reserved.
!  The source code contained  or  described herein and all documents related to
!  the source code ("Material") are owned by Intel Corporation or its suppliers
!  or licensors.  Title to the  Material remains with  Intel Corporation or its
!  suppliers and licensors. The Material contains trade secrets and proprietary
!  and  confidential  information of  Intel or its suppliers and licensors. The
!  Material  is  protected  by  worldwide  copyright  and trade secret laws and
!  treaty  provisions. No part of the Material may be used, copied, reproduced,
!  modified, published, uploaded, posted, transmitted, distributed or disclosed
!  in any way without Intel's prior express written permission.
!  No license  under any  patent, copyright, trade secret or other intellectual
!  property right is granted to or conferred upon you by disclosure or delivery
!  of the Materials,  either expressly, by implication, inducement, estoppel or
!  otherwise.  Any  license  under  such  intellectual property  rights must be
!  express and approved by Intel in writing.
!
!*******************************************************************************
!  Content:
!    Calculation of basic statistics Example Program Text
!******************************************************************************/


#include "mkl.h"
#include "vsl_ss.h"
#include "stdio.h"

#define SEED    777
#define BRNG    VSL_BRNG_MCG31
#define METHOD  0

#define DIM 4        /* dimension of the task */ 
#define N   1000     /* number of observations */

int main()
{
    VSLSSTaskPtr task=0;
    VSLStreamStatePtr stream;
    MKL_INT dim=DIM, n=N, storage=VSL_SS_MATRIX_COLUMNS_STORAGE;
    
    float *x;  /* matrix of observations */
    
    float min_estimator[DIM], max_estimator[DIM];
    float mean[DIM];
    float raw2[DIM], raw3[DIM], raw4[DIM];
    float central2[DIM], central3[DIM], central4[DIM];
    float skewness[DIM], kurtosis[DIM], variation[DIM];
    float *cov, *cor;
    MKL_INT cov_storage=VSL_SS_MATRIX_FULL_STORAGE, 
        cor_storage=VSL_SS_MATRIX_FULL_STORAGE;
    float a=5.0, sigma=1.0;
    int status, errcode;
    int i,j;
    vsl_ss_uint64_t to_compute = 0;

    x = (float *) malloc(DIM*N*sizeof(float));
    if(x==0) {printf("Error in malloc\n");return(0);}

    cov = (float *) malloc(DIM*DIM*sizeof(float));
    if(cov==0) {printf("Error in malloc\n");return(0);}

    cor = (float *) malloc(DIM*DIM*sizeof(float));
    if(cor==0) {printf("Error in malloc\n");return(0);}

    for(i=0;i<DIM;i++) 
    {
        min_estimator[i]=0.0f;
        max_estimator[i]=0.0f;
        mean[i]=0.0f;
        raw2[i]=0.0f;
        raw3[i]=0.0f;
        raw4[i]=0.0f;
        central2[i]=0.0f;
        central3[i]=0.0f;
        central4[i]=0.0f;
        skewness[i]=0.0f;
        kurtosis[i]=0.0f;
        variation[i]=0.0f;
    }

    for(i=0;i<DIM*DIM;i++) 
    {
        cov[i]=cor[i]=0.0f;
    }
    
    /***** Initialize *****/
    errcode = vslNewStream( &stream, BRNG,  SEED );
    if(errcode < 0) {printf("Error in VSL\n");return(0);}
    
    /***** Call RNG *****/
    errcode = vsRngGaussian( METHOD, stream, N*DIM, (float*)x, a, sigma );
    if(errcode < 0) {printf("Error in VSL\n");return(0);}

    /***** Deinitialize *****/
    errcode = vslDeleteStream( &stream );
    if(errcode < 0) {printf("Error in VSL\n");return(0);}
    

    
    status = vslsSSNewTask( &task, &dim, &n, &storage, (float*)x, 0, 0 );
    if(status!=0) printf("status=%i after vsldSSNewTask\n", status);
    
    status = vslsSSEditTask( task, VSL_SS_MIN_ARRAY, min_estimator );
    if(status!=0) printf("status=%i after vsldSSEditTask\n", status);
    
    status = vslsSSEditTask( task, VSL_SS_MAX_ARRAY, max_estimator );
    if(status!=0) printf("status=%i after vsldSSEditTask\n", status);

    to_compute |= VSL_SS_MIN|VSL_SS_MAX;

    status = vslsSSEditMoments( task, mean, raw2, raw3, raw4, 
                                central2, central3, central4 ); 
    if(status!=0) printf("status=%i after vsldSSEditMoments\n", status);

    to_compute |= VSL_SS_MEAN | 
        VSL_SS_2RAW_MOMENT | VSL_SS_3RAW_MOMENT | VSL_SS_4RAW_MOMENT |
        VSL_SS_2CENTRAL_MOMENT | VSL_SS_3CENTRAL_MOMENT | VSL_SS_4CENTRAL_MOMENT;
    
    status = vslsSSEditTask( task, VSL_SS_KURTOSIS_ARRAY, kurtosis );
    if(status!=0) printf("status=%i after vsldSSEditTask\n", status);
    
    status = vslsSSEditTask( task, VSL_SS_SKEWNESS_ARRAY, skewness );
    if(status!=0) printf("status=%i after vsldSSEditTask\n", status);

    status = vslsSSEditTask( task, VSL_SS_VARIATION_ARRAY, variation );
    if(status!=0) printf("status=%i after vsldSSEditTask\n", status);

    to_compute |= VSL_SS_KURTOSIS | VSL_SS_SKEWNESS | VSL_SS_VARIATION;

    status = vslsSSEditCovCor( task, mean, (float*)cov, &cov_storage , (float*)cor, &cor_storage );
    if(status!=0) printf("status=%i after vsldSSEditCovCor\n", status);

    to_compute |= VSL_SS_COVARIANCE_MATRIX|VSL_SS_CORRELATION_MATRIX;
    
    status = vslsSSCompute(task, to_compute, VSL_SS_FAST_METHOD );
    if(status!=0) printf("status=%i after vsldSSCompute\n", status);
    
    status = vslSSDeleteTask( &task );
    if(status!=0) printf("status=%i after vslSSDeleteTask\n", status);

    printf(" Dimension of the task: %d\n", DIM);
    printf("Number of observations: %d\n\n", N);

    printf("               Min        Max        Mean       2nd_raw    3rd_raw    4th_raw    2nd_cen    3rd_cen    4th_cen\n");
    for(i=0;i<DIM;i++)
    {
        printf("Variable #%i:  %+lf  %+lf  %+lf  %+lf  %+lf  %+lf  %+lf  %+lf  %+lf\n", i, 
            min_estimator[i], max_estimator[i], mean[i], 
            raw2[i], raw3[i], raw4[i], 
            central2[i], central3[i], central4[i]);
    }

    printf("               Kurtosis   Skewness   Variation\n");
    for(i=0;i<DIM;i++)
    {
        printf("Variable #%i:  %+lf  %+lf  %+lf\n", i, 
            kurtosis[i], skewness[i], variation[i]);
    }

    printf("\n Computed covariance matrix                 Computed correlation matrix\n");
    for(i=0;i<4;i++)
    {
        for(j=0;j<4;j++)
        {
            printf("%+9lf ", cov[i+j*DIM]);
        }

        printf("   ");

        for(j=0;j<4;j++)
        {
            printf("%+9lf ", cor[i+j*DIM]);
        }
        printf("\n");
    }

    printf("\n");

    free(x);
    free(cov);
    free(cor);

    printf("\nTEST PASSED\n");

    return 0;
}