• C

# Complex SVD algorithm

I'm implementing the SVD with complex class elements.
So far I've modified the svdcmp() from the book "Numerical Recipes for C" by Press et al. for complex elements as outlined in the LINPACK user manual.  The bidiagonal form now has real elements (w[i] and rv1[i], and its appropriate Householder transforms in U and V matrices).  But the accumulation of the real rotations (Diagonalization part) in the _complex_ matrices U and V does not give the right results.  "Numerical ..." does not give explanation of its svdcmp(), so I wonder if anyone can explain how the accummulation of the rotations are done in the U and V matrices.
Thanks a huge bunch.
PK

Looks like I've to upload the source.  Don't have the source right now; will upload it soon.  By the way, this question might be harder than its current worth.
Thanks again.
PK

Here's the _incomplete_ source:
1. I've modified it half-way only.
2. To use complex arithmetic in the Householder part is
beyond my capability.
3. Anyway, I've just implemented the SVD by breaking the
complex matrices down to real matrices and process them
just like normal.  Only thing is that the    number of operations have increased by two.  But I'm    working with small matrices, so it's okay.
4. Thanks for the help guys, how do I clear this question
anyway?

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <conio.h>
#define NR_END 1
#include <math.h>
#include <complex.h>    // I used the complex C++ class
// in BC 5.0
void nrerror(char error_text[])
{
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(1);
}

complex *cpvector(long nl, long nh)
{
complex *v;
long n=nh-nl+1;
v=new complex[n+NR_END];
return v;
}

void free_cpvector(complex *v, long nl, long nh)
{
delete[] v;
}

double *dvector(long nl, long nh)
{
double *v;
long n=nh-nl+1;
v=new double[n+NR_END];
return v;
}

void free_dvector(double *v, long nl, long nh)
{
delete[] v;
}

complex **cpmatrix(long nrl, long nrh, long ncl, long nch)
{
long i, j, nrow=nrh-nrl+1, ncol=nch-ncl+1;
complex **m;

/* allocate pointers to rows */
m = new complex*[nrow+NR_END];
for (i = 0; i < nrow; i++)
m[i] = new complex[ncol+NR_END];
return m;
}

void free_cpmatrix(complex **m, long nrl, long nrh, long ncl, long nch)
{
int i, nrow = nrh-nrl;

for (int i = 0; i < nrow+1;  i++)
delete[] m[i];
delete[] m;
}

void svdcmp(complex **a, int m, int n, complex w[], complex **v)
{
double pythag(complex a, complex b);
int flag, i, its, j, jj, k, l, nm;
double anorm;
complex c, f, g, h, s, scale, x, y, z, *rv1;
double  C, F, G, H, S, X, Y, Z, *W, *RV1;
rv1=cpvector(1,n);
RV1=dvector(1,n);
W=dvector(1,n);
B=cpmatrix(1,m,1,n);
g=scale=complex(0.0, 0.0);
anorm=0.0;

for (i=1;i<=n;i++) {
l=i+1;
rv1[i]=scale*g;
g=s=scale=complex(0.0, 0.0);
if (i <= m) {
for(k=i;k<=m;k++) scale += complex(abs(a[k][i]),0);
if (scale != complex(0.0,0.0)) {
for (k=i;k<=m;k++) {
a[k][i] /= scale;
s += a[k][i]*a[k][i];
}
f=a[i][i];
//g=-SIGN(sqrt(s),f);
if (f != complex(0.0,0.0))
g=-sqrt(s)*f/abs(f);
else g=-sqrt(s);
h=f*g-s;
a[i][i]=f-g;
for(j=l;j<=n;j++) {
for(s=complex(0.0,0.0),k=i;k<=m;k++) s += a[k][i]*a[k][j];
f=s/h;
for(k=i;k<=m;k++) a[k][j] += f*a[k][i];
}
for (k=i;k<=m;k++) a[k][i] *= scale;
}
}
w[i]=scale*g;
g=s=scale=complex(0.0,0.0);
if (i<=m && i!=n) {
for(k=l;k<=n;k++) scale += complex(abs(a[i][k]),0.0);
if (scale != complex(0.0,0.0)) {
for (k=l;k<=n;k++) {
a[i][k] /= scale;
s += a[i][k]*a[i][k];
}
f=a[i][l];
//g=-SIGN(sqrt(s),f);
if (f != complex(0.0,0.0))
g=-sqrt(s)*f/abs(f);
else g=-sqrt(s);
h=f*g-s;
a[i][l]=f-g;
for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
for (j=l;j<=m;j++) {
for(s=complex(0.0,0.0),k=l;k<=n;k++)
s +=  a[j][k]*a]i][k];
for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
}
for (k=l;k<=n;k++) a[i][k] *= scale;
}
}
anorm=FMAX(anorm, (abs(w[i])+abs(rv1[i])));
}

for(i=n;i>=1;i--){  // Accumulation of RH transformations
if (i<n) {
if (g != complex(0.0,0.0)) {
for(j=l;j<=n;j++)
v[j][i]=(a[i][j]/a[i][l])/g; //Double division
for(j=l;j<=n;j++) {
for(s=complex(0.0,0.0),k=l;k<=n;k++)
s += a[i][k]*v[k][j];
for(k=l;k<=n;k++) v[k][j] += s*v[k][i];
}
}
for(j=l;j<=n;j++) v[i][j]=v[j][i]=complex(0.0,0.0);
}
v[i][i]=complex(1.0,0.0);
g=rv1[i];
l=i;
}

for(i=IMIN(m,n);i>=1;i--){ //Accumulation of LH                                   transformations
l=i+1;
g=w[i];
for(j=l;j<=n;j++) a[i][j]=complex(0.0,0.0);
if (g != complex(0.0,0.0)) {
g=1.0/g;
for(j=l;j<=n;j++) {
for(s=complex(0.0,0.0),k=l;k<=m;k++)
s += a[k][i]*a[k][j];
f=(s/a[i][i])*g;
for(k=i;k<=m;k++) a[k][j] += f*a[k][i];
}
for(j=i;j<=m;j++) a[j][i] *= g;
} else for(j=i;j<=m;j++) a[j][i]=complex(0.0,0.0);
a[i][i]=1.0+a[i][i];
}

c=complex(1.0,0.0);
for(i=1;i<n;i++) {
if ((w[i] != complex(0.0,0.0))&&(rv1[i+1] !=                    complex(0.0,0.0))) {
for(k=1;k<=m;k++) {
a[k][i] *= (w[i]/abs(w[i]))*c;
v[k][i+1] *=                (rv1[i+1]/abs(rv1[i+1]))*(abs(w[i])/w[i])/c;
}
c=(rv1[i+1]/abs(rv1[i+1]))*(abs(w[i])/w[i])/c;
}
else if ((w[i] == complex(0.0,0.0))&&(rv1[i+1] !=                         complex(0.0,0.0))) {
for(k=1;k<=m;k++)
v[k][i+1] *= (rv1[i+1]/abs(rv1[i+1]));
c=(rv1[i+1]/abs(rv1[i+1]));
}
else if ((w[i] != complex(0.0,0.0))&&(rv1[i+1] ==                         complex(0.0,0.0))) {
for(k=1;k<=m;k++) a[k][i] *= (w[i]/abs(w[i]))*c;
c=complex(1.0,0.0);
}
else c=complex(1.0,0.0);
}
for(k=1;k<=m;k++) a[k][n] *= (w[n]/abs(w[n]))*c;

for(i=1;i<=n;i++) w[i] = complex(W[i]=abs(w[i]),0.0);
for(i=1;i<=n;i++)
rv1[i] = complex(RV1[i]=abs(rv1[i]),0.0);

printf("Matrix LH\n");
for (k=1;k<=m;k++) {
for (l=1;l<=n;l++)
printf("%12.6f %12.6f", real(a[k][l]),                                 imag(a[k][l]));
printf("\n");
}

printf("Matrix RH\n");
for (k=1;k<=m;k++) {
for (l=1;l<=n;l++)
printf("%12.6f %12.6f", real(v[k][l]),                                 imag(v[k][l]));
printf("\n");
}

for(k=n;k>=1;k--){  // Diagonalization of bidiagonal form
for(its=1;its<=30;its++){ // Loop over singular values
flag=1;
for(l=k;l>=1;l--) {  // Test for splitting
nm=l-1;           // Note rv1[1] is always zero
if ((double)(abs(RV1[l])+anorm) == anorm) {
flag=0;
break;
}
if ((double)(abs(W[nm])+anorm) == anorm) break;
}
if (flag) {
C=0.0;//c=complex(0.0,0.0);
S=0.0;//s=complex(1.0,0.0);
for(i=l;i<=k;i++) {
F=S*RV1[i];
RV1[i]=C*RV1[i];
if ((double)(abs(F)+anorm) == anorm) break;
G=W[i];
H=pythag(F,G);
W[i]=H;
H=1.0/H;
C=G*H;
S=-F*H;
for(j=1;j<=m;j++) {
y=a[j][nm]);
z=a[j][i]);
a[j][nm]=y*C+z*S;
a[j][i]=z*C-y*S;
}
}
}
Z=W[k];
if (l==k) {     // Convergence
if (Z<0.0) { // Singular value is made                                   non-negative
W[k] = -Z;
for(j=1;j<=n;j++) v[j][k] = -v[j][k];
}
break;
}
if (its == 30) nrerror("no convergence in 30 svdcmp                                      iterations");
X=W[l];       // Shift from bottom 2 by 2 minor
nm=k-1;
Y=W[nm];
G=RV1[nm];
H=RV1[k];
F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y);
G=pythag(F,1.0);
F=((X-Z)*(X+Z)+H*((Y/(F+SIGN(G,F)))-H))/X;
C=S=1.0;           // Next QR transformation
for (j=l;j<=nm;j++) {
i=j+1;
G=RV1[i];
Y=W[i];
H=S*G;
G=C*G;
Z=pythag(F,H);
RV1[j]=Z;
C=F/Z;
S=H/Z;
F=X*C+G*S;
G=G*C-X*S;
H=Y*S;
Y*=C;
for(jj=1;jj<=n;jj++) {
x=v[jj][j];
z=v[jj][i];
v[jj][j]=x*C+z*S;
v[jj][i]=z*C-x*S;
}
Z=pythag(F,H);
W[j]=Z;   // Rotation can be arbitrary if z=0
if (Z != 0.0) {
Z=1.0/Z;
C=F*Z;
S=H*Z;
}
F=C*G+S*Y;
X=C*Y-S*G;
for(jj=1;jj<=m;jj++) {
y=a[jj][j];
z=a[jj][i];
a[jj][j]=y*C+z*S;                                           a[jj][i]=z*C-y*S;
}
}
RV1[l]=0.0;
RV1[k]=F;
W[k]=X;
}
}
for(i=1;i<=n;i++) w[i] = complex(W[i],0.0);
for(i=1;i<=n;i++) rv1[i] = complex(RV1[i],0.0);
free_dvector(W,1,n);
free_dvector(RV1,1,n);
free_cpvector(rv1,1,n);
}

double pythag(complex a, complex b)
{
double absa, absb;
absa=abs(a);
absb=abs(b);
if (absa>absb) return absa*sqrt(1.0+SQR(absb/absa));
else return (absb == 0.0 ? 0.0 :                 absb*sqrt(1.0+SQR(absa/absb)));
}

void svbksb(complex **u, complex w[], complex **v, int m, int n, complex b[], complex x[])
{
int jj,j,i;
complex s, *tmp;

tmp=cpvector(1,n);
for (j=1;j<=n;j++) {
s=0.0;
if (w[j] != complex(0.0,0.0)) {
for (i=1;i<=m;i++) s += u[i][j]*b[i];
s /= w[j];
}
tmp[j]=s;
}
for (j=1;j<=n;j++) {
s=complex(0.0,0.0);
for (jj=1;jj<=n;jj++) s += v[j][jj]*tmp[jj];
x[j]=s;
}
free_cpvector(tmp,1,n);
}

void main(void)
{
complex **A, **B, *W, **U, **V, *X, *C;
int j,k,l,m=3,n=3, sm=2, sn=3;;
double wmax,wmin;

W=cpvector(1,5);
X=cpvector(1,5);
C=cpvector(1,5);
A=cpmatrix(1,5,1,5);
B=cpmatrix(1,5,1,5);
U=cpmatrix(1,5,1,5);
V=cpmatrix(1,5,1,5);
A[1][1] = complex(1.0 , 0.0);
A[1][2] = complex(2.0 , 0.0);
A[1][3] = complex(3.0 , 0.0);
A[2][1] = complex(2.0 , 0.0);
A[2][2] = complex(2.0 , 0.0);
A[2][3] = complex(3.0 , 0.0);
A[3][1] = complex(3.0 , 0.0);
A[3][2] = complex(3.0 , 0.0);
A[3][3] = complex(3.0 , 0.0);
B[1][1] = complex(1.0 , 0.0);
B[2][1] = complex(1.0 , 0.0);
B[3][1] = complex(1.0 , 0.0);
B[1][2] = complex(1.0 , 0.0);
B[2][2] = complex(2.0 , 0.0);
B[3][2] = complex(3.0 , 0.0);
for(k=1;k<=m;k++)
for(l=1;l<=n;l++) U[k][l]=A[k][l];
if (n>m) {
for (k=m+1;k<=n;k++) {
for (l=1;l<=n;l++) {
A[k][l] = complex(0.0,0.0);
U[k][l] = complex(0.0,0.0);
}
}
m=n;
}
svdcmp(U,m,n,W,V);
printf("Decomposition Matrices:\n");
printf("Matrix A\n");
for (k=1;k<=m;k++) {
for (l=1;l<=n;l++)
printf("%12.6f %12.6f", real(A[k][l]),                                  imag(A[k][l]));
printf("\n");
}
printf("Matrix U\n");
for (k=1;k<=m;k++) {
for (l=1;l<=n;l++)
printf("%12.6f %12.6f", real(U[k][l]),                                  imag(U[k][l]));
printf("\n");
}
printf("Diagonal of matrix w\n");
for (k=1;k<=n;k++)
printf("%12.6f %12.6f\n", real(W[k]), imag(W[k]));
printf("Matrix v-transpose\n");
for (k=1;k<=n;k++) {
for (l=1;l<=n; l++)
printf("%12.6f %12.6f  ", real(V[l][k]),                                    imag(V[l][k]));
printf("\n");
}

wmax=0.0;
for (k=1;k<=sn;k++)
if (abs(W[k]) > wmax) wmax=abs(W[k]);
wmin=wmax*(1.0e-6);
for (k=1;k<=sn;k++)
if (abs(W[k]) < wmin) W[k]=complex(0.0,0.0);

for (l=1;l<=sm;l++) {
printf("\nVector number %2d\n",l);
for (k=1;k<=sn;k++) C[k]=B[k][l];
svbksb(U,W,V,sn,sn,C,X);
printf(" solution vector is:\n");
for (k=1;k<=sn;k++) printf("%12.6f %12.6f",real(X[k]),                                                 imag(X[k]));
printf("\n original right-hand side vector:\n");
for (k=1;k<=sn;k++) printf("%12.6f %12.6f",real(C[k]),                                                 imag(C[k]));
printf("\n (matrix)*(sol'n vector):\n");
for (k=1;k<=sn;k++) {
C[k]=0.0;
for (j=1;j<=sn;j++)
C[k] += A[k][j]*X[j];
}
for (k=1;k<=sn;k++) printf("%12.6f %12.6f",real(C[k]),                                                 imag(C[k]));
printf("\n");
}

free_cpmatrix(V,1,5,1,5);
free_cpmatrix(U,1,5,1,5);
free_cpmatrix(B,1,5,1,5);
free_cpmatrix(A,1,5,1,5);
free_cpvector(C,1,5);
free_cpvector(X,1,5);
free_cpvector(W,1,5);
}
###### Who is Participating?
I wear a lot of hats...

"The solutions and answers provided on Experts Exchange have been extremely helpful to me over the last few years. I wear a lot of hats - Developer, Database Administrator, Help Desk, etc., so I know a lot of things but not a lot about one thing. Experts Exchange gives me answers from people who do know a lot about one thing, in a easy to use platform." -Todd S.

Commented:
You need to include your source or make it available, along with an explanation of the SVD algorithm and what you  are trying to accomplish by modifying it.

Mike
0
Author Commented:
Edited text of question
0
Commented:
I assume you are talking about Householder QU-decomposition.
Why do you want to use Householder decomposition. If you want
to solve ordinary lineair equations you should use Gauss Pivoting, its much simpler and faster.
If however you still want to use it, for least squares or something you should buy a good book about mathematics or pay (much) more points for source. I have to warn you: the source will be C++
0
Author Commented:
Edited text of question
0
Author Commented:
Edited text of question
0
Commented:
I believe the only way to clear a question is to receive an answer (like this one) and give it a passing grade. If points are valuable to you then you can adjust the points to something small then edit your question to be easily answered. I guess your complex SVD algorithm question has been answered?  :-)
0

Experts Exchange Solution brought to you by