examples about CGESVD
1 post
• Page 1 of 1
examples about CGESVD
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);
}
- warpv7aa
- Posts: 1
- Joined: Tue Feb 09, 2016 10:23 pm
1 post
• Page 1 of 1
Who is online
Users browsing this forum: No registered users and 2 guests