← Back to branch summary

~siesta-maint/siesta/trunk

« back to all changes in this revision

Viewing changes to Src/create_Sparsity_Union.F90

  • Committer: Nick Papior
  • Date: 2018-04-07 19:24:04 UTC
  • mfrom: (560.1.325 4.1)
  • Revision ID: nickpapior@gmail.com-20180407192404-z5414ljigmh07vd3
Merged r882-885, bug-fixes non-equilibrium contours in example and added tiling

Show diffs side-by-side

added added

removed removed

Lines of Context:
Src/create_Sparsity_Union.F90
31
31
 
32
32
  private
33
33
 
 
34
  interface crtSparsity_Union
 
35
    module procedure crtSparsity_Union_explicit
 
36
    module procedure crtSparsity_Union_region
 
37
  end interface crtSparsity_Union
34
38
  public :: crtSparsity_Union
35
 
  public :: crtSparsity_Union_region
36
39
 
37
40
contains
38
41
 
47
50
  ! Needless to say, this routine makes not much sense when the final sparsity is > 50%
48
51
  ! Also you should consider that this will only work properly on unit-cell sparse
49
52
  ! patterns.
50
 
  subroutine crtSparsity_Union(dit,in,in_i,in_j,ni,nj,out)
 
53
  subroutine crtSparsity_Union_explicit(dit,in,in_i,in_j,ni,nj,out)
51
54
    use class_Sparsity
52
55
    use class_OrbitalDistribution
53
56
    use parallel, only : Node
115
118
    end if
116
119
 
117
120
    ! Prepare creation of num and listptr arrays
118
 
    allocate(num    (n_rows))
 
121
    allocate(num(n_rows))
119
122
    allocate(listptr(n_rows))
120
123
 
121
124
    ! The list pointer for the first entry is always the "first"
235
238
    ! The newSparsity copies the values...
236
239
    deallocate(num,listptr,list)
237
240
 
238
 
  end subroutine crtSparsity_Union
 
241
  end subroutine crtSparsity_Union_explicit
239
242
 
240
243
  ! Subroutine will create a new pattern which matches will be the union of 
241
244
  ! the 'in' sparse pattern and a segmented block in the full matrix.
269
272
    integer, allocatable :: num(:), listptr(:), list(:)
270
273
 
271
274
    integer :: lir,ir,i,iu,ptr,l,j
272
 
    type(tRgn) :: sr
273
275
 
274
276
    ! Local variables...
275
277
    integer, pointer :: ncol(:), l_ptr(:), l_col(:)
 
278
    integer :: n_r
 
279
    logical, allocatable :: log_r(:)
276
280
    logical :: dense_inserted
277
281
 
278
282
    ! Save the rows ( this is the same for all cases)
283
287
    ncol     => n_col   (in)
284
288
    l_ptr    => list_ptr(in)
285
289
    l_col    => list_col(in)
 
290
    n_r = r%n
286
291
 
287
292
    ! We no a check of arguments before anything...
288
293
    if ( any(r%r < 1) .or. any(n_rows_g < r%r) ) then
295
300
 
296
301
    ! When creating the sparsity, we search for the region
297
302
    ! several times, this is slow if it isn't sorted.
298
 
    call rgn_copy(r,sr)
299
 
    call rgn_sort(sr)
 
303
    allocate(log_r(n_rows_g))
 
304
    log_r(:) = .false.
 
305
    call rgn_2logical(r, log_r)
300
306
 
301
307
    ! We first check whether the parts are contained in the sparsity
302
308
    ! pattern.
304
310
    do lir = 1 , n_rows
305
311
       ! Retrieve the global row-index
306
312
       ir = index_local_to_global(dit,lir,Node)
307
 
       if ( in_rgn(sr,ir) ) then
 
313
       if ( log_r(ir) ) then
308
314
          dense_inserted = .true.
309
315
          exit
310
316
       end if
312
318
 
313
319
    ! if the dense part is not present, we can easily return
314
320
    if ( .not. dense_inserted ) then
315
 
       out = in
316
 
       call rgn_delete(sr)
317
 
       return
 
321
      out = in
 
322
      deallocate(log_r)
 
323
      return
318
324
    end if
319
325
 
320
326
    ! Prepare creation of num and listptr arrays
321
 
    allocate(num    (n_rows))
 
327
    allocate(num(n_rows))
322
328
    allocate(listptr(n_rows))
323
329
 
324
330
    ! The list pointer for the first entry is always the "first"
332
338
       ir = index_local_to_global(dit,lir,Node)
333
339
 
334
340
       ! The easy part is if we are not in the row containing the dense part
335
 
       if ( .not. in_rgn(sr,ir) ) then
 
341
       if ( .not. log_r(ir) ) then
336
342
          num(lir) = ncol(lir)
337
343
       else
338
344
          ! We assume a block-cyclic distribution
339
345
          ! Which means that the current row has the dense part
340
346
          ! We know all the dense elements, then we only need to
341
347
          ! count the non-dense
342
 
          num(lir) = sr%n
 
348
          num(lir) = n_r
343
349
          ! Retrieve the pointer to the original sparsity
344
350
          do i = l_ptr(lir) + 1 , l_ptr(lir) + ncol(lir)
345
351
             ! we check whether the sparse element is
346
352
             ! outside the dense part (in that case we need
347
353
             ! to count it)
348
 
             if ( .not. in_rgn(sr,l_col(i)) ) then
 
354
             if ( .not. log_r(l_col(i)) ) then
349
355
                num(lir) = num(lir) + 1
350
356
             end if
351
357
          end do
391
397
          do i = l_ptr(lir) + 1 , l_ptr(lir) + ncol(lir)
392
398
             ! When we are on either side of the dense part
393
399
             ! we simply copy the column index
394
 
             if ( .not. in_rgn(sr,l_col(i)) ) then
 
400
             if ( .not. log_r(l_col(i)) ) then
395
401
                j = j + 1
396
402
                list(iu+j) = l_col(i)
397
403
             else if ( .not. dense_inserted ) then
399
405
                ! we are within the dense part of the array
400
406
                dense_inserted = .true.
401
407
                ! insert the dense format
402
 
                do l = 1 , sr%n
 
408
                do l = 1 , n_r
403
409
                   j = j + 1
404
 
                   list(iu+j) = sr%r(l)
 
410
                   list(iu+j) = r%r(l)
405
411
                end do
406
412
                
407
413
             end if
413
419
          ! back of the list arry with the dense matrix elements
414
420
          if ( .not. dense_inserted ) then
415
421
             ! insert the dense format
416
 
             do l = 1 , sr%n
 
422
             do l = 1 , n_r
417
423
                j = j + 1
418
 
                list(iu+j) = sr%r(l)
 
424
                list(iu+j) = r%r(l)
419
425
             end do
420
426
          end if
421
427
          
426
432
       end if
427
433
    end do
428
434
 
429
 
    ! Clean-up sorted region
430
 
    call rgn_delete(sr)
431
 
 
 
435
    deallocate(log_r)
 
436
    
432
437
    call newSparsity(out,n_rows,n_rows_g,n_nzs,num,listptr,list, &
433
438
         name='(DU of: '//name(in)//')', &
434
439
         ncols=ncols(in),ncols_g=ncols_g(in))