SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
st_csr.f90
Go to the documentation of this file.
1 !
2 !Authors: Jean-Luc Guermond, Lugi Quartapelle, Copyright 1994
3 !
4 MODULE st_matrix
5 
6  PUBLIC :: st_aij_csr_glob_block, extract, create_my_ghost, st_aij_csr, tri_jlg, st_aij_csr_loc_block, &
7  st_csr, st_csr_block, st_csr_bloc
8  PRIVATE
9 #include "petsc/finclude/petsc.h"
10 CONTAINS
11 
12  !=========================================================================
13 
14  SUBROUTINE create_my_ghost(mesh,LA,ifrom)
15  USE def_type_mesh
16  IMPLICIT NONE
17  TYPE(mesh_type) :: mesh
18  type(petsc_csr_la) :: la
19  INTEGER, DIMENSION(:), POINTER :: ifrom
20  INTEGER :: kmax, nifrom, start, fin, k
21  kmax = SIZE(la%loc_to_glob,1)
22  nifrom = mesh%np-mesh%dom_np
23  ALLOCATE(ifrom(kmax*nifrom))
24  IF (nifrom/=0) THEN
25  DO k = 1, kmax
26  start = (k-1)*nifrom+1
27  fin = start + nifrom - 1
28  ifrom(start:fin)= la%loc_to_glob(k,mesh%dom_np+1:)-1
29  END DO
30  END IF
31  END SUBROUTINE create_my_ghost
32 
33  SUBROUTINE extract(xghost,ks,ke,LA,phi)
34  USE def_type_mesh
35  IMPLICIT NONE
36  REAL(KIND=8), DIMENSION(:), INTENT(OUT):: phi
37  type(petsc_csr_la) :: la
38  INTEGER :: ks, ke
39  INTEGER :: k, start, fin, nbghost, s, f
40  !#include "petsc/finclude/petsc.h"
41 #include "petsc/finclude/petscvec.h90"
42  vec :: xghost
43  petscerrorcode :: ierr
44  petscscalar, POINTER :: x_loc(:)
45  CALL vecgetarrayf90(xghost, x_loc, ierr)
46  DO k = ks, ke
47  start = sum(la%dom_np(1:k-1)) + 1
48  fin = start + la%dom_np(k) - 1
49  s = sum(la%np(ks:k-1)) + 1
50  f = s + la%dom_np(k) - 1
51  phi(s:f) = x_loc(start:fin)
52  nbghost = la%np(k) - la%dom_np(k)
53  start = sum(la%dom_np) + sum(la%np(1:k-1)-la%dom_np(1:k-1)) + 1
54  fin = start + nbghost - 1
55  s = sum(la%np(ks:k-1))+la%dom_np(k)+1
56  f = s + nbghost - 1
57  phi(s:f) = x_loc(start:fin)
58  END DO
59  CALL vecrestorearrayf90(xghost, x_loc, ierr)
60  END SUBROUTINE extract
61 
62  SUBROUTINE block_index(communicator,kmax,mesh,loc_to_glob_LA)
63  USE def_type_mesh
64  USE my_util
65  IMPLICIT NONE
66  TYPE(mesh_type) :: mesh
67  INTEGER, INTENT(IN) :: kmax
68  INTEGER, DIMENSION(:,:), POINTER :: loc_to_glob_la
69  INTEGER, DIMENSION(:), POINTER :: dom_np, disp
70  INTEGER :: code, nb_procs, rank
71  INTEGER :: i, p, n, k, i_loc, proc, iglob
72  !#include "petsc/finclude/petsc.h"
73  mpi_comm :: communicator
74 
75  CALL mpi_comm_size(communicator,nb_procs,code)
76  CALL mpi_comm_rank(communicator,rank,code)
77  ALLOCATE(dom_np(nb_procs) ,disp(nb_procs+1))
78  CALL mpi_allgather(mesh%dom_np, 1, mpi_integer, dom_np, 1, &
79  mpi_integer, communicator ,code)
80  disp(1)=1
81  DO n = 1, nb_procs
82  disp(n+1) = disp(n) + dom_np(n)
83  END DO
84  IF (ASSOCIATED(mesh%disp)) THEN
85  nullify(mesh%disp)
86  END IF
87  IF (ASSOCIATED(mesh%domnp)) THEN
88  nullify(mesh%domnp)
89  END IF
90  ALLOCATE(mesh%disp(nb_procs+1))
91  ALLOCATE(mesh%domnp(nb_procs))
92  mesh%disp = disp
93  mesh%domnp=dom_np
94 
95  ALLOCATE(loc_to_glob_la(kmax,mesh%np))
96  proc = rank + 1
97 
98  DO i = 1, mesh%dom_np
99  DO k = 1, kmax
100  loc_to_glob_la(k,i) = kmax*(disp(proc)-1)+(k-1)*dom_np(proc)+i
101  END DO
102  END DO
103 
104 !!$!TEST
105 !!$ DO i = 1, mesh%dom_np
106 !!$ iglob = mesh%loc_to_glob(i)
107 !!$ DO p = 2, nb_procs+1
108 !!$ IF (disp(p) > iglob) THEN
109 !!$ proc = p - 1
110 !!$ EXIT
111 !!$ END IF
112 !!$ END DO
113 !!$ IF (rank+1/=proc) THEN
114 !!$ write(*,*) 'BUG2', rank+1, proc
115 !!$ STOP
116 !!$ END IF
117 !!$
118 !!$ DO k = 2, kmax
119 !!$ IF (loc_to_glob_LA(k,i) - dom_np(proc) /= loc_to_glob_LA(k-1,i)) THEN
120 !!$ write(*,*) ' BUG1 ', rank
121 !!$ stop
122 !!$ END IF
123 !!$ END DO
124 !!$ END DO
125 !!$!TEST
126 
127  DO i = mesh%dom_np + 1, mesh%np
128  iglob = mesh%loc_to_glob(i)
129  DO p = 2, nb_procs+1
130  IF (disp(p) > iglob) THEN
131  proc = p - 1
132  EXIT
133  END IF
134  END DO
135  i_loc = iglob - disp(proc) + 1
136  DO k = 1, kmax
137  loc_to_glob_la(k,i) = kmax*(disp(proc)-1)+(k-1)*dom_np(proc)+i_loc
138  END DO
139  END DO
140 
141 !!$!TEST
142 !!$ DO i = 1, mesh%np
143 !!$ iglob = mesh%loc_to_glob(i)
144 !!$ DO p = 2, nb_procs+1
145 !!$ IF (disp(p) > iglob) THEN
146 !!$ proc = p - 1
147 !!$ EXIT
148 !!$ END IF
149 !!$ END DO
150 !!$ DO k = 2, kmax
151 !!$ IF (loc_to_glob_LA(k,i) - dom_np(proc) /= loc_to_glob_LA(k-1,i)) THEN
152 !!$ write(*,*) ' BUG ', rank
153 !!$ stop
154 !!$ END IF
155 !!$ END DO
156 !!$ END DO
157 !!$!TEST
158 
159  DEALLOCATE(dom_np,disp)
160 
161  END SUBROUTINE block_index
162 
163  SUBROUTINE st_aij_csr_glob_block(communicator,kmax,mesh_glob,mesh,LA, opt_per)
164  ! input coefficient structure of the matrix and
165  ! perform the ordering of the unknowns
166  ! jj(nodes_per_element, number_of_elements)
167  ! ---> node number in the grid
168  USE def_type_mesh
169  USE my_util
170  IMPLICIT NONE
171  TYPE(mesh_type), INTENT(IN) :: mesh_glob,mesh
172  INTEGER, INTENT(IN) :: kmax
173  TYPE(periodic_type), OPTIONAL, INTENT(IN) :: opt_per
174  TYPE(petsc_csr_la), INTENT(OUT) :: la
175  INTEGER :: nparm=200
176  INTEGER :: me, nw, nmax, np, knp, ki, kj, k, njt, i1, i2
177  INTEGER :: m, ni, nj, i, j, n_a_d, iloc, jloc, jglob, nb_procs, p, proc=-1
178  INTEGER, DIMENSION(:,:), ALLOCATABLE :: ja_work
179  INTEGER, DIMENSION(:), ALLOCATABLE :: nja, a_d
180  INTEGER, DIMENSION(SIZE(mesh_glob%jj,1)) :: jj_loc
181  INTEGER, DIMENSION(:), ALLOCATABLE :: per_loc
182  LOGICAL :: out
183  !#include "petsc/finclude/petsc.h"
184  mpi_comm :: communicator
185 
186  CALL block_index(communicator,kmax,mesh,la%loc_to_glob)
187  nw = SIZE(mesh%jj, 1)
188  me = mesh%dom_me
189  np = mesh%dom_np
190  knp = kmax*np
191  nb_procs = SIZE(mesh%domnp)
192 
193  la%kmax = kmax
194  ALLOCATE(la%dom_np(kmax),la%np(kmax))
195  la%dom_np(:) = mesh%dom_np
196  la%np(:) = mesh%np
197 
198  IF (np==0) THEN
199  ALLOCATE(la%ia(0:0),la%ja(0))
200  la%ia(0) = 0
201  RETURN
202  END IF
203 
204  ALLOCATE (ja_work(knp,kmax*nparm), a_d(kmax*nparm), nja(knp))
205  ALLOCATE (per_loc(knp))
206  ja_work = 0
207  nja = 1
208  DO ki = 1, kmax
209  DO i = 1, np
210  ja_work((ki-1)*np+i,1) = la%loc_to_glob(ki,i)
211  END DO
212  END DO
213 
214  DO m = 1, mesh_glob%me
215  jj_loc = mesh_glob%jj(:,m)
216  IF (maxval(jj_loc)< mesh%loc_to_glob(1) .OR. minval(jj_loc)> mesh%loc_to_glob(1) + mesh%np -1) cycle
217  DO ni = 1, nw
218  iloc = jj_loc(ni)-mesh%loc_to_glob(1)+1
219  IF (iloc<1 .OR. iloc>np) cycle
220  DO ki = 1, kmax
221  i = iloc + (ki-1)*np
222  DO nj = 1, nw
223  jglob = jj_loc(nj)
224  IF (jglob< mesh%loc_to_glob(1) .OR. jglob> mesh%loc_to_glob(2)) THEN
225  DO p = 2, nb_procs+1
226  IF (mesh%disp(p) > jglob) THEN
227  proc = p - 1
228  EXIT
229  END IF
230  END DO
231  out = .true.
232  jloc = jglob - mesh%disp(proc) + 1
233  ELSE
234  out = .false.
235  jloc = jglob - mesh%loc_to_glob(1) + 1
236  END IF
237  DO kj = 1, kmax
238  IF (out) THEN
239  j = kmax*(mesh%disp(proc)-1)+(kj-1)*mesh%domnp(proc)+jloc
240  ELSE
241  j = la%loc_to_glob(kj,jloc)
242  END IF
243 
244  IF (minval(abs(ja_work(i,1:nja(i))-j)) /= 0) THEN
245  nja(i) = nja(i) + 1
246  ja_work(i,nja(i)) = j
247  END IF
248  END DO
249  END DO
250  END DO
251  END DO
252  END DO
253 
254  IF (present(opt_per)) THEN
255  DO k = 1, opt_per%n_bord
256  DO i = 1, SIZE(opt_per%list(k)%DIL)
257  per_loc=0
258  i1 = opt_per%list(k)%DIL(i)
259  i2 = opt_per%perlist(k)%DIL(i)
260  njt = nja(i1)+nja(i2)
261  IF (njt > kmax*nparm) THEN
262  CALL error_petsc('BUG in st_aij_glob_block, SIZE(ja) not large enough')
263  END IF
264  per_loc(1:nja(i1)) = ja_work(i1,1:nja(i1))
265  per_loc(nja(i1)+1:nja(i1)+nja(i2)) = ja_work(i2,1:nja(i2))
266  nja(i1) = njt
267  nja(i2) = njt
268  ja_work(i1,1:njt) = per_loc(1:njt)
269  ja_work(i2,1:njt) = per_loc(1:njt)
270  END DO
271  END DO
272  END IF
273 
274  IF (maxval(nja)>nparm) THEN
275  WRITE(*,*) 'ST_SPARSEKIT: dimension de ja doit etre >= ',nparm
276  stop
277  END IF
278 
279  nmax = 0
280  DO i = 1, knp
281  nmax = nmax + nja(i)
282  END DO
283  ALLOCATE(la%ia(0:knp),la%ja(0:nmax-1))
284  la%ia(0) = 0
285  DO i = 1, knp
286  CALL tri_jlg(ja_work(i,1:nja(i)), a_d, n_a_d)
287  IF (n_a_d /= nja(i)) THEN
288  WRITE(*,*) ' BUG : st_p1_CSR'
289  WRITE(*,*) 'n_a_d ', n_a_d, 'nja(i)', nja(i)
290  stop
291  END IF
292  la%ia(i) = la%ia(i-1) + nja(i)
293  la%ja(la%ia(i-1):la%ia(i)-1) = a_d(1:nja(i))-1
294  END DO
295  DEALLOCATE (ja_work, nja, a_d )
296  DEALLOCATE (per_loc)
297  END SUBROUTINE st_aij_csr_glob_block
298 
299 !!$ SUBROUTINE st_aij_csr_glob_block_without_per(communicator,kmax,mesh_glob,mesh,LA)
300 !!$ ! input coefficient structure of the matrix and
301 !!$ ! perform the ordering of the unknowns
302 !!$ ! jj(nodes_per_element, number_of_elements)
303 !!$ ! ---> node number in the grid
304 !!$ USE def_type_mesh
305 !!$ IMPLICIT NONE
306 !!$ TYPE(mesh_type), INTENT(IN) :: mesh_glob,mesh
307 !!$ INTEGER, INTENT(IN) :: kmax
308 !!$ TYPE(petsc_csr_LA), INTENT(OUT) :: LA
309 !!$ INTEGER :: nparm=200
310 !!$ INTEGER :: me, nw, nmax, np, knp, ki, kj
311 !!$ INTEGER :: m, ni, nj, i, j, n_a_d, iloc, jloc, jglob, nb_procs, p, proc=-1
312 !!$ INTEGER, DIMENSION(:,:), ALLOCATABLE :: ja_work
313 !!$ INTEGER, DIMENSION(:), ALLOCATABLE :: nja, a_d
314 !!$ INTEGER, DIMENSION(SIZE(mesh_glob%jj,1)) :: jj_loc
315 !!$ LOGICAL :: out
316 !!$ !#include "petsc/finclude/petsc.h"
317 !!$ MPI_Comm :: communicator
318 !!$
319 !!$ CALL block_index(communicator,kmax,mesh,LA%loc_to_glob)
320 !!$ nw = SIZE(mesh%jj, 1)
321 !!$ me = mesh%dom_me
322 !!$ np = mesh%dom_np
323 !!$ knp = kmax*np
324 !!$ nb_procs = SIZE(mesh%domnp)
325 !!$
326 !!$ LA%kmax = kmax
327 !!$ ALLOCATE(LA%dom_np(kmax),LA%np(kmax))
328 !!$ LA%dom_np(:) = mesh%dom_np
329 !!$ LA%np(:) = mesh%np
330 !!$
331 !!$ IF (np==0) THEN
332 !!$ ALLOCATE(LA%ia(0:0),LA%ja(0))
333 !!$ LA%ia(0) = 0
334 !!$ RETURN
335 !!$ END IF
336 !!$
337 !!$ ALLOCATE (ja_work(knp,kmax*nparm), a_d(kmax*nparm), nja(knp))
338 !!$ ja_work = 0
339 !!$ nja = 1
340 !!$ DO ki = 1, kmax
341 !!$ DO i = 1, np
342 !!$ ja_work((ki-1)*np+i,1) = LA%loc_to_glob(ki,i)
343 !!$ END DO
344 !!$ END DO
345 !!$
346 !!$ DO m = 1, mesh_glob%me
347 !!$ jj_loc = mesh_glob%jj(:,m)
348 !!$ IF (MAXVAL(jj_loc)< mesh%loc_to_glob(1) .OR. MINVAL(jj_loc)> mesh%loc_to_glob(1) + mesh%np -1) CYCLE
349 !!$ DO ni = 1, nw
350 !!$ iloc = jj_loc(ni)-mesh%loc_to_glob(1)+1
351 !!$ IF (iloc<1 .OR. iloc>np) CYCLE
352 !!$ DO ki = 1, kmax
353 !!$ i = iloc + (ki-1)*np
354 !!$ DO nj = 1, nw
355 !!$ jglob = jj_loc(nj)
356 !!$ IF (jglob< mesh%loc_to_glob(1) .OR. jglob> mesh%loc_to_glob(2)) THEN
357 !!$ DO p = 2, nb_procs+1
358 !!$ IF (mesh%disp(p) > jglob) THEN
359 !!$ proc = p - 1
360 !!$ EXIT
361 !!$ END IF
362 !!$ END DO
363 !!$ out = .TRUE.
364 !!$ jloc = jglob - mesh%disp(proc) + 1
365 !!$ ELSE
366 !!$ out = .FALSE.
367 !!$ jloc = jglob - mesh%loc_to_glob(1) + 1
368 !!$ END IF
369 !!$ DO kj = 1, kmax
370 !!$ IF (out) THEN
371 !!$ j = kmax*(mesh%disp(proc)-1)+(kj-1)*mesh%domnp(proc)+jloc
372 !!$ ELSE
373 !!$ j = LA%loc_to_glob(kj,jloc)
374 !!$ END IF
375 !!$
376 !!$ IF (MINVAL(ABS(ja_work(i,1:nja(i))-j)) /= 0) THEN
377 !!$ nja(i) = nja(i) + 1
378 !!$ ja_work(i,nja(i)) = j
379 !!$ END IF
380 !!$ END DO
381 !!$ END DO
382 !!$ END DO
383 !!$ END DO
384 !!$ END DO
385 !!$
386 !!$ IF (MAXVAL(nja)>nparm) THEN
387 !!$ WRITE(*,*) 'ST_SPARSEKIT: dimension de ja doit etre >= ',nparm
388 !!$ STOP
389 !!$ END IF
390 !!$
391 !!$ nmax = 0
392 !!$ DO i = 1, knp
393 !!$ nmax = nmax + nja(i)
394 !!$ END DO
395 !!$ ALLOCATE(LA%ia(0:knp),LA%ja(0:nmax-1))
396 !!$ LA%ia(0) = 0
397 !!$ DO i = 1, knp
398 !!$ CALL tri_jlg (ja_work(i,1:nja(i)), a_d, n_a_d)
399 !!$ IF (n_a_d /= nja(i)) THEN
400 !!$ WRITE(*,*) ' BUG : st_p1_CSR'
401 !!$ WRITE(*,*) 'n_a_d ', n_a_d, 'nja(i)', nja(i)
402 !!$ STOP
403 !!$ END IF
404 !!$ LA%ia(i) = LA%ia(i-1) + nja(i)
405 !!$ LA%ja(LA%ia(i-1):LA%ia(i)-1) = a_d(1:nja(i))-1
406 !!$ END DO
407 !!$ DEALLOCATE (ja_work, nja, a_d )
408 !!$ END SUBROUTINE st_aij_csr_glob_block_without_per
409 
410  SUBROUTINE st_aij_csr_loc_block(communicator,kmax,mesh,LA)
411  ! input coefficient structure of the matrix and`
412  ! perform the ordering of the unknowns
413  ! jj(nodes_per_element, number_of_elements)
414  ! ---> node number in the grid
415  USE def_type_mesh
416  IMPLICIT NONE
417  TYPE(mesh_type), INTENT(IN) :: mesh
418  INTEGER, INTENT(IN) :: kmax
419  TYPE(petsc_csr_la), INTENT(OUT) :: la
420  INTEGER :: nparm=90
421  INTEGER :: me, nw, nmax, np, knp, ki, kj
422  INTEGER :: m, ni, nj, i, j, n_a_d, iloc
423  INTEGER, DIMENSION(:,:), ALLOCATABLE :: ja_work
424  INTEGER, DIMENSION(:), ALLOCATABLE :: nja, a_d
425  INTEGER, DIMENSION(SIZE(mesh%jj,1)) :: j_loc
426  !#include "petsc/finclude/petsc.h"
427  mpi_comm :: communicator
428 
429  CALL block_index(communicator,kmax,mesh,la%loc_to_glob)
430  nw = SIZE(mesh%jj, 1)
431  me = mesh%dom_me
432  np = mesh%dom_np
433  knp = kmax*np
434  ALLOCATE (ja_work(knp,kmax*nparm), a_d(kmax*nparm), nja(knp))
435  ja_work = 0
436  nja = 1
437  DO ki = 1, kmax
438  DO i = 1, np
439  ja_work((ki-1)*np+i,1) = la%loc_to_glob(ki,i)
440  END DO
441  END DO
442 
443  DO m = 1, me
444  j_loc = mesh%jj(:,m)
445  DO ki = 1, kmax
446  DO ni = 1, nw
447  iloc = j_loc(ni)
448  IF (iloc>np) cycle
449  i = iloc + (ki-1)*np
450  DO kj = 1, kmax
451  DO nj = 1, nw
452  j = la%loc_to_glob(kj,j_loc(nj))
453  IF (minval(abs(ja_work(i,1:nja(i))-j)) /= 0) THEN
454  nja(i) = nja(i) + 1
455  ja_work(i,nja(i)) = j
456  END IF
457  END DO
458  END DO
459  END DO
460  END DO
461  END DO
462 
463  IF (maxval(nja)>nparm) THEN
464  WRITE(*,*) 'ST_SPARSEKIT: dimension de ja doit etre >= ',nparm
465  stop
466  END IF
467 
468  nmax = 0
469  DO i = 1, knp
470  nmax = nmax + nja(i)
471  END DO
472  ALLOCATE(la%ia(0:knp),la%ja(0:nmax-1))
473  la%ia(0) = 0
474  DO i = 1, knp
475  CALL tri_jlg(ja_work(i,1:nja(i)), a_d, n_a_d)
476  IF (n_a_d /= nja(i)) THEN
477  WRITE(*,*) ' BUG : st_p1_CSR'
478  WRITE(*,*) 'n_a_d ', n_a_d, 'nja(i)', nja(i)
479  stop
480  END IF
481  la%ia(i) = la%ia(i-1) + nja(i)
482  la%ja(la%ia(i-1):la%ia(i)-1) = a_d(1:nja(i))-1
483  END DO
484  DEALLOCATE (ja_work, nja, a_d )
485  END SUBROUTINE st_aij_csr_loc_block
486 
487  SUBROUTINE st_aij_csr(jj, aij)
488  ! input coefficient structure of the matrix and
489  ! perform the ordering of the unknowns
490  ! jj(nodes_per_element, number_of_elements)
491  ! ---> node number in the grid
492  USE def_type_mesh
493  IMPLICIT NONE
494  TYPE(aij_type), INTENT(OUT) :: aij
495  INTEGER, DIMENSION(:,:), INTENT(IN) :: jj
496  INTEGER :: nparm=90
497  INTEGER :: me, nw, nmax, np
498  INTEGER :: m, ni, nj, i, j, n_a_d
499  INTEGER, DIMENSION(:,:), ALLOCATABLE :: ja_work
500  INTEGER, DIMENSION(:), ALLOCATABLE :: nja, a_d
501 
502  nw = SIZE(jj, 1)
503  me = SIZE(jj, 2)
504  np = maxval(jj)
505  ALLOCATE (ja_work(np,nparm), aij%ia(np+1), a_d(nparm), nja(np))
506  ja_work = 0
507  nja = 1
508  DO i = 1, np
509  ja_work(i,1) = i
510  END DO
511 
512  DO m = 1, me
513  DO ni = 1, nw
514  i = jj(ni,m)
515  DO nj = 1, nw
516  j = jj(nj,m)
517  IF (minval(abs(ja_work(i,1:nja(i))-j)) /= 0) THEN
518  nja(i) = nja(i) + 1
519  ja_work(i,nja(i)) = j
520  END IF
521  END DO
522  END DO
523  END DO
524 
525  IF (maxval(nja)>nparm) THEN
526  WRITE(*,*) 'ST_SPARSEKIT: dimension de ja doit etre >= ',nparm
527  stop
528  END IF
529 
530  nmax = 0
531  DO i = 1, np
532  nmax = nmax + nja(i)
533  END DO
534  ALLOCATE(aij%ja(nmax))
535  aij%ia(1) = 1
536  DO i = 1, np
537  CALL tri_jlg(ja_work(i,1:nja(i)), a_d, n_a_d)
538  IF (n_a_d /= nja(i)) THEN
539  WRITE(*,*) ' BUG : st_p1_CSR'
540  WRITE(*,*) 'n_a_d ', n_a_d, 'nja(i)', nja(i)
541  stop
542  END IF
543  aij%ia(i+1) = aij%ia(i) + nja(i)
544  aij%ja(aij%ia(i):aij%ia(i+1)-1) = a_d(1:nja(i))
545  END DO
546  DEALLOCATE ( ja_work, nja, a_d )
547  END SUBROUTINE st_aij_csr
548 
549  SUBROUTINE st_csr_block(n_b,in,out)
550  !SUBROUTINE st_csr_block(in,vv_rt_loc_to_glob_LA,out)
551  !
552  ! Builds the CSR structure of a parallel block matrix
553  ! from its scalar counterpart.
554  !
555  USE def_type_mesh
556  IMPLICIT NONE
557  TYPE(aij_type), INTENT(IN) :: in
558  TYPE(aij_type), INTENT(OUT) :: out
559  INTEGER, INTENT(IN) :: n_b ! Number of blocs
560  INTEGER :: ki, kj, i, ib, jb, p, pb, bloc_size, np, nnz
561 
562  np = SIZE(in%ia) - 1
563  nnz = in%ia(np+1) - in%ia(1)
564  ALLOCATE (out%ia(n_b*np+1), out%ja(n_b*n_b*nnz))
565 
566  bloc_size = SIZE(in%ia) - 1
567 
568  out%ia(1) = in%ia(1)
569 
570  DO ki = 1, n_b
571  DO i = 2, bloc_size + 1
572  ib = i + (ki-1)*bloc_size
573  out%ia(ib) = out%ia(ib-1) + n_b*(in%ia(i) - in%ia(i-1))
574  END DO
575  END DO
576 
577  DO ki = 1, n_b
578  DO i = 1, bloc_size
579  ib = i + (ki-1)*bloc_size
580 
581  DO kj = 1, n_b
582  DO p = in%ia(i), in%ia(i+1) - 1
583  jb = in%ja(p) + (kj-1)*bloc_size
584 
585  pb = out%ia(ib) + p - in%ia(i) + (kj-1)*(in%ia(i+1)-in%ia(i))
586  out%ja(pb) = jb
587 
588  END DO
589  END DO
590 
591  END DO
592  END DO
593  END SUBROUTINE st_csr_block
594 
595  SUBROUTINE st_csr_bloc(ia,ja,ia_b,ja_b,n_b)
596  !
597  ! Builds the CSR structure of a 3D vector matrix
598  ! from its scalar counterpart.
599  !
600  IMPLICIT NONE
601  INTEGER, DIMENSION(:), INTENT(IN) :: ia, ja
602  INTEGER, DIMENSION(:), POINTER :: ia_b, ja_b
603  INTEGER, INTENT(IN) :: n_b ! Number of blocs
604  INTEGER :: ki, kj, i, ib, jb, p, pb, bloc_size, np, nnz
605 
606  np = SIZE(ia) - 1
607  nnz = ia(np+1) - ia(1)
608  ALLOCATE (ia_b(n_b*np+1), ja_b(n_b*n_b*nnz))
609 
610  bloc_size = SIZE(ia) - 1
611 
612  ia_b(1) = ia(1)
613 
614  DO ki = 1, n_b
615  DO i = 2, bloc_size + 1
616  ib = i + (ki-1)*bloc_size
617  ia_b(ib) = ia_b(ib-1) + n_b*(ia(i) - ia(i-1))
618  END DO
619  END DO
620 
621  DO ki = 1, n_b
622  DO i = 1, bloc_size
623  ib = i + (ki-1)*bloc_size
624 
625  DO kj = 1, n_b
626  DO p = ia(i), ia(i+1) - 1
627  jb = ja(p) + (kj-1)*bloc_size
628 
629  pb = ia_b(ib) + p - ia(i) + (kj-1)*(ia(i+1)-ia(i))
630  ja_b(pb) = jb
631 
632  END DO
633  END DO
634  END DO
635  END DO
636 
637  END SUBROUTINE st_csr_bloc
638 
639  SUBROUTINE st_csr(jj, ia, ja)
640  ! input coefficient structure of the matrix and
641  ! perform the ordering of the unknowns
642  ! jj(nodes_per_element, number_of_elements)
643  ! ---> node number in the grid
644  IMPLICIT NONE
645  INTEGER, DIMENSION(:,:), INTENT(IN) :: jj
646  INTEGER, DIMENSION(:), POINTER :: ia, ja
647  INTEGER :: nparm=90
648  INTEGER :: me, nw, nmax, np
649  INTEGER :: m, ni, nj, i, j, n_a_d
650  INTEGER, DIMENSION(:,:), ALLOCATABLE :: ja_work
651  INTEGER, DIMENSION(:), ALLOCATABLE :: nja, a_d
652 
653  nw = SIZE(jj, 1)
654  me = SIZE(jj, 2)
655  np = maxval(jj)
656 
657  ALLOCATE (ja_work(np,nparm), ia(np+1), a_d(nparm), nja(np))
658 
659  ja_work = 0
660  nja = 1
661  DO i = 1, np
662  ja_work(i,1) = i
663  END DO
664 
665  DO m = 1, me
666  DO ni = 1, nw
667  i = jj(ni,m)
668 
669  DO nj = 1, nw
670  j = jj(nj,m)
671  IF (minval(abs(ja_work(i,1:nja(i))-j)) /= 0) THEN
672  nja(i) = nja(i) + 1
673  ja_work(i,nja(i)) = j
674  END IF
675  END DO
676 
677 
678  END DO
679  END DO
680 
681  IF (maxval(nja)>nparm) THEN
682  WRITE(*,*) 'ST_SPARSEKIT: dimension de ja doit etre >= ',nparm
683  stop
684  END IF
685 
686  nmax = 0
687  DO i = 1, np
688  nmax = nmax + nja(i)
689  END DO
690  ALLOCATE(ja(nmax))
691 
692  ia(1) = 1
693  DO i = 1, np
694  CALL tri_jlg(ja_work(i,1:nja(i)), a_d, n_a_d)
695  IF (n_a_d /= nja(i)) THEN
696  WRITE(*,*) ' BUG : st_p1_CSR'
697  WRITE(*,*) 'n_a_d ', n_a_d, 'nja(i)', nja(i)
698  stop
699  END IF
700  ia(i+1) = ia(i) + nja(i)
701  ja(ia(i):ia(i+1)-1) = a_d(1:nja(i))
702  END DO
703 
704 
705  DEALLOCATE ( ja_work, nja, a_d )
706 
707  END SUBROUTINE st_csr
708 
709 !!$ SUBROUTINE st_csr_edge_stab(jji, ia, ja)
710 !!$ ! input coefficient structure of the matrix and
711 !!$ ! perform the ordering of the unknowns
712 !!$ ! jj(nodes_per_element, number_of_elements)
713 !!$ ! ---> node number in the grid
714 !!$ IMPLICIT NONE
715 !!$ INTEGER, DIMENSION(:,:,:), INTENT(IN) :: jji
716 !!$ INTEGER, DIMENSION(:), POINTER :: ia, ja
717 !!$ INTEGER :: nparm=180
718 !!$ INTEGER :: mi, nw, nmax, np
719 !!$ INTEGER :: ms, ni, nj, i, j, n_a_d, cotei, cotej
720 !!$ INTEGER, DIMENSION(:,:), ALLOCATABLE :: ja_work
721 !!$ INTEGER, DIMENSION(:), ALLOCATABLE :: nja, a_d
722 !!$
723 !!$ nw = SIZE(jji, 1)
724 !!$ mi = SIZE(jji, 3)
725 !!$ np = MAXVAL(jji)
726 !!$
727 !!$ ALLOCATE (ja_work(np,nparm), ia(np+1), a_d(nparm), nja(np))
728 !!$
729 !!$ ja_work = 0
730 !!$ nja = 1
731 !!$ DO i = 1, np
732 !!$ ja_work(i,1) = i
733 !!$ END DO
734 !!$
735 !!$ DO ms = 1, mi
736 !!$ DO cotei = 1, 2
737 !!$ DO ni = 1, nw
738 !!$ i = jji(ni,cotei,ms)
739 !!$ DO cotej = 1, 2
740 !!$ DO nj = 1, nw
741 !!$ j = jji(nj,cotej,ms)
742 !!$ IF (MINVAL(ABS(ja_work(i,1:nja(i))-j)) /= 0) THEN
743 !!$ nja(i) = nja(i) + 1
744 !!$ ja_work(i,nja(i)) = j
745 !!$ END IF
746 !!$ END DO
747 !!$ END DO
748 !!$ END DO
749 !!$ END DO
750 !!$ END DO
751 !!$
752 !!$ IF (MAXVAL(nja)>nparm) THEN
753 !!$ WRITE(*,*) 'ST_SPARSEKIT: dimension de ja doit etre >= ',nparm
754 !!$ STOP
755 !!$ END IF
756 !!$
757 !!$ nmax = 0
758 !!$ DO i = 1, np
759 !!$ nmax = nmax + nja(i)
760 !!$ END DO
761 !!$ ALLOCATE(ja(nmax))
762 !!$
763 !!$ ia(1) = 1
764 !!$ DO i = 1, np
765 !!$ CALL tri_jlg (ja_work(i,1:nja(i)), a_d, n_a_d)
766 !!$ IF (n_a_d /= nja(i)) THEN
767 !!$ WRITE(*,*) ' BUG : st_p1_CSR'
768 !!$ WRITE(*,*) 'n_a_d ', n_a_d, 'nja(i)', nja(i)
769 !!$ STOP
770 !!$ END IF
771 !!$ ia(i+1) = ia(i) + nja(i)
772 !!$ ja(ia(i):ia(i+1)-1) = a_d(1:nja(i))
773 !!$ END DO
774 !!$
775 !!$ DEALLOCATE ( ja_work, nja, a_d )
776 !!$
777 !!$ END SUBROUTINE st_csr_edge_stab
778 
779  SUBROUTINE tri_jlg (a, a_d, n_a_d)
780  ! sort in ascending order of the integer array a and generation
781  ! of the integer array a_d whose first n_a_d leading entries
782  ! contain different values in ascending order, while all the
783  ! remaining entries are set to zero
784  ! sorting by Shell's method.
785  IMPLICIT NONE
786  INTEGER, DIMENSION(:), INTENT(INOUT) :: a
787  INTEGER, DIMENSION(:), INTENT(OUT) :: a_d
788  INTEGER, INTENT(OUT) :: n_a_d
789  INTEGER :: n, na, inc, i, j, k, ia
790 
791  na = SIZE(a)
792 
793  ! sort phase
794 
795  IF (na == 0) THEN
796  n_a_d = 0
797  RETURN
798  ENDIF
799 
800  inc = 1
801  DO WHILE (inc <= na)
802  inc = inc * 3
803  inc = inc + 1
804  ENDDO
805 
806  DO WHILE (inc > 1)
807  inc = inc/3
808  DO i = inc + 1, na
809  ia = a(i)
810  j = i
811  DO WHILE (a(j-inc) > ia)
812  a(j) = a(j-inc)
813  j = j - inc
814  IF (j <= inc) EXIT
815  ENDDO
816  a(j) = ia
817  ENDDO
818  ENDDO
819 
820  ! compression phase
821 
822  n = 1
823  a_d(n) = a(1)
824  DO k = 2, na
825  IF (a(k) > a(k-1)) THEN
826  n = n + 1
827  a_d(n) = a(k)
828  !TEST JUIN 13 2008
829  ELSE
830  WRITE(*,*) 'We have a problem in the compression phase of tri_jlg', a(k), a(k-1)
831  !TEST JUIN 13 2008
832  ENDIF
833  ENDDO
834 
835  n_a_d = n
836 
837  a_d(n_a_d + 1 : na) = 0
838 
839  END SUBROUTINE tri_jlg
840 
841 END MODULE st_matrix
842 
843 
844 
845 
846 
847 
848 
849 
850 
851 
st_matrix::st_csr_bloc
subroutine, public st_csr_bloc(ia, ja, ia_b, ja_b, n_b)
Definition: st_csr.f90:595
def_type_mesh::petsc_csr_la
Definition: def_type_mesh.f90:18
st_matrix::create_my_ghost
subroutine, public create_my_ghost(mesh, LA, ifrom)
Definition: st_csr.f90:14
st_matrix::st_aij_csr_loc_block
subroutine, public st_aij_csr_loc_block(communicator, kmax, mesh, LA)
Definition: st_csr.f90:410
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:73
st_matrix
Definition: st_csr.f90:4
def_type_mesh::aij_type
Definition: def_type_mesh.f90:8
st_matrix::block_index
subroutine block_index(communicator, kmax, mesh, loc_to_glob_LA)
Definition: st_csr.f90:62
f
section doc_intro_frame_work_num_app Numerical approximation subsection doc_intro_fram_work_num_app_Fourier_FEM Fourier Finite element representation The SFEMaNS code uses a hybrid Fourier Finite element formulation The Fourier decomposition allows to approximate the problem’s solutions for each Fourier mode modulo nonlinear terms that are made explicit The variables are then approximated on a meridian section of the domain with a finite element method The numerical approximation of a function f $f f is written in the following generic f with f(r,\theta, z)\f $the cylindrical coordinates
st_matrix::st_aij_csr
subroutine, public st_aij_csr(jj, aij)
Definition: st_csr.f90:487
def_type_mesh
Definition: def_type_mesh.f90:4
my_util
Definition: my_util.f90:1
st_matrix::tri_jlg
subroutine, public tri_jlg(a, a_d, n_a_d)
Definition: st_csr.f90:779
def_type_mesh::periodic_type
Definition: def_type_mesh.f90:54
st_matrix::st_aij_csr_glob_block
subroutine, public st_aij_csr_glob_block(communicator, kmax, mesh_glob, mesh, LA, opt_per)
Definition: st_csr.f90:163
st_matrix::st_csr_block
subroutine, public st_csr_block(n_b, in, out)
Definition: st_csr.f90:549
my_util::error_petsc
subroutine error_petsc(string)
Definition: my_util.f90:15
st_matrix::st_csr
subroutine, public st_csr(jj, ia, ja)
Definition: st_csr.f90:639
st_matrix::extract
subroutine, public extract(xghost, ks, ke, LA, phi)
Definition: st_csr.f90:33