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

Last change on this file since 3359 was 3298, checked in by kanani, 6 years ago

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