svd.zip_SVD_zip

#include <stdio.h> #include "export.h" #include <cula.h> #define min(a,b) (((a) < (b)) ? (a) : (b)) /* DLM for single precision SVD on CUDA, called from IDL as CUDA_SVD, A, W, U, VT where A is an M rows by N cols IDL array (input) S is a min(M,N) array of singular values (output) U is an N by min(M,N) array (output) VT is an min(M,N) by M array (output) such that A = U diag(S) VT This DLM encapsulates the CULA basic function culaSgesvd It is programmed, unnecessarily since there is no device code (yet), as a CUDA C project. NOTE: The input is an n by m IDL array representing AT in column major form as required by culaSgesvd. Since AT = V diag(S) UT CULA will return S, V, UT. Therefore, switching the last two, the output IDL arrays are S, UT, V which represent the matrices S, U VT. Mort Canty (2010) */ void IDL_CDECL cuda_svd(int argc, IDL_VPTR argv[]) { // output array pointers float * w, * u, * vt; // get the input matrix float * a = (float *) argv[0]->value.arr->data; // get its dimensions long ndim_a = argv[0]->value.arr->n_dim; long * dim_a = argv[0]->value.arr->dim; int m = (int) dim_a[0]; int n = (int) dim_a[1]; // IDL output arrays IDL_VPTR ivWptr; long ndim_w = 1; long dim_w[] = {min(m,n)}; w = (float * ) IDL_MakeTempArray( (int) IDL_TYP_FLOAT, ndim_w, dim_w, IDL_ARR_INI_ZERO, &ivWptr); IDL_VPTR ivUptr; long dim_u[] = {m,min(m,n)}; u = (float * ) IDL_MakeTempArray( (int) IDL_TYP_FLOAT, ndim_a, dim_u, IDL_ARR_INI_ZERO, &ivUptr); IDL_VPTR ivVTptr; long dim_v[] = {min(m,n),n}; vt = (float * ) IDL_MakeTempArray( (int) IDL_TYP_FLOAT, ndim_a, dim_v, IDL_ARR_INI_ZERO, &ivVTptr); // CULA general single precision SVD with host pointers culaStatus s = culaInitialize(); if(s == culaNoError) { s = culaSgesvd('S','S',m,n,a,m,w,u,m,vt,min(m,n)); culaShutdown(); } // return results to IDL (all zeroes if CULA failed to initialize) IDL_VarCopy(ivWptr,argv[1]); IDL_VarCopy(ivVTptr,argv[2]); IDL_VarCopy(ivUptr,argv[3]); } // the entry point, which loads the routine into IDL int IDL_Load(void) { static IDL_SYSFUN_DEF2 procedure_addr[] = { { (IDL_SYSRTN_GENERIC) cuda_svd, "CUDA_SVD", 4, 4, 0, 0 } }; return IDL_SysRtnAdd(procedure_addr, IDL_FALSE, 1); }