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