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*V^T

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))
• V^T 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[0]=(double*) malloc((row*col)*sizeof(double));
	if (!arr[0]) {
		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[0]);
	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*V^T . So, You can compute output from this svd code in c and I get this output from this svd code in c is right.