examples about CGESVD

[warpv7aa]

Is there any example for CGESVD ? My code is listed below, but it does not work. Can you shed any light on how to use this function ? thanks a lot.

Code: Select all

/*#########################################################################
*svd 0.1
*mex -g -O svd01.c -lwsock32 -lKernel32 -DWIN32 COMPFLAGS="/openmp $COMPFLAGS"  LDFLAGS="\$LDFLAGS /openmp" cula_lapack_basic.lib
#########################################################################*/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#include <assert.h>
#include <complex.h>
#include <mex.h>

#include <C:\Program Files\CULA\R17\include\cula_lapack.h>
#include <C:\Program Files\CULA\R17\include\cula.h>
//#include <lapack.h>

#define imin(X, Y)  ((X) < (Y) ? (X) : (Y))
#define MAX_ITER 200

void checkStatus(culaStatus status);

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    /* Declare all the necessary variables */
    /* Dimensions of matrices */
//     int M  = 2 ;
//     int N  = 16;

   int M  = (int)mxGetM(prhs[0]);
   int N  = (int)mxGetN(prhs[0]);
//    int M  = (int) mxGetScalar(prhs[0]);
//    int N  = (int) mxGetScalar(prhs[1]);
   
    culaStatus status;
   
    /* Setup SVD Parameters */
    int LDA = M;
    int LDU = M;
    int LDVT = N;
   
//     double* A = NULL;
//     double* S = NULL;
//     double* U = NULL;
//     double* VT = NULL;

//     doublecomplex * A = NULL;
//     doublecomplex * S = NULL;
//     doublecomplex * U = NULL;
//     doublecomplex * VT = NULL;
   


    time_t begin_time;
    time_t end_time;
    int cula_time;

    char jobu = 'A';
    char jobvt = 'A';

//     A = (doublecomplex*)malloc(M*N*sizeof(doublecomplex));
//     S = (doublecomplex*)malloc(imin(M,N)*sizeof(doublecomplex));
//     U = (doublecomplex*)malloc(LDU*M*sizeof(doublecomplex));
//     VT =(doublecomplex*)malloc(LDVT*N*sizeof(doublecomplex));
//
//     if(!A || !S || !U || !VT) /* Memory allocation failed */
//     {
//         free(A);
//         free(U);
//         free(S);
//         free(VT);
//
//         return ;//EXIT_FAILURE;
//     }

   
    double *samples_real = mxGetPr(prhs[0]);
    double *samples_imag = mxGetPi(prhs[0]);
//     A = (float*)mxGetData(prhs[2]);
    int index;
   int index_M = 0;
   int index_N = 0;
    for(index_M = 0; index_M < M; index_M = index_M + 1)
    {
        printf("row: %d--",index_M);
        for(index_N = 0; index_N < N; index_N = index_N + 1)
        {
            index = index_N*M+index_M ;
//             real(A[index])=samples_real[index];
//             imag(A[index])=samples_imag[index];
//             printf("A_real[%d]=%f\n",index,real(A[index]));
//             printf("A_imag[%d]=%f\n",index,imag(A[index]));
            printf("%f+i%f;\t",samples_real[index],samples_imag[index]);
        }
        printf("\n");
    }
   
   
//  [U,S,V] = svd01(A);
    plhs[0] = mxCreateDoubleMatrix( (mwSize)LDU, (mwSize)M, mxCOMPLEX);
    plhs[1] = mxCreateDoubleMatrix( (mwSize)imin(M,N), (mwSize)imin(M,N), mxREAL);
    plhs[2] = mxCreateDoubleMatrix( (mwSize)LDVT, (mwSize)N, mxCOMPLEX);
   
    /* Initialize CULA */
    status = culaInitialize();
    checkStatus(status);

    /* Perform singular value decomposition CULA */
    printf("Performing singular value decomposition using CULA ... ");

    culaFloatComplex * A = (culaFloatComplex*)prhs[0];
    culaFloatComplex * U = (culaFloatComplex*)plhs[0];
    culaFloat * S = (culaFloat*)plhs[1];
    culaFloatComplex * VT= (culaFloatComplex*)plhs[2];
    time(&begin_time);
//     status = culaCgesvd(jobu, jobvt, M, N, A, LDA, S, U, LDU, VT, LDVT);
    status = culaCgesvd(jobu, jobvt, M, N, A, LDA, S, U, LDU, VT, LDVT);

    checkStatus(status);
    time(&end_time);

//     abc;
   
    cula_time = (int)difftime( end_time, begin_time);
    printf("done. (%d seconds)\n", cula_time);

   
   
    double *output_array_U[2]    = {NULL, NULL};
    output_array_U[0]  =mxGetPr(plhs[0]);
    output_array_U[1]  =mxGetPi(plhs[0]);

    for(index_M = 0; index_M < M; index_M = index_M + 1)
    {
        printf("row: %d--",index_M);
        for(index_N = 0; index_N < M; index_N = index_N + 1)
        {
            index = index_N*M+index_M ;
//             real(A[index])=samples_real[index];
//             imag(A[index])=samples_imag[index];
//             printf("A_real[%d]=%f\n",index,real(A[index]));
//             printf("A_imag[%d]=%f\n",index,imag(A[index]));
            printf("%f+i%f;\t",output_array_U[0][index],output_array_U[1][index]);
        }
        printf("\n");
    }
   
   
    culaShutdown();

    /* Clean up workspace, input, and output memory allocations */
//     free(A);
//     free(U);
//     free(S);
//     free(VT);

//     return ;//EXIT_SUCCESS;
}

/* Check for errors and exit if one occurred */
void checkStatus(culaStatus status)
{
    char buf[256];

    if(!status)
        return;

    culaGetErrorInfoString(status, culaGetErrorInfo(), buf, sizeof(buf));
    printf("%s\n", buf);

    culaShutdown();
    exit(EXIT_FAILURE);
}