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 @ 3494

Last change on this file since 3494 was 3458, checked in by kanani, 6 years ago

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