Solved

Complex SVD algorithm

Posted on 1997-03-09
6
2,273 Views
Last Modified: 2010-04-16
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:
Note before reading:
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);
}
0
Comment
Question by:PK030997
[X]
Welcome to Experts Exchange

Add your voice to the tech community where 5M+ people just like you are talking about what matters.

  • Help others & share knowledge
  • Earn cash & points
  • Learn & ask questions
6 Comments
 
LVL 2

Expert Comment

by:lockhart
ID: 1249553
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 Comment

by:PK030997
ID: 1249554
Edited text of question
0
 
LVL 3

Expert Comment

by:LucHoltkamp
ID: 1249555
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
Independent Software Vendors: We Want Your Opinion

We value your feedback.

Take our survey and automatically be enter to win anyone of the following:
Yeti Cooler, Amazon eGift Card, and Movie eGift Card!

 

Author Comment

by:PK030997
ID: 1249556
Edited text of question
0
 

Author Comment

by:PK030997
ID: 1249557
Edited text of question
0
 
LVL 1

Accepted Solution

by:
roy020697 earned 200 total points
ID: 1249558
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

Featured Post

Technology Partners: We Want Your Opinion!

We value your feedback.

Take our survey and automatically be enter to win anyone of the following:
Yeti Cooler, Amazon eGift Card, and Movie eGift Card!

Question has a verified solution.

If you are experiencing a similar issue, please ask a related question

Suggested Solutions

Title # Comments Views Activity
undefined reference to `bswap_128' 9 186
voltage to force translation ? 8 100
NEED HELP WITH VISUAL STUDIO 2017 (beginner) 6 83
Test the speeds on my PC Drives 12 71
Have you thought about creating an iPhone application (app), but didn't even know where to get started? Here's how: ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ Important pre-programming comments: I’ve never tri…
This is a short and sweet, but (hopefully) to the point article. There seems to be some fundamental misunderstanding about the function prototype for the "main" function in C and C++, more specifically what type this function should return. I see so…
The goal of this video is to provide viewers with basic examples to understand how to use strings and some functions related to them in the C programming language.
Video by: Grant
The goal of this video is to provide viewers with basic examples to understand and use for-loops in the C programming language.

740 members asked questions and received personalized solutions in the past 7 days.

Join the community of 500,000 technology professionals and ask your questions.

Join & Ask a Question