1
! @@LICENSE@@ see Copyright notice in the top-directory
3
! This module may be controlled via these
4
! preprocessor variables:
5
! - BUD_MOD_NAME name of the module
6
! - BUD_TYPE_NAME name of the public type
7
! - BUD_TYPE_NAME_ name of the internal data pointer
8
! - BUD_TYPE_NAME_STR the routine name in "string" format (for IO)
10
! Additionally the CSR matrix format may be
11
! controlled via these flags:
13
! case ( 0 ) ! default
14
! standard CSR format (MKL-SpBLAS)
16
! zero based pointer CSR format
17
! The row-pointer is zero based
18
! This will have an additional array with
19
! the number of non-zero elements per row
20
! (this is equivalent to the SIESTA sparsity)
21
! This does not allow C-interoperability (do NOT set
22
! BUD_SM_INTEROP_C in this case)
24
! - BUD_SM_INTEROP_C=0|1
25
! Use C-indexing with BUD_SM_INTEROP_C == 1
26
! All subsequent operations correspond to a
28
! Cannot be used together with BUD_SM_CSR == 1
34
! Retrieve the settings for the sparse matrix
39
# if BUD_SM_INTEROP_C == 1
40
# error "SM_INTEROP_C may not be defined when using SM_CSR == 1"
43
# error "SM_CSR *MUST* be either 0 or 1"
47
#include "bud_utils.inc"
49
use BUD_CC3(BUD_MOD,_,SM_common)
51
! This *MUST* be the first statement
52
! Common parameters and variables used
53
# include "bud_common_declarations.inc"
55
!> Sparse matrix type in the CSR format
58
!> @cond BUD_DEVELOPER
60
!> Stored pointer which contains the reference counting etc.
61
type(BUD_TYPE_NAME_), pointer :: D => null()
63
!> @endcond BUD_DEVELOPER
65
# include "bud_common_type.inc"
66
#if BUD_FORTRAN >= 2003
68
!> @name Private procedures
72
procedure, private :: new_dim_
73
procedure, private :: new_copy_
75
procedure, private :: column_p_
76
procedure, private :: column_rp_
78
procedure, private :: remove_row_el_
79
procedure, private :: remove_row_list_
81
procedure, private :: remove_col_el_
82
procedure, private :: remove_col_list_
84
procedure, private :: translate_row_el_
85
procedure, private :: translate_row_list_
87
procedure, private :: translate_col_el_
88
procedure, private :: translate_col_list_
93
generic, public :: new => new_dim_, new_copy_
96
procedure, public :: nonzeros => nonzeros_
98
!> @iSee #max_nonzeros
99
procedure, public :: max_nonzeros => max_nonzeros_
102
procedure, public :: rows => rows_
105
procedure, public :: columns => columns_
108
procedure, public :: size => size_
111
procedure, public :: index => index_
114
procedure, public :: offset_p => offset_p_
117
procedure, public :: nrow_p => nrow_p_
120
generic, public :: column_p => column_p_, column_rp_
123
!> @iSee #add_element
124
procedure, public :: add_element => add_element_
128
generic, public :: remove_row => remove_row_el_, &
131
!> @iSee #remove_column
132
generic, public :: remove_column => remove_col_el_, &
135
!> @iSee #translate_row
136
generic, public :: translate_row => translate_row_el_, &
139
!> @iSee #translate_column
140
generic, public :: translate_column => translate_col_el_, &
145
procedure, public :: sort => sort_
148
procedure, public :: equivalent => equivalent_
151
procedure, public :: finalize => finalize_
154
procedure, public :: attach => attach_
157
end type BUD_TYPE_NAME
160
!> @cond BUD_DEVELOPER
162
! These fields are used in the sparse matrix stuff.
163
integer(BUD_TYPE_VAR_PREC), parameter :: ONE = BUD_CC2(1_,BUD_TYPE_VAR_PREC)
164
integer(BUD_TYPE_VAR_PREC), parameter :: ZERO = BUD_CC2(0_,BUD_TYPE_VAR_PREC)
166
!> @bud container for BUD_TYPE_NAME
168
!! Contains the sparsity pattern for a CSR matrix.
171
!> Number of rows in the matrix
172
integer(BUD_TYPE_VAR_PREC) :: nr = ZERO
173
!> Number of columns in the matrix
174
integer(BUD_TYPE_VAR_PREC) :: nc = ZERO
175
!> Number of non-zero elements in the sparse matrix
176
integer(BUD_TYPE_VAR_PREC) :: nz = ZERO
177
!> Total number of elements possible in the sparse matrix
178
integer(BUD_TYPE_VAR_PREC) :: nt = ZERO
180
!> Whether the sparse matrix has been sorted, see #sort
181
logical :: sorted = .false.
183
!> Whether the sparse matrix has been finalized, see #finalize
184
logical :: finalized = .false.
186
!> Index of the equivalent row (size `nr+1`)
187
integer(BUD_TYPE_VAR_PREC), allocatable :: rptr(:)
189
!> Number of non-zero elements per row
190
integer(BUD_TYPE_VAR_PREC), allocatable :: nrow(:)
192
!> The column index of the equivalent sparse matrix (size `nz`)
193
integer(BUD_TYPE_VAR_PREC), allocatable :: col(:)
196
! Consistent data in the reference counted object
197
# include "bud_common_type_.inc"
199
end type BUD_TYPE_NAME_
201
!> @endcond BUD_DEVELOPER
204
!> Create a new sparse matrix
206
!! You may either create an empty sparse matrix with
207
!! a fixed size, or create a sparse matrix from
210
module procedure new_dim_
211
module procedure new_copy_
216
!> Retrieve a pointer to the row offsets/pointers
218
!! A pointer with the row offsets.
221
# if BUD_SM_INTEROP_C == 0
222
!! `M(ir,col(pointer(ir)))`
224
!! `M(ir,col(pointer(ir)+1))`
226
!! is the first sparse element in row `ir`.
227
#elif BUD_SM_CSR == 1
228
!! `M(ir,col(pointer(ir)+1))` is the first sparse
229
!! element in row `ir`.
232
module procedure offset_p_
237
!> Retrieve a pointer to the column indices
239
module procedure column_p_
240
module procedure column_rp_
244
!> Retrieve a pointer to the number of entries per row
246
module procedure nrow_p_
250
!> Retrieve the sparse index for the row and column (-1 if non existing)
252
module procedure index_
257
!> Remove rows from the sparse matrix
259
module procedure remove_row_el_
260
module procedure remove_row_list_
264
!> Remove columns from the sparse matrix
265
interface remove_column
266
module procedure remove_col_el_
267
module procedure remove_col_list_
269
public :: remove_column
272
!> Translate rows in the sparse matrix to new rows
274
!! This will essentially change all entries for given
275
!! rows into a new set of rows.
276
!! Furthermore, by denoting the new row with a negative
277
!! index the element will be deleted.
278
interface translate_row
279
module procedure translate_row_el_
280
module procedure translate_row_list_
282
public :: translate_row
285
!> Translate columns in the sparse matrix to new columns
287
!! This will essentially change all entries for given
288
!! columns into a new set of columns.
289
!! Furthermore, by denoting the new column with a negative
290
!! index the column will be deleted.
291
interface translate_column
292
module procedure translate_col_el_
293
module procedure translate_col_list_
295
public :: translate_column
298
!> Sorts columns in the sparse matrix
300
!! Sorts the sparse matrix such that the column
301
!! index is always increasing for each row.
302
!! This will generally allow faster access patterns.
304
!! One may query the error of the object to check
305
!! whether the input is correct.
307
module procedure sort_
312
!> Checks whether two sparse matrices have the same elements
314
!! This routine returns `.true.` if the two
315
!! sparse matrices have all the same elements in the
318
!! In case it is the same object it returns immediately
321
module procedure equivalent_
326
!> Finalize the sparsity pattern by removing non-used elements
328
!! After finalization the number of non-zero elements
329
!! and the total number of possible elements are the same.
330
!! I.e. this routine removes all non-used elements in the
333
!! @note this routine re-creates the sparse matrix,
334
!! so any pointers to the internal data-structure
335
!! must be re-instantiated before use.
337
module procedure finalize_
342
!> Adds a non-zero element to the sparse matrix
344
!! Add a non-zero element to the sparse matrix.
345
!! In case there is not enough space the sparse matrix
346
!! will be re-allocated and copied to it-self.
347
!! Adding an element will thus not necessarily preserve
348
!! the allocated elements and any pointers to the
349
!! data contained needs to be updated.
351
!! If the error is `0` there was no need to
352
!! extend the sparse matrix data.
353
!! If the error is `-1` the element has not been added
354
!! If the error is above `0` the sparse matrix
355
!! has been re-allocated and the element is contained.
356
!! The re-allocation ensures that all rows can
357
!! contain `maxval(nrow)+8` and can thus result in
358
!! a large memory increase.
359
interface add_element
360
module procedure add_element_
362
public :: add_element
365
! Include the common elements of a sparsity method
366
! This includes "bud_common.inc"
367
# include "SM_common.inc"
370
!> @cond BUD_DEVELOPER
372
!> Internal routine for cleaning up the data container.
375
!! This routine is only used internally to clean-up
376
!! any data in the type.
377
!! Should never be made public.
378
subroutine delete_(this)
379
type(BUD_TYPE_NAME), intent(inout) :: this
380
integer :: stat, istat
386
this%D%sorted = .false.
387
this%D%finalized = .false.
389
! Currently we do not allow external memory
391
if ( .not. allocated(this%D%col) ) return
393
deallocate(this%D%rptr, stat=istat)
394
if ( 0 /= istat ) stat = istat
395
deallocate(this%D%col, stat=stat)
396
if ( 0 /= istat ) stat = istat
397
deallocate(this%D%nrow, stat=stat)
398
if ( 0 /= istat ) stat = istat
400
end subroutine delete_
402
!> @endcond BUD_DEVELOPER
405
!> @param[in] from the original `bud` which is copied to `to`
406
!! @param[inout] to the output `bud` with the full copied data
407
!! @param[in] dealloc @opt=.true. whether `to` should be de-allocated before
408
subroutine copy_(from, to, dealloc)
409
BUD_CLASS(BUD_TYPE_NAME), intent(in) :: from
410
BUD_CLASS(BUD_TYPE_NAME), intent(inout) :: to
411
logical, intent(in), optional :: dealloc
413
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: fptr(:), nrow(:), fcol(:)
414
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: tptr(:)
415
integer(BUD_TYPE_VAR_PREC) :: i, ir, nr
417
call set_error(to, 0)
421
if ( present(dealloc) ) ldealloc = dealloc
422
if ( ldealloc ) call delete(to)
424
! quick return if there is nothing to copy
425
if ( .not. is_initd(from) ) return
427
if ( is_initd(to) ) then
429
! Check if to has room
430
if ( rows(from) /= rows(to) .or. &
431
columns(from) /= columns(to) ) then
433
! We do not have the same shape
434
call set_error(to, -1)
439
! They have the same dimensions
440
! Check there is room
441
call attach(from, nr=nr, rptr=fptr, nrow=nrow, col=fcol)
442
call attach(to, rptr=tptr)
446
if ( nrow(ir) > tptr(ir+ONE) - tptr(ir) ) then
447
call set_error(to, ir)
453
! quick return if there is not room
454
! The error message will return the column
455
! where there is not enough room
456
if ( error(to) /= 0 ) return
458
! Now reset all elements in `to`
459
call attach(to, nrow=nrow)
464
! Loop all elements in from and copy them to
468
fcol => column_p(from, ir)
471
call add_element(to, ir, fcol(i))
478
! We make a copy by creating it
479
call new(to, from%D%nr, from%D%nc, from%D%nz, &
480
from%D%rptr, from%D%col, from%D%nrow)
481
to%D%sorted = from%D%sorted
487
!> @param[inout] this the sparse matrix
488
!! @param[in] nr number of rows of the matrix
489
!! @param[in] nc number of columns of the matrix
490
!! @param[in] nt @opt total number of possible non-zero elements
491
!! @param[in] npr @opt number of non-zero elements per row (defaults to 10 if `nt` is not specified)
492
subroutine new_dim_(this, nr, nc, nt, npr)
493
BUD_CLASS(BUD_TYPE_NAME), intent(inout) :: this
494
integer(BUD_TYPE_VAR_PREC), intent(in) :: nr, nc
495
integer(BUD_TYPE_VAR_PREC), intent(in), optional :: nt, npr
496
integer(BUD_TYPE_VAR_PREC) :: r
498
call initialize(this)
502
if ( present(nt) ) then
504
else if ( present(npr) ) then
507
this%D%nt = nr * BUD_CC2(10_,BUD_TYPE_VAR_PREC)
510
this%D%sorted = .false.
511
this%D%finalized = .false.
513
allocate(this%D%rptr(nr+1))
514
allocate(this%D%nrow(nr))
515
allocate(this%D%col(this%D%nt))
517
! Initialize sparse matrix
518
if ( present(npr) ) then
519
this%D%rptr(1) = BUD_SM_PTR
521
this%D%rptr(r) = this%D%rptr(r-ONE) + npr
523
this%D%rptr(nr+1) = this%D%nt + BUD_SM_PTR
525
this%D%rptr(1) = BUD_SM_PTR
528
this%D%nrow(r) = ZERO
531
end subroutine new_dim_
534
!> @param[inout] this the new sparse matrix
535
!! @param[in] nr number of rows of the matrix
536
!! @param[in] nc number of columns of the matrix
537
!! @param[in] nz number of non-zero elements of the matrix
538
!! @param[in] rptr row pointers (at least size `nr`)
539
!! @param[in] col column indices for the sparse elements
540
!! @param[in] nrow @opt number of column elements per row
541
subroutine new_copy_(this, nr, nc, nz, rptr, col, nrow)
542
BUD_CLASS(BUD_TYPE_NAME), intent(inout) :: this
544
integer(BUD_TYPE_VAR_PREC), intent(in) :: nr, nc, nz
546
integer(BUD_TYPE_VAR_PREC), intent(in) :: rptr(nr), col(nz)
547
integer(BUD_TYPE_VAR_PREC), intent(in), optional :: nrow(nr)
549
integer(BUD_TYPE_VAR_PREC) :: r
551
! Copy over information
552
! It must be based on the starting index
553
if ( rptr(1) /= BUD_SM_PTR ) then
559
call new(this, nr, nc, nz)
564
this%D%rptr(r) = rptr(r)
567
! create last pointer (to one plus number of elements)
568
! This ensures simple loops without taking care of
570
this%D%rptr(this%D%nr+1) = this%D%nz + BUD_SM_PTR
572
if ( present(nrow) ) then
574
this%D%nrow(r) = nrow(r)
578
this%D%nrow(r) = this%D%rptr(r+1) - this%D%rptr(r)
583
this%D%col(r) = col(r)
586
end subroutine new_copy_
589
!> @param[in] this sparse matrix
590
!! @return a pointer to the row offsets for the sparse matrix (contiguous)
591
function offset_p_(this) result(rptr)
592
BUD_CLASS(BUD_TYPE_NAME), intent(in) :: this
593
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: rptr(:)
597
end function offset_p_
599
!> @param[in] this sparse matrix
600
!! @return a pointer to the column indices for the sparse matrix (contiguous)
601
function column_p_(this) result(col)
602
BUD_CLASS(BUD_TYPE_NAME), intent(in) :: this
603
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: col(:)
607
end function column_p_
609
!> @param[in] this sparse matrix
610
!> @param[in] r only retrieve the columns that reside in row `r`
611
!! @return a pointer to the column indices for the sparse matrix of row `r` (contiguous)
612
function column_rp_(this, r) result(col)
613
BUD_CLASS(BUD_TYPE_NAME), intent(in) :: this
614
integer(BUD_TYPE_VAR_PREC), intent(in) :: r
615
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: col(:)
617
col => this%D%col(this%D%rptr(r)BUD_SM_PTR_A:this%D%rptr(r)+this%D%nrow(r)BUD_SM_PTR_B)
619
end function column_rp_
621
!> @param[in] this sparse matrix
622
!! @return a pointer to the array that holds the number of entries per row
623
function nrow_p_(this) result(nrow)
624
BUD_CLASS(BUD_TYPE_NAME), intent(in) :: this
625
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: nrow(:)
632
!> @param[in] this sparse matrix @bud
633
!! @param[out] nr @opt number of rows in SM
634
!! @param[out] nc @opt number of columns in SM
635
!! @param[out] nz @opt number of non-zero elements in SM
636
!! @param[out] nt @opt total number of possible non-zero elements in SM
637
!! @param[out] rptr @opt row pointer (`rptr(2)BUD_SM_PTR_A` is starting index of `ir=2`)
638
!! @param[out] nrow @opt number of non-zero elements per row
639
subroutine attach_(this, D, nr, nc, nz, nt, rptr, col, nrow)
640
BUD_CLASS(BUD_TYPE_NAME), intent(in) :: this
641
! Having this forces the user to explicitly specify the
642
! wanted information.
643
logical, intent(in), optional :: D
644
integer(BUD_TYPE_VAR_PREC), intent(out), optional :: nr, nc, nz, nt
645
integer(BUD_TYPE_VAR_PREC), intent(out), pointer BUD_FORTRAN_CONTIGUOUS, optional :: rptr(:), col(:)
646
integer(BUD_TYPE_VAR_PREC), intent(out), pointer BUD_FORTRAN_CONTIGUOUS, optional :: nrow(:)
648
if ( is_initd(this) ) then
649
if ( present(nr) ) nr = this%D%nr
650
if ( present(nc) ) nc = this%D%nc
651
if ( present(nz) ) nz = this%D%nz
652
if ( present(nt) ) nt = this%D%nt
653
if ( present(rptr) ) rptr => this%D%rptr
654
if ( present(col) ) col => this%D%col
655
if ( present(nrow) ) nrow => this%D%nrow
657
if ( present(nr) ) nr = ZERO
658
if ( present(nc) ) nc = ZERO
659
if ( present(nz) ) nz = ZERO
660
if ( present(nt) ) nt = ZERO
661
if ( present(rptr) ) nullify(rptr)
662
if ( present(col) ) nullify(col)
663
if ( present(nrow) ) nullify(nrow)
666
end subroutine attach_
669
!> @param[in] this the sparse matrix (sorted, @isee #sp_sort)
670
#if BUD_SM_CSR == 0 && BUD_SM_INTEROP_C == 1
671
!! @param[in] ir row index (0-based)
672
!! @param[in] ic column index (0-based)
673
!! @return the sparse index of `M(ir,ic)`, `<0` if `M(ir,ic) = 0` (0-based)
675
!! @param[in] ir row index (1-based)
676
!! @param[in] ic column index (1-based)
677
!! @return the sparse index of `M(ir,ic)`, `<0` if `M(ir,ic) = 0` (1-based)
679
pure function index_(this, ir, ic) result(idx)
680
use BUD_CC2(BUD_MOD, _utils), only: find_bin
681
BUD_CLASS(BUD_TYPE_NAME), intent(in) :: this
682
integer(BUD_TYPE_VAR_PREC), intent(in) :: ir
683
integer(BUD_TYPE_VAR_PREC), intent(in) :: ic
684
integer(BUD_TYPE_VAR_PREC) :: idx
686
if ( .not. this%D%sorted ) then
688
do idx = this%D%rptr(ir) BUD_SM_PTR_A , this%D%rptr(ir) + this%D%nrow(ir) BUD_SM_PTR_B
689
if ( this%D%col(idx) == ic ) return
698
call find_bin(this%D%nrow(ir), &
699
this%D%col(this%D%rptr(ir) BUD_SM_PTR_A:), ic, idx)
700
if ( idx > ZERO ) idx = this%D%rptr(ir) BUD_SM_PTR_B + idx
706
!> @param[inout] this sparse matrix to sort (in-place)
707
!! @param[out] pvt @opt=@none if requested the pivoting array for the sorted sparsity pattern
708
subroutine sort_(this, pvt)
709
! We use the quick-sort algorithm in this general finalization
711
use BUD_CC2(BUD_MOD, _utils), only: sort_quick
713
BUD_CLASS(BUD_TYPE_NAME), intent(inout) :: this
714
integer(BUD_TYPE_VAR_PREC), intent(out), target, optional :: pvt(:)
717
integer(BUD_TYPE_VAR_PREC) :: ir, nr, nz, ptr, i
718
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: rptr(:), nrow(:), col(:), tvt(:)
721
call set_error(this, 0)
723
! Get sparse matrix array
724
call attach(this, nr=nr, nz=nz, rptr=rptr, col=col, nrow=nrow)
726
if ( present(pvt) ) then
728
if ( size(pvt) < nz ) then
729
call set_error(this, SM_INPUT + 3)
737
ptr = rptr(ir) - BUD_SM_PTR
742
call sort_quick(nz, col(ptr+1:), tvt)
744
tvt(i) = ptr + tvt(i)
747
col(ptr + i) = col(tvt(i))
756
col => column_p(this, ir)
759
call sort_quick(nz, col(1:))
765
this%D%sorted = .true.
770
!> @param[inout] a first sparse matrix
771
!! @param[inout] b second sparse matrix
772
!! @return whether they have the same entries
773
function equivalent_(a, b) result(is)
774
BUD_CLASS(BUD_TYPE_NAME), intent(inout) :: a, b
777
integer(BUD_TYPE_VAR_PREC) :: i, nr
779
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: arow(:), acol(:)
780
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: brow(:), bcol(:)
783
if ( same(a, b) ) return
785
is = rows(a) == rows(b) .and. &
786
columns(a) == columns(b)
787
if ( .not. is ) return
790
call attach(a, nr=nr, nrow=arow)
791
call attach(b, nrow=brow)
795
! Check number of elements per row
796
is = arow(i) == brow(i)
797
if ( .not. is ) return
799
! check that there are elements
800
if ( arow(i) == 0 ) cycle
802
acol => column_p(a, i)
803
bcol => column_p(b, i)
805
is = all(acol == bcol)
806
if ( .not. is ) return
810
end function equivalent_
813
!> @param[inout] this sparse matrix to finalize (in-place)
814
subroutine finalize_(this)
815
use BUD_CC3(BUD_MOD,_,BUD_LIST_NAME)
817
BUD_CLASS(BUD_TYPE_NAME), intent(inout) :: this
820
integer(BUD_TYPE_VAR_PREC) :: ir, nr, nc, nz
821
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: rptr(:), col(:)
823
type(BUD_LIST_NAME) :: lptr, lcol
824
integer(BUD_TYPE_VAR_PREC) :: p
828
call set_error(this, 0)
830
! Get sparse matrix array
831
call attach(this, nr=nr, nc=nc, nz=nz, col=col)
837
! Initialize the pointer
843
col => column_p(this, ir)
846
! Push the new pointer index
850
! Push the column-indices to the list
851
call push(lcol, nz, col)
855
! Determine whether we should sort it after-wards
856
! I.e. maintain state of object
857
sorted = is_sorted(this)
863
call new(this, nr, nc, p, rptr, col)
870
this%D%finalized = .true.
872
end subroutine finalize_
875
!> @param[inout] this the sparse matrix to add an element to
876
!! @param[in] ir the row index of the element
877
!! @param[in] ic the column index of the element
878
!! @param[in] dry_run @opt=.false., if `.true.` will do nothing but issue error messages as though it had runned.
879
recursive subroutine add_element_(this, ir, ic, dry_run)
880
BUD_CLASS(BUD_TYPE_NAME), intent(inout) :: this
881
integer(BUD_TYPE_VAR_PREC), intent(in) :: ir, ic
882
logical, intent(in), optional :: dry_run
885
type(BUD_TYPE_NAME) :: nthis
886
integer(BUD_TYPE_VAR_PREC) :: i, r, ix, nr, nc, npr
888
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: ptr(:), nrow(:), col(:)
889
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: ptr_(:), nrow_(:), col_(:)
892
if ( present(dry_run) ) ldry = dry_run
894
call attach(this, nr=nr, nc=nc, &
895
rptr=ptr, nrow=nrow, col=col)
897
! check if there is place in the sparse matrix
898
! for an entry without re-allocation
899
if ( ptr(ir+ONE) - ptr(ir) > nrow(ir) ) then
901
ix = ptr(ir) BUD_SM_PTR_A +nrow(ir)
902
if ( nrow(ir) == 0 ) then
907
if ( .not. ldry ) then
909
nrow(ir) = nrow(ir) + ONE
910
this%D%nz = this%D%nz + ONE
912
if ( this%D%sorted .and. ic < i ) then
914
! denote it non-sorted
915
this%D%sorted = .false.
920
call set_error(this, 0)
926
! Re-allocate the sparse matrix and copy it
927
! We add 8 elements per row.
928
! Perhaps this should be defined by the user.
929
npr = maxval(nrow) + BUD_CC2(8_,BUD_TYPE_VAR_PREC)
930
! Quick return if user request dry-run
932
call set_error(this, npr)
936
call new(nthis, nr, nc, npr=npr)
937
call attach(nthis, rptr=ptr_, nrow=nrow_, col=col_)
939
! Copy all elements correctly
944
ix = ptr_(r) BUD_SM_PTR_A
946
do i = ptr(r) BUD_SM_PTR_A , ptr(r) + nrow(r) BUD_SM_PTR_B
953
! add the element to the new sparse pattern
954
call add_element(nthis, ir, ic)
955
! copy the sparse pattern to this one
959
! specify the error message which then
960
! is the new size per row
961
call set_error(this, npr)
963
end subroutine add_element_
966
!> @param[in] sp1 the first sparse matrix
967
!! @param[in] sp2 the second sparse matrix
968
!! @param[inout] sp the union of `sp1` and `sp2`
969
subroutine union_(sp1, sp2, sp)
971
use BUD_CC3(BUD_MOD,_,BUD_LIST_NAME)
973
BUD_CLASS(BUD_TYPE_NAME), intent(in) :: sp1, sp2
974
BUD_CLASS(BUD_TYPE_NAME), intent(inout) :: sp
976
type(BUD_LIST_NAME) :: L1, L2, L
977
integer(BUD_TYPE_VAR_PREC) :: ir, n, nr, nc, nz
978
integer(BUD_TYPE_VAR_PREC), allocatable :: rptr(:)
979
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: col(:)
981
! A union removes duplicate entries...
982
! This may be undesired.. :(
983
call attach(sp1, nr=nr, nc=nc)
985
! Create a list with the pointers
991
col => column_p(sp1, ir)
994
col => column_p(sp2, ir)
998
! Sort both lists (makes unions faster)
1003
call union(L1, L2, L2)
1004
! We also ensure it is sorted.
1007
! push the rows to the new list
1010
rptr(ir+1) = rptr(ir) + size(L2)
1015
call new(sp, nr, nc, nz, rptr, col)
1021
end subroutine union_
1024
!> @param[in] from the originating sparse matrix
1025
!! @param[in] n number of elements in `remove`
1026
!! @param[in] remove columns that are removed
1027
!! @param[inout] to sparse matrix after deleting the columns
1028
subroutine remove_col_el_(from, n, remove, to)
1030
use BUD_CC3(BUD_MOD,_,BUD_LIST_NAME)
1032
BUD_CLASS(BUD_TYPE_NAME), intent(in) :: from
1033
integer(BUD_TYPE_VAR_PREC), intent(in) :: n, remove(n)
1034
BUD_CLASS(BUD_TYPE_NAME), intent(inout) :: to
1036
type(BUD_LIST_NAME) :: lrem
1039
call new(lrem, n, remove)
1042
call remove_col_list_(from, lrem, to)
1046
end subroutine remove_col_el_
1048
!> @param[in] from the originating sparse matrix
1049
!! @param[in] remove list with columns that are removed
1050
!! @param[inout] to sparse matrix after deleting the columns
1051
subroutine remove_col_list_(from, remove, to)
1053
use BUD_CC3(BUD_MOD,_,BUD_LIST_NAME)
1055
BUD_CLASS(BUD_TYPE_NAME), intent(in) :: from
1056
BUD_CLASS(BUD_LIST_NAME), intent(in) :: remove
1057
BUD_CLASS(BUD_TYPE_NAME), intent(inout) :: to
1059
type(BUD_LIST_NAME) :: ll
1060
integer(BUD_TYPE_VAR_PREC) :: i, n
1061
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: lst(:)
1071
call translate_col_list_(from, remove, ll, to)
1075
end subroutine remove_col_list_
1078
!> @param[in] from the originating sparse matrix
1079
!! @param[in] nin number of elements in `in_col`
1080
!! @param[in] in_col the set of columns that will be translated into `out_col` (preferentially this should be sorted)
1081
!! @param[in] nout number of elements in `out_col`
1082
!! @param[in] out_col the set of translation columns
1083
!! @param[inout] to the resulting sparse matrix after translating the columns
1084
subroutine translate_col_el_(from, nin, in_col, nout, out_col, to)
1086
use BUD_CC3(BUD_MOD,_,BUD_LIST_NAME)
1088
BUD_CLASS(BUD_TYPE_NAME), intent(in) :: from
1089
integer(BUD_TYPE_VAR_PREC), intent(in) :: nin, in_col(nin), nout, out_col(nout)
1090
BUD_CLASS(BUD_TYPE_NAME), intent(inout) :: to
1092
type(BUD_LIST_NAME) :: lin, lout
1094
call new(lin, nin, in_col)
1095
call new(lout, nout, out_col)
1097
call translate_col_list_(from, lin, lout, to)
1102
end subroutine translate_col_el_
1105
!> @param[in] from the originating sparse matrix
1106
!! @param[in] in_col a list set of columns that will be translated into `out_col` (preferentially this should be sorted)
1107
!! @param[in] out_col a list set of translation columns
1108
!! @param[inout] to the resulting sparse matrix after translating the columns
1109
subroutine translate_col_list_(from, in_col, out_col, to)
1111
use BUD_CC3(BUD_MOD,_,BUD_LIST_NAME)
1113
BUD_CLASS(BUD_TYPE_NAME), intent(in) :: from
1114
BUD_CLASS(BUD_LIST_NAME) :: in_col, out_col
1115
BUD_CLASS(BUD_TYPE_NAME), intent(inout) :: to
1117
type(BUD_LIST_NAME) :: lcol
1118
integer(BUD_TYPE_VAR_PREC) :: nr, nc, nz, ir, i, idx
1119
integer(BUD_TYPE_VAR_PREC), allocatable :: trptr(:)
1120
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: rptr(:), nrow(:), col(:)
1121
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: oc(:)
1123
call attach(from, nr=nr, nc=nc, nz=nz, rptr=rptr, col=col, nrow=nrow)
1126
oc => list_p(out_col)
1128
! Create a list with the columns
1130
allocate(trptr(nr+1))
1132
! Now start create the sparse matrix
1136
! figure out if they should be translated
1137
do i = rptr(ir) BUD_SM_PTR_A, rptr(ir) + nrow(ir) BUD_SM_PTR_B
1138
idx = index(in_col, col(i))
1139
if ( idx <= 0 ) then
1140
call push(lcol, col(i))
1141
else if ( oc(idx) > 0 ) then
1142
! we have a translation
1143
! check if the resulting translated column is
1144
! negative, if so, it means deletion
1146
! we retain the item
1147
call push(lcol, oc(idx))
1151
! update the following pointer
1152
trptr(ir+1) = size(lcol) + BUD_SM_PTR
1157
call new(to, nr, nc, size(lcol), trptr, col)
1161
end subroutine translate_col_list_
1164
!> @param[in] from the originating sparse matrix
1165
!! @param[in] n number of elements in `remove`
1166
!! @param[in] remove rows that are removed
1167
!! @param[inout] to sparse matrix after deleting the rows
1168
subroutine remove_row_el_(from, n, remove, to)
1170
use BUD_CC3(BUD_MOD,_,BUD_LIST_NAME)
1172
BUD_CLASS(BUD_TYPE_NAME), intent(in) :: from
1173
integer(BUD_TYPE_VAR_PREC), intent(in) :: n, remove(n)
1174
BUD_CLASS(BUD_TYPE_NAME), intent(inout) :: to
1176
type(BUD_LIST_NAME) :: lrem
1179
call new(lrem, n, remove)
1182
call remove_row_list_(from, lrem, to)
1186
end subroutine remove_row_el_
1188
!> @param[in] from the originating sparse matrix
1189
!! @param[in] remove list with rows that are removed
1190
!! @param[inout] to sparse matrix after deleting the rows
1191
subroutine remove_row_list_(from, remove, to)
1193
use BUD_CC3(BUD_MOD,_,BUD_LIST_NAME)
1195
BUD_CLASS(BUD_TYPE_NAME), intent(in) :: from
1196
BUD_CLASS(BUD_LIST_NAME), intent(in) :: remove
1197
BUD_CLASS(BUD_TYPE_NAME), intent(inout) :: to
1199
type(BUD_LIST_NAME) :: ll
1200
integer(BUD_TYPE_VAR_PREC) :: i, n
1201
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: lst(:)
1211
call translate_row_list_(from, remove, ll, to)
1215
end subroutine remove_row_list_
1219
!> @param[in] from the originating sparse matrix
1220
!! @param[in] nin number of elements in `in_row`
1221
!! @param[in] in_row the set of rows that will be translated into `out_row` (preferentially this should be sorted)
1222
!! @param[in] nout number of elements in `out_row`
1223
!! @param[in] out_row the set of translation rows
1224
!! @param[inout] to the resulting sparse matrix after translating the rows
1225
subroutine translate_row_el_(from, nin, in_row, nout, out_row, to)
1227
use BUD_CC3(BUD_MOD,_,BUD_LIST_NAME)
1229
BUD_CLASS(BUD_TYPE_NAME), intent(in) :: from
1230
integer(BUD_TYPE_VAR_PREC), intent(in) :: nin, in_row(nin), nout, out_row(nout)
1231
BUD_CLASS(BUD_TYPE_NAME), intent(inout) :: to
1233
type(BUD_LIST_NAME) :: lin, lout
1235
call new(lin, nin, in_row)
1236
call new(lout, nout, out_row)
1238
call translate_row_list_(from, lin, lout, to)
1243
end subroutine translate_row_el_
1245
!> @param[in] from the originating sparse matrix
1246
!! @param[in] in_row a list set of rows that will be translated into `out_row`
1247
!! @param[in] out_row a list set of translation rows (preferentially this should be sorted)
1248
!! @param[inout] to the resulting sparse matrix after translating the rows
1249
subroutine translate_row_list_(from, in_row, out_row, to)
1251
use BUD_CC3(BUD_MOD,_,BUD_LIST_NAME)
1253
BUD_CLASS(BUD_TYPE_NAME), intent(in) :: from
1254
BUD_CLASS(BUD_LIST_NAME) :: in_row, out_row
1255
BUD_CLASS(BUD_TYPE_NAME), intent(inout) :: to
1257
type(BUD_LIST_NAME), allocatable :: lrow(:)
1258
type(BUD_LIST_NAME) :: tmp
1259
integer(BUD_TYPE_VAR_PREC) :: nr, nc, nnr, r, ir, i, idx
1260
integer(BUD_TYPE_VAR_PREC), allocatable :: rptr(:)
1261
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: row(:), col(:)
1262
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: or(:)
1264
call attach(from, nr=nr, nc=nc)
1266
! This algorithm moves in_row to out_row, and possibly deletes
1270
! Initialize the sparse lists
1272
call new(lrow(ir), ONE)
1276
or => list_p(out_row)
1278
! Now actually populate the sparse rows
1281
! Get rows in current column
1282
col => column_p(from, ir)
1283
i = int(size(col), BUD_TYPE_VAR_PREC)
1285
! Figure out if this is in in_row (i.e. should it be translated)
1286
idx = index(in_row, ir)
1289
! The column is not translated
1290
call push(lrow(ir), i, col)
1292
else if ( or(idx) >= BUD_SM_IDX(1) ) then
1294
! Find where the new column has moved
1295
r = BUD_SM_IDXF(or(idx))
1296
call push(lrow(r), i, col)
1302
! Now re-create the new sparsity pattern
1303
! First figure out the new number of rows...
1304
row => list_p(in_row)
1306
do i = 1 , size(in_row)
1307
idx = index(out_row, row(i))
1308
if ( idx <= 0 ) then
1309
! We also delete the row, so we may recognize it
1310
call delete(lrow(row(i)))
1315
! Now re-create the pointers
1316
allocate(rptr(nnr+1))
1317
rptr(1) = BUD_SM_PTR
1322
! if the list has been deleted, simply
1324
if ( .not. is_initd(lrow(r)) ) cycle
1327
! Create full list of columns
1328
call push(tmp, lrow(r))
1330
! Create the pointer
1331
rptr(ir+1) = rptr(ir) + size(lrow(r))
1333
call delete(lrow(r))
1337
! Now we have, rptr, and column
1339
! finally we also know that all elements of lrow
1342
call new(to, nnr, nc, size(tmp), rptr, col)
1347
end subroutine translate_row_list_
1350
!> @param[inout] f `File` bud
1351
!! @param[in] this the sparse matrix bud
1352
subroutine write_(f, this)
1353
use BUD_CC2(BUD_MOD,_File)
1355
BUD_CLASS( BUD_CC2(BUD_TYPE,File) ), intent(inout) :: f
1356
BUD_CLASS(BUD_TYPE_NAME), intent(in) :: this
1358
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: col(:), nrow(:)
1360
logical :: formatted
1363
integer(BUD_TYPE_VAR_PREC) :: nz
1364
integer(BUD_TYPE_VAR_PREC) :: ir
1365
integer(BUD_TYPE_VAR_PREC) :: nr, nc
1367
! If file is not opened, return immediately
1368
if ( .not. is_open(f) ) return
1369
if ( .not. is_initd(this) ) return
1371
! First figure out if the file is an unformatted file
1372
formatted = is_formatted(f)
1375
call attach(this, nr=nr, nc=nc, nz=nz, nrow=nrow)
1377
! First we write the size of the sparse matrix
1378
if ( formatted ) then
1379
write(iu, '(i16)') nr, nc, nz
1380
write(iu, '(l16)') this%D%sorted
1381
write(iu, '(i16)') nrow(:)
1383
write(iu) nr, nc, nz
1384
write(iu) this%D%sorted
1388
! Write the sparse columns
1390
col => column_p(this, ir)
1391
if ( formatted ) then
1392
write(iu, '(i16)') col(:)
1398
end subroutine write_
1400
!> @param[inout] f `File` bud
1401
!! @param[inout] this the sparse matrix bud
1402
subroutine read_(f, this)
1403
use BUD_CC2(BUD_MOD,_File)
1405
BUD_CLASS( BUD_CC2(BUD_TYPE,File) ), intent(inout) :: f
1406
BUD_CLASS(BUD_TYPE_NAME), intent(inout) :: this
1408
integer(BUD_TYPE_VAR_PREC), pointer BUD_FORTRAN_CONTIGUOUS :: ptr(:), col(:), nrow(:)
1410
logical :: formatted
1414
integer(BUD_TYPE_VAR_PREC) :: nz, i
1415
integer(BUD_TYPE_VAR_PREC) :: ir
1416
integer(BUD_TYPE_VAR_PREC) :: nr, nc
1418
! If file is not opened, return immediately
1419
if ( .not. is_open(f) ) return
1421
! First figure out if the file is an unformatted file
1422
formatted = is_formatted(f)
1425
! First we need to read the array dimensions...
1426
if ( formatted ) then
1427
read(iu, '(i16)') nr, nc, nz
1428
read(iu, '(l16)') sorted
1434
call new(this, nr, nc, nz)
1435
call attach(this, rptr=ptr, nrow=nrow)
1438
if ( formatted ) then
1439
read(iu, '(i16)') nrow
1444
! Create the pointer array
1447
ptr(i) = ptr(i-1) + nrow(i-1)
1450
! Read the sparse columns
1452
col => column_p(this, ir)
1453
if ( formatted ) then
1454
read(iu, '(i16)') col(:)
1460
this%D%sorted = sorted
1461
this%D%finalized = .true.
1463
end subroutine read_
1466
! Local pre-processor variables that
1467
! undefine the variables that are not needed anymore.
1469
#undef BUD_LIST_NAME
1470
#undef BUD_TYPE_NAME
1471
#undef BUD_TYPE_NAME_
1472
#undef BUD_TYPE_NAME_STR
1474
#undef BUD_TYPE_VAR_PREC
1478
#undef BUD_SM_INTEROP_C
1483
#include "bud_cleanup.inc"
1486
! project-buds -- local file settings
1487
! Anything below this line may be overwritten by scripts
1488
! Below are non-editable settings
1493
! f90-type-indent: 2
1494
! f90-associate-indent: 2
1495
! f90-continuation-indent: 2
1496
! f90-structure-indent: 2
1497
! f90-critical-indent: 2
1498
! f90-program-indent: 2
added
removed
Src/buds/src/SM_CSR.inc
