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);
}