有人知道c++的稀疏SVD求解器吗?我的问题涉及一些条件恶劣的矩阵,这些矩阵可能将列/行归零。我的数据存储在uBLAS矩阵中,该矩阵是Harwell Boeing稀疏格式。
我在查找时遇到了一些问题
SVD求解器
- 可以对稀疏矩阵进行操作的SVD解算器。拉帕克好像做不到?我想要将稀疏矩阵传递给函数并输出稀疏矩阵
- 一种重新组合结果的方法。。。这样我就可以从x=b(A^-1(中读出xs。我希望这是x=(b((v.(d^-1(.(u^t((
我希望从GSL 中重新创建以下两个步骤
gsl_linalg_SV_decomp_jacobi (gsl_matrix * A, gsl_matrix * V, gsl_vector * S)
gsl_linalg_SV_solve (const gsl_matrix * U, const gsl_matrix * V, const gsl_vector * S, const gsl_vector * b, gsl_vector * x)
我也不知道如何用c++封装FORTRAN库。哪里/有任何PROPACK c/c++绑定吗?
编辑1:我在使用PROPACK时遇到了一些问题。PROPACK输出稀疏矩阵吗?它似乎将V输出为"V(LDV,KMAX(:DOUBLE PRECISION array"。这意味着它没有?
SVDLIBC是一个C库,部分支持Harwell Boeing格式。我对这个图书馆不熟悉,但从表面上看,它似乎符合你的要求。
您提到了PROPACK。Fortran是C兼容的,您只需要知道调用约定是如何工作的。我不确定,但我认为您想在PROPACK中调用的函数是dlansvd
(假设为双精度(,其文档如下:
subroutine dlansvd(jobu,jobv,m,n,k,kmax,aprod,U,ldu,Sigma,bnd,
c V,ldv,tolin,work,lwork,iwork,liwork,doption,ioption,info,
c dparm,iparm)
c DLANSVD: Compute the leading singular triplets of a large and
c sparse matrix by Lanczos bidiagonalization with partial
c reorthogonalization.
c
c Parameters:
c
c JOBU: CHARACTER*1. If JOBU.EQ.'Y' then compute the left singular vectors.
c JOBV: CHARACTER*1. If JOBV.EQ.'Y' then compute the right singular
c vectors.
c M: INTEGER. Number of rows of A.
c N: INTEGER. Number of columns of A.
c K: INTEGER. Number of desired singular triplets. K <= MIN(KMAX,M,N)
c KMAX: INTEGER. Maximal number of iterations = maximal dimension of
c the generated Krylov subspace.
c APROD: Subroutine defining the linear operator A.
c APROD should be of the form:
c
c SUBROUTINE DAPROD(TRANSA,M,N,X,Y,DPARM,IPARM)
c CHARACTER*1 TRANSA
c INTEGER M,N,IPARM(*)
c DOUBLE PRECISION X(*),Y(*),DPARM(*)
c
c If TRANSA.EQ.'N' then the function should compute the matrix-vector
c product Y = A * X.
c If TRANSA.EQ.'T' then the function should compute the matrix-vector
c product Y = A^T * X.
c The arrays IPARM and DPARM are a means to pass user supplied
c data to APROD without the use of common blocks.
c U(LDU,KMAX+1): DOUBLE PRECISION array. On return the first K columns of U
c will contain approximations to the left singular vectors
c corresponding to the K largest singular values of A.
c On entry the first column of U contains the starting vector
c for the Lanczos bidiagonalization. A random starting vector
c is used if U is zero.
c LDU: INTEGER. Leading dimension of the array U. LDU >= M.
c SIGMA(K): DOUBLE PRECISION array. On return Sigma contains approximation
c to the K largest singular values of A.
c BND(K) : DOUBLE PRECISION array. Error estimates on the computed
c singular values. The computed SIGMA(I) is within BND(I)
c of a singular value of A.
c V(LDV,KMAX): DOUBLE PRECISION array. On return the first K columns of V
c will contain approximations to the right singular vectors
c corresponding to the K largest singular values of A.
c LDV: INTEGER. Leading dimension of the array V. LDV >= N.
c TOLIN: DOUBLE PRECISION. Desired relative accuracy of computed singular
c values. The error of SIGMA(I) is approximately
c MAX( 16*EPS*SIGMA(1), TOLIN*SIGMA(I) )
c WORK(LWORK): DOUBLE PRECISION array. Workspace of dimension LWORK.
c LWORK: INTEGER. Dimension of WORK.
c If JOBU.EQ.'N' and JOBV.EQ.'N' then LWORK should be at least
c M + N + 9*KMAX + 2*KMAX**2 + 4 + MAX(M+N,4*KMAX+4).
c If JOBU.EQ.'Y' or JOBV.EQ.'Y' then LWORK should be at least
c M + N + 9*KMAX + 5*KMAX**2 + 4 +
c MAX(3*KMAX**2+4*KMAX+4, NB*MAX(M,N)), where NB>1 is a block
c size, which determines how large a fraction of the work in
c setting up the singular vectors is done using fast BLAS-3
c operation.
c IWORK: INTEGER array. Integer workspace of dimension LIWORK.
c LIWORK: INTEGER. Dimension of IWORK. Should be at least 8*KMAX if
c JOBU.EQ.'Y' or JOBV.EQ.'Y' and at least 2*KMAX+1 otherwise.
c DOPTION: DOUBLE PRECISION array. Parameters for LANBPRO.
c doption(1) = delta. Level of orthogonality to maintain among
c Lanczos vectors.
c doption(2) = eta. During reorthogonalization, all vectors with
c with components larger than eta along the latest Lanczos vector
c will be purged.
c doption(3) = anorm. Estimate of || A ||.
c IOPTION: INTEGER array. Parameters for LANBPRO.
c ioption(1) = CGS. If CGS.EQ.1 then reorthogonalization is done
c using iterated classical GRAM-SCHMIDT. IF CGS.EQ.0 then
c reorthogonalization is done using iterated modified Gram-Schmidt.
c ioption(2) = ELR. If ELR.EQ.1 then extended local orthogonality is
c enforced among u_{k}, u_{k+1} and v_{k} and v_{k+1} respectively.
c INFO: INTEGER.
c INFO = 0 : The K largest singular triplets were computed succesfully
c INFO = J>0, J<K: An invariant subspace of dimension J was found.
c INFO = -1 : K singular triplets did not converge within KMAX
c iterations.
c DPARM: DOUBLE PRECISION array. Array used for passing data to the APROD
c function.
c IPARM: INTEGER array. Array used for passing data to the APROD
c function.
c
c (C) Rasmus Munk Larsen, Stanford, 1999, 2004
c
在Fortran中,需要记住的重要事项是,所有参数都通过引用传递,并且非稀疏数组以列主格式存储。因此,这个函数在C++中的正确声明应该如下(未经测试(:
extern "C"
void dlansvd(const char *jobu,
const char *jobv,
int *m,
int *n,
int *k,
int *kmax,
void (*aprod)(const char *transa,
int *m,
int *n,
int *iparm,
double *x,
double *y,
double *dparm),
double *U,
int *ldu,
double *Sigma,
double *bnd,
double *V,
int *ldv,
double *tolin,
double *work,
int *lwork,
int *iwork,
int *liwork,
double *doption,
int *ioption,
int *info,
double *dparm,
int *iparm);
这真是一头野兽。祝你好运
不妨看看Tim Davis的稀疏线性代数软件:http://www.cise.ufl.edu/~戴维斯/
总的来说,我发现他的软件非常有用,通常非常高效和健壮。
他似乎一直在和一个学生一起做一个稀疏的SVD包,但我不确定这个项目处于什么阶段
希望这能有所帮助。