File: SparseMatrix.F90

package info (click to toggle)
paraview 5.1.2%2Bdfsg1-2
  • links: PTS, VCS
  • area: main
  • in suites: stretch
  • size: 221,108 kB
  • ctags: 236,092
  • sloc: cpp: 2,416,026; ansic: 190,891; python: 99,856; xml: 81,001; tcl: 46,915; yacc: 5,039; java: 4,413; perl: 3,108; sh: 1,974; lex: 1,926; f90: 748; asm: 471; pascal: 228; makefile: 198; objc: 83; fortran: 31
file content (177 lines) | stat: -rw-r--r-- 6,074 bytes parent folder | download | duplicates (6)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
module SparseMatrix
  ! this uses the MKL sparse matrix storage scheme. it assumes we have a
  ! node-based partitioning of the grid.
  implicit none
  public :: initialize, insert, finalize, allreducevector, matvec

  type SparseMatrixData
     real(kind=8), DIMENSION(:), allocatable :: values
     integer, DIMENSION(:), allocatable :: columns, rowindices, mapping
     integer :: globalsize
  end type SparseMatrixData

contains
  subroutine initialize(sm, ownedbox, dimensions)
    type(SparseMatrixData), intent(inout) :: sm
    integer, intent(in) :: ownedbox(6), dimensions(3)
    integer :: numpoints, i, j, k, numx, numy, numz, counter, rowcounter, globalindex
    logical :: negativey, positivey, negativez, positivez
    integer :: allocatestatus

    numpoints = 1
    numx = ownedbox(2)-ownedbox(1)+1
    if(numx < 1) numx = 1
    numy = ownedbox(4)-ownedbox(3)+1
    if(numy < 1) numy = 1
    numz = ownedbox(6)-ownedbox(5)+1
    if(numz < 1) numz = 1
    numpoints = numx*numy*numz

    ! over-allocate both values and columns
    allocate(sm%rowindices(numpoints+1), sm%columns(7*numpoints), sm%values(7*numpoints), STAT = allocatestatus)
    if (allocatestatus /= 0) STOP "*** SparseMatrix.F90: Not enough memory for arrays ***"

    sm%values(:) = 0.d0
    sm%columns(:) = -1

    ! a mapping from the global row index to the local row index
    sm%globalsize = dimensions(1)*dimensions(2)*dimensions(3)
    allocate(sm%mapping(sm%globalsize))
    sm%mapping(:) = -1 ! initialize with bad values

    rowcounter = 1
    counter = 1
    do k=1, numz
       negativez = .TRUE.
       if(k .eq. 1 .and. ownedbox(5) .eq. 1) negativez = .FALSE.
       positivez = .TRUE.
       if(k .eq. numz .and. ownedbox(6) .eq. dimensions(3)) positivez = .FALSE.
       do j=1, numy
          negativey = .TRUE.
          if(j .eq. 1 .and. ownedbox(3) .eq. 1) negativey = .FALSE.
          positivey = .TRUE.
          if(j .eq. numy .and. ownedbox(4) .eq. dimensions(2)) positivey = .FALSE.
          do i=1, numx
             globalindex = (ownedbox(5)+k-2)*dimensions(1)*dimensions(2)+(ownedbox(3)+j-2)*dimensions(1)+ownedbox(1)+i-1
             sm%rowindices(rowcounter) = counter
             sm%mapping(globalindex) = rowcounter
             rowcounter = rowcounter + 1
             if(negativez) then
                sm%columns(counter) = globalindex-dimensions(1)*dimensions(2)
                counter = counter + 1
             endif
             if(negativey) then
                sm%columns(counter) = globalindex-dimensions(1)
                counter = counter + 1
             endif
             if(i .ne. 1 .or. ownedbox(1) .ne. 1) then
                sm%columns(counter) = globalindex-1
                counter = counter + 1
             endif
             sm%columns(counter) = globalindex
             counter = counter + 1
             if(i .ne. numx .or. ownedbox(2) .ne. dimensions(1)) then
                sm%columns(counter) = globalindex+1
                counter = counter + 1
             endif
             if(positivey) then
                sm%columns(counter) = globalindex+dimensions(1)
                counter = counter + 1
             endif
             if(positivez) then
                sm%columns(counter) = globalindex+dimensions(1)*dimensions(2)
                counter = counter + 1
             endif
          end do
       end do
    end do
    ! the final value for ending the row
    sm%rowindices(rowcounter) = counter
  end subroutine initialize

  subroutine finalize(sm)
    type(SparseMatrixData), intent(inout) :: sm
    if(allocated(sm%values)) deallocate(sm%values)
    if(allocated(sm%columns)) deallocate(sm%columns)
    if(allocated(sm%rowindices)) deallocate(sm%rowindices)
    if(allocated(sm%mapping)) deallocate(sm%mapping)
  end subroutine finalize

  subroutine insert(sm, row, column, value)
    type(SparseMatrixData), intent(inout) :: sm
    integer, intent(in) :: row, column
    real(kind=8), intent(in) :: value
    integer :: i, rowindex

    rowindex = sm%mapping(row)

    do i=sm%rowindices(rowindex), sm%rowindices(rowindex+1)-1
       if(sm%columns(i) .eq. column) then
          sm%values(i) = value
          return
       endif
    end do

    write(*,*) 'SparseMatrix: bad indices for insert() ', row, column
  end subroutine insert

  subroutine get(sm, row, column, value)
    type(SparseMatrixData), intent(in) :: sm
    integer, intent(in) :: row, column
    real(kind=8), intent(out) :: value
    integer :: i, rowindex

    rowindex = sm%mapping(row)

    do i=sm%rowindices(rowindex), sm%rowindices(rowindex+1)-1
       if(sm%columns(i) .eq. column) then
          value = sm%values(i)
          return
       endif
    end do
    write(*,*) 'SparseMatrix: bad indices for get() ', row, column
  end subroutine get

  subroutine matvec(sm, x, b)
    !Ax=b
    type(SparseMatrixData), intent(inout) :: sm
    real(kind=8), intent(in) :: x(:)
    real(kind=8), intent(out) :: b(:)
    integer :: i, j

    do i=1, sm%globalsize
       b(i) = 0.d0
       if(sm%mapping(i) .gt. 0) then
          do j=sm%rowindices(sm%mapping(i)), sm%rowindices(sm%mapping(i)+1)-1
             b(i) = b(i) + sm%values(j)*x(sm%columns(j))
          end do
       endif
    end do

    call allreducevector(sm%globalsize, b)
  end subroutine matvec

  subroutine allreducevector(globalsize, vec)
    implicit none
    include 'mpif.h'
    integer, intent(in) :: globalsize
    real(kind=8), intent(inout) :: vec(:)
    real(kind=8), DIMENSION(:), allocatable :: temp
    integer :: ierr, i, numtasks, allocatestatus

    call mpi_comm_size(MPI_COMM_WORLD, numtasks, ierr)
    if(numtasks .eq. 1) then
       return
    endif

    allocate(temp(globalsize), STAT = allocatestatus)
    if (allocatestatus /= 0) STOP "*** SparseMatrix.F90: Not enough memory for temp array ***"
    do i=1, globalsize
       temp(i) = vec(i)
    enddo
    call mpi_allreduce(temp, vec, globalsize, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, ierr)
    deallocate(temp)

  end subroutine allreducevector

end module SparseMatrix