vsldcorrelationmatrixcol.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 correlation matrix  Example Program Text
!******************************************************************************/


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

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

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


int generate_input(int var, int n, double *x, double *a, double *C);

int main()
{
    int i, j, k1, k2, incx=1;
    int fail=0;
    int dim=DIM, n=N, dd=DIM*DIM, storage=VSL_SS_MATRIX_COLUMNS_STORAGE;
    VSLSSTaskPtr task=0;
    double x[N][DIM];  /* matrix of observations */

    double cov_target[DIM][DIM];

    double cov_f[DIM][DIM], cor_f[DIM][DIM];
    double cov_u[DIM*(DIM+1)/2], cor_u[DIM*(DIM+1)/2];
    double cov_l[DIM*(DIM+1)/2], cor_l[DIM*(DIM+1)/2];

    double check[DIM][DIM], norm;

    double mean[DIM], mean_target[DIM];
    MKL_INT cov_storage = VSL_SS_MATRIX_L_PACKED_STORAGE; // VSL_SS_MATRIX_FULL_STORAGE;
    MKL_INT cor_storage = VSL_SS_MATRIX_L_PACKED_STORAGE; // VSL_SS_MATRIX_FULL_STORAGE;
    MKL_INT errcode;    

    generate_input(DIM, N, (double *)x, mean_target, (double*)cov_target);

    for(i=0;i<DIM;i++)
    {
        for(j=0;j<DIM;j++)
        {
            cov_f[i][j]=0;
            cor_f[i][j]=0;
        }
    }

    for(i=0;i<DIM*(DIM+1)/2;i++)
    {
        cov_u[i]=0.0;
        cor_u[i]=0.0;
        cov_l[i]=0.0;
        cor_l[i]=0.0;
    }

    errcode = vsldSSNewTask( &task, &dim, &n, &storage, (double*)x, 0, 0 );

    cov_storage = VSL_SS_MATRIX_FULL_STORAGE;
    cor_storage = VSL_SS_MATRIX_FULL_STORAGE;

    errcode = vsldSSEditCovCor( task, mean, (double*)cov_f, &cov_storage , (double*)cor_f, &cor_storage );
    RETURN_ON_ERROR;

    for(i=0;i<DIM;i++) mean[i] = 0.0;

    errcode = vsldSSCompute( task, VSL_SS_COVARIANCE_MATRIX|VSL_SS_CORRELATION_MATRIX, VSL_SS_FAST_METHOD );
    RETURN_ON_ERROR;

    cov_storage = VSL_SS_MATRIX_U_PACKED_STORAGE;
    cor_storage = VSL_SS_MATRIX_U_PACKED_STORAGE;

    errcode = vsldSSEditCovCor( task, mean, (double*)cov_u, &cov_storage , (double*)cor_u, &cor_storage );
    RETURN_ON_ERROR;

    for(i=0;i<DIM;i++) mean[i] = 0.0;

    errcode = vsldSSCompute( task, VSL_SS_COVARIANCE_MATRIX|VSL_SS_CORRELATION_MATRIX, VSL_SS_FAST_METHOD );
    RETURN_ON_ERROR;

    cov_storage = VSL_SS_MATRIX_L_PACKED_STORAGE;
    cor_storage = VSL_SS_MATRIX_L_PACKED_STORAGE;

    errcode = vsldSSEditCovCor( task, mean, (double*)cov_l, &cov_storage , (double*)cor_l, &cor_storage );
    RETURN_ON_ERROR;

    for(i=0;i<DIM;i++) mean[i] = 0.0;

    errcode = vsldSSCompute( task, VSL_SS_COVARIANCE_MATRIX|VSL_SS_CORRELATION_MATRIX, VSL_SS_FAST_METHOD );
    RETURN_ON_ERROR;

    errcode = vslSSDeleteTask( &task );
    RETURN_ON_ERROR;


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

    printf(" Expected means\n");
    for(i=0;i<DIM;i++)
    {
        printf("%+lf ", mean_target[i]);
    }

    printf("\n");
    printf("\n");

    printf(" Expected covariance matrix       Expected correlation matrix\n");
    for(i=0;i<DIM;i++)
    {
        for(j=0;j<DIM;j++)
        {
            printf("%+lf ", cov_target[i][j]);
        }

        printf("   ");

        for(j=0;j<DIM;j++)
        {
            printf("%+lf ", 
                cov_target[i][j] /
                ((j!=i) ? sqrt(cov_target[i][i] * cov_target[j][j]) : 1.0));
        }
        printf("\n");
    }

    printf("\n");

    printf(" Computed means\n");
    for(i=0;i<DIM;i++)
    {
        printf("%+lf ", mean[i]);
    }

    printf("\n");
    printf("\n");

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

            check[i][j] = cov_f[i][j] - cov_target[i][j];
        }

        printf("   ");

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

    norm = dnrm2(&dd, (double*)check, &incx);
    if(norm > 1.0e-10) fail++;

    printf("\n");

    printf(" Packed (upper) covariance        Packed (upper) correlation\n");

    k1=0;
    k2=0;
    for(i=0;i<DIM;i++)
    {
        for(j=0;j<i;j++) printf("          ");

        for(j=i;j<DIM;j++)
        {
            printf("%+9lf ", cov_u[k1]);

            check[i][j] = cov_u[k1] - cov_target[i][j];
            check[j][i] = check[i][j];

            k1++;
        }

        printf("   ");

        for(j=0;j<i;j++) printf("          ");

        for(j=i;j<DIM;j++)
        {
            printf("%+9lf ", cor_u[k2]);
            k2++;
        }
        printf("\n");
    }

    norm = dnrm2(&dd, (double*)check, &incx);

    printf("\n");

    printf(" Packed (lower) covariance        Packed (lower) correlation\n");

    k1=0;
    k2=0;
    for(i=0;i<DIM;i++)
    {
        for(j=0;j<=i;j++)
        {
            printf("%+9lf ", cov_l[k1]);

            check[i][j] = cov_l[k1] - cov_target[i][j];
            check[j][i] = check[i][j];

            k1++;
        }

        for(j=i+1;j<DIM;j++) printf("          ");

        printf("   ");

        for(j=0;j<=i;j++)
        {
            printf("%+9lf ", cor_l[k2]);
            k2++;
        }
        printf("\n");
    }

    norm = dnrm2(&dd, (double*)check, &incx);

    printf("\nTEST PASSED\n");

    return 0; 
}

/* This function generates samples with known covarince and mean */
int generate_input(int var, int n, double *x, double *a, double *C)
{
    int i, k;

    for(k=0;k<var;k++)
    {
        for(i=0;i<n;i++)
        {
            x[i*var+k]=
                (double)(i+1)*(k+1) /
                ((double)n * (double)var); 
        }
    }

    for(k=0;k<var;k++)
    {
        for(i=0;i<var;i++)
        {
            C[i+k*var]=
                ((double)(i+1) * (k+1)) /
                ((double)n * var * var) *
                (n+1.0) / 12.0;
        }
    }

    for(k=0;k<var;k++)
    {
        a[k]= (double)(k+1) * (n+1) * 0.5 / ((double)var*n);
    }

    return 0;
}