矩阵填充中的SVT重构算法

/* * Stephen Becker, 11/10/08 * Computes A(omega), where A = U*V' is never explicitly computed * srbecker@caltech.edu * * Partial Changelog: * * 3/24/09 * adding in "mxIsComplex" and "mxGetPi" routines to handle complex data seamlessly * Not fully implemented yet, though; need to define .' vs ' convention (i.e. transpose vs adjoint) * 5/13/09 * Finished implementing the complex stuff. * The convention is A = U*V', not A = U*V.' * Having problems with complex data on 64-bit linux though: zdot isn't working * see * http://matlab.izmiran.ru/help/techdoc/matlab_external/ch04cr17.html * and the fort.h file in this directory, for alternative methods * 5/19/09 * try mex -largeArrayDims -lblas -DWINDOWS XonOmega.c ?? * it only affects zdotu, not the level-3 blas * tried -fPIC (as opposed to -fpic) in LDFLAGS, no luck * So, my current solution is to implement the blas function myself for 64-bit processors * * 10/29/09 * see XonOmegaTranspose.c for a much faster implementation * The idea is that we pass in U' and V' (instead of U and V) * Very little modification of the code, but benefit is enormous: * instead of stride lengths of M and N, stride lengths are now 1 * (because C and MATLAB use column-major order). * Also, adding in new level-3 BLAS feature: if omega isn't sorted, * and even if length(omega) isn't quite N*M, we still do level-3 BLAS * call, and THEN sort it or subsample it. Also, extend to omegaX,omegaY case * Also, why were omega, omegaI and omegaJ all doubles, not ints, or similar?? * Because MATLAB passes them in as doubles. Do NOT cast as integers!! * 11/9/09 * Consider including blas.h and using ptrdiff_t as the types, not ints * Added _WIN32 definition (implies WINDOWS) definition * Switched malloc and free to mxMalloc and mxFree * * * Note to the user: * to really take advantage of this function, you should recompile it on * your computer so that you can take advantage of any good BLAS libraries you might have * However, it is unlikely that this will be the speed bottleneck for most applications... * * There are 3 ways to call this program: * * (1) b = XonOmega(U,V,omega) * implicity forms the matrix A = U*V', then returns b = A(omega) * where omega are linear indices. * This is perhaps 20% faster if omega is sorted * * (2) b = XonOmega(U,V,omegaX,omegaY) * does the same, but omegaX and omegaY are subscripts * (i.e. omega = ind2sub(... OmegaX, OmegaY), or [OmegaX,OmegaY] = sub2ind(...,omega) ) * This is sometimes faster than method (1). If in addition, omegaX and omegaY * are created from the indices of a sorted linear index * (so that omegaY is sorted), then it will be about 20% faster * * (3) b = XonOmega(U,V,Y) * gets information about omega from the sparsity pattern of the sparse matrix Y * This is is actually the fastest method usually! * * (note: the relative speed of the different calling sequences * will depend on your computer) * * Update: all of the above methods to call the program can be modified by adding * an extra argument, cuttoff, to the end of the input sequence, e.g. * (1) b = XonOmega(U,V,omega,cutoff) * (2) b = XonOmega(U,V,omegaX,omegaY,cutoff) * (3) b = XonOmega(U,V,Y,cutoff) * where cutoff is a number between 0 and 1 that affects how the internal * computation is done; specifically, if length(omega) > cutoff*size(U,1)*size(V,1), * then the intermediate matrix A = U*V' is actually computed (using level-3 BLAS, * which is usually very optimized). The default value of cutoff is .25, * which is roughly the level where you begin to see benefit from not explicitly constructing A. * * Compilation notes: * For efficieny, this program makes BLAS calls whenever possible; * it uses level-3 BLAS if the entire matrix A = U*V' is needed, * and if only a small subset is needed, then it makes many level-1 BLAS calls. * * On Windows computers, BLAS functions have no underscore, but for linux/unix/darwin, they do * So, for Windows, either uncomment the line that says #define WINDOWS, or alternatively, * when compiling, use mex -L(location) -lmwblas -largeArrayDims -DWINDOWS XonOmega.c * (the -DWINDOWS is the same as putting "#define WINDOWS" in the code; * "location<S-Del>) is fullfile(matlabroot,'extern','lib','win32',cc) * where cc is usually 'lcc' (default) or 'microsoft' (for visual studio); * The location string should be enclosed in quotes) * ex. on my Windows machine, I do this mex -O -L"C:\Program Files\MATLAB\R2008b\extern\lib\win32\microsoft" ... -largeArrayDims -lmwblas -DWINDOWS XonOmega.c * * * * For linux, etc., you can specifically undefine the WINDOWS symbol * by passing in the -UWINDOWS option to mex, but that's not * necessary unless you uncommented the #define WINDOWS line below. * You can use the mathwork's blas, -lmwblas, or, recommended, * use a good blas that's already installed, and probably available * with just -lblas. So, try this: * mex -O -lblas -largeArrayDims updateSparse.c * the -O is for optimized * the -lblas may actually not even be necessary * * For both windows and linux, the -largeArrayDims flag is important if you have * a 64-bit system, otherwise you can skip it. * * If you've never compiled with mex before, please see the MATHWORKS help page. * You probably want to run "mex -setup" the first time. * * * Dec '09 notes: * zdotu -- if using libmwblas and matlab/extern/.../blas.h exists, use this for prototype and use blas call on 64-bit * also, use fortran_wrapper stuff * Source: zdotu.f #define zdotu FORTRAN_WRAPPER(zdotu) #ifndef FORTRAN_COMPLEX_FUNCTIONS_RETURN_VOID extern doublecomplex zdotu( #else extern void zdotu( doublecomplex* retval, #endif ptrdiff_t *n, double *zx, ptrdiff_t *incx, double *zy, ptrdiff_t *incy ); * extern complex16 zdotu_( int *K, complex16 *x, int *incx, complex16 *y, int *incy ); * */ #include "mex.h" /* #include "limits.h" in mex.h already */ /* #include "blas.h" */ /* might want to try */ /* #define WINDOWS */ /* I use the "WINDOWS" symbol to change between underscores and no underscores * _WIN32 *should* be automatically defined on WINDOWS machines, for most compilers */ #if defined(_WIN32) #define WINDOWS #endif /* to use "mwSize", need matrix.h, but doesn't work for R2006a */ #include "matrix.h" /* So, use the following definitions instead: */ #ifndef mwSize #define mwSize size_t #endif #ifndef mwIndex #define mwIndex size_t /* should make it compatible w/ 64-bit systems */ #endif /* R2009a: BLAS is now 64-bit on 64-bit systems, so use mwSignedIndex * For R2006a and earlier, when matrix.h doesn't exist, the following * definitions can be used: */ #ifndef mwSignedIndex #define mwSignedIndex ptrdiff_t #endif /* However, I'm not changing this just yet... * http://www.mathworks.com/access/helpdesk/help/techdoc/index.html?/access/helpdesk/help/techdoc/apiref/mwsignedindex.html&http://www.mathworks.com/access/helpdesk/help/techdoc/rn/brvak9c-1.html * *http://www.mathworks.com/access/helpdesk/help/techdoc/index.html?/access/helpdesk/help/techdoc/rn/brvak9c-1.html&http://www.mathworks.com/access/helpdesk/help/techdoc/rn/bryg9vd-1.html#br02jub-1 * */ #ifndef true #define true 1 #endif #ifndef false #define false 0 #endif typedef struct{ double re; double im; } complex16; /* a hack */ /* this is what mathworks does in lapack.h: */ /* #ifndef