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);
}
与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…