Looking for SVD example.

[eeisemann]

Hi,

I try to calculate the kernel of am(rows) and n(columns) matrix. (n>m)

Therefore I wanted to decompose it with svd.

Unfortunately I do not really understand the description in the manual.

Has anyone a sample code?

Isn't it sufficient to call or sgebrd and then sorgbr?

The matrix P should then contain the kernel vectors in the last two columns.

Thank you very much.

I apologize asking this question, but I am trying to figure this out for quite some time now...

This is my sample code:

int rows=2;

int cols=3;

//matrix a little bigger to save result

float * a=new float[cols*cols];

a[0]=1; a[1]=0; a[2]=0;

a[3]=0; a[4]=1; a[5]=0;

//any values for the rest(not used for matrix)

a[6]=-2; a[7]=-2; a[8]=-2;

//create values as described in manual

//first dimension of a

int lda= rows;

//diagonal B

float * d=new float [min(rows,cols)];

//details of Q and P

float * tauq= new float[min(rows,cols)];

float * taup= new float[min(rows,cols)];

//off diagonal of B

float * e=new float[max(1,min(rows,cols)-1)];

//temp

float * work=new float[64*(rows+cols)];

int worksize=64*(rows+cols);

int info;

//void sgebrd(int *m,int *n,float *a,int *lda,float *d,float *e,float *tauq,float *taup,float *work,int *lwork,int *info);

sgebrd(&rows,&cols,a,&lda,d,e,tauq,taup,work,&worksize,&info);

if (info<0)

return info;

//lda=rows did not work either

lda=cols;

//read P

sorgbr ( "P", &cols, &cols, &rows, a, &lda, taup , work, &worksize, &info );

if (info<0)

return info;

//the last column is not (0,0,1)

[Sergey_K_Intel1]

There is a very good site
http://www.netlib.org/lapack/
where you can find samples and all the information.

In fact, your code doesn't solve your problem and the call

sorgbr ( "P", &cols, &cols, &rows, a, &lda, taup , work, &worksize, &info );

is wrong. This is because SORGBR generatesoneof the real orthogonal matrices Q or P**T determined by SGEBRD when reducing a real matrix A to bidiagonal form: A = Q * B * P**T. Q and P**T are defined as products of
elementary reflectors H(i) or G(i) respectively. And if you call

sgebrd(&rows, &cols,....)

then the correct call to sorgbr MUST be the following

sorgbr ( "P", &rows, &cols, &k,...

where k< = rows.

It's qite important that you will not obtain the kernel vectors by you using sorgbr becausem < n and according to the sorgbr manual

" If VECT = 'P', A is assumed to have been a K-by-N matrix, and P**T
is of order N:
if k < n, P**T = G(k) . . . G(2) G(1) and SORGBR returns the first m
rows of P**T, where n >= m >= k. "

This means you can obtain the first m rows of the matrix P**T but the whole P**T can not be obtained, and these m rows form the basis which is the orthogonal compliment to the kernel vectors if rank A=m.
So your code doesn't produce the kernel vectors.

There exist several ways to compute the kernel vectors using LAPACK routines. Let's us give an example. For example, it can be done through the LQ decomposition of the initial matrix A. If we have m by n matrix A andm < n, rank A = m, then Q**T Z give us the kernel vectors where Z is cols by n-m matrix which non-zero elements are defined as follows
a(j+m, j)=1. j=1,... n-m

Then A Z= L Q * (Q**T *Z)= 0. So if you replacethecalls to sgebrd and sorgbr with

vt[0]=0.;
vt[1]=0.;
vt[2]=1.;

sgelqf(&rows, &cols, a, &rows, taup, work, &worksize,&info);

sormlq("L", "T", &rows, &cols, &rows, a, &rows, taup, vt, &cols, work,&worksize,&info);

for (i=0; i<3; i++) {
printf("%f ", vt);
};

you obtain what you are looking for.