QR_MUMPS
 All Classes Files Functions Variables Enumerations Enumerator Pages
dqrm_detect_singletons.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 
38 
81 subroutine dqrm_detect_singletons(graph, scol, srow, mrperm, mcperm, nrsing, ncsing)
82 
83  use dqrm_spmat_mod
84  use qrm_common_mod
85  use qrm_mem_mod
86 
87  implicit none
88 
89  type(dqrm_spmat_type) :: graph
90  integer, allocatable :: scol(:), srow(:)
91  integer, allocatable :: mcperm(:), mrperm(:)
92  integer :: ncsing, nrsing
93 
94 
95  type(dqrm_spmat_type) :: graph_t
96 
97  integer, allocatable :: nel(:), isrow(:), tmp(:)
98  integer :: j, row, col, jj, ii, i, curr, ncol_p, nrow_p, nz_p, nz_save, idx
99  ! error management
100  integer :: err_act
101  character(len=*), parameter :: name='qrm_detect_singletons'
102 
103  call qrm_err_act_save(err_act)
104 
105 
106  ! nel contains the number of nnz in each column
107  ! isrow is the row corresponding to a singleton column (singleton row)
108  ! squeue is the singletons queue.
109  ! srow may not be big enough for internal work so we must allocate isrow
110  ! and then copy/compress its content back into srow
111  call qrm_arealloc(isrow , graph%n)
112  call qrm_arealloc(mcperm, graph%n)
113  call qrm_arealloc(mrperm, graph%m)
114  call qrm_arealloc(scol , graph%n)
115  call qrm_arealloc(srow , graph%m)
116  call move_alloc(mcperm, nel)
117  __qrm_check_ret(name,'qrm_arealloc',9999)
118 
119 
120  ! the singleton detection will proceed like this:
121  ! a singleton queue will be built with a preliminary
122  ! pass over the columns of graph.
123  ! Then, for each singleton in the queue, the corresponding row
124  ! is removed (this is equivalent to updating nel on the
125  ! columns corresponding to nnzs in the singleton row). Once
126  ! the row is removed, if a new singleton is found, it is added
127  ! at the end of the queue.
128 
129  ! first pass to fill the queue and compute the degrees
130  ncsing = 0
131  do j=1, graph%n
132  nel(j) = graph%jptr(j+1)-graph%jptr(j)
133  if (nel(j) .eq. 0) then
134  ! j is a "dead" singleton
135  ncsing = ncsing+1
136  scol(ncsing) = j
137  isrow(ncsing) = 0
138  else if (nel(j) .eq. 1) then
139  ! j is a singleton
140  ncsing = ncsing+1
141  scol(ncsing) = j
142  isrow(ncsing) = graph%irn(graph%jptr(j))
143  end if
144  end do
145 
146  if((ncsing .eq. 0) .or. (ncsing .eq. graph%n)) then
147  ! shortcut!
148  call qrm_adealloc(isrow)
149  __qrm_check_ret(name,'qrm_dealloc',9999)
150  call move_alloc(nel, mcperm)
151  return
152  end if
153 
154  ! we need both a copy by columns and a copy by rows of the graph
155  call dqrm_spmat_convert(graph, graph_t, 'csr', values=.false.)
156  __qrm_check_ret(name,'qrm_spmat_convert',9999)
157 
158  ! loop on all the elements in the queue
159  curr = 1
160  do
161  if(curr .gt. ncsing) exit
162  col = scol(curr)
163  row = isrow(curr)
164  if ((row .eq. 0)) then
165  curr = curr+1
166  cycle
167  end if
168  if(graph_t%iptr(row) .lt.0) then
169  ! row was already eliminated so col
170  ! has become a dead singleton
171  isrow(curr) = 0
172  curr = curr+1
173  cycle
174  end if
175  graph_t%iptr(row) = flip(graph_t%iptr(row))
176  ! loop over all the elements j in row i
177  do jj= unflip(graph_t%iptr(row)), unflip(graph_t%iptr(row+1)) -1
178  j = graph_t%jcn(jj)
179  ! update the degree
180  nel(j) = nel(j)-1
181  if (nel(j) .eq. 1) then
182  ! new singleton found
183  ncsing = ncsing+1
184  scol(ncsing) = j
185  ! now we must find the corresponding singleton row
186  ! we loop over all the elements in the column and look
187  ! for that whose row pointer was not flipped in graph_t
188  do ii=graph%jptr(j), graph%jptr(j+1)-1
189  i = graph%irn(ii)
190  if(graph_t%iptr(i) .gt. 0) exit
191  end do
192  isrow(ncsing) = i
193  end if
194  end do
195  curr = curr+1
196  end do
197 
198  ! remove all dead singletons from isrow
199  curr=0
200  do i=1, ncsing
201  if(isrow(i) .ne. 0) then
202  curr = curr+1
203  srow(curr) = isrow(i)
204  end if
205  end do
206 
207  ! at this point we must prune the graph by removing all the singleton columns
208  ! and the corresponding singleton rows.
209 
210  ! (non dead) singleton rows are those for which graph_t%iptr is negative
211  ! to identify (non dead) singleton columns we must built an inverse permutation array
212  ! which (hopefully) will come handy later on
213  call move_alloc(nel, mcperm)
214 
215  ! purge
216  mcperm = 0
217  do j = 1,ncsing
218  idx = scol(j)
219  mcperm(idx) = -j
220  end do
221 
222  mrperm = 0
223  nrow_p = 0
224  do i=1, graph_t%m
225  if(graph_t%iptr(i) .lt. 0) cycle
226  nrow_p = nrow_p+1
227  mrperm(i) = nrow_p
228  end do
229 
230  ncol_p = 0
231  nz_p = 1
232  nz_save = 1
233  do j=1, graph%n
234  ! if mcperm < 0 j is a column singleton. just skip it
235  if (mcperm(j) .lt. 0) cycle
236  ncol_p = ncol_p+1
237  ! this next instruction works because j>=ncol_p
238  mcperm(ncol_p) = j
239  do ii= graph%jptr(j), graph%jptr(j+1)-1
240  i = graph%irn(ii)
241  ! go down the column and skip all the elements belonging to singleton rows
242  if (graph_t%iptr(i) .lt. 0) cycle
243  graph%irn(nz_p) = mrperm(i)
244  nz_p = nz_p+1
245  end do
246  graph%jptr(ncol_p) = nz_save
247  nz_save = nz_p
248  end do
249  graph%jptr(ncol_p+1) = nz_save
250 
251 #if defined (debug)
252  if((ncsing + ncol_p) .ne. graph%n) then
253  __qrm_prnt_dbg('("Inconsistency in singleton detection")')
254  end if
255 #endif
256 
257  nrow_p=0
258  do i=1, graph_t%m
259  if(mrperm(i) .gt. 0) then
260  ! this next instruction works because i>=nrow_p
261  nrow_p = nrow_p+1
262  mrperm(nrow_p) = i
263  end if
264  end do
265  mrperm(nrow_p+1:graph_t%m) = 0
266 
267 
268  nrsing = graph%m-nrow_p
269  graph%m = nrow_p
270  graph%n = ncol_p
271  graph%nz = nz_p-1
272 
273  call dqrm_spmat_destroy(graph_t, all=.true.)
274  call qrm_adealloc(isrow)
275  __qrm_check_ret(name,'spmat_destroy/dealloc',9999)
276 
277  call qrm_aalloc(tmp, ncsing)
278  tmp(1:ncsing) = scol(1:ncsing)
279  call qrm_adealloc(scol)
280  call move_alloc(tmp, scol)
281 
282  call qrm_aalloc(tmp, nrsing)
283  tmp(1:nrsing) = srow(1:nrsing)
284  call qrm_adealloc(srow)
285  call move_alloc(tmp, srow)
286 
287 
288  call qrm_err_act_restore(err_act)
289  return
290 
291 9999 continue ! error management
292  call qrm_err_act_restore(err_act)
293  if(err_act .eq. qrm_abort_) then
294  call qrm_err_check()
295  end if
296 
297  return
298 
299 
300 contains
301 
302  function flip(i)
303  integer :: flip
304  integer :: i
305  flip = -i
306  return
307  end function flip
308 
309  function unflip(i)
310  integer :: unflip
311  integer :: i
312  if (i .lt. 0) then
313  unflip = -i
314  else
315  unflip = i
316  end if
317  return
318  end function unflip
319 
320 
321 end subroutine dqrm_detect_singletons
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.
subroutine dqrm_detect_singletons(graph, scol, srow, mrperm, mcperm, nrsing, ncsing)
This subroutine detects singletons in a matrix.
This module contains the interfaces of all non-typed routines.
subroutine qrm_err_act_restore(err_act)
Restores the value of the qrm_err_act variable.
integer function unflip(i)
This module contains the definition of the basic sparse matrix type and of the associated methods...
subroutine qrm_err_check()
This subroutine checks the errors stack. If something is found all the entries in the stack are poppe...
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 dqrm_spmat_destroy(qrm_spmat, all)
This subroutine destroyes a qrm_spmat instance.
subroutine dqrm_spmat_convert(in_mat, out_mat, fmt, values)
This subroutine converts an input matrix into a different storage format. Optionally the values may b...
This type defines the data structure used to store a matrix.
Generic interface for the qrm_arealloc_i qrm_arealloc_s qrm_arealloc_d qrm_arealloc_c qrm_arealloc_z...
This module implements the memory handling routines. Pretty mucch allocations and deallocations...
Definition: qrm_mem_mod.F90:38
int i
Definition: secs.c:40
integer function flip(i)
subroutine qrm_err_act_save(err_act)
Saves a copy of the qrm_err_act variable.