【问题标题】:SVD to solve harwell-boeing sparse a.x=b system in C/C++?SVD 在 C/C++ 中解决 harwell-boeing 稀疏 a.x=b 系统?
【发布时间】:2011-09-28 11:07:47
【问题描述】:

有人知道用于 c++ 的稀疏 SVD 求解器吗?我的问题涉及一些可能将列/行归零的条件不佳的矩阵。我的数据存储在 uBLAS 矩阵中,该矩阵是 Harwell-Boeing 稀疏格式。

我很难找到:

SVD 求解器

  1. 可以对稀疏矩阵进行运算的 SVD 求解器。 Lapack似乎无法做到这一点?我想将稀疏矩阵传递给函数并输出稀疏矩阵。
  2. 一种重新组合结果的方法...这样我就可以从 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”。这意味着它没有?

【问题讨论】:

    标签: c++ linear-algebra sparse-matrix numerical-methods svd


    【解决方案1】:

    SVDLIBC 是一个partial support for the Harwell-Boeing 格式的C 库。我对这个库不熟悉,但表面上它似乎符合您的要求。

    【讨论】:

    • 最接近我想要的东西。虽然我不知道它是否真的有效!
    【解决方案2】:

    你提到了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 中,要记住的重要一点是所有参数都通过引用传递,非稀疏数组以 column-major 格式存储。所以,这个函数在 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);
    

    真是一头野兽。祝你好运!

    【讨论】:

    • 我目前正在调查这条路线:虽然它根本不是我想去的方向:-)。有一个小错字:应该说 double *V, int *ldv... (v not u)
    • 这输出密集(非稀疏)数组?有没有支持稀疏数组的?我的矩阵 A 至少有 10000x10000+ 单位宽...
    【解决方案3】:

    也许值得看看 Tim Davis 的稀疏线性代数软件:http://www.cise.ufl.edu/~davis/

    总的来说,我发现他的软件非常有用,通常非常高效和强大。

    他似乎一直在和一个学生一起研究一个稀疏的 SVD 包,但我不确定该项目处于什么阶段。

    希望这会有所帮助。

    【讨论】:

    • @Misha:如果您点击他网站的链接,请找到“当前学生”并查看有关稀疏 SVD 的项目。根据我的快速浏览,如果您发送电子邮件,您似乎可以访问代码......
    猜你喜欢
    • 1970-01-01
    • 2017-03-26
    • 1970-01-01
    • 1970-01-01
    • 2017-06-26
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2018-10-07
    相关资源
    最近更新 更多