Changeset 26891 for trunk/Ohana/src/opihi/cmd.data/svd.c

Timestamp:

Feb 10, 2010, 7:24:46 PM (16 years ago)

Author:

eugene

Message:

updates from eam_branches/20091201

File:

• trunk/Ohana/src/opihi/cmd.data/svd.c (modified) (1 diff)

trunk/Ohana/src/opihi/cmd.data/svd.c

@@ -57 +57 @@
   /* use a bcopy instead? */
 
+  // try C.R. Bond's version -- requires matrices in the form A[row][col] not A[row*Ncol + col]
+  if (1) { 
+      int j;
+      double **a, **u, **v, *q;
+      ALLOCATE (a, double *, Ny);
+      ALLOCATE (u, double *, Ny);
+      ALLOCATE (v, double *, Ny);
+      for (i = 0; i < Ny; i++) {
+          ALLOCATE (a[i], double, Nx);
+          ALLOCATE (u[i], double, Nx);
+          ALLOCATE (v[i], double, Nx);
+      }   
+      ALLOCATE (q, double, Nx);
+
+      for (j = 0; j < Ny; j++) {
+          for (i = 0; i < Nx; i++) {
+              a[j][i] = A[j*Nx + i];
+              u[j][i] = 0;
+              v[j][i] = 0;
+          }
+      }
+      for (i = 0; i < Nx; i++) {
+          q[i] = 0;
+      }
+
+      status = svdcmp_bond_raw (Ny, Nx, 1, 1, FLT_EPSILON, 1e-6, a, q, u, v);
+      fprintf (stderr, "status: %d\n", status);
+
+      // copy u q v back to U W V:
+      for (j = 0; j < Ny; j++) {
+          for (i = 0; i < Nx; i++) {
+              U[j*Nx + i] = u[j][i];
+              V[j*Nx + i] = v[j][i];
+          }
+      }
+      for (i = 0; i < Nx; i++) {
+          W[i] = q[i];
+      }
+
+      for (j = 0; j < Ny; j++) {
+          free(a[j]);
+          free(u[j]);
+          free(v[j]);
+      }
+      free (a);
+      free (u);
+      free (v);
+      free (q);
+
+      return TRUE;
+  }
+
   status = svdcmp (U, W, V, Nx, Ny);
   if (!status) {