LAPACK SVD(奇异值分解)

10

你知道如何使用LAPACK计算SVD的示例吗?

1个回答

14

例程dgesdd计算双精度矩阵的奇异值分解(SVD)。您只需要一个如何使用它的示例吗?您尝试过阅读文档吗?

以下是使用C LAPACK绑定的示例(请注意,我刚刚编写了此示例,并且尚未进行测试。还要注意,用于clapack参数的确切类型在平台之间有所不同,因此您可能需要将int更改为其他内容):

#include <clapack.h>

void SingularValueDecomposition(int m,     // number of rows in matrix
                                int n,     // number of columns in matrix
                                int lda,   // leading dimension of matrix
                                double *a) // pointer to top-left corner
{
    // Setup a buffer to hold the singular values:
    int numberOfSingularValues = m < n ? m : n;
    double *s = malloc(numberOfSingularValues * sizeof s[0]);

    // Setup buffers to hold the matrices U and Vt:
    double *u = malloc(m*m * sizeof u[0]);
    double *vt = malloc(n*n * sizeof vt[0]);

    // Workspace and status variables:
    double workSize;
    double *work = &workSize;
    int lwork = -1;
    int *iwork = malloc(8 * numberOfSingularValues * sizeof iwork[0]);
    int info = 0;

    // Call dgesdd_ with lwork = -1 to query optimal workspace size:
    dgesdd_("A", &m, &n, a, &lda, s, u, &m, vt, &n, work, &lwork, iwork, &info);
    if (info) // handle error conditions here

    // Optimal workspace size is returned in work[0].
    lwork = workSize;
    work = malloc(lwork * sizeof work[0]);

    // Call dgesdd_ to do the actual computation:
    dgesdd_("A", &m, &n, a, &lda, s, u, &m, vt, &n, work, &lwork, iwork, &info);
    if (info) // handle error conditions here

    // Cleanup workspace:
    free(work);
    free(iwork);

    // do something useful with U, S, Vt ...

    // and then clean them up too:
    free(s);
    free(u);
    free(vt);
}

请问您能否提供一个使用dgesdd的例子? - edgarmtze
@darkcminor:是Fortran还是C?再问一遍:你看过文档了吗? - Stephen Canon
1
还可以参考英特尔的示例:http://software.intel.com/sites/products/documentation/hpc/mkl/lapack/mkl_lapack_examples/dgesdd_ex.c.htm - Alec Jacobson
1
double *work = workSize 应该改为 double *work = &workSize - Warren Weckesser

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接