mumps_fortran_solver.F
Go to the documentation of this file.
1 C========================================================================
2 C Broken up fortran interface to mumps solver, callable from C++.
3 C Based on demo driver code in
4 C
5 C MUMPS 4.8.4, built on Mon Dec 15 15:31:38 UTC 2008
6 C
7 C========================================================================
8 
9 
10 c======================================================================
11 c Module as common block replacement
12 c======================================================================
13  module mumps_module
14 
15  include 'dmumps_struc.h'
16 
17  INTEGER n_mumps
18  INTEGER, dimension(:), allocatable :: workspace_scaling_factor
19  TYPE (dmumps_struc), dimension(:), allocatable :: mumps_par_pool
20 
21  end module mumps_module
22 
23 
24 c======================================================================
25 c Setup pool of mumps solvers
26 c======================================================================
27  subroutine mumps_setup_solver_pool(max_n_solvers)
28 
29  use mumps_module
30 
31  INTEGER max_n_solvers
32 
33  n_mumps=max_n_solvers
34  ALLOCATE( mumps_par_pool(max_n_solvers) )
35  ALLOCATE( workspace_scaling_factor(max_n_solvers) )
36 
37  end
38 
39 
40 c======================================================================
41 c Catch wrong ID
42 c======================================================================
43  subroutine mumps_catch_wrong_id(i_pool)
44 
45  use mumps_module
46 
47  INTEGER i_pool
48 
49  if (i_pool .gt. n_mumps) then
50  write(*,*) "i_pool ", i_pool," exceeds max. number of"
51  write(*,*) "mumps solvers in pool (n_mumps=", n_mumps,")"
52  stop
53  endif
54 
55  end
56 
57 c======================================================================
58 c Setup routine
59 c======================================================================
60  subroutine mumps_setup(i_pool,initial_workspace_scaling_factor)
61 
62  use mumps_module
63 
64  integer i_pool
65 
66  INTEGER ierr, i, nproc, nrows_per_proc
67  integer initial_workspace_scaling_factor
68  double precision error, max_error
69 
70 c.....Range checking
71  call mumps_catch_wrong_id(i_pool)
72 
73 C.....Initialize an instance of the package
74 C.....for L U factorization (sym = 0, with working host)
75 
76 C.....Initialise
77  mumps_par_pool(i_pool)%JOB = -1
78 
79 C.....Matrix is not symmetric
80  mumps_par_pool(i_pool)%SYM = 0
81 
82 C.....We're computing on the host too
83  mumps_par_pool(i_pool)%PAR = 1
84 
85 C.....Do it...
86  CALL dmumps(mumps_par_pool(i_pool))
87 
88 C.....Output stream for global info on host. Negative value suppresses printing
89  mumps_par_pool(i_pool)%icntl(3)=-1
90 
91 C.....Only show error messages and stats
92  mumps_par_pool(i_pool)%icntl(4)=2
93 
94 C.....Assembled matrix (rather than element-by_element)
95  mumps_par_pool(i_pool)%icntl(5)=0
96 
97 C.....Distributed problem with user-specified distribution
98  mumps_par_pool(i_pool)%icntl(18)=3
99 
100 C.....Dense RHS
101  mumps_par_pool(i_pool)%icntl(20)=0
102 
103 C.....Non-distributed solution
104  mumps_par_pool(i_pool)%icntl(21)=0
105 
106 C.....Default scaling factor for workspace
107  workspace_scaling_factor(i_pool)=initial_workspace_scaling_factor
108 
109  return
110  end
111 
112 
113 
114 
115 c======================================================================
116 c Specify multiplier for workspace scaling factor. Default is 2.
117 c======================================================================
119 
120  use mumps_module
121 
122  integer i_pool
123  double precision s
124 
125 c.....Range checking
126  call mumps_catch_wrong_id(i_pool)
127 
128 C.....Assign...
129  workspace_scaling_factor(i_pool)=s
130 
131  return
132  end
133 
134 
135 
136 c======================================================================
137 c Switch on doc of stats etc. Redirects output to "channel 6"
138 c======================================================================
139  subroutine mumps_switch_on_doc(i_pool)
140 
141  use mumps_module
142 
143  integer i_pool
144 
145 c.....Range checking
146  call mumps_catch_wrong_id(i_pool)
147 
148 C.....Output stream for global info on host. Negative value suppresses printing
149  mumps_par_pool(i_pool)%icntl(3)=6
150 
151  return
152  end
153 
154 
155 
156 c======================================================================
157 c Switch off doc of stats etc.
158 c======================================================================
159  subroutine mumps_switch_off_doc(i_pool)
160 
161 
162  use mumps_module
163 
164  integer i_pool
165 
166 c.....Range checking
167  call mumps_catch_wrong_id(i_pool)
168 
169 C.....Output stream for global info on host. Negative value suppresses printing
170  mumps_par_pool(i_pool)%icntl(3)=-1
171 
172  return
173  end
174 
175 
176 
177 
178 c======================================================================
179 c Solver routine
180 c======================================================================
181  subroutine mumps_solve(i_pool,n,nz_loc,irn_loc,jcn_loc,a_loc,rhs)
182 
183  use mumps_module
184 
185  integer i_pool
186 
187  integer n, nz_loc
188  integer irn_loc(nz_loc), jcn_loc(nz_loc)
189  integer initial_workspace_scaling_factor
190  double precision a_loc(nz_loc), rhs(n)
191 
192 c.....Range checking
193  call mumps_catch_wrong_id(i_pool)
194 
195 C.....call mumps_setup
196  initial_workspace_scaling_factor=2
197  call mumps_setup(i_pool,initial_workspace_scaling_factor)
198 
199 C.....analyse/factorise
200  call mumps_factorise(i_pool,n,nz_loc,irn_loc,jcn_loc,a_loc)
201 
202 C.....back-substitute
203  call mumps_backsub(i_pool,n,rhs)
204 
205 C.....Shutdown
206  CALL mumps_shutdown(i_pool)
207 
208  return
209  end
210 
211 
212 c======================================================================
213 c Cleanup memory routine
214 c======================================================================
215  subroutine mumps_cleanup_memory(i_pool)
216 
217  use mumps_module
218 
219  integer i_pool
220 
221 c.....Range checking
222  call mumps_catch_wrong_id(i_pool)
223 
224 C.....Deallocate user data
225  IF ( associated(mumps_par_pool(i_pool)%IRN_loc) )THEN
226  DEALLOCATE( mumps_par_pool(i_pool)%IRN_loc)
227  end if
228  IF ( associated(mumps_par_pool(i_pool)%JCN_loc) )THEN
229  DEALLOCATE( mumps_par_pool(i_pool)%JCN_loc)
230  end if
231  IF ( associated(mumps_par_pool(i_pool)%A_loc) )THEN
232  DEALLOCATE( mumps_par_pool(i_pool)%A_loc)
233  end if
234  IF ( associated(mumps_par_pool(i_pool)%RHS) )THEN
235  DEALLOCATE( mumps_par_pool(i_pool)%RHS)
236  END IF
237 
238  RETURN
239 
240  END
241 
242 
243 
244 
245 c======================================================================
246 c Shutdown mumps
247 c======================================================================
248  subroutine mumps_shutdown(i_pool)
249 
250  use mumps_module
251 
252  integer i_pool
253 
254 c.....Range checking
255  call mumps_catch_wrong_id(i_pool)
256 
257 C.....Cleanup memory just to be on the safe side
258  CALL mumps_cleanup_memory(i_pool)
259 
260 C.....Destroy the instance (deallocate internal data structures)
261  mumps_par_pool(i_pool)%JOB = -2
262  CALL dmumps(mumps_par_pool(i_pool))
263 
264  RETURN
265 
266  END
267 
268 
269 c======================================================================
270 c Analyse and factorise
271 c======================================================================
272  subroutine mumps_factorise(i_pool,n,nz_loc,irn_loc,jcn_loc,a_loc)
273 
274  use mumps_module
275 
276  integer i_pool
277 
278  integer n, nz_loc
279  integer irn_loc(nz_loc),jcn_loc(nz_loc)
280  double precision a_loc(nz_loc)
281 
282  INTEGER ierr, i, nproc, nrows_per_proc
283  double precision error, max_error
284 
285 c.....Range checking
286  call mumps_catch_wrong_id(i_pool)
287 
288 C.....Specify size of system
289  mumps_par_pool(i_pool)%N=n
290  mumps_par_pool(i_pool)%NZ_loc=nz_loc
291 
292 C.....Allocate storage for entries
293  ALLOCATE( mumps_par_pool(i_pool)%IRN_loc(
294  & mumps_par_pool(i_pool)%NZ_loc ) )
295  ALLOCATE( mumps_par_pool(i_pool)%JCN_loc(
296  & mumps_par_pool(i_pool)%NZ_loc ) )
297  ALLOCATE( mumps_par_pool(i_pool)%A_loc(
298  & mumps_par_pool(i_pool)%NZ_loc ) )
299 
300  DO i = 1, mumps_par_pool(i_pool)%NZ_loc
301  mumps_par_pool(i_pool)%IRN_loc(i)=irn_loc(i)
302  mumps_par_pool(i_pool)%JCN_loc(i)=jcn_loc(i)
303  mumps_par_pool(i_pool)%A_loc(i)=a_loc(i)
304  END DO
305 
306 
307 C.....Do analysis
308  mumps_par_pool(i_pool)%JOB = 1
309  CALL dmumps(mumps_par_pool(i_pool))
310 
311 C.....Document estimate for size of workspace
312  if (mumps_par_pool(i_pool)%icntl(3).ge.0) then
313  if (mumps_par_pool(i_pool)%myid.eq.0) then
314  write(*,*) 'Estimated max. workspace in MB: ',
315  & mumps_par_pool(i_pool)%infog(26)
316  end if
317  end if
318 
319 c.....Jump address for re-factorisation with larger workspace
320  10 continue
321 
322 C.....Set workspace to multiple of that -- ought to be "significantly
323 C.....larger than infog(26)", according to the manual :(
324  mumps_par_pool(i_pool)%icntl(23)=
325  & workspace_scaling_factor(i_pool)*
326  & mumps_par_pool(i_pool)%infog(26)
327 
328 C.....Do factorisation
329  mumps_par_pool(i_pool)%JOB = 2
330  CALL dmumps(mumps_par_pool(i_pool))
331 
332 C.....check for error
333  if (mumps_par_pool(i_pool)%infog(1).ne.0) then
334 C........Increase scaling factor for workspace and run again
335  if (mumps_par_pool(i_pool)%myid.eq.0) then
336  write(*,*) 'Processor', mumps_par_pool(i_pool)%myid,
337  & 'has error during mumps factorisation!'
338  write(*,*) 'Processor', mumps_par_pool(i_pool)%myid,
339  & 'error codes: ',
340  & mumps_par_pool(i_pool)%info(1), " ",
341  & mumps_par_pool(i_pool)%info(2)
342  end if
343  workspace_scaling_factor(i_pool)=
344  & workspace_scaling_factor(i_pool)*2
345  if (mumps_par_pool(i_pool)%myid.eq.0) then
346  write(*,*) 'Increasing workspace_scaling_factor to ',
347  & workspace_scaling_factor(i_pool)
348  end if
349  go to 10
350  else
351  if (mumps_par_pool(i_pool)%icntl(3).ge.0) then
352  if (mumps_par_pool(i_pool)%myid.eq.0) then
353  write(*,*)
354  & 'Successfully completed factorisation in mumps'
355  end if
356  end if
357  end if
358 
359 
360  RETURN
361 
362  END
363 
364 
365 c======================================================================
366 c Back-subsititute/solve
367 c======================================================================
368  subroutine mumps_backsub(i_pool,n,rhs)
369 
370  use mumps_module
371 
372  integer i_pool
373 
374  INTEGER ierr, i, nproc, nrows_per_proc
375  double precision rhs(n)
376 
377 
378 c.....Range checking
379  call mumps_catch_wrong_id(i_pool)
380 
381  ALLOCATE( mumps_par_pool(i_pool)%RHS( n ) )
382  DO i=1,n
383  mumps_par_pool(i_pool)%RHS(i)=rhs(i)
384  END DO
385 
386 C.....Do solve/back-substitution
387  mumps_par_pool(i_pool)%JOB = 3
388  CALL dmumps(mumps_par_pool(i_pool))
389 
390 C.....copy it back
391  DO i=1,n
392  rhs(i)=mumps_par_pool(i_pool)%RHS(i)
393  END DO
394 
395  RETURN
396 
397  END
398 
subroutine mumps_set_workspace_scaling_factor(i_pool, s)
double & max_error()
Access fct for max. actual error in present solution (i.e. before re-solve on adapted mesh) ...
cstr elem_len * i
Definition: cfortran.h:607
subroutine mumps_catch_wrong_id(i_pool)
mumps_factorise
Definition: mumps.h:53
subroutine mumps_switch_off_doc(i_pool)
mumps_setup
Definition: mumps.h:37
subroutine mumps_solve(i_pool, n, nz_loc, irn_loc, jcn_loc, a_loc, rhs)
subroutine mumps_cleanup_memory(i_pool)
static char t char * s
Definition: cfortran.h:572
mumps_switch_on_doc
Definition: mumps.h:45
subroutine mumps_setup_solver_pool(max_n_solvers)
mumps_shutdown
Definition: mumps.h:68
mumps_backsub
Definition: mumps.h:60