update

system/seismic/Mkirchinv.c

@@ -24,15 +24,17 @@
 #include <rsf.h>
 
 #include "kirchnew.h"
+#include "tfweight.h"
 #include "invert.h"
 
 int main(int argc, char* argv[])
 {
-    int n12, n1, n2, n3, i3, sw, niter;
+    int n12, n1, n2, n3, i3, sw, nw, niter;
     bool hd, ps;
-    float *data, *modl, *modl0, *vrms, *error=NULL, o1,d1,o2,d2;
+    float *data, *modl, *modl0, *vrms, *error=NULL, *ww, *ff;
+    float o1,d1,o2,d2;
     char *errfile;
-    sf_file in, out, vel, in0, err=NULL;
+    sf_file in, out, vel, in0, err=NULL, fwght, wght;
 
     sf_init (argc,argv);
     in = sf_input("in");
@@ -42,6 +44,25 @@ int main(int argc, char* argv[])
     if (!sf_histint(in,"n2",&n2)) sf_error("No n2= in input");
     n3 = sf_leftsize(in,2);
 
+    if (NULL != sf_getstring("fweight")) {
+	fwght = sf_input("fweight");
+	nw = kiss_fft_next_fast_size((n1+1)/2)+1;
+	ff = sf_floatalloc(nw*n2);
+	sf_floatread(ff,nw*n2,fwght);
+	sf_fileclose(fwght);
+
+	wght = sf_input("weight");
+	ww = sf_floatalloc(n1*n2);
+	sf_floatread(ww,n1*n2,wght);
+	sf_fileclose(wght);
+
+	tfweight_init(n1,nw,n2,ww,ff);
+    } else {
+	fwght = NULL;
+	wght = NULL;
+    }
+
+
     if (!sf_getbool("hd",&hd)) hd=true;
     /* if y, apply half-derivative filter */
     if (!sf_getbool("ps",&ps)) ps=true;
@@ -95,7 +116,12 @@ int main(int argc, char* argv[])
     for (i3=0; i3 < n3; i3++) {
 	sf_floatread (data,n12,in);
 
-	invert(kirchnew_lop,niter,niter,n12,n12,modl,modl0,data,error);
+	if (NULL != fwght) {
+	    invert(kirchnew_lop,tfweight_lop,
+		   niter,n12,n12,modl,modl0,data,error);
+	} else {
+	    invert(kirchnew_lop,NULL,niter,n12,n12,modl,modl0,data,error);
+	}
 
 	sf_floatwrite (modl,n12,out);
 	if (NULL != err) sf_floatwrite(error,niter,err);

system/seismic/Mkirchnew.c

@@ -40,7 +40,7 @@ int main(int argc, char* argv[])
     if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input");
     if (!sf_histint(in,"n2",&n2)) sf_error("No n2= in input");
     n3 = sf_leftsize(in,2);
-    
+
     if (!sf_getbool("adj",&adj)) adj=true;
     /* yes: migration, no: modeling */
     if (!sf_getbool("hd",&hd)) hd=true;

system/seismic/invert.c

@@ -19,9 +19,9 @@
 
 #include <rsf.h>
 
-void invert(sf_operator oper /* linear operator */, 
+void invert(sf_operator oper /* linear operator */,
+	    sf_operator prec /* preconditioning */,
 	    int niter        /* number of iterations */,
-	    int miter        /* memory */, 
 	    int nx, int ny   /* model and data size */,
 	    float *x         /* model */,
             float *x0        /* inital model */, 
@@ -32,10 +32,14 @@ void invert(sf_operator oper /* linear operator */,
     int iy, iter;
     float norm;
 
-    sf_cdstep_init();
-    sf_solver(oper,sf_cdstep,nx,ny,x,y,niter,
-	      "x0",x0,"nmem",0,"nfreq",miter,"err",error,"verb",true,"end");
-    sf_cdstep_close();
+    if (NULL != prec) {
+	sf_solver_prec(oper,sf_cgstep,prec,nx,nx,ny,x,y,niter,0.0f,
+		  "x0",x0,"err",error,"verb",true,"end");
+    } else {
+	sf_solver(oper,sf_cgstep,nx,ny,x,y,niter,
+		  "x0",x0,"err",error,"verb",true,"end");
+    }
+    sf_cgstep_close();
 
     norm = 0.;
     for (iy=0; iy < ny; iy++) {

system/seismic/tfweight.c

@@ -0,0 +1,99 @@
+/* time and frequency weights */
+/*
+  Copyright (C) 2004 University of Texas at Austin
+  
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2 of the License, or
+  (at your option) any later version.
+  
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+  
+  You should have received a copy of the GNU General Public License
+  along with this program; if not, write to the Free Software
+  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+*/
+
+#include <rsf.h>
+
+static int nt, nw, nf, nx;
+static float *w, *f, *tmp;
+static kiss_fftr_cfg forw, invs;
+static kiss_fft_cpx *cdata;
+
+void tfweight_init(int n1, int nw1, int n2,
+		   float* ww /* [n1*n2] time weight */,
+		   float* ff /* [nw1*n2] frequency weight */)
+/*< initialoze >*/
+{
+    nt = n1;
+    nw = nw1;
+    nx = n2;
+    w = ww;
+    f = ff;
+
+    nf = 2*(nw-1);
+
+    forw = kiss_fftr_alloc(nf,0,NULL,NULL);
+    invs = kiss_fftr_alloc(nf,1,NULL,NULL);
+
+    tmp = sf_floatalloc(nf);    
+    cdata =  (kiss_fft_cpx*) sf_alloc(nw,sizeof(kiss_fft_cpx));
+}
+
+void tfweight_close(void)
+/*< free allocated storage >*/
+{
+    free(tmp);
+    free(cdata);
+    free(forw);
+    free(invs);
+}
+
+void tfweight_lop(bool adj, bool add, int nxx, int nyy, float* x, float* y) 
+/*< linear operator >*/
+{
+    int iw, i1, i2;
+
+    if (nxx != nt*nx || nyy != nt*nx) sf_error("%s: Wrong size",__FILE__);
+
+    sf_adjnull(adj,add,nxx,nyy,x,y);
+
+    for (i2=0; i2 < nx; i2++) {
+	if (adj) {
+	    for (i1=0; i1 < nt; i1++) {
+		tmp[i1] = w[i1+i2*nt]*y[i1+i2*nt];
+	    }
+	    for (i1=nt; i1 < nf; i1++) {
+		tmp[i1] = 0.0f;
+	    }
+	    kiss_fftr(forw, tmp, cdata);
+	    for (iw=0; iw < nw; iw++) {
+		cdata[iw]=sf_crmul(cdata[iw],f[iw+i2*nw]/nf);
+	    }
+	    kiss_fftri(invs, cdata, tmp);
+	    for (i1=0; i1 < nt; i1++) {
+		x[i1+i2*nt] += tmp[i1];
+	    }
+	} else {
+	    for (i1=0; i1 < nt; i1++) {
+		tmp[i1] = x[i1+i2*nt];
+	    }
+	    for (i1=nt; i1 < nf; i1++) {
+		tmp[i1] = 0.0f;
+	    }
+	    kiss_fftr(forw, tmp, cdata);
+	    for (iw=0; iw < nw; iw++) {
+		cdata[iw]=sf_crmul(cdata[iw],f[iw+i2*nw]/nf);
+	    }
+	    kiss_fftri(invs, cdata, tmp);
+	    for (i1=0; i1 < nt; i1++) {
+		y[i1+i2*nt] += w[i1+i2*nt]*tmp[i1];
+	    }
+	}
+    }
+}
+

user/fomels/Mpchain.c

@@ -24,7 +24,8 @@
 
 int main(int argc, char* argv[])
 {
-    int i, ic, n, nc, n2, iter, niter, liter, rect;
+    bool verb;
+    int i, ic, n, nc, n2, iter, niter, liter, rect, it, nt;
     float f;
     sf_complex *xn, *x1, *y1, *dx, *r, *p;
     sf_file inp, out, pef;
@@ -36,10 +37,12 @@ int main(int argc, char* argv[])
 
     if (SF_COMPLEX != sf_gettype(inp)) sf_error("Need complex input");
     if (!sf_histint(inp,"n1",&n)) sf_error("No n1= in input");
+    nt = sf_leftsize(inp,1);
 
     if (!sf_getint("nc",&nc)) nc=1; /* number of components */
 
     sf_putint(pef,"n2",nc);
+    sf_shiftdim(inp, pef, 2);
 
     n2 = (2*nc-1)*n;
 
@@ -48,31 +51,14 @@ int main(int argc, char* argv[])
     p  = sf_complexalloc(n2);
 
     x1 = sf_complexalloc(n);
-    sf_complexread(x1,n,inp);
 
     y1 = sf_complexalloc(n);
-    for (i=0; i < n; i++) {
-	y1[i] = sf_cmplx(0.0f,0.0f);
-    }
 
     pchain_init(n,nc,x1,xn,xn+n*nc);
-
-    /* initialize */
-    for (ic=0; ic < nc; ic++) {
-	f = SF_PI*ic/nc;
-	for (i=0; i < n; i++) {
-	    xn[ic*n+i] =  cexpf(sf_cmplx(0.0f,f)); 
-            /* distribute around the unit circle */
-	}
-    }
-    for (ic=0; ic < nc-1; ic++) {
-	for (i=0; i < n; i++) {
-	    xn[(nc+ic)*n+i] = sf_cmplx(0.0f,0.0f); 
-	}
-    }
-
     r = sf_complexalloc(n*nc);
 
+    if (!sf_getbool("verb",&verb)) verb=(bool) (1 == nt);
+    /* verbosity flag */
     if (!sf_getint("niter",&niter)) niter=0;
     /* number of iterations */
     if (!sf_getint("liter",&liter)) liter=50;
@@ -83,27 +69,53 @@ int main(int argc, char* argv[])
 
     csmooth1_init(n,nc,rect);
 
-    sf_cconjgrad_init(n2, n2, n*nc, n*nc, 1., 1.e-6, true, false);
+    for (it=0; it < nt; it++) {
+	sf_warning("trace %d of %d;",it+1,nt);
+	sf_complexread(x1,n,inp);
+
+	/* initialize */
+	for (i=0; i < n; i++) {
+	    y1[i] = sf_cmplx(0.0f,0.0f);
+	}
+
+	for (ic=0; ic < nc; ic++) {
+	    f = SF_PI*ic/nc;
+	    for (i=0; i < n; i++) {
+		xn[ic*n+i] =  cexpf(sf_cmplx(0.0f,f)); 
+		/* distribute around the unit circle */
+	    }
+	}
+	for (ic=0; ic < nc-1; ic++) {
+	    for (i=0; i < n; i++) {
+		xn[(nc+ic)*n+i] = sf_cmplx(0.0f,0.0f); 
+	    }
+	}
 
-    for (iter=0; iter < niter; iter++) {
-	pchain_apply(y1,r);
+	sf_cconjgrad_init(n2, n2, n*nc, n*nc, 1., 1.e-6, verb, false);
 
-	for (i=0; i < n*nc; i++) {
-	    r[i] = -r[i];
-	}
+	for (iter=0; iter < niter; iter++) {
+	    pchain_apply(y1,r);
+
+	    for (i=0; i < n*nc; i++) {
+		r[i] = -r[i];
+	    }
 
-	sf_cconjgrad(NULL, pchain_lop, csmooth1_lop,p,dx,r,liter);
+	    sf_cconjgrad(NULL, pchain_lop, csmooth1_lop,p,dx,r,liter);
 
-	for (i=0; i < n2; i++) {
-	    xn[i] += dx[i];
+	    for (i=0; i < n2; i++) {
+		xn[i] += dx[i];
+	    }
 	}
-    }
 
-    sf_complexwrite(xn,n*nc,pef);
+	sf_cconjgrad_close();
 
-    pchain(y1);
+	sf_complexwrite(xn,n*nc,pef);
 
-    sf_complexwrite(y1,n,out);
+	pchain(y1);
+
+	sf_complexwrite(y1,n,out);
+    }
+    sf_warning(".");
 
     exit(0);
 }