# svd code in c

This is my new post about how to create svd code in c with LAPACK routines. I get a problem in programming about gridding and I want to use svd code in c (I create svd code with LAPACK routines and call that from C) to solve this problem computation. In linear algebra, the singular value decomposition (SVD) is a factorization of a real or complex matrix, with many useful applications in signal processing and statistics (from wikipedia). We can use svd in Matlab with command [u, s, v] = svd(a). When you read in Matlab documentation, Matlab svd code uses the LAPACK routines to compute the svd (singular value decomposition). From this post, I try to create a simple svd code to solve svd computation with LAPACK routines. In my code, I use the last post about create 2d Array in C.

Singular value decomposition takes a rectangular matrix of gene expression data (defined as A, where A is a n x p matrix) in which the n rows represents the genes, and the p columns represents the experimental conditions. The SVD theorem states:

## A = U*S*VT

Result from this svd code computation is three matrix/array :

• U (Where the columns of U are the left singular vectors (gene coefficient vectors))
• S (S (the same dimensions as A) has singular values and is diagonal (mode amplitudes))
• VT has rows that are the right singular vectors (expression level vectors)

This is my simple svd code in c with LAPACK routine, you can save this code with name svdCLapack.c

```/* svdCLapack.c
* compute svd with lapack in C
* create by toto_plg@yahoo.com
* 08/03/2011
* http://toto-share.com
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

void svdlapack(double **a, int row, int col, double **u, double **s, double **v);
void printmatrix(char *var, double **a, int row, int col);
double *array1_create(int len);
double **array2_create(int row, int col);
void array1_free(double *arr);
void array2_free(double **arr);

int main(int argc, char **argv)
{
int row, col;
double idx;
int i, j;
double **a, **u, **s, **v;

row = 4;
col = 5;

a = array2_create(row, col);
u = array2_create(row, row);
s = array2_create(row, col);
v = array2_create(col, col);

idx= 1.0;
for(i=0; i<row; i++)
{
for(j=0; j<col; j++)
{
a[i][j] = idx;
idx = idx+1.0;
}
}

svdlapack(a, row, col, u, s, v);
printmatrix("A", a, row, col);
printmatrix("U", u, row, row);
printmatrix("S Transp", s, row, col);
printmatrix("V", v, col, col);

array2_free(a);
array2_free(u);
array2_free(v);
array2_free(s);

return(1);
}

void svdlapack(double **a, int row, int col, double **u, double **s, double **v)
{
/* Locals  */
int info, lwork;
double *tmpa, *tmpu, *tmpv, *tmps;
double wkopt;
double* work;
int i, j, ioff, iss;

tmpa = array1_create(row*col);
tmpu = array1_create(row*row);
tmpv = array1_create(col*col);
tmps = array1_create(col);

/* convert input to matrix 1D */
ioff = 0;
for(i=0; i<col; i++)
{
for(j=0; j<row; j++)
{
tmpa[ioff] = a[j][i];
ioff = ioff+1;
}
}

/* Query and allocate the optimal workspace */
lwork = -1;
dgesvd_( "All", "All", &row, &col, tmpa, &row, tmps, tmpu, &row, tmpv, &col, &wkopt,
&lwork, &info );

lwork = (int)wkopt;
work = (double*)malloc( lwork*sizeof(double) );

/* Compute SVD */
dgesvd_( "All", "All", &row, &col, tmpa, &row, tmps, tmpu, &row, tmpv, &col, work,
&lwork, &info );

/* Check for convergence */
if( info > 0 ) {
printf( "The algorithm computing SVD failed to converge.\n" );
exit( 1 );
}

/* Convert from tmpu (matrix 1D) to u (matrix 2D) */
ioff = 0;
for(i=0; i<row; i++)
{
for(j=0; j<row; j++)
{
u[j][i] = tmpu[ioff];
ioff = ioff+1;
}
}

/* Convert from tmpv (matrix 1D) to v (matrix 2D) */
ioff = 0;
for(i=0; i<col; i++)
{
for(j=0; j<col; j++)
{
v[j][i] = tmpv[ioff];
ioff = ioff+1;
}
}

/* get minimum size from row and columns */
if(row<col) iss = row;
else iss=col;

/* Convert from tmps (matrix 1D) to s (matrix 2D) */
for(i=0; i<iss; i++)
s[i][i] = tmps[i];

/* free allocated memory */
free( (void*)work );
array1_free(tmpa);
array1_free(tmpu);
array1_free(tmpv);
array1_free(tmps);
}

void printmatrix(char *var, double **a, int row, int col)
{
int i, j;
fprintf(stderr, "\nMatrix %s, row=%i col=%i \n", var, row, col);
for(i=0; i<row; i++)
{
for(j=0; j<col; j++)
{
fprintf(stderr, "%f   ", a[i][j]);
}
fprintf(stderr, "\n");
}
}

double *array1_create(int len)
{
double *arr=NULL;

/* allocate pointers to rows */
arr = (double *) malloc((len*sizeof(double)));
if (!arr) {
fprintf(stderr, "array1_create : error allocation array");
exit(0);
}
return(arr);
}

double **array2_create(int row, int col)
{
int i;
double **arr=NULL;

/* allocate pointers to rows */
arr = (double **) malloc((row*sizeof(double*)));
if (!arr) {
fprintf(stderr, "array2_create : error allocation row");
exit(0);
}

/* allocate rows and set pointers to them */
arr=(double*) malloc((row*col)*sizeof(double));
if (!arr) {
fprintf(stderr, "array2_create : error allocation column");
exit(0);
}

for(i=1; i<row; i++)
arr[i]=arr[i-1] + col;

/* return pointer to array of pointers to rows */
return arr;
}

void array1_free(double *arr)
{
free(arr);
}

void array2_free(double **arr)
{
free (arr);
free (arr);
}```

You can compile this svd code in c (Singular Value Decomposition code) with command :

gcc svdCLapack.c -lm -llapack -o svd

Maybe you want to ask, why output from this svd LAPACK is different with Matlab. I dont know why this happen. But I think Matlab modified this LAPACK code :). But if you want to check output from this svd LAPACK code is right or wrong, You can compute output from U, S, and V. From above equation, A=U*S*VT . So, You can compute output from this svd code in c and I get this output from this svd code in c is right.

Source:
http://en.wikipedia.org/wiki/Singular_value_decomposition
http://web.mit.edu/be.400/www/SVD/Singular_Value_Decomposition.htm
1. 2. 