Index: /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/Makefile =================================================================== --- /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/Makefile (revision 26417) +++ /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/Makefile (revision 26418) @@ -69,4 +69,5 @@ $(SRC)/load.$(ARCH).o \ $(SRC)/lookup.$(ARCH).o \ +$(SRC)/matrix.$(ARCH).o \ $(SRC)/mkrgb.$(ARCH).o \ $(SRC)/mcreate.$(ARCH).o \ Index: /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/init.c =================================================================== --- /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/init.c (revision 26417) +++ /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/init.c (revision 26418) @@ -58,4 +58,5 @@ int load PROTO((int, char **)); int lookup PROTO((int, char **)); +int matrix PROTO((int, char **)); int mkrgb PROTO((int, char **)); int mcreate PROTO((int, char **)); @@ -195,4 +196,5 @@ {1, "minterp", minterp, "interpolate image pixels"}, {1, "iminterp", minterp, "interpolate image pixels"}, + {1, "matrix", matrix, "matrix math operations"}, {1, "mkrgb", mkrgb, "convert 3 images to rgb jpeg (use Kapa for better control)"}, {1, "mset", mset, "insert a vector in an image"}, Index: /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/matrix.c =================================================================== --- /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/matrix.c (revision 26418) +++ /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/matrix.c (revision 26418) @@ -0,0 +1,106 @@ +# include "data.h" + +int matrix (int argc, char **argv) { + + int ix, iy, ik, Nx, Ny, Nk; + float *Va, *Vb, *Vc, sum; + Buffer *bufA, *bufB, *bufC; + + if ((argc != 5) && (argc != 7)) { + gprint (GP_ERR, "USAGE: matrix set A = B * C : matrix multiplication\n"); + gprint (GP_ERR, "USAGE: matrix transpose A to B\n"); + return (FALSE); + } + + if (!strcasecmp (argv[1], "set")) { + + if ((argc != 7) || strcasecmp (argv[3], "=") || strcasecmp (argv[5], "*")) { + gprint (GP_ERR, "USAGE: matrix set A = B * C : matrix multiplication\n"); + return (FALSE); + } + if ((bufB = SelectBuffer (argv[4], OLDBUFFER, TRUE)) == NULL) return (FALSE); + if ((bufC = SelectBuffer (argv[6], OLDBUFFER, TRUE)) == NULL) return (FALSE); + + if (bufB[0].matrix.Naxis[0] != bufC[0].matrix.Naxis[1]) { + gprint (GP_ERR, "size mis-match in matrices: (%d x %d) * (%d x %d)\n", + bufB[0].matrix.Naxis[0], bufB[0].matrix.Naxis[1], + bufC[0].matrix.Naxis[0], bufC[0].matrix.Naxis[1]); + return (FALSE); + } + if ((bufA = SelectBuffer (argv[2], ANYBUFFER, TRUE)) == NULL) return (FALSE); + + Nx = bufC[0].matrix.Naxis[0]; + Ny = bufB[0].matrix.Naxis[1]; + Nk = bufB[0].matrix.Naxis[0]; + + gfits_free_matrix (&bufA[0].matrix); + gfits_free_header (&bufA[0].header); + + bufA[0].bitpix = bufB[0].bitpix; + bufA[0].unsign = bufB[0].unsign; + bufA[0].bscale = bufB[0].bscale; + bufA[0].bzero = bufB[0].bzero; + gfits_copy_header (&bufB[0].header, &bufA[0].header); + gfits_modify (&bufA[0].header, "NAXIS1", "%d", 1, Nx); + gfits_modify (&bufA[0].header, "NAXIS2", "%d", 1, Ny); + bufA[0].header.Naxis[0] = Nx; + bufA[0].header.Naxis[1] = Ny; + gfits_create_matrix (&bufA[0].header, &bufA[0].matrix); + + Va = (float *)bufA[0].matrix.buffer; + Vb = (float *)bufB[0].matrix.buffer; + Vc = (float *)bufC[0].matrix.buffer; + + for (iy = 0; iy < Ny; iy++) { + for (ix = 0; ix < Nx; ix++) { + sum = 0.0; + for (ik = 0; ik < Nk; ik++) { + sum += Vb[iy*Nk + ik] * Vc[ik*Nx + ix]; + } + Va[iy*Nx + ix] = sum; + } + } + return (TRUE); + } + + if (!strcasecmp (argv[1], "transpose")) { + if ((argc != 5) || strcasecmp (argv[3], "to")) { + gprint (GP_ERR, "USAGE: matrix transpose A to B\n"); + return (FALSE); + } + if ((bufA = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE); + if ((bufB = SelectBuffer (argv[4], ANYBUFFER, TRUE)) == NULL) return (FALSE); + + Nx = bufA[0].matrix.Naxis[0]; + Ny = bufA[0].matrix.Naxis[1]; + + gfits_free_matrix (&bufB[0].matrix); + gfits_free_header (&bufB[0].header); + + bufB[0].bitpix = bufA[0].bitpix; + bufB[0].unsign = bufA[0].unsign; + bufB[0].bscale = bufA[0].bscale; + bufB[0].bzero = bufA[0].bzero; + gfits_copy_header (&bufA[0].header, &bufB[0].header); + gfits_modify (&bufB[0].header, "NAXIS1", "%d", 1, Ny); + gfits_modify (&bufB[0].header, "NAXIS2", "%d", 1, Nx); + bufB[0].header.Naxis[0] = Ny; + bufB[0].header.Naxis[1] = Nx; + gfits_create_matrix (&bufB[0].header, &bufB[0].matrix); + + Va = (float *)bufA[0].matrix.buffer; + Vb = (float *)bufB[0].matrix.buffer; + + for (iy = 0; iy < Ny; iy++) { + for (ix = 0; ix < Nx; ix++) { + Vb[ix*Ny + iy] = Va[iy*Nx + ix]; + } + } + + return (TRUE); + } + + gprint (GP_ERR, "USAGE: matrix set A = B * C : matrix and/or vector multiplication\n"); + gprint (GP_ERR, "USAGE: matrix transpose A to B\n"); + return (FALSE); +} Index: /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/svd.c =================================================================== --- /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/svd.c (revision 26417) +++ /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/svd.c (revision 26418) @@ -57,4 +57,56 @@ /* 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) { Index: /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/test/matrix.sh =================================================================== --- /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/test/matrix.sh (revision 26418) +++ /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/test/matrix.sh (revision 26418) @@ -0,0 +1,78 @@ + +list tests + test1 + test2 +end + +# a very simple diagonal matrix equation +macro test1 + $PASS = 1 + break -auto off + + delete -q A B C + mcreate B 2 2 + mcreate C 2 2 + + B[0][0] = 2 + B[1][0] = -1 + B[0][1] = -1 + B[1][1] = 4 + + C[0][0] = 2 + C[1][0] = 1 + C[0][1] = 1 + C[1][1] = 1 + + matrix set A = B * C + + # A[0][0] = B[0][0]*C[0][0] + B[1][0]*C[0][1] = 4 - 1 = 3 + # A[1][0] = B[0][0]*C[1][0] + B[1][0]*C[1][1] = 2 - 1 = 1 + # A[0][1] = B[0][1]*C[0][0] + B[1][1]*C[0][1] = -2 + 4 = 2 + # A[1][1] = B[0][1]*C[1][0] + B[1][1]*C[1][1] = -1 + 4 = 3 + + if (abs(A[0][0] - 3.0) > 0.0001) + $PASS = 0 + end + if (abs(A[0][1] - 2.0) > 0.0001) + $PASS = 0 + end + if (abs(A[1][0] - 1.0) > 0.0001) + $PASS = 0 + end + if (abs(A[1][1] - 3.0) > 0.0001) + $PASS = 0 + end +end + +# a very simple off-diagonal matrix equation +macro test2 + $PASS = 1 + break -auto off + + delete -q A B + + mcreate A 2 2 + + A[0][0] = 2 + A[1][0] = +1 + A[0][1] = -1 + A[1][1] = 4 + + matrix transpose A to B + if (not($STATUS)) + $PASS = 0 + end + + if (abs(B[0][0] - 2.0) > 0.0001) + $PASS = 0 + end + if (abs(B[0][1] - 1.0) > 0.0001) + $PASS = 0 + end + if (abs(B[1][0] + 1.0) > 0.0001) + $PASS = 0 + end + if (abs(B[1][1] - 4.0) > 0.0001) + $PASS = 0 + end +end Index: /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/test/svd.sh =================================================================== --- /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/test/svd.sh (revision 26418) +++ /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/test/svd.sh (revision 26418) @@ -0,0 +1,91 @@ + +list tests + test1 + test2 +end + +# a very simple diagonal matrix equation +macro test1 + $PASS = 1 + break -auto off + + delete -q A U w Vt W A1 + + mcreate A 2 2 + + A[0][0] = 1.0 + A[1][0] = 0.5 + + A[0][1] = 0.5 + A[1][1] = 2.0 + + # svd give A = U w Vt, but this function returns V, not Vt: + svd A = U w V + if (not($STATUS)) + $PASS = 0 + end + + # can we recreate A? + mcreate W w[] w[] + for i 0 w[] + W[$i][$i] = w[$i] + end + + matrix transpose V to Vt + matrix set t1 = W * Vt + matrix set A1 = U * t1 + set dA = A1 - A + + stats -q dA + if ($MIN < -0.0001) + $PASS = 0 + end + if ($MAX > 0.0001) + $PASS = 0 + end +end + +# a very simple diagonal matrix equation +macro test2 + $PASS = 1 + break -auto off + + delete -q A U w Vt W A1 + + mcreate A 2 3 + + A[0][0] = 1.0 + A[1][0] = 0.0 + + A[0][1] = 0.0 + A[1][1] = 2.0 + + A[0][2] = 0.0 + A[1][2] = 1.0 + + # svd give A = U w Vt, but this function returns V, not Vt: + svd A = U w V + if (not($STATUS)) + $PASS = 0 + end + + # can we recreate A? + mcreate W w[] w[] + for i 0 w[] + W[$i][$i] = w[$i] + end + + matrix transpose V to Vt + matrix set t1 = W * Vt + matrix set A1 = U * t1 + set dA = A1 - A + + stats -q dA + if ($MIN < -0.0001) + $PASS = 0 + end + if ($MAX > 0.0001) + $PASS = 0 + end +end + Index: /branches/eam_branches/20091201/Ohana/src/opihi/doc/svdcmp_eispack.f =================================================================== --- /branches/eam_branches/20091201/Ohana/src/opihi/doc/svdcmp_eispack.f (revision 26418) +++ /branches/eam_branches/20091201/Ohana/src/opihi/doc/svdcmp_eispack.f (revision 26418) @@ -0,0 +1,358 @@ + subroutine svd(nm,m,n,a,w,matu,u,matv,v,ierr,rv1) +c + integer i,j,k,l,m,n,ii,i1,kk,k1,ll,l1,mn,nm,its,ierr + double precision a(nm,n),w(n),u(nm,n),v(nm,n),rv1(n) + double precision c,f,g,h,s,x,y,z,tst1,tst2,scale,pythag + logical matu,matv +c +c this subroutine is a translation of the algol procedure svd, +c num. math. 14, 403-420(1970) by golub and reinsch. +c handbook for auto. comp., vol ii-linear algebra, 134-151(1971). +c +c this subroutine determines the singular value decomposition +c t +c a=usv of a real m by n rectangular matrix. householder +c bidiagonalization and a variant of the qr algorithm are used. +c +c on input +c +c nm must be set to the row dimension of two-dimensional +c array parameters as declared in the calling program +c dimension statement. note that nm must be at least +c as large as the maximum of m and n. +c +c m is the number of rows of a (and u). +c +c n is the number of columns of a (and u) and the order of v. +c +c a contains the rectangular input matrix to be decomposed. +c +c matu should be set to .true. if the u matrix in the +c decomposition is desired, and to .false. otherwise. +c +c matv should be set to .true. if the v matrix in the +c decomposition is desired, and to .false. otherwise. +c +c on output +c +c a is unaltered (unless overwritten by u or v). +c +c w contains the n (non-negative) singular values of a (the +c diagonal elements of s). they are unordered. if an +c error exit is made, the singular values should be correct +c for indices ierr+1,ierr+2,...,n. +c +c u contains the matrix u (orthogonal column vectors) of the +c decomposition if matu has been set to .true. otherwise +c u is used as a temporary array. u may coincide with a. +c if an error exit is made, the columns of u corresponding +c to indices of correct singular values should be correct. +c +c v contains the matrix v (orthogonal) of the decomposition if +c matv has been set to .true. otherwise v is not referenced. +c v may also coincide with a if u is not needed. if an error +c exit is made, the columns of v corresponding to indices of +c correct singular values should be correct. +c +c ierr is set to +c zero for normal return, +c k if the k-th singular value has not been +c determined after 30 iterations. +c +c rv1 is a temporary storage array. +c +c calls pythag for dsqrt(a*a + b*b) . +c +c questions and comments should be directed to burton s. garbow, +c mathematics and computer science div, argonne national laboratory +c +c this version dated august 1983. +c +c ------------------------------------------------------------------ +c + ierr = 0 +c + do 100 i = 1, m +c + do 100 j = 1, n + u(i,j) = a(i,j) + 100 continue +c .......... householder reduction to bidiagonal form .......... + g = 0.0d0 + scale = 0.0d0 + x = 0.0d0 +c + do 300 i = 1, n + l = i + 1 + rv1(i) = scale * g + g = 0.0d0 + s = 0.0d0 + scale = 0.0d0 + if (i .gt. m) go to 210 +c + do 120 k = i, m + 120 scale = scale + dabs(u(k,i)) +c + if (scale .eq. 0.0d0) go to 210 +c + do 130 k = i, m + u(k,i) = u(k,i) / scale + s = s + u(k,i)**2 + 130 continue +c + f = u(i,i) + g = -dsign(dsqrt(s),f) + h = f * g - s + u(i,i) = f - g + if (i .eq. n) go to 190 +c + do 150 j = l, n + s = 0.0d0 +c + do 140 k = i, m + 140 s = s + u(k,i) * u(k,j) +c + f = s / h +c + do 150 k = i, m + u(k,j) = u(k,j) + f * u(k,i) + 150 continue +c + 190 do 200 k = i, m + 200 u(k,i) = scale * u(k,i) +c + 210 w(i) = scale * g + g = 0.0d0 + s = 0.0d0 + scale = 0.0d0 + if (i .gt. m .or. i .eq. n) go to 290 +c + do 220 k = l, n + 220 scale = scale + dabs(u(i,k)) +c + if (scale .eq. 0.0d0) go to 290 +c + do 230 k = l, n + u(i,k) = u(i,k) / scale + s = s + u(i,k)**2 + 230 continue +c + f = u(i,l) + g = -dsign(dsqrt(s),f) + h = f * g - s + u(i,l) = f - g +c + do 240 k = l, n + 240 rv1(k) = u(i,k) / h +c + if (i .eq. m) go to 270 +c + do 260 j = l, m + s = 0.0d0 +c + do 250 k = l, n + 250 s = s + u(j,k) * u(i,k) +c + do 260 k = l, n + u(j,k) = u(j,k) + s * rv1(k) + 260 continue +c + 270 do 280 k = l, n + 280 u(i,k) = scale * u(i,k) +c + 290 x = dmax1(x,dabs(w(i))+dabs(rv1(i))) + 300 continue +c .......... accumulation of right-hand transformations .......... + if (.not. matv) go to 410 +c .......... for i=n step -1 until 1 do -- .......... + do 400 ii = 1, n + i = n + 1 - ii + if (i .eq. n) go to 390 + if (g .eq. 0.0d0) go to 360 +c + do 320 j = l, n +c .......... double division avoids possible underflow .......... + 320 v(j,i) = (u(i,j) / u(i,l)) / g +c + do 350 j = l, n + s = 0.0d0 +c + do 340 k = l, n + 340 s = s + u(i,k) * v(k,j) +c + do 350 k = l, n + v(k,j) = v(k,j) + s * v(k,i) + 350 continue +c + 360 do 380 j = l, n + v(i,j) = 0.0d0 + v(j,i) = 0.0d0 + 380 continue +c + 390 v(i,i) = 1.0d0 + g = rv1(i) + l = i + 400 continue +c .......... accumulation of left-hand transformations .......... + 410 if (.not. matu) go to 510 +c ..........for i=min(m,n) step -1 until 1 do -- .......... + mn = n + if (m .lt. n) mn = m +c + do 500 ii = 1, mn + i = mn + 1 - ii + l = i + 1 + g = w(i) + if (i .eq. n) go to 430 +c + do 420 j = l, n + 420 u(i,j) = 0.0d0 +c + 430 if (g .eq. 0.0d0) go to 475 + if (i .eq. mn) go to 460 +c + do 450 j = l, n + s = 0.0d0 +c + do 440 k = l, m + 440 s = s + u(k,i) * u(k,j) +c .......... double division avoids possible underflow .......... + f = (s / u(i,i)) / g +c + do 450 k = i, m + u(k,j) = u(k,j) + f * u(k,i) + 450 continue +c + 460 do 470 j = i, m + 470 u(j,i) = u(j,i) / g +c + go to 490 +c + 475 do 480 j = i, m + 480 u(j,i) = 0.0d0 +c + 490 u(i,i) = u(i,i) + 1.0d0 + 500 continue +c .......... diagonalization of the bidiagonal form .......... + 510 tst1 = x +c .......... for k=n step -1 until 1 do -- .......... + do 700 kk = 1, n + k1 = n - kk + k = k1 + 1 + its = 0 +c .......... test for splitting. +c for l=k step -1 until 1 do -- .......... + 520 do 530 ll = 1, k + l1 = k - ll + l = l1 + 1 + tst2 = tst1 + dabs(rv1(l)) + if (tst2 .eq. tst1) go to 565 +c .......... rv1(1) is always zero, so there is no exit +c through the bottom of the loop .......... + tst2 = tst1 + dabs(w(l1)) + if (tst2 .eq. tst1) go to 540 + 530 continue +c .......... cancellation of rv1(l) if l greater than 1 .......... + 540 c = 0.0d0 + s = 1.0d0 +c + do 560 i = l, k + f = s * rv1(i) + rv1(i) = c * rv1(i) + tst2 = tst1 + dabs(f) + if (tst2 .eq. tst1) go to 565 + g = w(i) + h = pythag(f,g) + w(i) = h + c = g / h + s = -f / h + if (.not. matu) go to 560 +c + do 550 j = 1, m + y = u(j,l1) + z = u(j,i) + u(j,l1) = y * c + z * s + u(j,i) = -y * s + z * c + 550 continue +c + 560 continue +c .......... test for convergence .......... + 565 z = w(k) + if (l .eq. k) go to 650 +c .......... shift from bottom 2 by 2 minor .......... + if (its .eq. 30) go to 1000 + its = its + 1 + x = w(l) + y = w(k1) + g = rv1(k1) + h = rv1(k) + f = 0.5d0 * (((g + z) / h) * ((g - z) / y) + y / h - h / y) + g = pythag(f,1.0d0) + f = x - (z / x) * z + (h / x) * (y / (f + dsign(g,f)) - h) +c .......... next qr transformation .......... + c = 1.0d0 + s = 1.0d0 +c + do 600 i1 = l, k1 + i = i1 + 1 + g = rv1(i) + y = w(i) + h = s * g + g = c * g + z = pythag(f,h) + rv1(i1) = z + c = f / z + s = h / z + f = x * c + g * s + g = -x * s + g * c + h = y * s + y = y * c + if (.not. matv) go to 575 +c + do 570 j = 1, n + x = v(j,i1) + z = v(j,i) + v(j,i1) = x * c + z * s + v(j,i) = -x * s + z * c + 570 continue +c + 575 z = pythag(f,h) + w(i1) = z +c .......... rotation can be arbitrary if z is zero .......... + if (z .eq. 0.0d0) go to 580 + c = f / z + s = h / z + 580 f = c * g + s * y + x = -s * g + c * y + if (.not. matu) go to 600 +c + do 590 j = 1, m + y = u(j,i1) + z = u(j,i) + u(j,i1) = y * c + z * s + u(j,i) = -y * s + z * c + 590 continue +c + 600 continue +c + rv1(l) = 0.0d0 + rv1(k) = f + w(k) = x + go to 520 +c .......... convergence .......... + 650 if (z .ge. 0.0d0) go to 700 +c .......... w(k) is made non-negative .......... + w(k) = -z + if (.not. matv) go to 700 +c + do 690 j = 1, n + 690 v(j,k) = -v(j,k) +c + 700 continue +c + go to 1001 +c .......... set error -- no convergence to a +c singular value after 30 iterations .......... + 1000 ierr = k + 1001 return + end Index: /branches/eam_branches/20091201/Ohana/src/opihi/include/data.h =================================================================== --- /branches/eam_branches/20091201/Ohana/src/opihi/include/data.h (revision 26417) +++ /branches/eam_branches/20091201/Ohana/src/opihi/include/data.h (revision 26418) @@ -89,4 +89,7 @@ int svdcmp (float *a, opihi_flt *w, float *v, int Nx, int Ny); +/* in svdcmp_bond_raw.c */ +int svdcmp_bond_raw(int m, int n, int withu, int withv, double eps, double tol, double **a, double *q, double **u, double **v); + /* mrqmin.c */ opihi_flt mrqcof (opihi_flt *x, opihi_flt *y, opihi_flt *dy, int Npts,