Index: trunk/Ohana/src/opihi/cmd.astro/mkgauss.c =================================================================== --- trunk/Ohana/src/opihi/cmd.astro/mkgauss.c (revision 26884) +++ trunk/Ohana/src/opihi/cmd.astro/mkgauss.c (revision 26891) @@ -11,5 +11,6 @@ Buffer *buf; - Xo = Yo = 0; + // this should be Nx/2, Ny/2 if not set + Xo = Yo = NAN; if ((N = get_argument (argc, argv, "-c"))) { remove_argument (N, &argc, argv); @@ -30,4 +31,7 @@ Ny = buf[0].header.Naxis[1]; + if (isnan(Xo)) Xo = Nx / 2; + if (isnan(Yo)) Yo = Ny / 2; + /* gaussian parameters */ Sig_x = atof (argv[2]); Index: trunk/Ohana/src/opihi/cmd.astro/profile.c =================================================================== --- trunk/Ohana/src/opihi/cmd.astro/profile.c (revision 26884) +++ trunk/Ohana/src/opihi/cmd.astro/profile.c (revision 26891) @@ -31,12 +31,12 @@ ResetVector (yvec, OPIHI_FLT, (int)SQ(2*N+1)); - bzero (yvec[0].elements.Flt, (int)SQ(2*N+1)*sizeof(opihi_flt)+1); + bzero (yvec[0].elements.Flt, (int)SQ(2*N+1)*sizeof(opihi_flt)); V = (float *)(buf[0].matrix.buffer); Npt = 0; Nx = buf[0].matrix.Naxis[0]; - for (i = sx - N; i < sx + N; i++) { - for (j = sy - N; j < sy + N; j++, Npt++) { + for (i = sx - N; i <= sx + N; i++) { + for (j = sy - N; j <= sy + N; j++, Npt++) { yvec[0].elements.Flt[Npt] = V[i + j*Nx]; - xvec[0].elements.Flt[Npt] = hypot (i - sx, j - sy); + xvec[0].elements.Flt[Npt] = hypot (i + 0.5 - sx, j + 0.5 - sy); } } Index: trunk/Ohana/src/opihi/cmd.astro/star.c =================================================================== --- trunk/Ohana/src/opihi/cmd.astro/star.c (revision 26884) +++ trunk/Ohana/src/opihi/cmd.astro/star.c (revision 26891) @@ -24,4 +24,5 @@ if ((argc != 4) && (argc != 5)) { gprint (GP_ERR, "USAGE: star (buffer) x y [dx] [-border N] [-sat cnts]\n"); + gprint (GP_ERR, " dx is the aperture diameter, but is adjusted up to the next odd number\n"); return (FALSE); } Index: trunk/Ohana/src/opihi/cmd.basic/Makefile =================================================================== --- trunk/Ohana/src/opihi/cmd.basic/Makefile (revision 26884) +++ trunk/Ohana/src/opihi/cmd.basic/Makefile (revision 26891) @@ -39,4 +39,5 @@ $(SRC)/mkdir.$(ARCH).o \ $(SRC)/module.$(ARCH).o \ +$(SRC)/nop.$(ARCH).o \ $(SRC)/output.$(ARCH).o \ $(SRC)/quit.$(ARCH).o \ Index: trunk/Ohana/src/opihi/cmd.basic/init.c =================================================================== --- trunk/Ohana/src/opihi/cmd.basic/init.c (revision 26884) +++ trunk/Ohana/src/opihi/cmd.basic/init.c (revision 26891) @@ -25,4 +25,5 @@ int mkdir_opihi PROTO((int, char **)); int module PROTO((int, char **)); +int nop PROTO((int, char **)); int output PROTO((int, char **)); int pwd PROTO((int, char **)); @@ -66,4 +67,5 @@ {1, "?", list_help, "list commands *"}, {1, "??", list_vars, "list variables *"}, + {1, "#", nop, "a NOP function"}, {1, "local", local, "define local variables"}, {1, "macro", macro, "deal with the macros *"}, @@ -71,4 +73,5 @@ {1, "mkdir", mkdir_opihi, "built-in mkdir command"}, {1, "module", module, "load script file from the modules directories"}, + {1, "nop", nop, "a NOP function"}, {1, "output", output, "redirect output to file"}, {1, "pwd", pwd, "print current working directory"}, Index: trunk/Ohana/src/opihi/cmd.basic/nop.c =================================================================== --- trunk/Ohana/src/opihi/cmd.basic/nop.c (revision 26891) +++ trunk/Ohana/src/opihi/cmd.basic/nop.c (revision 26891) @@ -0,0 +1,5 @@ +# include "basic.h" + +int nop (int argc, char **argv) { + return (TRUE); +} Index: trunk/Ohana/src/opihi/cmd.data/Makefile =================================================================== --- trunk/Ohana/src/opihi/cmd.data/Makefile (revision 26884) +++ trunk/Ohana/src/opihi/cmd.data/Makefile (revision 26891) @@ -53,4 +53,6 @@ $(SRC)/ungridify.$(ARCH).o \ $(SRC)/histogram.$(ARCH).o \ +$(SRC)/hermitian1d.$(ARCH).o \ +$(SRC)/hermitian2d.$(ARCH).o \ $(SRC)/imcut.$(ARCH).o \ $(SRC)/imhist.$(ARCH).o \ @@ -69,4 +71,5 @@ $(SRC)/load.$(ARCH).o \ $(SRC)/lookup.$(ARCH).o \ +$(SRC)/matrix.$(ARCH).o \ $(SRC)/mkrgb.$(ARCH).o \ $(SRC)/mcreate.$(ARCH).o \ Index: trunk/Ohana/src/opihi/cmd.data/hermitian1d.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/hermitian1d.c (revision 26891) +++ trunk/Ohana/src/opihi/cmd.data/hermitian1d.c (revision 26891) @@ -0,0 +1,61 @@ +# include "data.h" + +int hermitian1d (int argc, char **argv) { + + int i, N, order, Nvec, Poly; + opihi_int *inI; + opihi_flt *out, *inF; + opihi_flt mean, sigma, x; + double nf, norm; + Vector *xvec, *yvec; + + Poly = FALSE; + if ((N = get_argument (argc, argv, "-poly"))) { + Poly = TRUE; + remove_argument (N, &argc, argv); + } + + if (argc != 6) goto usage; + + /* select input / output buffers */ + if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((yvec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE); + Nvec = xvec[0].Nelements; + ResetVector (yvec, OPIHI_FLT, Nvec); + + /* gaussian parameters */ + mean = atof (argv[3]); + sigma = atof (argv[4]); + order = atoi (argv[5]); + if (order < 0) goto usage; + if (order > 10) goto usage; + + out = yvec[0].elements.Flt; + inF = xvec[0].elements.Flt; + inI = xvec[0].elements.Int; + + nf = exp(lgamma(order + 1)); + norm = 1.0 / sqrt(nf*sqrt(2*M_PI)*sigma); + + // a little sub-optimal : split this up with macros? + for (i = 0; i < Nvec; i++, out++) { + if (xvec[0].type == OPIHI_FLT) { + x = (inF[i] - mean) / sigma; + } else { + x = (inI[i] - mean) / sigma; + } + *out = hermitian_polynomial (x, order); + if (!Poly) { + *out *= norm*exp (-0.25*x*x); + } + } + return (TRUE); + +usage: + gprint (GP_ERR, "USAGE: hermitian1d (x) (y) (mean) (sigma) (order)\n"); + gprint (GP_ERR, " Note : only orders up to 10 are supported\n"); + return (FALSE); + +} + +// {\psi}_n(x) = \frac{1}{\sqrt{n! \, 2^n\sqrt{\pi}}}\, \mathrm{e}^{-x^2/2}H_n(x).\,\! Index: trunk/Ohana/src/opihi/cmd.data/hermitian2d.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/hermitian2d.c (revision 26891) +++ trunk/Ohana/src/opihi/cmd.data/hermitian2d.c (revision 26891) @@ -0,0 +1,79 @@ +# include "astro.h" + +int hermitian2d (int argc, char **argv) { + + int i, j, N, Nx, Ny, Xorder, Yorder, Poly; + float *in; + double Xo, Yo, sigma; + double x, y, xf, yf, value, norm, weight; + double nf, mf; + Buffer *buf; + + // optionally return the Hermitian Polynomial or the Hermitian Function + Poly = FALSE; + if ((N = get_argument (argc, argv, "-poly"))) { + Poly = TRUE; + remove_argument (N, &argc, argv); + } + + Xo = Yo = NAN; + if ((N = get_argument (argc, argv, "-c"))) { + remove_argument (N, &argc, argv); + Xo = atof (argv[N]); + remove_argument (N, &argc, argv); + Yo = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + if (argc != 5) goto usage; + + /* select input / output buffers */ + if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); + Nx = buf[0].header.Naxis[0]; + Ny = buf[0].header.Naxis[1]; + + if (isnan(Xo)) Xo = Nx / 2; + if (isnan(Yo)) Yo = Ny / 2; + + /* gaussian parameters */ + sigma = atof (argv[2]); + Xorder = atof (argv[3]); + Yorder = atof (argv[4]); + if (Xorder < 0) goto usage; + if (Yorder < 0) goto usage; + if (Xorder > 10) goto usage; + if (Yorder > 10) goto usage; + + nf = exp(lgamma(Xorder + 1)); + mf = exp(lgamma(Yorder + 1)); + norm = 1.0 / sqrt(nf*mf*2*M_PI) / sigma; + + in = (float *) buf[0].matrix.buffer; + for (j = 0; j < Ny; j++) { + for (i = 0; i < Nx; i++, in++) { + + x = (i - Xo) / sigma; + y = (j - Yo) / sigma; + + // not sure what a 2D hermitian looks like + // is it H_i(x) * H_j(y)? + + xf = hermitian_polynomial (x, Xorder); + yf = hermitian_polynomial (y, Yorder); + value = xf*yf; + if (!Poly) { + weight = norm*exp(-0.25*(x*x + y*y)); + value *= weight; + } + + *in += value; + } + } + + return (TRUE); + +usage: + gprint (GP_ERR, "USAGE: hermitian2d (buffer) (sigma) Xorder Yorder\n"); + gprint (GP_ERR, " Note : only orders up to 10 are supported\n"); + return (FALSE); +} Index: trunk/Ohana/src/opihi/cmd.data/init.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/init.c (revision 26884) +++ trunk/Ohana/src/opihi/cmd.data/init.c (revision 26891) @@ -42,4 +42,6 @@ int ungridify PROTO((int, char **)); int histogram PROTO((int, char **)); +int hermitian1d PROTO((int, char **)); +int hermitian2d PROTO((int, char **)); int imcut PROTO((int, char **)); int imhist PROTO((int, char **)); @@ -58,4 +60,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 **)); @@ -172,4 +175,6 @@ {1, "header", header, "print image header"}, {1, "histogram", histogram, "generate histogram from vector"}, + {1, "hermitian1d", hermitian1d, "generate 1-D Hermitian Polynomial"}, + {1, "hermitian2d", hermitian2d, "generate 2-D Hermitian Polynomial"}, {1, "imbin", rebin, "rebin image data by factor of N"}, {1, "imclip", imclip, "clip values in an image to be within a range"}, @@ -195,4 +200,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: trunk/Ohana/src/opihi/cmd.data/matrix.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/matrix.c (revision 26891) +++ trunk/Ohana/src/opihi/cmd.data/matrix.c (revision 26891) @@ -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: trunk/Ohana/src/opihi/cmd.data/svd.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/svd.c (revision 26884) +++ trunk/Ohana/src/opihi/cmd.data/svd.c (revision 26891) @@ -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: trunk/Ohana/src/opihi/cmd.data/test/matrix.sh =================================================================== --- trunk/Ohana/src/opihi/cmd.data/test/matrix.sh (revision 26891) +++ trunk/Ohana/src/opihi/cmd.data/test/matrix.sh (revision 26891) @@ -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: trunk/Ohana/src/opihi/cmd.data/test/svd.sh =================================================================== --- trunk/Ohana/src/opihi/cmd.data/test/svd.sh (revision 26891) +++ trunk/Ohana/src/opihi/cmd.data/test/svd.sh (revision 26891) @@ -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: trunk/Ohana/src/opihi/doc/svdcmp_eispack.f =================================================================== --- trunk/Ohana/src/opihi/doc/svdcmp_eispack.f (revision 26891) +++ trunk/Ohana/src/opihi/doc/svdcmp_eispack.f (revision 26891) @@ -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: trunk/Ohana/src/opihi/include/data.h =================================================================== --- trunk/Ohana/src/opihi/include/data.h (revision 26884) +++ trunk/Ohana/src/opihi/include/data.h (revision 26891) @@ -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, @@ -162,3 +165,17 @@ void FreeBooks (void); +/* hermitian functions */ +double hermitian_polynomial (double x, int order); +double hermitian_00(double x); +double hermitian_01(double x); +double hermitian_02(double x); +double hermitian_03(double x); +double hermitian_04(double x); +double hermitian_05(double x); +double hermitian_06(double x); +double hermitian_07(double x); +double hermitian_08(double x); +double hermitian_09(double x); +double hermitian_10(double x); + # endif Index: trunk/Ohana/src/opihi/include/pclient.h =================================================================== --- trunk/Ohana/src/opihi/include/pclient.h (revision 26884) +++ trunk/Ohana/src/opihi/include/pclient.h (revision 26891) @@ -52,4 +52,6 @@ void FreePclient (void); void gotsignal (int signum); +void pipe_signal (int signum); +void pipe_signal_clear(void); # define DTIME(A,B) ((A.tv_sec - B.tv_sec) + 1e-6*(A.tv_usec - B.tv_usec)) Index: trunk/Ohana/src/opihi/lib.data/Makefile =================================================================== --- trunk/Ohana/src/opihi/lib.data/Makefile (revision 26884) +++ trunk/Ohana/src/opihi/lib.data/Makefile (revision 26891) @@ -20,4 +20,5 @@ $(SDIR)/fft.$(ARCH).o \ $(SDIR)/svdcmp.$(ARCH).o \ +$(SDIR)/svdcmp_bond_raw.$(ARCH).o \ $(SDIR)/convert.$(ARCH).o \ $(SDIR)/bracket.$(ARCH).o \ @@ -27,4 +28,5 @@ $(SDIR)/precess.$(ARCH).o \ $(SDIR)/starfuncs.$(ARCH).o \ +$(SDIR)/hermitian.$(ARCH).o \ $(SDIR)/gaussian.$(ARCH).o \ $(SDIR)/graphtools.$(ARCH).o \ Index: trunk/Ohana/src/opihi/lib.data/hermitian.c =================================================================== --- trunk/Ohana/src/opihi/lib.data/hermitian.c (revision 26891) +++ trunk/Ohana/src/opihi/lib.data/hermitian.c (revision 26891) @@ -0,0 +1,120 @@ +# include "data.h" + +double hermitian_polynomial (double x, int order) { + double value; + switch (order) { + case 0: + value = hermitian_00(x); + break; + case 1: + value = hermitian_01(x); + break; + case 2: + value = hermitian_02(x); + break; + case 3: + value = hermitian_03(x); + break; + case 4: + value = hermitian_04(x); + break; + case 5: + value = hermitian_05(x); + break; + case 6: + value = hermitian_06(x); + break; + case 7: + value = hermitian_07(x); + break; + case 8: + value = hermitian_08(x); + break; + case 9: + value = hermitian_09(x); + break; + case 10: + value = hermitian_10(x); + break; + default: + value = NAN; + break; + } + return value; +} + +double hermitian_00(double x) { + double value; + // H_0(x) = 1 + value = 1; + return value; +} +double hermitian_01(double x) { + double value; + // H_1(x) = x + value = x; + return value; +} +double hermitian_02(double x) { + double value, x2; + // H_2(x) = x^2-1 + x2 = x*x; + value = x2 - 1.0; + return value; +} +double hermitian_03(double x) { + double value, x2; + // H_3(x) = x^3-3x + x2 = x*x; + value = x*(x2 - 3.0); + return value; +} +double hermitian_04(double x) { + double value, x2; + // H_4(x) = x^4-6x^2+3 + x2 = x*x; + value = (x2 - 6.0)*x2 + 3.0; + return value; +} +double hermitian_05(double x) { + double value, x2; + // H_5(x) = x^5-10x^3+15x + x2 = x*x; + value = ((x2 - 10.0)*x2 + 15.0)*x; + return value; +} +double hermitian_06(double x) { + double value, x2; + // H_6(x) = x^6-15x^4+45x^2-15 + x2 = x*x; + value = (((x2 - 15.0)*x2 + 45.0)*x2) - 15.0; + return value; +} +double hermitian_07(double x) { + double value, x2; + // H_7(x) = x^7-21x^5+105x^3-105x + x2 = x*x; + value = (((x2 - 21.0)*x2+105.0)*x2 - 105.0)*x; + return value; +} +double hermitian_08(double x) { + double value, x2; + // H_8(x) = x^8-28x^6+210x^4-420x^2+105 + x2 = x*x; + value = ((((x2 - 28.0)*x2 + 210.0)*x2 - 420.0)*x2 + 105.0); + return value; +} +double hermitian_09(double x) { + double value, x2; + // H_9(x) = x^9-36x^7+378x^5-1260x^3+945x + x2 = x*x; + value = ((((x2 - 36.0)*x2 + 378.0)*x2 - 1260.0)*x2 + 945.0)*x; + return value; +} +double hermitian_10(double x) { + double value, x2; + // H_{10}(x) = x^{10}-45x^8+630x^6-3150x^4+4725x^2-945 + x2 = x*x; + value = (((((x2 - 45.0)*x2 + 630.0)*x2 - 3150.0)*x2 + 4725.0)*x2 - 945.0); + return value; +} Index: trunk/Ohana/src/opihi/lib.data/svdcmp_bond_raw.c =================================================================== --- trunk/Ohana/src/opihi/lib.data/svdcmp_bond_raw.c (revision 26891) +++ trunk/Ohana/src/opihi/lib.data/svdcmp_bond_raw.c (revision 26891) @@ -0,0 +1,244 @@ +/* svd.c -- Singular value decomposition. Translated to 'C' from the + * original Algol code in "Handbook for Automatic Computation, + * vol. II, Linear Algebra", Springer-Verlag. + * + * (C) 2000, C. Bond. All rights reserved. + * + * This is almost an exact translation from the original, except that + * an iteration counter is added to prevent stalls. This corresponds + * to similar changes in other translations. + * + * Returns an error code = 0, if no errors and 'k' if a failure to + * converge at the 'kth' singular value. + * + */ +#include /* for array allocation */ +#include /* for 'fabs' */ + +int svdcmp_bond_raw(int m, int n, int withu, int withv, double eps, double tol, double **a, double *q, double **u, double **v) +{ + int i,j,k,l,l1,iter,retval; + double c,f,g,h,s,x,y,z; + double *e; + + e = (double *)calloc(n,sizeof(double)); + retval = 0; + +/* Copy 'a' to 'u' */ + for (i=0;i x) + x = y; + } /* end i */ + +/* accumulation of right-hand transformations */ + if (withv) { + for (i=n-1;i>=0;i--) { + if (g != 0.0) { + h = u[i][i+1] * g; + for (j=l;j=0;i--) { + l = i + 1; + g = q[i]; + for (j=l;j=0;k--) { + iter = 0; +test_f_splitting: + for (l=k;l>=0;l--) { + if (fabs(e[l]) <= eps) goto test_f_convergence; + if (fabs(q[l-1]) <= eps) goto cancellation; + } /* end l */ + +/* cancellation of e[l] if l > 0 */ +cancellation: + c = 0.0; + s = 1.0; + l1 = l - 1; + for (i=l;i<=k;i++) { + f = s * e[i]; + e[i] *= c; + if (fabs(f) <= eps) goto test_f_convergence; + g = q[i]; + h = q[i] = sqrt(f*f + g*g); + c = g / h; + s = -f / h; + if (withu) { + for (j=0;j 30) { + retval = k; + break; + } + x = q[l]; + y = q[k-1]; + g = e[k-1]; + h = e[k]; + f = ((y-z)*(y+z) + (g-h)*(g+h)) / (2*h*y); + g = sqrt(f*f + 1.0); + f = ((x-z)*(x+z) + h*(y/((f<0)?(f-g):(f+g))-h))/x; +/* next QR transformation */ + c = s = 1.0; + for (i=l+1;i<=k;i++) { + g = e[i]; + y = q[i]; + h = s * g; + g *= c; + e[i-1] = z = sqrt(f*f+h*h); + c = f / z; + s = h / z; + f = x * c + g * s; + g = -x * s + g * c; + h = y * s; + y *= c; + if (withv) { + for (j=0;jtype == 's') { \ + if (V2->type == 'S') { \ opihi_flt M2 = V2[0].FltValue; \ for (i = 0; i < Nx*Ny; i++, out++, M1++) { \ @@ -635,5 +635,5 @@ break; \ } \ - if (V2->type == 'S') { \ + if (V2->type == 's') { \ opihi_int M2 = V2[0].IntValue; \ for (i = 0; i < Nx*Ny; i++, out++, M1++) { \ @@ -704,5 +704,5 @@ # define SM_FUNC(OP) { \ - if (V1->type == 's') { \ + if (V1->type == 'S') { \ opihi_flt M1 = V1[0].FltValue; \ for (i = 0; i < Nx*Ny; i++, out++, M2++) { \ @@ -711,5 +711,5 @@ break; \ } \ - if (V1->type == 'S') { \ + if (V1->type == 's') { \ opihi_int M1 = V1[0].IntValue; \ for (i = 0; i < Nx*Ny; i++, out++, M2++) { \ Index: trunk/Ohana/src/opihi/pclient/ChildOps.c =================================================================== --- trunk/Ohana/src/opihi/pclient/ChildOps.c (revision 26884) +++ trunk/Ohana/src/opihi/pclient/ChildOps.c (revision 26891) @@ -1,6 +1,6 @@ # include "pclient.h" -#include -#include -#include +# include +# include +# include // #include @@ -67,4 +67,5 @@ case 0: /* pipe is closed */ /** change child state? **/ + pipe_signal_clear(); break; default: /* data in pipe */ @@ -81,4 +82,5 @@ case 0: /* pipe is closed */ /** change child state? **/ + pipe_signal_clear(); break; default: /* data in pipe */ Index: trunk/Ohana/src/opihi/pclient/pclient.c.in =================================================================== --- trunk/Ohana/src/opihi/pclient/pclient.c.in (revision 26884) +++ trunk/Ohana/src/opihi/pclient/pclient.c.in (revision 26891) @@ -46,5 +46,15 @@ /* ignore the history file. to change this, see, eg, mana.c */ - signal (SIGPIPE, gotsignal); + + /* We do not want pclient to exit just because the connection to the + child process is broken. However, if the connection to the parent + process is broken, we should exit. To handle this, we define + pipe_signal() which tracks the number of times it has been called. + Every time an I/O operation on the child file descriptor raises this signal, + we need to manage the crashed child and then clear the count. Signals raised by + I/O to the parent do not result in clearing the count. If the count passes a maximum, + we exit pclient */ + + signal (SIGPIPE, pipe_signal); signal (SIGTSTP, gotsignal); signal (SIGTTIN, gotsignal); @@ -67,4 +77,19 @@ } +static int Npipe = 0; + +void pipe_signal (int signum) { + Npipe ++; + if (Npipe > 100) { + exit (5); + } + gprint (GP_ERR, "closed pipe : did child process exit unexpectedly?\n"); + return; +} + +void pipe_signal_clear(void) { + Npipe = 0; +} + void gotsignal (int signum) { gprint (GP_ERR, "got signal : %d\n", signum);