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

Last change on this file since 4654 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
31void 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
175void 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.