Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
493 views
in Technique[技术] by (71.8m points)

linear algebra - LAPACK SVD (Singular Value Decomposition)

Do yo know any example to use LAPACK To calculate SVD?

See Question&Answers more detail:os

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Reply

0 votes
by (71.8m points)

The routine dgesdd computes the SVD for a double precision matrix. Do you just need an example of how to use it? Have you tried reading the documentation?

An example using the C LAPACK bindings (note that I wrote this just now, and haven't actually tested it. Also note that the exact types for arguments to clapack vary somewhat between platforms so you may need to change int to something else):

#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);
}

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
OGeek|极客中国-欢迎来到极客的世界,一个免费开放的程序员编程交流平台!开放,进步,分享!让技术改变生活,让极客改变未来! Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

...