vslsgrouppooledcovariance.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 group/pooled covariance matrices  Example Program Text
!******************************************************************************/


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

#define RETURN_ON_ERROR                 \
    if(errcode<0)                        \
    {                                   \
        printf("Error: %i\n", errcode);  \
        printf("\nTEST FAILED\n");		\
        return 0;                       \
    }

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


#define DIM 3      /* dimension of the task */ 
#define N   1000   /* number of observations */
#define G   2      /* number of groups */  
#define GN  2      /* number of group covariance matrices */

int main()
{
    int i, j, dim=DIM, n=N, storage=VSL_SS_MATRIX_ROWS_STORAGE;
    VSLSSTaskPtr task=0;
    VSLStreamStatePtr stream;
    float a=0.0,sigma=1.0;
    float x[DIM][N];  /* matrix of observations */
    int group_indices[N];   /* indices of the groups */
    int group_matrix_indices[G]={1,1}; /* 1st and 2nd group matrix to be returned */

    float cov[DIM*DIM]; /* array to hold covariance */
    float mean[DIM]; /* array to hold covariance */

    float pooled_cov[DIM*DIM]; /* array to hold pooled covariance */
    float pooled_mean[DIM];    /* array to hold pooled mean */ 

    float group_cov[DIM*DIM*GN]; /* array to hold group covariance matrices */
    float group_mean[DIM*GN];    /* array to hold group means */
    int errcode;

    MKL_INT cov_storage, pooled_cov_storage, group_cov_storage; 
    unsigned long long task_params;

    float groups_check[DIM*DIM];
    int fail=0, n_odd, n_even, incx=1;
    float C_n, C_m, C_nm;
    int dd = DIM*DIM;
    float norm;

    cov_storage = VSL_SS_MATRIX_FULL_STORAGE;
    pooled_cov_storage = VSL_SS_MATRIX_FULL_STORAGE;
    group_cov_storage = VSL_SS_MATRIX_FULL_STORAGE;

    for(i=0;i<DIM;i++)          mean[i]=0.0;
    for(i=0;i<DIM*DIM;i++)      cov[i]=0.0;
    for(i=0;i<DIM*DIM;i++)      pooled_cov[i]=0.0;
    for(i=0;i<DIM;i++)          pooled_mean[i]=0.0;
    for(i=0;i<DIM*DIM*GN;i++)   group_cov[i]=0.0;
    for(i=0;i<DIM*GN;i++)       group_mean[i]=0.0;

    /***** Initialize *****/
    errcode = vslNewStream( &stream, BRNG,  SEED );

    /***** Call RNG *****/
    errcode = vsRngGaussian( METHOD, stream, N*DIM, (float*)x, a, sigma );

    /***** Deinitialize *****/
    errcode = vslDeleteStream( &stream );

    /* Dividing elements into odd and even */
    for ( i = 0; i < N;   i++) group_indices[i] = (i%2); 

    errcode = vslsSSNewTask( &task, &dim, &n, &storage, (float*)x, 0, 0 );
    RETURN_ON_ERROR;

    errcode = vsliSSEditTask( task, VSL_SS_POOLED_COV_MATRIX_STORAGE,
        &pooled_cov_storage );
    RETURN_ON_ERROR;

    errcode = vsliSSEditTask( task, VSL_SS_GROUP_COV_MATRIX_STORAGE,
        &group_cov_storage );
    RETURN_ON_ERROR;

    errcode = vslsSSEditCovCor( task, mean, (float*)cov, &cov_storage, NULL, NULL );
    RETURN_ON_ERROR;

    errcode = vslsSSEditPooledCovariance( task, group_indices,
        pooled_mean, pooled_cov, group_matrix_indices, group_mean, group_cov );
    RETURN_ON_ERROR;

    errcode = vslsSSCompute( task,
        VSL_SS_COVARIANCE_MATRIX | 
        VSL_SS_POOLED_COVARIANCE_MATRIX | VSL_SS_GROUP_COVARIANCE_MATRIX,
        VSL_SS_FAST_METHOD );
    RETURN_ON_ERROR;

    errcode = vslSSDeleteTask( &task );
    RETURN_ON_ERROR;



    printf("Covariance:\n");
    for(i=0;i<DIM;i++)
    {
        for(j=0;j<DIM;j++)
        {
            printf("%+lf ", cov[i*DIM+j]);
        }
        printf("\n");
    }

    printf("\n");

    printf("Group covariance matrices:\n");
    for(i=0;i<DIM;i++)
    {
        for(j=0;j<DIM;j++)
        {
            printf("%+lf ", group_cov[i*DIM+j]);
        }

        printf("     ");

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

    printf("\n");

    printf("Pooled covariance matrix:\n");
    for(i=0;i<DIM;i++)
    {
        for(j=0;j<DIM;j++)
        {
            printf("%+lf ", pooled_cov[i*DIM+j]);
        }
        printf("\n");
    }

    n_odd  = (N/2);         // n
    n_even = (N/2) + (N&1); // m

    C_n = (double)(n_odd*n_odd - n_odd) * n / ((double)n_odd * (n*n - n));

    C_m = (double)(n_even*n_even - n_even) * n / ((double)n_even * (n*n - n));

    C_nm = (double)n_even * n_odd / (double)(n*n-n);

    for(i=0;i<DIM;i++)
    {
        for(j=0;j<DIM;j++)
        {
            groups_check[i*DIM+j] =
                group_cov[i*DIM+j        ] * C_n +
                group_cov[i*DIM+j+DIM*DIM] * C_m +
                (group_mean[i] - group_mean[i+DIM]) *
                (group_mean[j] - group_mean[j+DIM]) * C_nm;

            groups_check[i*DIM+j] -= cov[i*DIM+j];
        }
    }

    norm = snrm2(&dd, (double*)groups_check, &incx);
    if(norm > 1.0e-5) fail++;

	printf("\nTEST PASSED\n");

    return 0;
}