lbfgs.c 4.53 KB
Newer Older
Simone Butzer's avatar
Simone Butzer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
/*------------------------------------------------------------------------
 * Copyright (C) 2015 For the list of authors, see file AUTHORS.
 *
 * This file is part of IFOS3D.
 * 
 * IFOS3D 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, version 2.0 of the License only.
 * 
 * IFOS3D 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 IFOS3D. See file COPYING and/or 
 * <http://www.gnu.org/licenses/gpl-2.0.html>.
--------------------------------------------------------------------------*/

/*-------------------------------------------------------------------------
 * calculatipon of L-BFGS update
 * S. Butzer 2014
 --------------------------------------------------------------------------*/

#include "fd.h"
void lbfgs(float ***grad1, float ***grad2, float ***grad3, float *bfgsscale, float **bfgsmod, float **bfgsgrad, int iteration){

28
	int m=0,v=0,w=0;
Simone Butzer's avatar
Simone Butzer committed
29
30
31
32
33
34
	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;
35
	extern int BFGSNUM, NUMPAR;
Simone Butzer's avatar
Simone Butzer committed
36
37

	
38
39
	dummy1 = vector(1,BFGSNUM);
	q = vector(1,NUMPAR*NX*NY*NZ);
Simone Butzer's avatar
Simone Butzer committed
40
41
42
	dummy[0]=0.0; dummy[1]=0.0;
	buf[0]=0.0; buf[1]=0.0;
		
43
	m=iteration-BFGSNUM;
Simone Butzer's avatar
Simone Butzer committed
44
45
46
47
48
	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*/
49
50
	w=(iteration-1)%BFGSNUM;
	if(w==0) w=BFGSNUM;
Simone Butzer's avatar
Simone Butzer committed
51
52
53
54
55
56
57
58
59
60
	l=0;
	for (j=1;j<=NY;j++){
		for (i=1;i<=NX;i++){
			for (k=1;k<=NZ;k++){
				l++;
				bfgsgrad[w][l]+=-grad1[j][i][k];
				q[l]=grad1[j][i][k];
			}
		}
	}
61
	if(NUMPAR==2){
Simone Butzer's avatar
Simone Butzer committed
62
63
64
65
66
67
68
69
70
71
72
73
74
75
		l=NX*NY*NZ;
		for (j=1;j<=NY;j++){
			for (i=1;i<=NX;i++){
				for (k=1;k<=NZ;k++){
					l++;
					bfgsgrad[w][l]+=-grad2[j][i][k];
					q[l]=grad2[j][i][k];
				}
			}
		}
	}
	
	
	/*calculate bfgsscale and H_0*/
76
77
	w=(iteration-1)%BFGSNUM; if(w==0) w=BFGSNUM;
	for (l=1;l<=NUMPAR*NX*NY*NZ;l++){
Simone Butzer's avatar
Simone Butzer committed
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
		dummy[0]+=bfgsgrad[w][l]*bfgsmod[w][l];
		dummy[1]+=bfgsgrad[w][l]*bfgsgrad[w][l];
	}
	
	MPI_Allreduce(&dummy,&buf,2,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
		
	bfgsscale[w]=1/buf[0];
	/*printf("bfgsscale(%d)=%e\n",w,bfgsscale[w]);*/
	h0=buf[0]/buf[1];
	/*printf("h0=%e\n",h0);*/
	/*---------------------------------------*/
	
		
	/*L-BFGS two-loop recursion (numerical optimisation Algorithm 9.1)*/
		  
	l=0;			  
	for(v=iteration-1; v>=m;v--){
	 /* printf("dummy10[%d]=%e",w,dummy1[w]);
	  fprintf(FP,"v1=%d\n",v);*/
97
98
		 w=v%BFGSNUM;
		 if(w==0) w=BFGSNUM;
Simone Butzer's avatar
Simone Butzer committed
99
		 fprintf(FP,"w1=%d\n",w);
100
		 for (l=1;l<=NUMPAR*NX*NY*NZ;l++){
Simone Butzer's avatar
Simone Butzer committed
101
102
103
104
105
106
107
108
			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]);*/
		 }
		/* printf("dummy11[%d]=%e",w,dummy1[w]);*/
		 buf[0]=0.0;
		 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]);*/
109
		 for (l=1;l<=NUMPAR*NX*NY*NZ;l++){
Simone Butzer's avatar
Simone Butzer committed
110
111
112
113
114
115
			q[l]+=-dummy1[w]*bfgsgrad[w][l];
			
		 }
	}
	
	
116
		 for (l=1;l<=NUMPAR*NX*NY*NZ;l++){
Simone Butzer's avatar
Simone Butzer committed
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
			dummy2=q[l];
			q[l]=dummy2*h0;
		 }
	
	/* if(HESS){l=0;
	   	for (j=1;j<=NY;j++){
			for (i=1;i<=NX;i++){
				for (k=1;k<=NZ;k++){ 
					l++;
					dummy2=q[l];
					q[l]=dummy2*hess[j][i][k];
				}
			}
		}
	 }*/
	 
	for(v=m; v<=iteration-1;v++){
	  /*fprintf(FP,"v2=%d\n",v);*/
135
136
		w=v%BFGSNUM;
		if(w==0) w=BFGSNUM;
Simone Butzer's avatar
Simone Butzer committed
137
138
139
		/*fprintf(FP,"w2=%d\n",w);*/
		dummy2=0.0;
		
140
		for (l=1;l<=NUMPAR*NX*NY*NZ;l++){
Simone Butzer's avatar
Simone Butzer committed
141
142
143
144
145
146
147
			dummy2+=bfgsscale[w]*bfgsgrad[w][l]*q[l];
		}
		
		buf[0]=0.0;
		MPI_Allreduce(&dummy2,&buf[0],1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
		dummy2=buf[0];
		
148
		for (l=1;l<=NUMPAR*NX*NY*NZ;l++){
Simone Butzer's avatar
Simone Butzer committed
149
150
151
152
153
			q[l]+=bfgsmod[w][l]*(dummy1[w]-dummy2);	
		}
	}
	
	/*save gradients of this iteration and update gradient*/
154
155
	w=iteration%BFGSNUM;
		if(w==0) w=BFGSNUM;
Simone Butzer's avatar
Simone Butzer committed
156
157
158
159
160
161
162
163
164
165
166
	
	l=0;
	for (j=1;j<=NY;j++){
		for (i=1;i<=NX;i++){
			for (k=1;k<=NZ;k++){  
				l++;  /*w=(j-1)*Nx*Nz+(i-1)*Nz+k*/  
				bfgsgrad[w][l]=grad1[j][i][k];
				grad1[j][i][k]=q[l];
			}
		}
	}
167
	if(NUMPAR==2){
Simone Butzer's avatar
Simone Butzer committed
168
169
170
171
172
173
174
175
176
177
178
		l=NX*NY*NZ;
		for (j=1;j<=NY;j++){
			for (i=1;i<=NX;i++){
				for (k=1;k<=NZ;k++){  
					l++;  /*w=(j-1)*Nx*Nz+(i-1)*Nz+k*/  
					bfgsgrad[w][l]=grad2[j][i][k];
					grad2[j][i][k]=q[l];
				}
			}
		}
	}
179
free_vector(q,1,NUMPAR*NX*NY*NZ);
Simone Butzer's avatar
Simone Butzer committed
180
181
free_vector(dummy1,1,m);
}