Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

LAPACK SVD (Singular Value Decomposition)

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

like image 409
edgarmtze Avatar asked Feb 18 '11 23:02

edgarmtze


1 Answers

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);
}
like image 95
Stephen Canon Avatar answered Sep 25 '22 20:09

Stephen Canon