QR_MUMPS
qrm_init_front.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 
115 
116 subroutine _qrm_init_front(qrm_mat, fnum, par, work)
118  use _qrm_spmat_mod
119  use _qrm_fdata_mod
120  use _qrm_utils_mod
121  use qrm_common_mod
122  use qrm_sort_mod
123  implicit none
124 
125  type(_qrm_spmat_type), target :: qrm_mat
126  integer :: fnum
127  logical :: par
128  integer, optional, target :: work(:)
129 
130  integer :: n_rows_orig, i, j, row, col, c, child, roff, f
131  integer :: m, n, npiv, k, p, ne, cbr, cbc
132  integer :: nb, b, father
133 
134  integer, allocatable :: first(:)
135  integer, pointer :: gcolmap(:)
136 
137  type(_qrm_front_type), pointer :: front, cfront
138  type(qrm_adata_type), pointer :: adata
139  type(_qrm_fdata_type), pointer :: fdata
140  logical :: map, sfront
141 
142  ! error management
143  integer :: err_act
144  character(len=*), parameter :: name='qrm_init_front'
145 
146  call qrm_err_act_save(err_act)
147 
148  front => null()
149  cfront => null()
150  nullify(gcolmap)
151 
152  ! to make things easier
153  adata => qrm_mat%adata
154  fdata => qrm_mat%fdata
155  front => fdata%front_list(fnum)
156  front%m = adata%nfrows(fnum)
157  front%n = adata%rc(fnum)
158  if(par) then
159  front%nb = qrm_mat%icntl(qrm_nb_)
160  else
161  front%nb = qrm_mat%icntl(qrm_ib_)
162  end if
163  front%ib = qrm_mat%icntl(qrm_ib_)
164  father = adata%parent(fnum)
165 
166  if( (front%n .le. 0) .or. (front%m .le. 0)) then
167  ! nothing to do here. Mark the front as done
168  front%status = qrm_done_
169  goto 10
170  end if
171 
172  if (fnum .eq. 1) then
173  roff = 0
174  else
175  roff = adata%stair(fnum-1)
176  end if
177 
178  ! set up a few variables that come handy in the activation and
179  ! factorization of the front
180  m = front%m
181  n = front%n
182  ! The number of eliminations to be performed on the front f
183  front%ne = min(m,n)
184  ! The number of panels in f
185  front%np = max((front%ne-1)/front%nb + 1,0)
186  ! The number of block-columns in f
187  front%nc = (front%n-1)/front%nb + 1
188  ! The number of pivots in f
189  front%npiv = min(adata%cp_ptr(fnum+1)-adata%cp_ptr(fnum),front%ne)
190  cbr = front%m-front%npiv ! The number of rows in R
191  cbc = front%n-front%npiv ! The number of rows in the cb
192 
193  ! if the size of the front is smaller than ib, then there's no need
194  ! to compute front%stair
195  sfront = front%ne .lt. front%ib
196  ! if(sfront) write(*,*)'====> sfront!!!'
197  ! allocate front%cols and copy the list f column indices in it
198  call qrm_aalloc(front%cols, front%n)
199  __qrm_check_ret(name,'qrm_aalloc',9999)
200 
201  front%cols(1:front%n) = adata%fcol(adata%fcol_ptr(fnum): &
202  & adata%fcol_ptr(fnum+1)-1)
203 
204  if(par) then
205  ! front%ptable is the progress table for keeping track of the work
206  ! done on this front. ptable(i) is initialized to the number of
207  ! contributions that have to be assembled into block-column i
208  call qrm_aalloc(front%ptable, front%nc )
209  __qrm_check_ret(name,'qrm_aalloc',9999)
210  front%ptable = 0
211  end if
212 
213  ! build the gcolmap for front f. gcolmap(k)=j means that global
214  ! column k is column j inside front f
215  if(present(work)) then
216  gcolmap => work
217  else
218  call qrm_palloc(gcolmap, qrm_mat%n)
219  __qrm_check_ret(name,'qrm_palloc',9999)
220  end if
221 
222  do j=1, front%n
223  k = front%cols(j)
224  gcolmap(k)=j
225  end do
226 
227  call qrm_aalloc(front%stair, front%n+1)
228  __qrm_check_ret(name,'qrm_aalloc',9999)
229 
230  if(.not. sfront) then
231  front%stair = 0
232  call qrm_aalloc(first, front%anrows)
233  __qrm_check_ret(name,'qrm_aalloc',9999)
234 
235  first = front%n+1
236  ! count in the rows from the original matrix
237  do i=1, front%anrows
238  ! sweep this row and determine its first coefficient
239  do p=front%aiptr(i), front%aiptr(i+1)-1
240  j = gcolmap(front%ajcn(p))
241  ! TODO: can be optimized id ajcn is sorted
242  if (j .lt. first(i)) first(i)=j
243  end do
244  front%stair(first(i)+1) = front%stair(first(i)+1)+1
245  end do
246 
247  ! count in rows coming from CBs
248  do p = adata%childptr(fnum), adata%childptr(fnum+1)-1
249  c = adata%child(p)
250  cfront => fdata%front_list(c)
251 
252  ! ne is the number of Householder vectors computed on the
253  ! child c. npiv is the number of fully assembled pivots in c
254  ne = cfront%ne
255  npiv = cfront%npiv
256  if(npiv.eq.ne) cycle
257 
258  ! count in all the rows on the CB of c
259  do i=1, ne-npiv
260  f = gcolmap(cfront%cols(npiv+i)) !cfront%colmap(i)
261  front%stair(f+1) = front%stair(f+1)+1
262  end do
263  end do
264 
265  ! finalize stair
266  do i=2, front%n+1
267  front%stair(i) = front%stair(i)+front%stair(i-1)
268  end do
269  else
270  front%stair = front%m
271  end if
272 
273  call qrm_aalloc(front%rows, front%m)
274  call qrm_aalloc(front%front, front%m, front%n)
275  __qrm_check_ret(name,'qrm_aalloc',9999)
276 
277  front%front = _qrm_zero
278  row = 0
279  ! at this point we're ready to assemble the rows from the original matrix
280  ! FIXME: deal with the presence of duplicates
281  do i=1, front%anrows
282  if(sfront) then
283  row = i
284  else
285  f = first(i) ! f is the front-local column index of the first
286  ! coefficient in this row
287  front%stair(f) = front%stair(f)+1
288  row = front%stair(f)
289  end if
290 
291  front%rows(row) = adata%rperm(roff+i)
292  !sweep this row and assemble it
293  do p=front%aiptr(i), front%aiptr(i+1)-1
294  col = gcolmap(front%ajcn(p))
295  front%front(row,col) = front%front(row,col)+front%aval(p)
296  end do
297  end do
298 
299 
300  ! now we can build the row-mappings for assemblying the children
301  ! later or, if required, assemble the children directly
302  do p = adata%childptr(fnum), adata%childptr(fnum+1)-1
303  c = adata%child(p)
304  cfront => fdata%front_list(c)
305 
306  map = par .and. (adata%small(c) .ne. 1)
307 
308  ! ne is the number of Householder vectors computed on the
309  ! child c. npiv is the number of fully assembled pivots in c
310  ne = cfront%ne
311  npiv = cfront%npiv
312  if(npiv.eq.ne) cycle
313 
314  ! this is the row mapping on the child. cfront%rowmap(k)=i
315  ! means that the k-th row of cfront will be assembled into the
316  ! i-th row of front
317  do i=npiv+1, ne
318  if(sfront) then
319  row = row+1
320  else
321  f = gcolmap(cfront%cols(i))
322  front%stair(f) = front%stair(f)+1
323  row = front%stair(f)
324  end if
325  front%rows(row) = cfront%rows(i)
326  cfront%rowmap(i-npiv) = row
327  end do
328 
329  if(.not. map) then
330  do j=npiv+1, cfront%n
331  ! this is the column mapping on the child. cfront%colmap(k)=j
332  ! means that the k-th column of cfront will be assembled into the
333  ! j-th column of front
334  cfront%colmap(j-npiv) = gcolmap(cfront%cols(j))
335 
336  ! fill in the front with coefficients front chil front cfront
337  do i=npiv+1, min(j,ne)
338  row = cfront%rowmap(i-npiv)
339  col = cfront%colmap(j-npiv)
340  front%front(row, col) = cfront%front(i,j)
341  end do
342  end do
343  else
344  ! initialize ptable
345  do j=1, cfront%n - npiv
346  ! this is the column mapping on the child. cfront%colmap(k)=j
347  ! means that the k-th column of cfront will be assembled into the
348  ! j-th column of front
349  cfront%colmap(j) = gcolmap(cfront%cols(npiv+j))
350 
351  ! count the number of contributions that have to be
352  ! assembled inside the corresponding block-column
353  b = (cfront%colmap(j)-1)/front%nb+1 ! j in cfront goes into b-column b of front
354  front%ptable(b) = front%ptable(b)-1
355  end do
356  end if
357  end do
358 
359  call qrm_adealloc(first)
360  if(present(work)) then
361  nullify(gcolmap)
362  else
363  call qrm_pdealloc(gcolmap)
364  nullify(gcolmap)
365  end if
366  __qrm_check_ret(name,'qrm_adealloc',9999)
367 
368  ! allocate tau and t
369  call qrm_aalloc(front%tau, front%ne)
370  call qrm_aalloc(front%t, front%ib, front%ne)
371  __qrm_check_ret(name,'qrm_aalloc',9999)
372 
373 
374  ! FIXME: check if this is really necessary
375  front%tau = _qrm_zero
376  front%t = _qrm_zero
377 
378  if (father .ne. 0) then
379  call qrm_aalloc(front%rowmap, cbr)
380  call qrm_aalloc(front%colmap, cbc)
381  __qrm_check_ret(name,'qrm_aalloc',9999)
382  end if
383 
384  ! no more needed, can be deallocated
385  call qrm_adealloc(front%aiptr)
386  call qrm_adealloc(front%ajcn)
387  call qrm_adealloc(front%aval)
388  __qrm_check_ret(name,'qrm_adealloc',9999)
389 
390 
391 10 continue
392 
393  if (father .ne. 0) then
394  front => fdata%front_list(father)
395  !$ if(par) call omp_set_lock(front%lock)
396  front%status = front%status+1
397  !$ if(par) call omp_unset_lock(front%lock)
398  end if
399 
400 
401  call qrm_err_act_restore(err_act)
402  return
403 
404 9999 continue
405  ! We get to this point if some allocation failed. Deallocate
406  ! everything and return
407  call qrm_adealloc(first)
408  if(present(work)) then
409  nullify(gcolmap)
410  else
411  call qrm_pdealloc(gcolmap)
412  end if
413  call qrm_adealloc(front%cols)
414  call qrm_adealloc(front%ptable)
415  call qrm_adealloc(front%stair)
416  call qrm_adealloc(front%rows)
417  call qrm_adealloc(front%front)
418  call qrm_adealloc(front%tau)
419  call qrm_adealloc(front%t)
420  call qrm_adealloc(front%colmap)
421  call qrm_adealloc(front%rowmap)
422 
423  call qrm_err_act_restore(err_act)
424  if(err_act .eq. qrm_abort_) then
425  call qrm_err_check()
426  end if
427 
428  return
429 
430 end subroutine _qrm_init_front
This module contains generic interfaces for a number of auxiliary tools.
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 _qrm_init_front(qrm_mat, fnum, par, work)
This routine initializes a front.
This module contains routines for sorting.
This module contains the interfaces of all non-typed routines.
integer, parameter qrm_done_
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
Generic interface for the qrm_pdealloc_i, qrm_pdealloc_2i, qrm_pdealloc_s, qrm_pdealloc_2s, qrm_pdealloc_d, qrm_pdealloc_2d, qrm_pdealloc_c, qrm_pdealloc_2c, qrm_pdealloc_z, qrm_pdealloc_2z, routines.
Definition: qrm_mem_mod.F90:98
This type defines the data structure used to store a matrix.
The data structure meant to store all the results of the factorization phase.
This module contains the definition of the basic sparse matrix type and of the associated methods...
This module contains the definition of all the data related to the factorization phase.
Generic interface for the qrm_palloc_i, qrm_palloc_2i, qrm_palloc_s, qrm_palloc_2s, qrm_palloc_d, qrm_palloc_2d, qrm_palloc_c, qrm_palloc_2c, qrm_palloc_z, qrm_palloc_2z, routines.
Definition: qrm_mem_mod.F90:57
This module implements the memory handling routines. Pretty mucch allocations and deallocations...
Definition: qrm_mem_mod.F90:38
This type defines a data structure containing all the data related to a front.