vslsparametrizedmatrix_lu_pack.c

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

/*******************************************************************************
!                             INTEL CONFIDENTIAL
!  Copyright(C) 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:
!  Parametrization of correlation matrix  Example Program Text
!******************************************************************************/

#include <stdio.h>
#include "vsl_ss.h"

/* parameters of the task */
#define DIM 3      /* dimension of the task */

#define RETURN_ON_ERROR                 \
    if(status<0)                        \
    {                                   \
        printf("Error: %i\n", status);  \
        return 0;                       \
    }
/* Distorted correlation matrix */
static float Cin[DIM*(DIM+1)/2] = { 1.0,
                                      0.9, 1.0,
                                      0.7, 0.29, 1.0  };

int main()
{
    VSLSSTaskPtr task=0;

    MKL_INT i, j, k;
    MKL_INT dim = DIM;
    MKL_INT cor_storage, parametrized_cor_storage, lwork = 25*DIM;
    int status;

    float Cout[DIM*DIM];
    float B[DIM*DIM], Work[25*DIM];
    float EigenVect[DIM];
    char uplo, jobz;

    printf("\nDisturbed correlation matrix\n\n");
    k = 0;
    for( i = 0 ; i < DIM ; i++ )
    {
        for( j = 0; j <= i; j++)
        {
            printf("%+1.5f ", Cin[k]);
            k++;
        }
        for( j = 0; j < DIM-i; j++ )
        {
            printf("         ");
        }
        printf("\n");
    }
    printf("\n");
    
    for( j = 0 ; j < DIM*(DIM+1)/2; j++ )
    {
        B[j] = Cin[j];
    }

    uplo = 'U';
    jobz = 'N';
    sspev( &jobz, &uplo, &dim, B, EigenVect, 0, &dim, Work, &status);
    RETURN_ON_ERROR;

    printf("Eigenvalue of the disturbed correlation matrix\n\n");
    for( j = 0 ; j < DIM ; j++ )
    {
        printf("%+1.5f ", EigenVect[j]);
    }
    printf("\n");

    /* Create task */
    status = vslsSSNewTask( &task, &dim, 0, 0, 0, 0, 0 );
    RETURN_ON_ERROR;

    cor_storage = VSL_SS_MATRIX_L_PACKED_STORAGE;
    parametrized_cor_storage = VSL_SS_MATRIX_U_PACKED_STORAGE;

    /* Register arrays for parameterization of the correlation matrix */
    status = vslsSSEditMatrixParameterization( task,
             Cin, &cor_storage,
             Cout, &parametrized_cor_storage );
    RETURN_ON_ERROR;

    /* Execute task */
    status = vslsSSCompute( task, VSL_SS_MATRIX_PARAMETERIZATION,
                            VSL_SS_SPECTRAL_DECOMPOSITION_METHOD );
    if( status != 0 )
    {
        printf("Error/Warning: %i\n", status);
    }

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

    printf("\n\nParametrized correlation matrix\n\n");

    k = 0;
    for( i = 0 ; i < DIM ; i++ )
    {
        for( j = 0; j < i; j++ )
        {
            printf("         ");
        }
        for( j = 0; j < DIM-i; j++ )
        {
            printf("%+1.5f ", Cout[k]);
            k++;
        }
        printf("\n");
    }

    for( j = 0 ; j < DIM*(DIM+1)/2; j++ )
    {
        B[j] = Cout[j];
    }

    uplo = 'L';
    jobz = 'N';
    sspev( &jobz, &uplo, &dim, B, EigenVect, 0, &dim, Work, &status);
    RETURN_ON_ERROR;

    printf("\nEigenvalue of the parametrized correlation matrix\n\n");
    for( j = 0 ; j < DIM ; j++ )
    {
        printf("%+1.5f ", EigenVect[j]);
    }
    printf("\n");

    if( EigenVect[0] < -1e-6 )
    {
        printf("\nTEST FAILED\n");
    }
    else
    {
        printf("\nTEST PASSED\n");
    }

    return 0;
}