Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

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

Last change on this file since 3789 was 3789, checked in by forkel, 5 years ago

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