Choral
csr_mod.F90
Go to the documentation of this file.
1 
11 
12 
13 module csr_mod
14 
15  use real_type
16  use algebra_set, only: shellsort_dec
17  use basic_tools
18  !$ use OMP_LIB
19 
20  implicit none
21  private
22 
23 
24  ! %----------------------------------------%
25  ! | |
26  ! | PUBLIC DATA |
27  ! | |
28  ! %----------------------------------------%
29  !! TYPE
30  public :: csr
31 
32  !! GENERIC
33  public :: clear
34  public :: write
35  public :: print
36  public :: valid !! TESTED
37  public :: transpose !! TESTED
38  public :: matvecprod !! TESTED
39  public :: setrow !! TESTED
40  public :: getrow !! TESTED
41  public :: addentry !! TESTED
42  public :: equal
43  public :: copy
44  public :: add
45  public :: scale
46  public :: reorder
47  public :: sortrows
48  public :: invllt
49 
50  !! NON GENERIC
51  public :: addentry_2 !! TESTED
52  public :: getdiag
53  public :: add_submatrix
54 
55  ! %----------------------------------------%
56  ! | |
57  ! | DERIVED TYPE |
58  ! | |
59  ! %----------------------------------------%
60 
62  type csr
63  integer :: nl = 0
64  integer :: nnz = 0
65  logical :: sorted = .false.
66  integer , dimension(:), allocatable :: row
67  integer , dimension(:), allocatable :: col
68  integer :: width = 0
69  real(RP), dimension(:), allocatable :: a
70  integer , dimension(:), allocatable :: diag
71 
72  contains
73 
75  final :: csr_clear
76 
77  end type csr
78 
79 
80  ! %----------------------------------------%
81  ! | |
82  ! | GENERIc SUBROUTINES |
83  ! | |
84  ! %----------------------------------------%
86  interface clear
87  module procedure csr_clear
88  end interface clear
89 
90  !! constructor
91  interface csr
92  module procedure csr_create
93  end interface csr
94 
96  interface print
97  module procedure csr_print
98  end interface print
99 
101  interface valid
102  module procedure csr_valid !! TESTED
103  end interface valid
104 
105  interface getrow
106  module procedure csr_getrow !! TESTED
107  end interface getrow
108 
109  interface setrow
110  module procedure csr_setrow !! TESTED
111  end interface setrow
112 
113  interface sortrows
114  module procedure csr_sortrows
115  end interface sortrows
116 
117  interface addentry
118  module procedure csr_addentry !! TESTED
119  end interface addentry
120 
122  interface equal
123  module procedure csr_equal
124  end interface equal
125 
127  interface write
128  module procedure csr_write
129  end interface write
130 
132  interface matvecprod
133  module procedure csr_matvecprod !! TESTED
134  end interface matvecprod
135 
138  interface add
139  module procedure csr_add_1, csr_add_2
140  end interface add
141 
145  interface add_submatrix
146  module procedure add_submatrix_1, add_submatrix_2
147  end interface add_submatrix
148 
150  interface scale
151  module procedure csr_scale
152  end interface scale
153 
156  interface reorder
157  module procedure csr_reorder
158  end interface reorder
159 
161  interface copy
162  module procedure csr_copy
163  end interface copy
164 
166  interface transpose
167  module procedure csr_transpose !! TESTED
168  end interface transpose
169 
172  interface invllt
173  module procedure csr_invllt
174  end interface invllt
175 
176 contains
177 
179  subroutine csr_clear(m)
180  type(csr), intent(inout) :: m
181 
182  call freemem(m%row)
183  call freemem(m%col)
184  call freemem(m%a )
185  call freemem(m%diag)
186  m%nl = 0
187  m%nnz = 0
188  m%sorted = .false.
189  m%width = 0
190 
191  end subroutine csr_clear
192 
203  function csr_create(nnz) result(m)
204  type(csr) :: m
205  integer, dimension(:), intent(in) :: nnz
206 
207  integer :: nz, ii, jj, wdt, length
208 
209  call clear(m)
210  m%nl = size(nnz,1)
211  nz = sum(nnz)
212  m%nnz = nz
213  wdt = 0
214  call allocmem(m%row,m%nl+1)
215  call allocmem(m%col, nz )
216  call allocmem(m%a , nz )
217 
218  m%row = 0
219  m%col = 0
220  m%a = 0._rp
221 
222  ! Construction of m%row
223  jj = 1
224  do ii = 1, m%nl
225  m%row(ii) = jj
226  length = nnz(ii)
227  if (length>wdt) wdt = length
228  jj = jj + length
229  end do
230  m%row(m%nl+1) = jj
231 
232  m%width = wdt
233  if (.NOT.valid(m)) call quit( &
234  & 'csr_mod: csr_create')
235 
236  end function csr_create
237 
238 
241  function csr_valid(m) result(b)
242  type(csr), intent(in) :: m
243  logical :: b
244 
245  integer :: ii, wdt, length
246 
247  b = .false.
248  if (.NOT. (allocated(m%row)) ) return
249  if (.NOT. (allocated(m%row)) ) return
250  if (.NOT. (allocated(m%a)) ) return
251 
252  b = ( m%nl>0 ) .AND. ( m%nnz>0 )
253  b = b .AND. ( size(m%row,1) == m%nl+1 )
254  b = b .AND. ( size(m%col,1) == m%nnz )
255  b = b .AND. ( size(m%a ,1) == m%nnz )
256  b = b .AND. ( m%row(m%nl+1) == size(m%col)+1 )
257 
258  if (allocated(m%diag)) then
259  b = b.AND.(size(m%diag,1) == m%nl)
260  b = b.AND.( maxval(m%diag) <= m%nnz )
261  end if
262 
263  wdt = 0
264  do ii = 1, m%nl
265  length = m%row(ii+1) - m%row(ii)
266  if (length>wdt) wdt = length
267  end do
268  b = ( wdt == m%width )
269 
270  end function csr_valid
271 
272 
273 
285  subroutine csr_getrow(sz, v, p, n, g, ii)
286  real(RP) , dimension(n), intent(out):: v
287  integer , dimension(n), intent(out):: p
288  integer , intent(out):: sz
289  type(csr) , intent(in) :: g
290  integer , intent(in) :: ii, n
291 
292  integer :: j1, j2
293 
294 #if DBG>0
295  if ( .NOT. valid(g) ) call quit( &
296  & "csr_mod: csr_getRow_1: CSR matrix 'g' not valid")
297 
298  if ( (11<=0) .OR. (ii>g%nl) ) call quit( &
299  & "csr_mod: csr_getRow_1: wrong line index 'ii'")
300 #endif
301 
302  j1 = g%row(ii)
303  j2 = g%row(ii+1)
304  sz = j2 - j1
305 
306 #if DBG>0
307  if (sz>n ) call quit( &
308  & "csr_mod: csr_getRow_1: unsufficient size&
309  & for the array 'p'")
310 #endif
311 
312  p(1:sz) = g%col(j1:j2-1)
313  v(1:sz) = g%a(j1:j2-1)
314 
315  end subroutine csr_getrow
316 
318  subroutine csr_setrow(m, val, col, ii)
319  type(csr) , intent(inout):: m
320  real(RP), dimension(:), intent(in) :: val
321  integer , dimension(:), intent(in) :: col
322  integer , intent(in) :: ii
323 
324  integer :: j1, j2
325 
326  j1 = m%row(ii)
327  j2 = m%row(ii+1)
328  m%col(j1:j2-1) = col
329  m%a( j1:j2-1) = val
330 
331  end subroutine csr_setrow
332 
333 
337  function csr_equal(g1, g2) result(b)
338  logical :: b
339  type(csr), intent(in) :: g1, g2
340 
341  integer , dimension(g1%width) :: p1, p2, q1, q2
342  real(RP), dimension(g1%width) :: a1, a2
343  integer :: ii, jj, sz, wdt
344 
345  b = valid(g1) .AND. valid(g2)
346  b = b .AND. ( g1%nl == g2%nl ) .AND. ( g1%nnz == g2%nnz )
347  if(.NOT.b) return
348 
349  b = b .AND. ( g1%width == g2%width )
350  if(.NOT.b) return
351 
352  b = b .AND. (all( g1%row == g2%row ) )
353  if(.NOT.b) return
354 
355  wdt = g1%width
356  do ii=1, g1%nl
357  call getrow(sz, a1, p1, wdt, g1, ii)
358  call getrow(sz, a2, p2, wdt, g2, ii)
359 
360  if (sz==0) cycle
361 
362  call shellsort_dec(p1(1:sz), q1(1:sz), sz)
363  call shellsort_dec(p2(1:sz), q2(1:sz), sz)
364 
365  do jj=1, sz
366  b = b .AND. ( p1(q1(jj)) == p2(q2(jj)) )
367  b = b .AND. equal(a1(q1(jj)), a2(q2(jj)))
368  end do
369 
370  if (.NOT. b) return
371 
372  end do
373 
374  end function csr_equal
375 
376 
384  subroutine csr_addentry(m, rw, cl, val)
385  type(csr), intent(inout) :: m
386  integer , intent(in) :: rw, cl
387  real(RP), intent(in) :: val
388 
389  integer :: ii
390 
391  do ii = m%row(rw),m%row(rw+1)-1
392 
393  if (m%col(ii)==0) then
394  ! the entry m(rw,cl) does not exist yet
395  m%col(ii)= cl
396  m%a(ii) = val
397  goto 10
398  else if (m%col(ii)==cl) then
399  ! the entry m(rw,cl) already exists
400  m%a(ii) = m%a(ii)+val
401  goto 10
402  end if
403  end do
404  ! the line rw is already filled
405  ! write(*,*) m%col(m%row(rw):m%row(rw+1)-1)
406  ! write(*,*) cl
407  call quit("csr_mod: csr_addEntry: 1")
408 
409 10 continue
410  end subroutine csr_addentry
411 
412 
417  subroutine addentry_2(m, rw, cl, val)
418  type(csr), intent(inout) :: m
419  integer , intent(in) :: rw, cl
420  real(RP), intent(in) :: val
421 
422  integer :: ii
423 
424  ii = m%row(rw)
425 
426  do while(.true.)
427  if (m%col(ii)==cl) then
428  m%a(ii) = m%a(ii)+val
429  return
430  end if
431  ii = ii + 1
432  end do
433 
434  end subroutine addentry_2
435 
436 
437 
439  subroutine csr_write (m, fileName, fileFormat)
440  type(csr) , intent(in) :: m
441  character(len=*), intent(in) :: fileName, fileFormat
442 
443  write(*,*)"csr_mod : csr_write = "//trim(filename)
444 
445  select case(trim(fileformat))
446  case("matlab")
447  call matlab_write (m, filename)
448 
449  case default
450  call quit("csr_mod: csr_write: unknown fileFormat")
451  end select
452 
453  end subroutine csr_write
454 
456  subroutine matlab_write (m, fileName)
457  type(csr) , intent(in) :: m
458  character(len=*), intent(in) :: fileName
459 
460  integer :: ii, jj, col, nnz
461 
462 #if DBG>0
463  if (.NOT.valid(m)) call quit( &
464  & 'csr_mod: matlab_write: 1')
465 #endif
466 
467  write(*,*)" Marlab ASCII format"
468 
469  nnz = size(m%a,1)
470  open(unit = 999,file = filename)
471  do ii = 1,m%nl
472  do jj=m%row(ii), m%row(ii+1)-1
473  col = m%col(jj)
474  if (col==0) exit
475  write(999,*) ii, col, m%a(jj)
476  end do
477  end do
478  close(999)
479 
480  end subroutine matlab_write
481 
482 
483 
486  subroutine csr_copy(s2,s1)
487  type(csr), intent(inout) :: s2
488  type(csr), intent(in) :: s1
489 
490 #if DBG>0
491  if (.NOT.valid(s1)) call quit( &
492  & 'csr_mod: csr_copy: 1')
493 #endif
494 
495  call clear(s2)
496 
497  s2%nl = s1%nl
498  s2%nnz = s1%nnz
499  s2%sorted = s1%sorted
500 
501  call allocmem(s2%row, s2%nl+1)
502  call allocmem(s2%col, s2%nnz )
503  call allocmem(s2%a , s2%nnz )
504  s2%row = s1%row
505  s2%col = s1%col
506  s2%a = s1%a
507  s2%width = s1%width
508 
509  if (allocated(s1%diag)) then
510  call allocmem(s2%diag, s2%nl)
511  s2%diag = s1%diag
512  end if
513 
514  end subroutine csr_copy
515 
516 
517 
521  subroutine csr_matvecprod(v,m,u)
522  real(RP), dimension(:), intent(out) :: v
523  type(csr) , intent(in) :: m
524  real(RP), dimension(:), intent(in) :: u
525 
526  real(RP) :: entry
527  integer :: ii, jj
528 
529 #if DBG>0
530  if (.NOT.valid(m)) call quit( &
531  & 'csr_mod: csr_matVecProd: CSR matrix not valid')
532  if (size(v,1) /= m%nl) call quit(&
533  & "csr_mod: csr_matVecProd: output vectr size uncorrect")
534  ! if (size(u,1) /= m%nc ) call quit(&
535  ! & "csr_mod: csr_matVecProd: intput vectr size uncorrect2")
536 #endif
537 
538  !$OMP PARALLEL PRIVATE(entry)
539  !$OMP DO
540  do ii = 1,m%nl
541  ! m%row(ii) = line ii storage strarting index
542  ! m%row(ii+1)-1 = line ii storage ending index
543 
544  entry = 0._rp
545 
546  do jj = m%row(ii), m%row(ii+1)-1
547  entry = entry + m%a(jj)*u(m%col(jj))
548  end do
549 
550  v(ii) = entry
551 
552  end do
553  !$OMP END DO
554  !$OMP END PARALLEL
555 
556  end subroutine csr_matvecprod
557 
558 
562  subroutine csr_add_2(m, m1, a, m2, b)
563  type(csr), intent(inout) :: m
564  type(csr), intent(in) :: m1, m2
565  real(RP) , intent(in) :: a, b
566 
567 #if DBG>0
568 
569  !! checks if patterns are identical
570 
571  logical :: b1, b2
572 
573  if (.NOT.valid(m1)) call quit( &
574  & "csr_mod: csr_add_2: 'm1' not valid")
575  if (.NOT.valid(m2)) call quit( &
576  & "csr_mod: csr_add_2: m2' not valid")
577 
578  b1 = (size(m1%row,1) == size(m2%row,1))
579  b2 = (size(m1%col,1) == size(m2%col,1))
580  if (b1) b1 = (all(m1%row==m2%row))
581  if (b2) b2 = (all(m1%col==m2%col))
582  b1 = (b1.AND.b2)
583  if (.NOT.b1) call quit("csr_mod: csr_add_2:&
584  & 'm1' and 'm2' have different patternsx")
585 #endif
586 
587  call copy(m, m1)
588  m%a = a*m1%a + b*m2%a
589 
590  end subroutine csr_add_2
591 
592 
593 
597  subroutine csr_add_1(m1, a, m2, b)
598  type(csr) , intent(inout) :: m1
599  type(csr) , intent(in) :: m2
600  real(RP), intent(in) :: a, b
601 
602 #if DBG>0
603 
604  !! checks if patterns are identical
605 
606  logical :: b1, b2
607 
608  if (.NOT.valid(m1)) call quit( &
609  & 'csr_mod: csr_add_1: 1')
610  if (.NOT.valid(m2)) call quit( &
611  & 'csr_mod: csr_add_1: 2')
612 
613  b1 = (size(m1%row,1) == size(m2%row,1))
614  b2 = (size(m1%col,1) == size(m2%col,1))
615  if (b1) b1 = (all(m1%row==m2%row))
616  if (b2) b2 = (all(m1%col==m2%col))
617  b1 = (b1.AND.b2)
618  if (.NOT.b1) call quit('csr_mod : csr_add_3: 1')
619 #endif
620 
621  m1%a = a*m1%a + b*m2%a
622 
623  end subroutine csr_add_1
624 
625 
628  subroutine csr_scale(m1, a)
629  type(csr), intent(inout) :: m1
630  real(RP) , intent(in) :: a
631 
632 #if DBG>0
633  if (.NOT.valid(m1)) call quit( &
634  & 'csr_mod: csr_scale: 1')
635 #endif
636 
637  m1%a = m1%a * a
638 
639  end subroutine csr_scale
640 
641 
642  ! %----------------------------------------%
643  ! | |
644  ! | PRECONDITIONNING |
645  ! | |
646  ! %----------------------------------------%
647 
648 
651  subroutine csr_reorder(s, perm, permInv)
652  type(csr) , intent(inout) :: s
653  integer, dimension(:), intent(in) :: perm, permInv
654 
655  type(csr) :: s2
656  integer :: ii, jj, jDeb, jFin, kk, cpt
657 
658 #if DBG>0
659  if (.NOT.valid(s)) call quit( &
660  & 'csr_mod: csr_reorder: 1')
661  if (size(perm, 1)/=s%nl) call quit( &
662  & 'csr_mod: csr_reorder: 2')
663  if (size(perminv, 1)/=s%nl) call quit( &
664  & 'csr_mod: csr_reorder: 3')
665 #endif
666 
667  call copy(s2, s)
668 
669  ! s2 <= reordering of s
670  cpt = 1
671  do ii=1, s%nl
672 
673  kk = perm(ii)
674  s2%row(ii) = cpt
675 
676  jdeb = s%row(kk)
677  jfin = s%row(kk+1) - 1
678 
679  do jj = jdeb,jfin
680 
681  kk = s%col(jj)
682 
683  s2%col(cpt) = perminv(kk)
684  s2%a(cpt) = s%a(jj)
685 
686  cpt = cpt + 1
687  end do
688  end do
689 
690  call copy(s, s2)
691 
692  end subroutine csr_reorder
693 
694 
695 
698  subroutine csr_sortrows(s)
699  type(csr), intent(inout) :: s
700 
701  integer :: ll, j1, j2, sz, ii
702  integer , dimension(s%width) :: new_col,old_col
703  real(RP), dimension(s%width) :: old_a
704 
705  if (.NOT.valid(s)) call quit( &
706  & "csr_mod: csr_sortRows: csr matrix 's' not valid")
707 
708  do ll = 1,s%nl
709  j1 = s%row(ll)
710  j2 = s%row(ll+1)
711  sz = j2 - j1
712  if (sz<=1) cycle
713 
714  j2 = j2 - 1
715  old_col(1:sz) = s%col(j1:j2)
716  old_a(1:sz) = s%a(j1:j2)
717  call shellsort_dec(old_col(1:sz), new_col(1:sz), sz)
718 
719  do ii = 0, sz-1
720  s%col(j1+ii) = old_col(new_col(sz-ii))
721  s%a (j1+ii) = old_a(new_col(sz-ii))
722  end do
723 
724  end do
725 
726  s%sorted = .true.
727 
728  end subroutine csr_sortrows
729 
730 
733  subroutine csr_invllt(x, LLT, b)
734  type(csr) , intent(in) :: LLT
735  real(RP), dimension(:), intent(in) :: b
736  real(RP), dimension(:), intent(out) :: x
737 
738  real(RP), dimension(:), allocatable :: z
739 
740 #if DBG>0
741  if (.NOT.valid(llt)) call quit( &
742  & 'csr_mod: csr_invllt: 1')
743  if (size(b,1)/=llt%nl) call quit( &
744  & 'csr_mod: csr_invllt: 2')
745  if (size(x,1)/=llt%nl) call quit( &
746  & 'csr_mod: csr_invllt: 3')
747 #endif
748 
749  call allocmem(z, llt%nl)
750 
751  call invl(z, llt, b)
752  call invlt(x, llt, z)
753 
754  end subroutine csr_invllt
755 
756 
760  subroutine invlt(x, LT, b)
761  type(csr) , intent(in) :: LT
762  real(RP), dimension(:), intent(in) :: b
763  real(RP), dimension(:), intent(out) :: x
764 
765  integer :: ii,jj, jDeb,jFin
766 
767 #if DBG>0
768  if (.NOT.valid(lt)) call quit( &
769  & 'csr_mod: invLT: 1')
770  if (size(b,1)/=lt%nl) call quit( &
771  & 'csr_mod: invLT: 2')
772  if (size(x,1)/=lt%nl) call quit( &
773  & 'csr_mod: invLT: 3')
774 #endif
775 
776  ! solve LT * x = b
777  x(lt%nl) = b(lt%nl) / lt%a(lt%nnz)
778  do ii = lt%nl-1,1,-1
779 
780  jdeb = lt%diag(ii)
781  jfin = lt%row(ii+1)-1
782 
783  x(ii) = b(ii)
784  do jj = jdeb+1, jfin
785  x(ii) = x(ii) - lt%a(jj) * x(lt%col(jj))
786  end do
787 
788  x(ii) = x(ii) / lt%a(jdeb)
789  end do
790 
791  end subroutine invlt
792 
795  subroutine invl(x,L,b)
796  type(csr) , intent(in) :: L
797  real(RP), dimension(:), intent(in) :: b
798  real(RP), dimension(:), intent(out) :: x
799 
800  integer :: ii,jj, jDeb,jFin
801 
802 #if DBG>0
803  if (.NOT.valid(l)) call quit( &
804  & 'csr_mod: invL: 1')
805  if (size(b,1)/=l%nl) call quit( &
806  & 'csr_mod: invL: 2')
807  if (size(x,1)/=l%nl) call quit( &
808  & 'csr_mod: invL: 3')
809 #endif
810 
811  ! solve L * x = b
812  x(1) = b(1) / l%a(1)
813  do ii = 2, l%nl
814  x(ii) = b(ii)
815  jdeb = l%row(ii)
816  jfin = l%diag(ii)
817 
818  do jj = jdeb, jfin-1
819  x(ii) = x(ii) - l%a(jj) * x(l%col(jj))
820  end do
821 
822  x(ii) = x(ii) / l%a(jfin)
823 
824  end do
825 
826  end subroutine invl
827 
828 
834  subroutine getdiag(A)
835  type(csr), intent(inout) :: A
836 
837  integer :: i, k, kDeb, kFin
838 
839 #if DBG>0
840  if (.NOT.valid(a)) call quit( &
841  & 'csr_mod: getDiag: 1')
842 #endif
843 
844  !! check if the diagonal has already been set:
845  !! If so exit
846  if (allocated( a%diag) ) then
847  if (size(a%diag,1) == a%nl ) return
848  end if
849 
850  !! initialisation
851  call allocmem(a%diag, a%nl)
852  a%diag = -1
853 
854  do i = 1,a%nl
855  kdeb = a%row(i)
856  kfin = a%row(i+1)-1
857  do k = kdeb,kfin
858  if (a%col(k)==i) then
859  a%diag(i) = k
860  exit
861  end if
862  end do
863  end do
864 
865  end subroutine getdiag
866 
867 
869  subroutine csr_transpose(g2, g1)
870  type(csr), intent(inout) :: g2
871  type(csr), intent(in) :: g1
872 
873  integer, dimension(:), allocatable :: nnz
874  integer :: nc, jj, ii, j1, j2
875 
876 #if DBG>0
877  if (.NOT.valid(g1)) call quit( &
878  & 'csr_mod: csr_transpose: 1')
879 #endif
880 
881  call clear(g2)
882 
883  nc = maxval(g1%col)
884 
885  call allocmem(nnz, nc)
886  nnz = 0
887 
888  do ii=1, g1%nnz
889  jj = g1%col(ii)
890  nnz(jj) = nnz(jj) + 1
891  end do
892  g2 = csr(nnz)
893  call freemem(nnz)
894 
895  do ii=1, g1%nl
896 
897  j1 = g1%row(ii)
898  j2 = g1%row(ii+1)-1
899 
900  do jj=j1, j2
901  call addentry(g2, g1%col(jj), ii, g1%a(jj))
902  end do
903 
904  end do
905 
906  end subroutine csr_transpose
907 
908 
912  subroutine add_submatrix_1(m, m2)
913  type(csr) , intent(inout) :: m
914  type(csr) , intent(in) :: m2
915 
916  integer :: j1, j2, ii, jj
917 
918 #if DBG>0
919  if (.NOT.valid(m)) call quit( &
920  & 'csr_mod: add_subMatrix1: 1')
921  if (.NOT.valid(m2)) call quit( &
922  & 'csr_mod: add_subMatrix1: 2')
923 #endif
924 
925  do ii=1, m2%nl
926 
927  j1 = m2%row(ii)
928  j2 = m2%row(ii+1)
929 
930  do jj=j1, j2-1
931  call addentry(m, ii, m2%col(jj), m2%a(jj))
932  end do
933 
934  end do
935 
936  end subroutine add_submatrix_1
937 
938 
942  subroutine add_submatrix_2(m, m2, a)
943  type(csr) , intent(inout) :: m
944  type(csr) , intent(in) :: m2
945  real(RP) , intent(in) :: a
946 
947  integer :: j1, j2, ii, jj
948 
949 #if DBG>0
950  if (.NOT.valid(m)) call quit( &
951  & 'csr_mod: add_subMatrix1: 1')
952  if (.NOT.valid(m2)) call quit( &
953  & 'csr_mod: add_subMatrix1: 2')
954 #endif
955 
956  do ii=1, m2%nl
957 
958  j1 = m2%row(ii)
959  j2 = m2%row(ii+1)
960 
961  do jj=j1, j2-1
962  call addentry(m, ii, m2%col(jj), a*m2%a(jj))
963  end do
964 
965  end do
966 
967  end subroutine add_submatrix_2
968 
969 
970 
973  subroutine csr_print(g)
974  type(csr), intent(in) :: g
975 
976  write(*,*) 'csr_mod : print'
977 
978  write(*,*) ' Valid = ', valid(g)
979  write(*,*) ' Nb lines =' , g%nl
980  write(*,*) ' Nb non-zero entries =' , g%nnz
981  write(*,*) ' Sorted lines = ', g%sorted
982  write(*,*) ' Maximal line width =' , g%width
983 
984  end subroutine csr_print
985 
986 end module csr_mod
987 
subroutine csr_add_2(m, m1, a, m2, b)
m <= a*m1 + b*m2 m1, m2 with identical patterns
Definition: csr_mod.F90:563
subroutine csr_write(m, fileName, fileFormat)
Write matrix m to file &#39;fileName&#39; with format &#39;fileFormat&#39;.
Definition: csr_mod.F90:440
is the csr matrix well defined ?
Definition: csr_mod.F90:101
scale the matrix entries
Definition: csr_mod.F90:150
subroutine csr_addentry(m, rw, cl, val)
Filling matrix m entries.
Definition: csr_mod.F90:385
M2 = M1.
Definition: csr_mod.F90:161
BASIC TOOLS
Definition: basic_tools.F90:8
subroutine csr_sortrows(s)
Sort column indexing per line (increasingly)
Definition: csr_mod.F90:699
deallocate memory for real(RP) arrays
Definition: real_type.F90:163
logical function csr_equal(g1, g2)
test equality
Definition: csr_mod.F90:338
type(csr) function csr_create(nnz)
Constructor for the type csr
Definition: csr_mod.F90:204
write the matrix to a file for a given format
Definition: csr_mod.F90:127
subroutine csr_setrow(m, val, col, ii)
Set line ii.
Definition: csr_mod.F90:319
subroutine csr_invllt(x, LLT, b)
x solution of LLT x = b, L lower diagonal LT = transpose(L)
Definition: csr_mod.F90:734
x solution of LLT x = b, L lower diagonal LT = transpose(L)
Definition: csr_mod.F90:172
subroutine, public quit(message)
Stop program execution, display an error messahe.
subroutine add_submatrix_2(m, m2, a)
m <= m + a*m2 m2(i,j) /=0 ==> m(i,j)/=0
Definition: csr_mod.F90:943
addition of csr matrices (with the same pattern)
Definition: csr_mod.F90:138
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
Definition: real_type.F90:33
transpose matrix of a sparse csr matrix
Definition: csr_mod.F90:166
subroutine csr_clear(m)
Destructor for csr type.
Definition: csr_mod.F90:180
print a short description
Definition: csr_mod.F90:96
subroutine invlt(x, LT, b)
x solution of LT x = b, L lower diagonal LT = transpose(L)
Definition: csr_mod.F90:761
subroutine csr_reorder(s, perm, permInv)
reorder s according to the permutation perm, permInv = perm^{-1}
Definition: csr_mod.F90:652
test csr matrices equality
Definition: csr_mod.F90:122
subroutine csr_transpose(g2, g1)
Transpose csr.
Definition: csr_mod.F90:870
subroutine csr_scale(m1, a)
m1 <= a*m1
Definition: csr_mod.F90:629
Derived type for sparse matrices in CSR format.
Definition: csr_mod.F90:62
ALGEBRAIC OPERATIONS ON SETS
Definition: algebra_set.F90:14
subroutine invl(x, L, b)
x solution of L x = b
Definition: csr_mod.F90:796
subroutine matlab_write(m, fileName)
Write matrix m to file &#39;fileName&#39; matlab ASCII format.
Definition: csr_mod.F90:457
allocate memory for real(RP) arrays
Definition: real_type.F90:158
subroutine, public getdiag(A)
Computes the array Adiag of the indices in Acol for the diagonal coefficients.
Definition: csr_mod.F90:835
reorder rows and cols of a matrix following a permutation of the indexes
Definition: csr_mod.F90:156
add a CSR matrix M1 with pattern P1 to a CSR matrix M with pattern P such that P1 is included in P ...
Definition: csr_mod.F90:145
subroutine csr_add_1(m1, a, m2, b)
m1 = a*m1 + b*m2 m1, m2 with identical patterns
Definition: csr_mod.F90:598
subroutine csr_print(g)
Print a short description.
Definition: csr_mod.F90:974
matrix vector product
Definition: csr_mod.F90:132
subroutine csr_getrow(sz, v, p, n, g, ii)
Extract row ii of the csr matrix g.
Definition: csr_mod.F90:286
subroutine csr_copy(s2, s1)
s2 = s1, copy : s2 is nullified, reallocated and set to s1
Definition: csr_mod.F90:487
subroutine csr_matvecprod(v, m, u)
Matrix Vector Product v <= m*u.
Definition: csr_mod.F90:522
subroutine, public shellsort_dec(t, new_i, n)
Sort integer array : shell sort.
DERIVED TYPE csr for sparse matrices
Definition: csr_mod.F90:13
subroutine, public addentry_2(m, rw, cl, val)
Filling matrix m entries.
Definition: csr_mod.F90:418
destructor
Definition: csr_mod.F90:86
logical function csr_valid(m)
Check if structure contents is correct.
Definition: csr_mod.F90:242
subroutine add_submatrix_1(m, m2)
m <= m + m2 m2(i,j) /=0 ==> m(i,j)/=0
Definition: csr_mod.F90:913
subroutine a(y, x)
Definition: krylov_mod.f90:310