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

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 |
/* 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.

My name is Toto Sugito. This is my notes when I try something. Maybe, this is **NOT THE BEST SOLUTION** for your problems. If you have other method or idea that better with my post, please share in this blog.** Thank for visiting my site**.

Terima kasih atas informasi yang dipaparkan. Ini tentu semakin menambah pengetahuan saya tentang banyak hal

“BEYOND SVD”

“Advanced Eigenvalue Vector Decomposition”

We have made a significant improvement to the Singular Value

Decomposition methodology, This is actually an understatement.

We have discovered that the current Singular Value Decomposition mathematical techniques and resulting FORTRAN and C code is quite slow and inaccurate and what we have accomplished is to speed up computer execution by more than a factor of 1000 and to improve numerical accuracy by about 12 bits out of 48 bits for the standard double mantissa That is to say that there is no more than 1 bit different between the exact answer and my new solution method .Previous or current methods can barely achieve 24 bit accuracy (out of48).This improvement will make it possible to recognize red, green , blue, black, grey and white infrared images in real time as the first application .

The range of applications for this new technology can go well

beyond this.

Of course, we expect skeptics about these claims, but we have demo

programs that will read your data and produce results that they can compare and we can even executed this new code on their computers if they desire.

How AEVD Improves SVD Performance

AEVD Technology, LLC offers a fully developed, advanced form of the Singular Value Decomposition (SVD) algorithm, which offers a generational advance in speed and accuracy. Our Advanced Eigenvalue Vector Decomposition (or AEVD) was built upon a recognition of the shortcomings in how computers manipulate numbers, data and calculations, and reflects a painstaking analysis of the incremental steps in SVD processing to provide a faster process with fewer steps and improved inaccuracy.

The SVD mathematical proposition is linearly independent, in an algebraic sense, of the similarity theorem and as such provides a variety of available solution paths. One such path is to first reduce the input array to a real bidiagonal matrix with a sequence of intermediate left and right unitary transformations. This reduction to a real bidiagonal matrix is usually chosen to be a real diagonal and having one real super diagonal. All of the remaining matrix elements are numerically considered as zero. It is possible to choose other reduced forms of the input matrix, but the use of a real bidiagonal array provides for the most numerically accurate and computationally rapid solution. Additional numerical stability and computer accuracy is obtained by AEVD by choosing unitary transformations that place the larger bidiagonal elements in the upper left and the smaller elements in the lower right positions. This is true since the final determination of the left and right transformations and the SVD weights are always an iterative process for matrix sizes greater than four. Even for matrix sizes of four and less iterative methods are usually simpler and require computational steps comparable to closed form algebraic solutions. Also, when a real bidiagonal array format is chosen as the final iterate, the left and right transformations employed during iteration are orthogonal matrices. Other SVD iterative solution methods available employ orthogonal variations such as Jacobi rotation, Givens, and Householder reduction algorithms. Consistent among the more efficient and accurate SVD algorithms is the choice of the real

Regards,

John Rawlins

678-776-1343