Commit 70308326 authored by Simone Butzer's avatar Simone Butzer

add parameters BFGSNUM and NUMPAR

(Input file and source code)
parent 8293f938
......@@ -184,3 +184,5 @@ Water_level_Hessian_for_vp/vs/rho_(WATER_HESS) = 0.0192, 0.0192, 1.0e-14
#
#------------------------L-BFGS ----------------------------------------
Apply_L_BFGS_(LBFGS) = 0
Number_of_inverted_parameters_(NUMPAR) = 2
Number_iterations_used_for_LBFGS_(BFGSNUM) = 5
\ No newline at end of file
......@@ -184,3 +184,5 @@ Water_level_Hessian_for_vp/vs/rho_(WATER_HESS) = 0.0192, 0.0192, 1.0e-14
#
#------------------------L-BFGS ----------------------------------------
Apply_L_BFGS_(LBFGS) = 0
Number_of_inverted_parameters_(NUMPAR) = 2
Number_iterations_used_for_LBFGS_(BFGSNUM) = 5
\ No newline at end of file
......@@ -42,10 +42,11 @@ void exchange_par(void){
extern char FILEINP[STRING_SIZE];
extern int LITTLEBIG, ASCIIEBCDIC, IEEEIBM;
extern int METHOD;
extern int ITMIN, ITMAX, FILT, NFMAX, TAST, NSHOTS_STEP, DAMPTYPE, HESS, READ_HESS, REC_HESS,LBFGS,EXTOBS;
extern int ITMIN, ITMAX, FILT, NFMAX, TAST, NSHOTS_STEP, DAMPTYPE, HESS, READ_HESS, REC_HESS, EXTOBS;
extern int BFGSNUM, NUMPAR, LBFGS;
/*extern float F_INV;*/
extern float TESTSTEP, WATER_HESS[3], WEIGHT[3], VP0, VS0, RHO0;
int idum[55];
int idum[57];
float fdum[42];
......@@ -152,12 +153,14 @@ void exchange_par(void){
idum[52] = READ_HESS;
idum[53] = REC_HESS;;
idum[54] = LBFGS;
idum[55] = NUMPAR;
idum[56] = BFGSNUM;
}
if (MYID != 0) FL=vector(1,L);
MPI_Barrier(MPI_COMM_WORLD);
MPI_Bcast(&idum,55,MPI_INT,0,MPI_COMM_WORLD);
MPI_Bcast(&idum,57,MPI_INT,0,MPI_COMM_WORLD);
MPI_Bcast(&fdum,42,MPI_FLOAT,0,MPI_COMM_WORLD);
MPI_Bcast(&SOURCE_FILE,STRING_SIZE,MPI_CHAR,0,MPI_COMM_WORLD);
......@@ -279,6 +282,8 @@ void exchange_par(void){
READ_HESS = idum[52];
REC_HESS =idum[53];
LBFGS = idum[54];
NUMPAR = idum[55];
BFGSNUM = idum[56];
MPI_Bcast(&FL[1],L,MPI_FLOAT,0,MPI_COMM_WORLD);
......
......@@ -56,6 +56,6 @@ int METHOD;
float F_INV=0.0, TESTSTEP=0.0;
char MOD_OUT_FILE[STRING_SIZE], HESS_FILE[STRING_SIZE], GRAD_FILE[STRING_SIZE], SEIS_OBS_FILE[STRING_SIZE];
int DAMPTYPE, ITMIN, ITMAX, FILT, NFMAX, NSHOTS_STEP, TAST,EXTOBS;
int HESS, READ_HESS, REC_HESS, LBFGS;
int HESS, READ_HESS, REC_HESS, LBFGS, NUMPAR, BFGSNUM;
float WATER_HESS[3]={0.0, 0.0, 0.0}, WEIGHT[3]={0.0, 0.0, 0.0};
float VP0, VS0, RHO0;
\ No newline at end of file
......@@ -30,7 +30,7 @@ void hess_apply(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, float ***g
extern float WATER_HESS[3], VP0, VS0, RHO0;
extern int NXG, NYG, NZG, FW, NX, NY, NZ;
extern FILE *FP;
int LBFGS=0;
extern int LBFGS;
float wl[3],buf[3];
int i,j,k,buf1;
......
......@@ -80,7 +80,6 @@ int main(int argc, char **argv){
float dummy=0.0;
int hloop=0,pshot=0,pshot1=0, pshot_loc=0;
int LBFGS=0, bfgsnum=5;
float **bfgsmod1=NULL,**bfgsgrad1=NULL, *bfgsscale1=NULL; /* *bfgsscale2,*bfgsscale3;**bfgsmod2,**bfgsmod3,**bfgsgrad2, **bfgsgrad3,*/
......@@ -234,7 +233,7 @@ if(METHOD) nseismograms+=4;
memgrad=12*fac2*fac3;
memdynf=12*NFMAX*(ntr_hess+1)*fac1*fac2;
if(HESS) memgrad+=3*fac2*fac3;
if(LBFGS) membfgs=2*bfgsnum*3*fac3;
if(LBFGS) membfgs=NUMPAR*BFGSNUM*3*fac3;
}
memtotal=memdyn+memmodel+memseismograms+membuffer+(buffsize*pow(2.0,-20.0))+memgrad+memdynf+membfgs;
......@@ -361,13 +360,13 @@ MPI_Barrier(MPI_COMM_WORLD);
}
if(LBFGS){
bfgsmod1=fmatrix(1,bfgsnum,1,2*NX*NY*NZ);
bfgsmod1=fmatrix(1,BFGSNUM,1,NUMPAR*NX*NY*NZ);
/*bfgsmod2=fmatrix(1,bfgsnum,1,NX*NY*NZ);
bfgsmod3=fmatrix(1,bfgsnum,1,NX*NY*NZ);*/
bfgsgrad1=fmatrix(1,bfgsnum,1,2*NX*NY*NZ);
bfgsgrad1=fmatrix(1,BFGSNUM,1,NUMPAR*NX*NY*NZ);
/*bfgsgrad2=fmatrix(1,bfgsnum,1,NX*NY*NZ);
bfgsgrad3=fmatrix(1,bfgsnum,1,NX*NY*NZ);*/
bfgsscale1=vector(1,bfgsnum);
bfgsscale1=vector(1,BFGSNUM);
/*bfgsscale2=vector(1,bfgsnum);
bfgsscale3=vector(1,bfgsnum);*/
}
......@@ -1204,9 +1203,9 @@ CPML_coeff(K_x,alpha_prime_x,a_x,b_x,K_x_half,alpha_prime_x_half,a_x_half,b_x_ha
}
if(LBFGS){
free_matrix(bfgsmod1,1,bfgsnum,1,2*NX*NY*NZ);
free_matrix(bfgsgrad1,1,bfgsnum,1,2*NX*NY*NZ);
free_vector(bfgsscale1,1,bfgsnum);
free_matrix(bfgsmod1,1,BFGSNUM,1,NUMPAR*NX*NY*NZ);
free_matrix(bfgsgrad1,1,BFGSNUM,1,NUMPAR*NX*NY*NZ);
free_vector(bfgsscale1,1,BFGSNUM);
}
if(ABS_TYPE==1){
......
......@@ -25,28 +25,29 @@
#include "fd.h"
void lbfgs(float ***grad1, float ***grad2, float ***grad3, float *bfgsscale, float **bfgsmod, float **bfgsgrad, int iteration){
int m=0,v=0,w=0,bfgsnum=5, npar=2;
int m=0,v=0,w=0;
int i,j,k,l;
float dummy[2], *dummy1, dummy2=0.0, buf[2];
float *q;
float h0;
extern int NX,NY,NZ; /*,HESS*/
extern FILE *FP;
extern int BFGSNUM, NUMPAR;
dummy1 = vector(1,bfgsnum);
q = vector(1,npar*NX*NY*NZ);
dummy1 = vector(1,BFGSNUM);
q = vector(1,NUMPAR*NX*NY*NZ);
dummy[0]=0.0; dummy[1]=0.0;
buf[0]=0.0; buf[1]=0.0;
m=iteration-bfgsnum;
m=iteration-BFGSNUM;
if(m<1) m=1;
fprintf(FP,"Start calculation L-BFGS update");
/*save gradient for bfgsgrad(iteration-1) and set q=gradient, attention grad is -gradient of misfit function*/
w=(iteration-1)%bfgsnum;
if(w==0) w=bfgsnum;
w=(iteration-1)%BFGSNUM;
if(w==0) w=BFGSNUM;
l=0;
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
......@@ -57,7 +58,7 @@ void lbfgs(float ***grad1, float ***grad2, float ***grad3, float *bfgsscale, flo
}
}
}
if(npar==2){
if(NUMPAR==2){
l=NX*NY*NZ;
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
......@@ -72,8 +73,8 @@ void lbfgs(float ***grad1, float ***grad2, float ***grad3, float *bfgsscale, flo
/*calculate bfgsscale and H_0*/
w=(iteration-1)%bfgsnum; if(w==0) w=bfgsnum;
for (l=1;l<=npar*NX*NY*NZ;l++){
w=(iteration-1)%BFGSNUM; if(w==0) w=BFGSNUM;
for (l=1;l<=NUMPAR*NX*NY*NZ;l++){
dummy[0]+=bfgsgrad[w][l]*bfgsmod[w][l];
dummy[1]+=bfgsgrad[w][l]*bfgsgrad[w][l];
}
......@@ -93,10 +94,10 @@ void lbfgs(float ***grad1, float ***grad2, float ***grad3, float *bfgsscale, flo
for(v=iteration-1; v>=m;v--){
/* printf("dummy10[%d]=%e",w,dummy1[w]);
fprintf(FP,"v1=%d\n",v);*/
w=v%bfgsnum;
if(w==0) w=bfgsnum;
w=v%BFGSNUM;
if(w==0) w=BFGSNUM;
fprintf(FP,"w1=%d\n",w);
for (l=1;l<=npar*NX*NY*NZ;l++){
for (l=1;l<=NUMPAR*NX*NY*NZ;l++){
dummy1[w]+=bfgsscale[w]*bfgsmod[w][l]*q[l];
/*if(iteration==3)fprintf(FP,"bfgsgrad[%d][%d]=%e,bfgsscale[%d]=%e,bfgsmod[%d][%d]=%e,q[%d]=%e\n", w,l,bfgsgrad[w][l],w,bfgsscale[w],w,l,bfgsmod[w][l],l,q[l]);*/
}
......@@ -105,14 +106,14 @@ void lbfgs(float ***grad1, float ***grad2, float ***grad3, float *bfgsscale, flo
MPI_Allreduce(&dummy1[w],&buf[0],1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
dummy1[w]=buf[0];
/*printf("dummy1(%d)=%e\n",w,dummy1[w]);*/
for (l=1;l<=npar*NX*NY*NZ;l++){
for (l=1;l<=NUMPAR*NX*NY*NZ;l++){
q[l]+=-dummy1[w]*bfgsgrad[w][l];
}
}
for (l=1;l<=npar*NX*NY*NZ;l++){
for (l=1;l<=NUMPAR*NX*NY*NZ;l++){
dummy2=q[l];
q[l]=dummy2*h0;
}
......@@ -131,12 +132,12 @@ void lbfgs(float ***grad1, float ***grad2, float ***grad3, float *bfgsscale, flo
for(v=m; v<=iteration-1;v++){
/*fprintf(FP,"v2=%d\n",v);*/
w=v%bfgsnum;
if(w==0) w=bfgsnum;
w=v%BFGSNUM;
if(w==0) w=BFGSNUM;
/*fprintf(FP,"w2=%d\n",w);*/
dummy2=0.0;
for (l=1;l<=npar*NX*NY*NZ;l++){
for (l=1;l<=NUMPAR*NX*NY*NZ;l++){
dummy2+=bfgsscale[w]*bfgsgrad[w][l]*q[l];
}
......@@ -144,14 +145,14 @@ void lbfgs(float ***grad1, float ***grad2, float ***grad3, float *bfgsscale, flo
MPI_Allreduce(&dummy2,&buf[0],1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
dummy2=buf[0];
for (l=1;l<=npar*NX*NY*NZ;l++){
for (l=1;l<=NUMPAR*NX*NY*NZ;l++){
q[l]+=bfgsmod[w][l]*(dummy1[w]-dummy2);
}
}
/*save gradients of this iteration and update gradient*/
w=iteration%bfgsnum;
if(w==0) w=bfgsnum;
w=iteration%BFGSNUM;
if(w==0) w=BFGSNUM;
l=0;
for (j=1;j<=NY;j++){
......@@ -163,7 +164,7 @@ void lbfgs(float ***grad1, float ***grad2, float ***grad3, float *bfgsscale, flo
}
}
}
if(npar==2){
if(NUMPAR==2){
l=NX*NY*NZ;
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
......@@ -175,6 +176,6 @@ void lbfgs(float ***grad1, float ***grad2, float ***grad3, float *bfgsscale, flo
}
}
}
free_vector(q,1,npar*NX*NY*NZ);
free_vector(q,1,NUMPAR*NX*NY*NZ);
free_vector(dummy1,1,m);
}
......@@ -25,7 +25,7 @@
void lbfgs_savegrad(float ***grad1, float ***grad2, float ***grad3,float **bfgsgrad1){
int d,i,j,k;
extern int NX,NY, NZ;
extern int NX,NY, NZ, NUMPAR;
d=0;
......@@ -39,16 +39,17 @@ void lbfgs_savegrad(float ***grad1, float ***grad2, float ***grad3,float **bfgsg
}
}
}
d=NX*NY*NZ;
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
for (k=1;k<=NZ;k++){
d++;
bfgsgrad1[1][d]=+grad2[j][i][k];
if(NUMPAR>1){
d=NX*NY*NZ;
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
for (k=1;k<=NZ;k++){
d++;
bfgsgrad1[1][d]=+grad2[j][i][k];
}
}
}
}
......
......@@ -26,7 +26,7 @@
void modelupdate(int nx, int ny, int nz, float ***gradvp, float ***gradvs, float ***gradrho, float *** rho, float *** pi, float *** u, float **bfgsmod1, float step, float *beta, int it_group){
extern FILE *FP;
extern int NX, NY, NZ, LBFGS;
extern int NX, NY, NZ, LBFGS, BFGSNUM, NUMPAR;
extern float VP0, VS0, RHO0, WEIGHT[3];
......@@ -34,7 +34,7 @@ void modelupdate(int nx, int ny, int nz, float ***gradvp, float ***gradvs, float
float vpnew,vsnew,rhonew;
float max[3],buf[3],max1[3],dummy[3];
float vp,vs;
int bfgsnum=5,w;
int w;
float scale1=0.0, scale2=0.0, scale3=0.0;
buf[0]=0.0;buf[1]=0.0;buf[2]=0.0;
......@@ -49,8 +49,8 @@ void modelupdate(int nx, int ny, int nz, float ***gradvp, float ***gradvs, float
if(iteration>40)scale2=0.7;
if(iteration>60)scale2=0.4;*/
if(LBFGS) {w=it_group%bfgsnum;
if(w==0) w=bfgsnum;}
if(LBFGS) {w=it_group%BFGSNUM;
if(w==0) w=BFGSNUM;}
fprintf(FP,"Modelupdatefunction \n");
......@@ -143,8 +143,9 @@ void modelupdate(int nx, int ny, int nz, float ***gradvp, float ***gradvs, float
/*save normalised model differences for LBFGS*/
l++; l1++;
bfgsmod1[w][l]=(vpnew-sqrt(pi[j][i][k]/rho[j][i][k]))/VP0;
if(NUMPAR>1){
bfgsmod1[w][l1]=(vsnew-sqrt(u[j][i][k]/rho[j][i][k]))/VS0;
/*bfgsmod3[w][l]=(rhonew-rho[j][i][k])/rho0;*/
}
}
u[j][i][k]=rhonew*vsnew*vsnew;
......
......@@ -68,11 +68,11 @@ void outgrad(int nx,int ny,int nz,float ***grad1, float ***grad2,float ***grad3,
MPI_Barrier(MPI_COMM_WORLD);
if(MYID==0){
sprintf(gradfile4,"%s.gradvp_%4.2fHz_it%d",outfile,finv,iteration);
sprintf(gradfile4,"%s.vp_%4.2fHz_it%d",outfile,finv,iteration);
mergemod(gradfile4,3);
sprintf(gradfile5,"%s.gradvs_%4.2fHz_it%d",outfile,finv,iteration);
sprintf(gradfile5,"%s.vs_%4.2fHz_it%d",outfile,finv,iteration);
mergemod(gradfile5,3);
sprintf(gradfile6,"%s.gradrho_%4.2fHz_it%d",outfile,finv,iteration);
sprintf(gradfile6,"%s.rho_%4.2fHz_it%d",outfile,finv,iteration);
mergemod(gradfile6,3);
}
MPI_Barrier(MPI_COMM_WORLD);
......
......@@ -49,6 +49,7 @@ int read_par(FILE *fp_in){
extern int ITMIN, ITMAX, FILT, NFMAX, TAST, NSHOTS_STEP, DAMPTYPE, HESS, READ_HESS, REC_HESS,EXTOBS,LBFGS;
/*extern float F_INV;*/
extern float TESTSTEP, WATER_HESS[3], WEIGHT[3], VP0, VS0, RHO0;
extern int BFGSNUM, NUMPAR;
/* definition of local variables */
char s[256], cline[256]="";
......@@ -431,6 +432,12 @@ int read_par(FILE *fp_in){
sscanf(cline,"%s =%i",s,&LBFGS);
break;
case 83 :
sscanf(cline,"%s =%i",s,&NUMPAR);
break;
case 84 :
sscanf(cline,"%s =%i",s,&BFGSNUM);
break;
case 85 :
nvarin=sscanf(cline,"%s =%i ,%i ,%i",s,&ASCIIEBCDIC,&LITTLEBIG,&IEEEIBM);
switch(nvarin){
case 0: ;
......@@ -445,7 +452,7 @@ int read_par(FILE *fp_in){
}
}
LOG=0;
if(lineno<83) fprintf(stderr," Warning: only %d non-commentary lines of input parameters read (expected 83).\n",lineno);
if(lineno<85) fprintf(stderr," Warning: only %d non-commentary lines of input parameters read (expected 83).\n",lineno);
/* else if (lineno>71) fprintf(stderr," Warning: %d non-commentary lines of input parameters read \n \t(expected and interpreted: 67).\n",lineno); */
/* else fprintf(stderr," %d non-commentary lines of input parameters read.\n",lineno); */
fclose(fp_in);
......
......@@ -42,6 +42,8 @@ void writepar(FILE *fp, int ns){
extern int METHOD;
extern int NP, NPROCX, NPROCY, NPROCZ, MYID;
extern int ITMIN, ITMAX, FILT, NFMAX, TAST, NSHOTS_STEP, DAMPTYPE, HESS, READ_HESS, REC_HESS, LBFGS,EXTOBS;
extern int NUMPAR, BFGSNUM;
extern float TESTSTEP,WATER_HESS[3], WEIGHT[3], VP0, VS0, RHO0;
/* definition of local variables */
char th1[3], file_ext[8];
......@@ -446,8 +448,9 @@ fprintf(fp,"Read_Hessian_from_file %d\n",READ_HESS);
fprintf(fp,"Part_of_receivers_used_for_Hessian %d",REC_HESS);
fprintf(fp," \nWater_level_Hessian_for_vp/vs/rho %e \n",WATER_HESS[0]);
fprintf(fp," \n------------------------LBFGS-----------------------------------------\n");
fprintf(fp," \nLBFGS %i \n",LBFGS);
fprintf(fp," \nLBFGS: %i \n",LBFGS);
fprintf(fp," Number_of_inverted_parameters_(NUMPAR): %i \n",NUMPAR);
fprintf(fp," Number_iterations_used_for_LBFGS: %i \n",BFGSNUM);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment