Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_sparse_matrix_mp.F90
Go to the documentation of this file.
1 !**************************************************************
2 ! Copyright INRIA
3 ! Authors :
4 ! CALVI project team
5 !
6 ! This code SeLaLib (for Semi-Lagrangian-Library)
7 ! is a parallel library for simulating the plasma turbulence
8 ! in a tokamak.
9 !
10 ! This software is governed by the CeCILL-B license
11 ! under French law and abiding by the rules of distribution
12 ! of free software. You can use, modify and redistribute
13 ! the software under the terms of the CeCILL-B license as
14 ! circulated by CEA, CNRS and INRIA at the following URL
15 ! "http://www.cecill.info".
16 !**************************************************************
17 
18 module sll_m_sparse_matrix_mp
19 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
20 #include "sll_working_precision.h"
21 #include "sll_memory.h"
22 #include "sll_assert.h"
23 
24  use sll_m_qsort_partition, only: &
25  sll_s_qsortc
26 
27  use sll_m_sparse_matrix, only: &
28  sll_t_csr_matrix
29 
30  implicit none
31 
32  public :: &
33  sll_s_csr_matrix_mp_init, &
34  sll_f_new_csr_matrix_mp
35 
36  private
37 
38 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
39 contains
40 
56  function sll_f_new_csr_matrix_mp( &
57  num_rows, &
58  num_cols, &
59  num_patch, &
60  num_elements, &
61  local_to_global_row, &
62  num_local_dof_row, &
63  local_to_global_col, &
64  num_local_dof_col) &
65  result(mat)
66  type(sll_t_csr_matrix), pointer :: mat
67  sll_int32, intent(in) :: num_rows
68  sll_int32, intent(in) :: num_cols
69  sll_int32, intent(in) :: num_patch
70  sll_int32, dimension(:), intent(in) :: num_elements
71  sll_int32, dimension(:, :, :), intent(in) :: local_to_global_row
72  sll_int32, dimension(:), intent(in) :: num_local_dof_row
73  sll_int32, dimension(:, :, :), intent(in) :: local_to_global_col
74  sll_int32, dimension(:), intent(in) :: num_local_dof_col
75  sll_int32 :: ierr
76  sll_allocate(mat, ierr)
77  call sll_s_csr_matrix_mp_init( &
78  mat, &
79  num_rows, &
80  num_cols, &
81  num_patch, &
82  num_elements, &
83  local_to_global_row, &
84  num_local_dof_row, &
85  local_to_global_col, &
86  num_local_dof_col)
87 
88  end function sll_f_new_csr_matrix_mp
89 
104  subroutine sll_s_csr_matrix_mp_init( &
105  mat, &
106  num_rows, &
107  num_cols, &
108  num_patch, &
109  num_elements, &
110  local_to_global_row, &
111  num_local_dof_row, &
112  local_to_global_col, &
113  num_local_dof_col)
114  type(sll_t_csr_matrix), intent(inout) :: mat
115  sll_int32, intent(in) :: num_rows
116  sll_int32, intent(in) :: num_cols
117  sll_int32, intent(in) :: num_patch
118  sll_int32, dimension(:), intent(in) :: num_elements
119  sll_int32, dimension(:, :, :), intent(in) :: local_to_global_row
120  sll_int32, dimension(:), intent(in) :: num_local_dof_row
121  sll_int32, dimension(:, :, :), intent(in) :: local_to_global_col
122  sll_int32, dimension(:), intent(in) :: num_local_dof_col
123  !local variables
124  sll_int32 :: num_nz
125  sll_int32, dimension(:, :), pointer :: lpi_columns
126  sll_int32, dimension(:), pointer :: lpi_occ
127  sll_int32 :: li_coef
128  sll_int32 :: ierr
129  sll_int32 :: max_nen_c
130  !print *,'#num_rows=',num_rows
131  !print *,'#num_nz=',num_nz
132 
133  li_coef = 20
134  ! for we take the maximun of number of element basis in all patches
135  max_nen_c = maxval(num_local_dof_col(:))
136 
137  sll_allocate(lpi_columns(num_rows, 0:li_coef*max_nen_c), ierr)
138  sll_allocate(lpi_occ(num_rows + 1), ierr)
139  lpi_columns(:, :) = 0
140  lpi_occ(:) = 0
141 
142  ! COUNTING NON ZERO ELEMENTS
143  num_nz = sll_count_non_zero_elts_mp( &
144  num_rows, &
145  num_cols, &
146  num_patch, &
147  num_elements, &
148  local_to_global_col, &
149  num_local_dof_col, &
150  local_to_global_row, &
151  num_local_dof_row, &
152  lpi_columns, &
153  lpi_occ)
154 
155  mat%num_rows = num_rows
156  mat%num_cols = num_cols
157  mat%num_nz = num_nz
158 
159  sll_allocate(mat%row_ptr(num_rows + 1), ierr)
160  sll_allocate(mat%col_ind(num_nz), ierr)
161  sll_allocate(mat%val(num_nz), ierr)
162 
163  print *, '#num_rows=', num_rows
164  print *, '#num_nz=', num_nz
165 
166  call sll_init_sparsematrix_mp( &
167  mat, &
168  num_rows, &
169  num_cols, &
170  num_patch, &
171  num_elements, &
172  local_to_global_col, &
173  num_local_dof_col, &
174  local_to_global_row, &
175  num_local_dof_row, &
176  lpi_columns, &
177  lpi_occ)
178 
179  mat%val(:) = 0.0_f64
180  sll_deallocate_array(lpi_columns, ierr)
181  sll_deallocate_array(lpi_occ, ierr)
182 
183  end subroutine sll_s_csr_matrix_mp_init
184 
185  integer function sll_count_non_zero_elts_mp( &
186  ai_nR, &
187  ai_nC, &
188  ai_npatch, &
189  api_nel, &
190  LM_Columns, &
191  nen_C, &
192  LM_Rows, &
193  nen_R, &
194  api_columns, &
195  api_occ)
196  ! _C FOR ROWS
197  ! _R FOR COLUMNS
198  implicit none
199  sll_int32 :: ai_nc
200  sll_int32 :: ai_nr
201  sll_int32 :: ai_npatch
202  sll_int32, dimension(:) :: api_nel
203  sll_int32, dimension(:, :, :), intent(in) :: lm_columns, lm_rows
204  sll_int32, dimension(:), intent(in) :: nen_c, nen_r
205  sll_int32, dimension(:, :), pointer :: api_columns
206  sll_int32, dimension(:), pointer :: api_occ
207  !local var
208  sll_int32 :: li_e
209  sll_int32 :: li_nel
210  sll_int32 :: li_id
211  sll_int32 :: li_nen_c
212  sll_int32 :: li_nen_r
213  sll_int32 :: li_b_c
214  sll_int32 :: li_a_c
215  sll_int32 :: li_b_r
216  sll_int32 :: li_a_r
217  sll_int32 :: li_i
218  sll_int32 :: li_result
219  sll_int32 :: ierr
220  logical :: ll_done
221  sll_int32, dimension(2) :: lpi_size
222  sll_int32, dimension(:, :), pointer :: lpi_columns
223 
224  do li_id = 1, ai_npatch
225 
226  li_nel = api_nel(li_id)
227 
228  ! WE FIRST COMPUTE, FOR EACH ROW, THE NUMBER OF COLUMNS THAT WILL BE USED
229  do li_e = 1, li_nel
230  ! print *,"li_id, li_e=",li_id, li_e
231  li_nen_c = nen_c(li_id) ! ao_con_C % opi_nen(li_id, li_e) ????
232  do li_b_c = 1, li_nen_c
233 
234  li_a_c = lm_columns(li_id, li_b_c, li_e)
235  if (li_a_c == 0) then
236  cycle
237  end if
238 
239  li_nen_r = nen_r(li_id)!ao_con_R % opi_nen(li_id, li_e)
240  do li_b_r = 1, li_nen_r
241 
242  li_a_r = lm_rows(li_id, li_b_r, li_e)
243  if (li_a_r == 0) then
244  cycle
245  end if
246 
247  ll_done = .false.
248  ! WE CHECK IF IT IS THE FIRST OCCURANCE OF THE COUPLE (li_A_C, li_A_R)
249  do li_i = 1, api_columns(li_a_c, 0)
250 
251  if (api_columns(li_a_c, li_i) /= li_a_r) then
252  cycle
253  end if
254 
255  ll_done = .true.
256  exit
257 
258  end do
259 
260  if (.not. ll_done) then
261 
262  api_occ(li_a_c) = api_occ(li_a_c) + 1
263 
264  ! li_A_C IS THE ROW NUM, li_A_R THE COLUMN NUM
265  ! INITIALIZATION OF THE SPARSE MATRIX
266  api_columns(li_a_c, 0) = api_columns(li_a_c, 0) + 1
267  api_columns(li_a_c, api_columns(li_a_c, 0)) = li_a_r
268 
269  ! resizing the array
270  lpi_size(1) = SIZE(api_columns, 1)
271  lpi_size(2) = SIZE(api_columns, 2)
272  if (lpi_size(2) < api_columns(li_a_c, 0)) then
273  sll_allocate(lpi_columns(lpi_size(1), lpi_size(2)), ierr)
274  lpi_columns = api_columns
275 
276  sll_deallocate(api_columns, ierr)
277 
278  sll_allocate(api_columns(lpi_size(1), 2*lpi_size(2)), ierr)
279  api_columns(1:lpi_size(1), 1:lpi_size(2)) = lpi_columns(1:lpi_size(1), 1:lpi_size(2))
280 
281  sll_deallocate(lpi_columns, ierr)
282  end if
283 
284  end if
285 
286  end do
287 
288  end do
289 
290  end do
291 
292  END DO
293 
294  ! COUNT NON ZERO ELEMENTS
295  li_result = sum(api_occ(1:ai_nr))
296 
297  sll_count_non_zero_elts_mp = li_result
298 
299  return
300  sll_assert(ai_nc > 0)
301  end function sll_count_non_zero_elts_mp
302 
303  subroutine sll_init_sparsematrix_mp( &
304  self, &
305  ai_nR, &
306  ai_nC, &
307  ai_npatch, &
308  api_nel, &
309  LM_Columns, &
310  nen_C, &
311  LM_Rows, &
312  nen_R, &
313  api_columns, &
314  api_occ)
315  ! _C FOR ROWS
316  ! _R FOR COLUMNS
317  implicit none
318  type(sll_t_csr_matrix) :: self
319  sll_int32 :: ai_nc
320  sll_int32 :: ai_nr
321  sll_int32 :: ai_npatch
322  sll_int32, dimension(:) :: api_nel
323  sll_int32, dimension(:), intent(in) :: nen_c, nen_r
324  sll_int32, dimension(:, :, :), intent(in) :: lm_columns, lm_rows
325  sll_int32, dimension(:, :), pointer :: api_columns
326  sll_int32, dimension(:), pointer :: api_occ
327  !local var
328  sll_int32 :: li_nel
329  sll_int32 :: li_id
330  sll_int32 :: li_e
331  sll_int32 :: li_b_c
332  sll_int32 :: li_a_c
333  !sll_int32 :: li_b_R
334  !sll_int32 :: li_A_R
335  !sll_int32 :: li_index
336  sll_int32 :: li_i
337  sll_int32 :: li_size
338  sll_int32 :: li_nen_c
339  sll_int32 :: li_err
340  sll_int32 :: li_flag
341  sll_int32, dimension(:), pointer :: lpr_tmp
342 
343  ! INITIALIZING ia
344  self%row_ptr(1) = 1
345 
346  do li_i = 1, self%num_rows !self % oi_nR
347 
348  self%row_ptr(li_i + 1) = self%row_ptr(1) + sum(api_occ(1:li_i))
349 
350  end do
351 
352  ! INITIALIZING ja
353  DO li_id = 1, ai_npatch
354 
355  li_nel = api_nel(li_id)
356 
357  ! WE FIRST COMPUTE, FOR EACH ROW, THE NUMBER OF COLUMNS THAT WILL BE USED
358  do li_e = 1, li_nel
359 
360  li_nen_c = nen_c(li_id)!ao_con_C % opi_nen(li_id, li_e)
361  do li_b_c = 1, li_nen_c
362 
363  li_a_c = lm_columns(li_id, li_b_c, li_e)
364  !ao_con_C % opi_LM(li_id, li_b_C, li_e)
365 
366  if (li_a_c == 0) then
367  cycle
368  end if
369 
370  if (api_columns(li_a_c, 0) == 0) then
371  cycle
372  end if
373 
374  li_size = api_columns(li_a_c, 0)
375 
376  allocate (lpr_tmp(li_size), stat=li_err)
377  if (li_err .ne. 0) li_flag = 10
378 
379  lpr_tmp(1:li_size) = api_columns(li_a_c, 1:li_size)
380 
381  call sll_s_qsortc(lpr_tmp)
382 
383  do li_i = 1, li_size
384 
385  self%col_ind(self%row_ptr(li_a_c) + li_i - 1) = lpr_tmp(li_i)
386 
387  end do
388 
389  api_columns(li_a_c, 0) = 0
390  deallocate (lpr_tmp)
391 
392  end do
393 
394  end do
395 
396  end do
397 
398  return
399  sll_assert(ai_nc > 0)
400  sll_assert(ai_nr > 0)
401  print *, size(lm_rows, 1)
402  print *, nen_r
403 
404  end subroutine sll_init_sparsematrix_mp
405 
406 end module sll_m_sparse_matrix_mp
sll_m_qsort_partition
Definition: sll_m_qsort_partition.F90:1
sll_m_qsort_partition::sll_s_qsortc
recursive subroutine, public sll_s_qsortc(A)
Definition: sll_m_qsort_partition.F90:14
sll_m_sparse_matrix_mp
Definition: sll_m_sparse_matrix_mp.F90:18
sll_m_sparse_matrix_mp::sll_count_non_zero_elts_mp
integer function sll_count_non_zero_elts_mp(ai_nR, ai_nC, ai_npatch, api_nel, LM_Columns, nen_C, LM_Rows, nen_R, api_columns, api_occ)
Definition: sll_m_sparse_matrix_mp.F90:196
sll_m_sparse_matrix_mp::sll_f_new_csr_matrix_mp
type(sll_t_csr_matrix) function, pointer, public sll_f_new_csr_matrix_mp(num_rows, num_cols, num_patch, num_elements, local_to_global_row, num_local_dof_row, local_to_global_col, num_local_dof_col)
allocates the memory space for a new CSR type on the heap,
Definition: sll_m_sparse_matrix_mp.F90:66
sll_m_sparse_matrix_mp::sll_init_sparsematrix_mp
subroutine sll_init_sparsematrix_mp(self, ai_nR, ai_nC, ai_npatch, api_nel, LM_Columns, nen_C, LM_Rows, nen_R, api_columns, api_occ)
Definition: sll_m_sparse_matrix_mp.F90:315
sll_m_sparse_matrix_mp::sll_s_csr_matrix_mp_init
subroutine, public sll_s_csr_matrix_mp_init(mat, num_rows, num_cols, num_patch, num_elements, local_to_global_row, num_local_dof_row, local_to_global_col, num_local_dof_col)
initialization of CSR matrix type
Definition: sll_m_sparse_matrix_mp.F90:114
sll_m_sparse_matrix
Sparse matrix linear solver utilities.
Definition: sll_m_sparse_matrix.F90:25
sll_m_sparse_matrix::sll_t_csr_matrix
type for CSR format
Definition: sll_m_sparse_matrix.F90:68
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:40 for Semi-Lagrangian Library by   doxygen 1.9.1