QR_MUMPS
qrm_rowcount.F90
Go to the documentation of this file.
1 !! ##############################################################################################
2 !!
3 !! Copyright 2012 CNRS, INPT
4 !!
5 !! This file is part of qr_mumps.
6 !!
7 !! qr_mumps is free software: you can redistribute it and/or modify
8 !! it under the terms of the GNU Lesser General Public License as
9 !! published by the Free Software Foundation, either version 3 of
10 !! the License, or (at your option) any later version.
11 !!
12 !! qr_mumps is distributed in the hope that it will be useful,
13 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
14 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 !! GNU Lesser General Public License for more details.
16 !!
17 !! You can find a copy of the GNU Lesser General Public License
18 !! in the qr_mumps/doc directory.
19 !!
20 !! ##############################################################################################
21 
22 
23 !! ##############################################################################################
33 
34 
35 #include "qrm_common.h"
36 
37 
39 
103 
104 subroutine _qrm_rowcount(graph, parent, porder, rc)
106 
107  use _qrm_spmat_mod
108  use qrm_mem_mod
109  use qrm_common_mod
110  implicit none
111 
112  type(_qrm_spmat_type) :: graph
113  integer :: parent(:), porder(:), rc(:)
114 
115  ! level : it is the level of node i, i.e., the distance from its root
116  ! fst_desc : it is the first descendent of node i according to the postorder
117 
118  integer, allocatable :: fst_desc(:), stack(:)
119  integer, allocatable :: first(:), hptr(:), hjcn(:), rcnt(:), last(:), ipord(:)
120  integer, allocatable :: prev_f(:), prev_nbr(:), setpath(:)
121  integer :: i, curr, fd, hp, rlev, j, f, ptr, k, ii, u, ref, p_leaf, q, jj
122  ! error management
123  integer :: err_act
124  character(len=*), parameter :: name='qrm_rowcount'
125 
126  call qrm_err_act_save(err_act)
127 
128  call qrm_aalloc(fst_desc, graph%n)
129  call qrm_aalloc(stack , graph%n)
130  __qrm_check_ret(name,'qrm_aalloc',9999)
131 
132  rc = 0
133  fst_desc=-1
134  do i=1, graph%n
135  curr = porder(i)
136  fd = curr
137  if(fst_desc(curr) .eq. -1) then
138  ! this is a leaf node. Initialize its row_count to 1
139  ! see [1], Fig 3.2
140  rc(curr) = 1
141  end if
142 
143  ! now we start going up thre tree seeting up fd as first
144  ! descendent of all its ancestors
145  do
146  if (fst_desc(curr) .gt. 0) exit
147  fst_desc(curr) = fd
148  if (parent(curr) .eq. 0) exit
149  curr = parent(curr)
150  end do
151  end do
152 
153  ! now we need to compute the "higher adjacency sets" (see [1] pag. 704)
154  ! Assume first(i)=j means that the first nonzero in row i has column index j.
155  ! The higher adjaceny set of column j is defined as the union of the column indices
156  ! of all the rows i for which first(i)=j.
157  ! The computation of the higher adjacency sets may be expensive in memory or in perf.
158  ! Here we choose (for the moment) to spend some more time on it but have a better
159  ! memory efficiency.
160  call qrm_aalloc(ipord, graph%n)
161  call qrm_aalloc(hptr , graph%n+1)
162  call qrm_aalloc(rcnt , graph%n)
163  call qrm_aalloc(first, graph%m)
164  call move_alloc(stack, last)
165  __qrm_check_ret(name,'qrm_aalloc',9999)
166 
167  ! compute inverse postorder
168  do j=1, graph%n
169  ipord(porder(j)) = j
170  end do
171 
172  first = 0
173  rcnt = 0
174  last = 0
175  ! we make a first pass to count the number of elements in each higher adjacency set
176  do k=1, graph%n
177  j = porder(k)
178  do ii=graph%jptr(j), graph%jptr(j+1)-1
179  i = graph%irn(ii)
180  f = first(i)
181  if(f .eq. 0) then
182  first(i) = j
183  else if (k .gt. last(f)) then
184  rcnt(f) = rcnt(f)+1
185  last(f) = k
186  end if
187  end do
188  end do
189 
190  ! transform the counts into pointers
191  hptr(1) = 1
192  do k=1, graph%n
193  hptr(k+1) = hptr(k)+rcnt(k)
194  end do
195 
196  ! allocate hjcn
197  call qrm_aalloc(hjcn, hptr(graph%n+1))
198  __qrm_check_ret(name,'qrm_aalloc',9999)
199 
200  rcnt = 0
201  last = 0
202  ! second loop to fill
203  do k=1, graph%n
204  j = porder(k)
205  do ii=graph%jptr(j), graph%jptr(j+1)-1
206  i = graph%irn(ii)
207  f = first(i)
208  if ( (k .gt. ipord(f)) .and. (k .gt. last(f)) ) then
209  ptr = hptr(f) + rcnt(f)
210  hjcn(ptr) = j
211  rcnt(f) = rcnt(f)+1
212  last(f) = k
213  end if
214  end do
215  end do
216 
217 
218  ! at this point we have everything we need to compute the rowcount.
219  ! This is the algorithm on Figure 3.2 in [1]
220  call qrm_adealloc(first)
221  call move_alloc(last, prev_f)
222  call move_alloc(rcnt, prev_nbr)
223  ! setpath is used for doing path compression
224  call qrm_aalloc(setpath, graph%n)
225  __qrm_check_ret(name,'aalloc/dealloc',9999)
226 
227  do j=1, graph%n
228  setpath(j)=j
229  end do
230  prev_f = 0
231  prev_nbr = 0
232 
233  ! during the symbolic facto we flag as negative all the first
234  ! variables of a supernode. At the end, we make a simple loop
235  ! to modify the tree so that it reflects the supernodal structure
236 
237  do jj=1, graph%n
238  j = porder(jj)
239  if(parent(j) .ne. 0) then
240  rc(parent(j)) = rc(parent(j))-1
241  else
242  ! if porder(jj) is a root, porder(jj+1) is the beginning of a supernode
243  if(jj .ne. graph%n) porder(jj+1) = porder(jj+1)
244  end if
245 
246  do k=hptr(j), hptr(j+1)-1
247  u = hjcn(k)
248  if(prev_nbr(u) .eq. 0) then
249  ref = 0
250  else
251  ref = ipord(prev_nbr(u))
252  end if
253 
254  if(ipord(fst_desc(j)) .gt. ref) then
255  rc(j) = rc(j)+1
256  p_leaf = prev_f(u)
257  ! porder(jj) is a leaf or a row subtree, thus the first variable in a
258  ! supernode. flag it
259  if (p_leaf .ne. 0) then
260  q = setfind(setpath, p_leaf)
261  ! q is a leaf or a row subtree, thus the first variable in a
262  ! supernode. flag it
263  rc(q) = rc(q) -1
264  end if
265  prev_f(u) = j
266  end if
267  prev_nbr(u) = j
268  end do
269  call setunion(setpath, j, parent(j))
270 
271  end do
272 
273  do jj=1, graph%n-1
274  j=porder(jj)
275  if(parent(j) .ne. 0) rc(parent(j)) = rc(parent(j)) + rc(j)
276  end do
277 
278  ! WARNING: this commented code is meant to detect supernodes in the tree based
279  ! on how the nodes where flagged above. It is currently replaced by the small
280  ! loop below as supernodes are detected in the amalgamation routine. As a result,
281  ! on output all the entries of parent will be positive.
282 
283  ! ! now we have to update the tree in order to reflect the supernodal structure.
284  ! ! we simply traverse the postorder looking for the variables we flagged above
285  ! i=1
286  ! do while (i .lt. graph%n)
287  ! porder(i) = unflip(porder(i))
288  ! j = i+1
289  ! do while( (porder(j) .gt. 0) .and. (j .le. graph%n))
290  ! j = j+1
291  ! end do
292  ! ! at this point j is the first variable of the next supernode and
293  ! ! j-1 is the principal variable of the supernode starting at i.
294  ! ! All the variables i:j-2 are subordinate to j-1
295  ! parent(porder(i)) = parent(porder(j-1))
296  ! do k=i+1, j-1
297  ! parent(porder(k)) = -porder(i)
298  ! rc(porder(k)) = 0
299  ! end do
300  ! ! we want the supervariable to be in first position in the porder
301  ! i = j
302  ! end do
303 
304  ! unflip the last in case it's a supervariable
305  ! porder(graph%n) = unflip(porder(graph%n))
306 
307  do i=1, graph%n
308  f = parent(i)
309  if (f .gt. 0) then
310  if(parent(f) .lt. 0) parent(i) = -parent(f)
311  end if
312  end do
313 
314  call qrm_adealloc(setpath)
315  call qrm_adealloc(hptr)
316  call qrm_adealloc(hjcn)
317  call qrm_adealloc(prev_nbr)
318  call qrm_adealloc(prev_f)
319  call qrm_adealloc(fst_desc)
320  call qrm_adealloc(ipord)
321  __qrm_check_ret(name,'qrm_adealloc',9999)
322 
323 
324  call qrm_err_act_restore(err_act)
325  return
326 
327 9999 continue ! error management
328  call qrm_err_act_restore(err_act)
329  if(err_act .eq. qrm_abort_) then
330  call qrm_err_check()
331  end if
332 
333  return
334 
335 
336 contains
337 
338  function setfind(setpath, p_leaf)
339  implicit none
340  integer :: setpath(:), p_leaf, setfind
341 
342  integer :: q, c, tmp
343 
344  q=p_leaf
345 
346  do while (setpath(q) .ne.q)
347  q = setpath(q)
348  end do
349 
350  ! path compression
351  c = p_leaf
352  do while (c .ne.q)
353  tmp = setpath(c)
354  setpath(c) = q
355  c = tmp
356  end do
357 
358  setfind = q
359  return
360  end function setfind
361 
362 
363  subroutine setunion(setpath, j, pj)
364  implicit none
365  integer :: setpath(:), j, pj
366  if(pj .ne. 0) setpath(j) = pj
367  return
368  end subroutine setunion
369 
370 
371  function flip(p)
372  integer :: p, flip
373  if(p .gt. 0) then
374  flip = -p
375  else
376  flip = p
377  end if
378  end function flip
379 
380  function unflip(p)
381  integer :: p, unflip
382  if(p .lt. 0) then
383  unflip = -p
384  else
385  unflip = p
386  end if
387  end function unflip
388 
389 end subroutine _qrm_rowcount
Generic interface for the qrm_adealloc_i, qrm_adealloc_2i, qrm_adealloc_s, qrm_adealloc_2s, qrm_adealloc_3s, qrm_adealloc_d, qrm_adealloc_2d, qrm_adealloc_3d, qrm_adealloc_c, qrm_adealloc_2c, qrm_adealloc_3c, qrm_adealloc_z, qrm_adealloc_2z, qrm_adealloc_3z, routines.
This module contains the interfaces of all non-typed routines.
subroutine setunion(setpath, j, pj)
integer function setfind(setpath, p_leaf)
integer function unflip(i)
Generic interface for the qrm_aalloc_i, qrm_aalloc_2i, qrm_aalloc_s, qrm_aalloc_2s, qrm_aalloc_3s, qrm_aalloc_d, qrm_aalloc_2d, qrm_aalloc_3d, qrm_aalloc_c, qrm_aalloc_2c, qrm_aalloc_3c, qrm_aalloc_z, qrm_aalloc_2z, qrm_aalloc_3z, routines.
Definition: qrm_mem_mod.F90:78
subroutine _qrm_rowcount(graph, parent, porder, rc)
This subroutine computes the rowcount of the R factor.
This type defines the data structure used to store a matrix.
This module contains the definition of the basic sparse matrix type and of the associated methods...
This module implements the memory handling routines. Pretty mucch allocations and deallocations...
Definition: qrm_mem_mod.F90:38
integer function flip(i)