Home

Context Navigation

expand_decomp.C @ 4416

Last change on this file since 4416 was 3799, checked in by forkel, 6 years ago
editing in kpp4palm: add statements for avoiding unused variables, remove $Id
File size: 12.9 KB

Line
1
2	// ############################################################################
3	//
4	// create_mz_kpp_module
5	//
6	// create scalar code from .f90 sources created by KPP to be used in MESSy
7	//
8	// COPYRIGHT Klaus Ketelsen and MPI-CH April 2007
9	//
10	// ############################################################################
11	//
12	//
13	// Former revisions:
14	// -----------------------
15	// Deleted $Id since document_changes does not work for C and C++ (15.03.2019, forkel)
16	//
17	// Added switch for de-indexing == 2 (W is not needed then) (08.03.2019, forkel)
18	// commented add_line for 'INTEGER :: k, kk, j, jj' which
19	// only needed for de-indexing == 0
20	//
21	// initial version (Nov. 2016, ketelsen)
22
23
24	#include <fstream>
25	#include <stdio.h>
26
27
28	#include "expand_decomp.h"
29	#include "utils.h"
30
31	void expand_decomp::create_sparse_info (vector<fortran_file> & include_list, string module_name) {
32	vector<fortran_file>::iterator it;
33	vector<program_line>::iterator ip;
34	bool done;
35	int number;
36
37	// Generate a dummy vector entry [0] to be compatable with FORTRAN indexing
38	number = 0;
39	LU_IROW.push_back(number);
40	LU_ICOL.push_back(number);
41	LU_CROW.push_back(number);
42	LU_DIAG.push_back(number);
43
44	cout <<endl;
45	cout << "Sparse matrix: " ;
46
47	// Get NVAR
48
49	for(it=include_list.begin();it!=include_list.end();it++) {
50	if(it->get_name() == module_name + "_Parameters") {
51	for (ip=it->pline.begin(); ip != it->pline.end(); ip++) {
52	if(ip->get_token(0).substr(0,7) == "INTEGER" && ip->get_token(3) == "NVAR") {
53	sscanf (ip->get_token(5).c_str(),"%i",&NVAR);
54	cout << " NVAR = "<< NVAR <<endl;
55	}
56	}
57	}
58	}
59
60	// Get LU_IROW
61
62	for(it=include_list.begin();it!=include_list.end();it++) {
63	if(it->get_name() == module_name + "_JacobianSP") {
64	for (ip=it->pline.begin(); ip != it->pline.end(); ip++) {
65	if(ip->get_token(0).substr(0,7) == "INTEGER" &&
66	ip->get_token(7).substr(0,7) == "LU_IROW") {
67	done = false;
68	while(1) {
69	ip++;
70	for(int i=0; i<ip->get_token_size (); i++) {
71	if(ip->get_token(i).substr(0,1) == "/") {
72	done = true;
73	break;
74	} else if(ip->get_token(i).substr(0,1) != "&") {
75	if ( sscanf (ip->get_token(i).c_str(),"%i",&number) == 1)
76	LU_IROW.push_back(number);
77	}
78	}
79	if(done) break;
80	}
81	}
82	}
83	}
84	}
85	cout << "extracted LU_IROW: " << LU_IROW.size()-1 << " values " <<endl;
86
87	// Get LU_ICOL
88
89	for(it=include_list.begin();it!=include_list.end();it++) {
90	if(it->get_name() == module_name + "_JacobianSP") {
91	for (ip=it->pline.begin(); ip != it->pline.end(); ip++) {
92	if(ip->get_token(0).substr(0,7) == "INTEGER" &&
93	ip->get_token(7).substr(0,7) == "LU_ICOL") {
94	done = false;
95	while(1) {
96	ip++;
97	for(int i=0; i<ip->get_token_size (); i++) {
98	if(ip->get_token(i).substr(0,1) == "/") {
99	done = true;
100	break;
101	} else if(ip->get_token(i).substr(0,1) != "&") {
102	if (sscanf (ip->get_token(i).c_str(),"%i",&number) == 1)
103	LU_ICOL.push_back(number);
104	}
105	}
106	if(done) break;
107	}
108	}
109	}
110	}
111	}
112	cout << "extracted LU_ICOL: " << LU_ICOL.size()-1 << " values " <<endl;
113
114	// Get LU_CROW
115	// In case of large system, for LU_CROW and LU_DIAG might be the same logic
116	// necessary like LU_IROW and LU_ICOL
117
118	for(it=include_list.begin();it!=include_list.end();it++) {
119	if(it->get_name() == module_name + "_JacobianSP") {
120	for (ip=it->pline.begin(); ip != it->pline.end(); ip++) {
121	if(ip->get_token(0).substr(0,7) == "INTEGER" &&
122	ip->get_token(7).substr(0,7) == "LU_CROW") {
123	done = false;
124	while(1) {
125	ip++;
126	for(int i=0; i<ip->get_token_size (); i++) {
127	if(ip->get_token(i).substr(0,1) == "/") {
128	done = true;
129	break;
130	} else if(ip->get_token(i).substr(0,1) != "&") {
131	if (sscanf (ip->get_token(i).c_str(),"%i",&number) == 1)
132	LU_CROW.push_back(number);
133	}
134	}
135	if(done) break;
136	}
137	}
138	}
139	}
140	}
141	cout << "extracted LU_CROW: " << LU_CROW.size()-1 << " values " <<endl;
142
143	// Get LU_DIAG
144
145	for(it=include_list.begin();it!=include_list.end();it++) {
146	if(it->get_name() == module_name + "_JacobianSP") {
147	for (ip=it->pline.begin(); ip != it->pline.end(); ip++) {
148	if(ip->get_token(0).substr(0,7) == "INTEGER" &&
149	ip->get_token(7).substr(0,7) == "LU_DIAG") {
150	done = false;
151	while(1) {
152	ip++;
153	for(int i=0; i<ip->get_token_size (); i++) {
154	if(ip->get_token(i).substr(0,1) == "/") {
155	done = true;
156	break;
157	} else if(ip->get_token(i).substr(0,1) != "&") {
158	if (sscanf (ip->get_token(i).c_str(),"%i",&number) == 1)
159	LU_DIAG.push_back(number);
160	}
161	}
162	if(done) break;
163	}
164	}
165	}
166	}
167	}
168	cout << "extracted LU_DIAG: " << LU_DIAG.size()-1 << " values " <<endl;
169
170	cout <<endl;
171
172	return;
173	}
174
175	void expand_decomp::create_routine (vector<fortran_file> & routine_list) {
176
177	vector<fortran_file>::iterator it;
178	fortran_file de;
179	int k,kk,j,jj,i;
180	string line;
181	char cline[80];
182
183	// Delete original KppDecomp from subroutine list
184
185	for(it=routine_list.begin();it!=routine_list.end();it++) {
186	if(it->get_name() == "KppDecomp") {
187	routine_list.erase(it);
188	}
189	}
190
191	de.set_name("KppDecomp");
192
193	de.add_line (" SUBROUTINE KppDecomp( JVS, IER) ");
194	de.add_line (" ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ");
195	de.add_line (" ! Sparse LU factorization ");
196	de.add_line (" ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ");
197	de.add_line (" ! Loop expansion generated by kp4 ");
198	de.add_line (" ");
199	de.add_line (" INTEGER :: IER ");
200	// de.add_line (" REAL ( kind=dp ) :: JVS (:), W ( NVAR ), a ");
201	if( kpp_switches.de_indexing () == 2 ) {
202	de.add_line (" REAL ( kind=dp ) :: JVS ( LU_NONZERO ), a ");
203	} else {
204	de.add_line (" REAL ( kind=dp ) :: JVS ( LU_NONZERO ), W ( NVAR ), a ");
205	}
206	// Commented, k, kk, j, jj are only required if de_indexing == 0, which is not handled here
207	// de.add_line (" INTEGER :: k, kk, j, jj ");
208	de.add_line (" ");
209	de.add_line (" a = 0. ");
210	de.add_line (" IER = 0 ");
211	de.add_line (" ");
212
213	if( kpp_switches.de_indexing () == 1 ) {
214	cout << " *** de-indexing of deomposition algorithm: 1 " <<endl;
215	for (k=1; k<= NVAR; k++) {
216	if(LU_CROW[k] <= LU_DIAG[k]-1) {
217	sprintf(cline,"! k = %i",k);
218	line = cline; de.add_line(line);
219	for (kk=LU_CROW[k]; kk<=LU_CROW[k+1]-1; kk++) {
220	sprintf(cline," W ( %i )= JVS ( %i )",LU_ICOL[kk],kk);
221	line = cline; de.add_line(line);
222	}
223	for (kk=LU_CROW[k]; kk<=LU_DIAG[k]-1; kk++) {
224	j = LU_ICOL[kk];
225	sprintf(cline," a = - W ( %i ) / JVS ( %i )",j,LU_DIAG[j]);
226	line = cline; de.add_line(line);
227	sprintf(cline," W ( %i ) = -a",j);
228	line = cline; de.add_line(line);
229	for (jj=LU_DIAG[j]+1; jj<=LU_CROW[j+1]-1; jj++) {
230	sprintf(cline," W ( %i ) = W ( %i ) + a * JVS ( %i )",LU_ICOL[jj],LU_ICOL[jj],jj);
231	line = cline; de.add_line(line);
232	}
233	}
234	for (kk=LU_CROW[k]; kk<=LU_CROW[k+1]-1; kk++) {
235	sprintf(cline," JVS ( %i )= W ( %i )",kk,LU_ICOL[kk]);
236	line = cline; de.add_line(line);
237	}
238	} else {
239	sprintf(cline,"! k = %i Nothing to do LU_CROW[k] > LU_DIAG[k]-1 (%i %i)",
240	k,LU_CROW[k],LU_DIAG[k]-1);
241	line = cline; de.add_line(line);
242	}
243	}
244	}
245
246	if( kpp_switches.de_indexing () == 2 ) {
247
248	// No W array required
249	// No data copying
250	// better data reuse
251
252	int flag [NVAR+1] [NVAR+1];
253	int index [NVAR+1] [NVAR+1];
254	int ij,ji,jk,ki;
255	int count;
256	bool term;
257
258	cout << " *** de-indexing of decomposition algorithm: 2 (fast) " <<endl;
259
260	for (i=1; i<= NVAR; i++) {
261	for (j=1; j<= NVAR; j++) {
262	flag[i][j] = 0;
263	}
264	}
265	for (k=1; k<= LU_ICOL.size()-1; k++) {
266	i = LU_IROW[k];
267	j = LU_ICOL[k];
268	flag [j][i] = 1;
269	index[j][i] = k;
270	}
271
272	// for (i=1; i<= NVAR; i++) {
273	// cout << " " <<i << " ";
274	// for (j=1; j<= NVAR; j++) {
275	// cout << flag[i][j] << " ";
276	// }
277	// cout << endl;
278	// }
279
280	for (i=1; i<= NVAR; i++) {
281	sprintf(cline,"! i = %i",i);
282	line = cline; de.add_line(line);
283	for (j=1; j<= i-1; j++) {
284	if(flag[j][i] == 1) {
285	sprintf(cline," a = 0.0; a = a ");
286	line = cline;
287	term = false;
288	count=0;
289	for (k=1; k<= j-1; k++) {
290	if(flag[j][k] == 1 && flag[k][i] == 1) {
291	term = true;
292	count++;
293	if(count == 5) {
294	line += "&";
295	de.add_line(line);
296	count = 0;
297	line=" ";
298	}
299	jk = index [j][k];
300	ki = index [k][i];
301	sprintf(cline," - JVS ( %i ) * JVS ( %i )",jk,ki);
302	line += cline;
303	}
304	}
305	ji = index [j][i];
306	jj = index [j][j];
307	if(term) {
308	de.add_line(line);
309	sprintf(cline," JVS ( %i ) = ( JVS ( %i ) +a ) / JVS ( %i ) ",ji,ji,jj);
310	} else {
311	sprintf(cline," JVS ( %i ) = ( JVS ( %i ) ) / JVS ( %i ) ",ji,ji,jj);
312	}
313	line = cline; de.add_line(line);
314	}
315	}
316	for (j=i; j<= NVAR; j++) {
317	if(flag[j][i] == 1) {
318	ij = index [j][i];
319	sprintf(cline," JVS ( %i )= JVS ( %i )",ij,ij);
320	line = cline;
321	term = false;
322	count=0;
323	for (k=1; k<= i-1; k++) {
324	if(flag[j][k] == 1 && flag[k][i] == 1) {
325	term = true;
326	count++;
327	if(count == 5) {
328	line += "&";
329	de.add_line(line);
330	count = 0;
331	line=" ";
332	}
333	jk = index [j][k];
334	ki = index [k][i];
335	sprintf(cline," - JVS ( %i ) * JVS ( %i )",jk,ki);
336	line += cline;
337	}
338	}
339	if(term) {
340	de.add_line(line);
341	}
342	}
343	}
344	}
345	}
346
347	if( kpp_switches.de_indexing () == 3 ) {
348	cout << " *** de-indexing of decomposition algorithm: 3 (old method)"
349	<<endl;
350	for (k=1; k<= NVAR; k++) {
351	for (kk=LU_CROW[k]; kk<=LU_CROW[k+1]-1; kk++) {
352	sprintf(cline," W ( %i ) = JVS ( %i )", LU_ICOL[kk], kk);
353	line = cline; de.add_line(line);
354	}
355
356	for (kk=LU_CROW[k]; kk<=LU_DIAG[k]-1; kk++) {
357	j = LU_ICOL[kk] ;
358	sprintf(cline, " a = - W ( %i ) / JVS ( %i )", j, LU_DIAG[j]);
359	line = cline; de.add_line(line);
360
361	sprintf(cline, " W ( %i ) = -a ", j);
362	line = cline; de.add_line(line);
363
364	for (jj=LU_DIAG[j]+1; jj<=LU_CROW[j+1]-1; jj++) {
365	sprintf(cline," W ( %i ) = W ( %i ) + a * JVS ( %i )",
366	LU_ICOL[jj], LU_ICOL[jj], jj);
367	line = cline; de.add_line(line);
368	}
369	}
370
371	if ( LU_DIAG[k]-1 >= LU_CROW[k] ) {
372	// removes unnecessary statements
373	for (kk=LU_CROW[k]; kk<=LU_CROW[k+1]-1; kk++) {
374	sprintf(cline, " JVS ( %i ) = W ( %i )",kk, LU_ICOL[kk]);
375	line = cline; de.add_line(line);
376
377	}
378	}
379	}
380	}
381
382	de.add_line (" return ");
383	de.add_line (" ");
384	de.add_line (" End SUBROUTINE KppDecomp ");
385	routine_list.push_back(de);
386
387	return;
388	}

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/expand_decomp.C @ 4416

Download in other formats: