1

C ++のスパースSVDソルバーを知っている人はいますか? 私の問題には、列/行がゼロになっている可能性のある条件の悪い行列が含まれています。私のデータは、Harwell-Boeing スパース形式である uBLAS マトリックスに格納されています。

見つけるのに苦労しています：

SVD ソルバー

  1. スパース行列を操作できる SVD ソルバー。Lapackはこれを行うことができないようですか？関数に疎行列を渡し、疎行列を出力したいと考えています。
  2. 結果を再結合する方法... x=b(A^-1) から xs を読み取ることができるように。これは x=(b)(v.(d^-1).(u^t)) になると思います

GSLから次の2つのステップを再現したいと考えています

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(LDV,KMAX): DOUBLE PRECISION array」として V を出力するようです。そうでないことを意味するのはどれですか？

4

3 に答える 3

2

SVDLIBCは、Harwell-Boeing形式を部分的にサポートするCライブラリです。私はライブラリに精通していませんが、表面的にはあなたの要件に合っているようです。

于 2011-07-06T01:38:22.273 に答える
1

あなたは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);

なかなかの野獣です。幸運を！

于 2011-07-05T04:34:42.147 に答える
0

It might be worthwhile checking out Tim Davis's sparse linear algebra software: http://www.cise.ufl.edu/~davis/

Generally speaking I've found his software to be really useful, typically very efficient and robust.

It seems that he's been working on a sparse SVD package with a student, but I'm not sure what stage the project is at.

Hope this helps.

于 2011-07-05T01:34:37.227 に答える