vslsdataprocessinginchunks.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:
!    Data processing in chunks  Example Program Text
!******************************************************************************/
#include "mkl.h"
#include "vsl_ss.h"
#include "stdio.h"
#include "math.h"

#define BRNG    VSL_BRNG_MCG31     /* VSL basic generator to be used */
#define SEED    7777777      /* Initial value for stream initialization */
#define METHOD VSL_METHOD_DGAUSSIANMV_ICDF

#define RETURN_ON_ERROR                 \
    if(errcode < 0)                     \
    {                                   \
        printf("Error: %i\n", errcode); \
        return errcode;                 \
    }
int generate_input(int, int, float *, float *, float *, float *);

#define P         5          /* dimension of the task */ 
#define N         100000     /* number of observations */
static float X[N*P];        /* Input matrix of observations */
static float S[P*P];        /* Covariance matrix of input data */
static float Cov[P][P];     /* Estimated covariance matrix in chunks*/
static float CovFull[P][P];  /* Estimated covariance matrix*/


int main(void)
{
    MKL_INT p = P, n = N, n_portion, i, j, counter;
    float mean[P]     = {-1.0, 1.0, 2.0, -3.0, 4.0};/* True mean vector */
    float Variance[P] = { 1.0, 1.0, 1.0,  1.0, 1.0};/* True variance vector */
    float W[2] = {0.0, 0.0}; /* array of weights
                                 used in covariance estimation */
    float weights[N];
    int indices[P] = {0, 1, 1, 0, 1};
    float Xmean[P];         /* Mean vector  */
    float Raw2Mom[P];       /* raw 2 moment vector  */
    float Central2Mom[P];   /* Central 2 moment vector  */
    float Maximum[P];       /* vector of max values */
    float Minimum[P];       /* vector of min values */
    float MaximumFull[P];   /* vector of max values */
    float MinimumFull[P];   /* vector of min values */
    float *X_portion,dif,temp1,temp2;
    MKL_INT Cov_storage;
    int errcode, test=0;

    VSLSSTaskPtr task = 0;
    MKL_INT storage_format_x;

    n_portion = n/100;

   /* Generation of observations matrix */
    errcode = generate_input(n, p, X, mean, Variance, weights);
    RETURN_ON_ERROR;

    X_portion = X;

    storage_format_x = VSL_SS_MATRIX_COLUMNS_STORAGE;

    /* Create task */
    errcode = vslsSSNewTask( &task, &p, &n_portion,  &storage_format_x,
                             X_portion, weights, indices );
    RETURN_ON_ERROR;

    /* Register array of weights in the task */
    errcode = vslsSSEditTask( task, VSL_SS_ACCUMULATED_WEIGHT, W );
    RETURN_ON_ERROR;

    /* Register array for moments estimate */  
    errcode = vslsSSEditMoments( task, Xmean, Raw2Mom,0,0,Central2Mom, 0, 0);
    RETURN_ON_ERROR;

    Cov_storage = VSL_SS_MATRIX_FULL_STORAGE;

    /* Register array for covariance estimate */  
    errcode = vslsSSEditCovCor( task, Xmean, (float*)Cov, &Cov_storage , 0, 0);
    RETURN_ON_ERROR;

    /* Inialization of the arrays */
    for( i = 0; i < p; i++ )
    {
        Maximum[i]=0.;
        Minimum[i]=0.;
        Xmean[i] = 0.0;
        Central2Mom[i] = 0.0;
        Raw2Mom[i] = 0.0;
        for( j = 0; j < p; j++ )
        {
            Cov[i][j] = 0.0;
        }
    }

    for( counter = 0; counter < 100; counter++ )
    {
       X_portion = X + counter*n_portion*p;

    /* Register new pointer in the task */
       errcode = vslsSSEditTask( task, VSL_SS_OBSERVATIONS, X_portion );
       RETURN_ON_ERROR;

    errcode = vslsSSEditTask( task, VSL_SS_MAX_ARRAY, Maximum);
    RETURN_ON_ERROR;
    errcode = vslsSSEditTask( task, VSL_SS_MIN_ARRAY, Minimum );
    RETURN_ON_ERROR;

    /* Compute covariance matrix, variances and means estimates */
       errcode = vslsSSCompute(task,VSL_SS_MEAN| 
                                VSL_SS_2CENTRAL_MOMENT|
                                VSL_SS_COVARIANCE_MATRIX|
                                VSL_SS_MIN|VSL_SS_MAX,
                                VSL_SS_FAST_METHOD );
        RETURN_ON_ERROR;
    }

    printf("\n    Mins are: ");
    for(j = 0; j < p; j++)
    {
        if(indices[j] == 0) continue;
        printf("%lf ", Minimum[j]);
    }
    printf("\n");

    printf("\n    Maxs are: ");
    for(j = 0; j < p; j++)
    {
        if(indices[j] == 0) continue;
        printf("%lf ", Maximum[j]);
    }
    printf("\n");

    printf("\n    Means are: ");
    for(j = 0; j < p; j++)
    {
        if(indices[j] == 0) continue;
        printf("%f ", Xmean[j]);
    }
    printf("\n");

    printf("\n Variances are: ");
    for(j = 0; j < p; j++)
    {
        if(indices[j] == 0) continue;
        printf("%f ", Central2Mom[j]);
    }
    printf("\n");


    printf("\nCovariance matrix\n\n");
    for(i = 0; i < p; i++)
    {
        for(j = 0; j < p; j++)
        {
            printf("%f ", Cov[i][j]);
        }
        printf("\n");
    }
    printf("\n");

    /* Validation */

    errcode = vslsSSEditTask( task, VSL_SS_MAX_ARRAY, MaximumFull);
    RETURN_ON_ERROR;
    errcode = vslsSSEditTask( task, VSL_SS_MIN_ARRAY, MinimumFull);
    RETURN_ON_ERROR;
    W[0] = 0.0;
    W[1] = 0.0;
    for( i = 0; i < p; i++ )
    {
        MaximumFull[i]=0.;
        MinimumFull[i]=0.;
        Xmean[i] = 0.0;
        Central2Mom[i] = 0.0;
        Raw2Mom[i] = 0.0;
        for( j = 0; j < p; j++ )
        {
            CovFull[i][j] = 0.0;
        }
    }
    /* Matrix norm*/

    errcode = vslsSSEditTask( task, VSL_SS_OBSERVATIONS, X );
    RETURN_ON_ERROR;
    errcode = vsliSSEditTask( task, VSL_SS_OBSERVATIONS_NUMBER, &n );
    RETURN_ON_ERROR;
    errcode = vslsSSEditTask( task,VSL_SS_COVARIANCE_ARRAY,(float*)CovFull );
    RETURN_ON_ERROR;

    /* Compute covariance matrix, mins and maxs estimates */
    errcode = vslsSSCompute(task,VSL_SS_COVARIANCE_MATRIX|
                        VSL_SS_MIN|VSL_SS_MAX, VSL_SS_FAST_METHOD );

    dif = 1.0;
    temp1 = 0.0;
    temp2 = 0.0;
    for ( i = 0; i < p; i++ )
    {
        for ( j = 0; j < p; j++ )
        {
            temp1 += (Cov[i][j]-CovFull[i][j])*(Cov[i][j]-CovFull[i][j]);
            temp2 += Cov[i][j]*Cov[i][j];
        }
    }
    dif = sqrt(temp1/temp2);
    if(dif > 1e-1)
    {
       test = -1;
    }

    for(j = 0; j < p; j++)
    {
        if(indices[j] == 0) continue;
        if(!MinimumFull[j] == Minimum[j]) test = -1;
        if(!MaximumFull[j] == Maximum[j]) test = -1;
    }

    if( test != 0 )
    {
        printf("\nTEST FAILED\n");
    }
    else
    {
        printf("\nTEST PASSED\n");
    }

    /* Delete task */
    errcode = vslSSDeleteTask( &task );
    RETURN_ON_ERROR;

    return 0;
}

int generate_input(int n, int p,  float *X, float *mean, float *Variance,
                   float *weights)
{
    MKL_INT i, j;
    MKL_INT errcode;
    float Ro = 1. / ( p * 2 );

    /* Following variables are used in Cholesky factorization subroutine */
    char uplo;

    VSLStreamStatePtr stream_good;

    /* Definition of parameters (covariance matrix, means)
           for input data */
    for ( i = 0; i < n; i++)
    {
        weights[i] = 1.;
    }
    for ( i = 0; i < p; i++)
    {
        for( j = 0; j < i; j++ )
        {
            S[i*p+j] = Ro;
            S[j*p+i] = Ro;
        }
        S[i*p+i] = Variance [i];
    }
   /* Here start generation of observations matrix
      (multivarite Gaussian random numbers) */
    uplo = 'U';
    /* MKL Choelsky factorization routine call */
    spotrf( &uplo, &p, S, &p, &errcode );
    RETURN_ON_ERROR;

    /* Stream initialization */
    errcode = vslNewStream( &stream_good, BRNG, SEED );
    RETURN_ON_ERROR;

    /* Generating random numbers from multivariate normal distribution */
    errcode = vsRngGaussianMV( METHOD, stream_good, n, X, p,
        VSL_MATRIX_STORAGE_FULL, mean, S );
    RETURN_ON_ERROR;

    /* Stream  deallocation */
    errcode = vslDeleteStream( &stream_good );
    RETURN_ON_ERROR;

    return 0;
}