65 logical :: sorted = .false.
66 integer ,
dimension(:),
allocatable :: row
67 integer ,
dimension(:),
allocatable :: col
69 real(RP),
dimension(:),
allocatable ::
a 70 integer ,
dimension(:),
allocatable :: diag
180 type(
csr),
intent(inout) :: m
205 integer,
dimension(:),
intent(in) :: nnz
207 integer :: nz, ii, jj, wdt, length
227 if (length>wdt) wdt = length
234 &
'csr_mod: csr_create')
242 type(
csr),
intent(in) :: m
245 integer :: ii, wdt, length
248 if (.NOT. (
allocated(m%row)) )
return 249 if (.NOT. (
allocated(m%row)) )
return 250 if (.NOT. (
allocated(m%a)) )
return 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 )
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 )
265 length = m%row(ii+1) - m%row(ii)
266 if (length>wdt) wdt = length
268 b = ( wdt == m%width )
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
296 &
"csr_mod: csr_getRow_1: CSR matrix 'g' not valid")
298 if ( (11<=0) .OR. (ii>g%nl) )
call quit( &
299 &
"csr_mod: csr_getRow_1: wrong line index 'ii'")
307 if (sz>n )
call quit( &
308 &
"csr_mod: csr_getRow_1: unsufficient size& 309 & for the array 'p'")
312 p(1:sz) = g%col(j1:j2-1)
313 v(1:sz) = g%a(j1:j2-1)
319 type(
csr) ,
intent(inout):: m
320 real(RP),
dimension(:),
intent(in) :: val
321 integer ,
dimension(:),
intent(in) :: col
322 integer ,
intent(in) :: ii
339 type(
csr),
intent(in) :: g1, g2
341 integer ,
dimension(g1%width) :: p1, p2, q1, q2
342 real(RP),
dimension(g1%width) :: a1, a2
343 integer :: ii, jj, sz, wdt
346 b = b .AND. ( g1%nl == g2%nl ) .AND. ( g1%nnz == g2%nnz )
349 b = b .AND. ( g1%width == g2%width )
352 b = b .AND. (all( g1%row == g2%row ) )
357 call getrow(sz, a1, p1, wdt, g1, ii)
358 call getrow(sz, a2, p2, wdt, g2, ii)
366 b = b .AND. ( p1(q1(jj)) == p2(q2(jj)) )
367 b = b .AND.
equal(a1(q1(jj)), a2(q2(jj)))
385 type(
csr),
intent(inout) :: m
386 integer ,
intent(in) :: rw, cl
387 real(RP),
intent(in) :: val
391 do ii = m%row(rw),m%row(rw+1)-1
393 if (m%col(ii)==0)
then 398 else if (m%col(ii)==cl)
then 400 m%a(ii) = m%a(ii)+val
407 call quit(
"csr_mod: csr_addEntry: 1")
418 type(
csr),
intent(inout) :: m
419 integer ,
intent(in) :: rw, cl
420 real(RP),
intent(in) :: val
427 if (m%col(ii)==cl)
then 428 m%a(ii) = m%a(ii)+val
439 subroutine csr_write (m, fileName, fileFormat)
440 type(
csr) ,
intent(in) :: m
441 character(len=*),
intent(in) :: fileName, fileFormat
443 write(*,*)
"csr_mod : csr_write = "//trim(filename)
445 select case(trim(fileformat))
450 call quit(
"csr_mod: csr_write: unknown fileFormat")
457 type(
csr) ,
intent(in) :: m
458 character(len=*),
intent(in) :: fileName
460 integer :: ii, jj, col, nnz
464 &
'csr_mod: matlab_write: 1')
467 write(*,*)
" Marlab ASCII format" 470 open(unit = 999,file = filename)
472 do jj=m%row(ii), m%row(ii+1)-1
475 write(999,*) ii, col, m%a(jj)
487 type(
csr),
intent(inout) :: s2
488 type(
csr),
intent(in) :: s1
492 &
'csr_mod: csr_copy: 1')
499 s2%sorted = s1%sorted
509 if (
allocated(s1%diag))
then 522 real(RP),
dimension(:),
intent(out) :: v
523 type(
csr) ,
intent(in) :: m
524 real(RP),
dimension(:),
intent(in) :: u
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")
546 do jj = m%row(ii), m%row(ii+1)-1
547 entry = entry + m%a(jj)*u(m%col(jj))
563 type(
csr),
intent(inout) :: m
564 type(
csr),
intent(in) :: m1, m2
565 real(RP) ,
intent(in) :: a, b
574 &
"csr_mod: csr_add_2: 'm1' not valid")
576 &
"csr_mod: csr_add_2: m2' not valid")
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))
583 if (.NOT.b1)
call quit(
"csr_mod: csr_add_2:& 584 & 'm1' and 'm2' have different patternsx")
588 m%a = a*m1%a + b*m2%a
598 type(
csr) ,
intent(inout) :: m1
599 type(
csr) ,
intent(in) :: m2
600 real(RP),
intent(in) :: a, b
609 &
'csr_mod: csr_add_1: 1')
611 &
'csr_mod: csr_add_1: 2')
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))
618 if (.NOT.b1)
call quit(
'csr_mod : csr_add_3: 1')
621 m1%a = a*m1%a + b*m2%a
629 type(
csr),
intent(inout) :: m1
630 real(RP) ,
intent(in) :: a
634 &
'csr_mod: csr_scale: 1')
652 type(
csr) ,
intent(inout) :: s
653 integer,
dimension(:),
intent(in) :: perm, permInv
656 integer :: ii, jj, jDeb, jFin, kk, cpt
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')
677 jfin = s%row(kk+1) - 1
683 s2%col(cpt) = perminv(kk)
699 type(
csr),
intent(inout) :: s
701 integer :: ll, j1, j2, sz, ii
702 integer ,
dimension(s%width) :: new_col,old_col
703 real(RP),
dimension(s%width) :: old_a
706 &
"csr_mod: csr_sortRows: csr matrix 's' not valid")
715 old_col(1:sz) = s%col(j1:j2)
716 old_a(1:sz) = s%a(j1:j2)
720 s%col(j1+ii) = old_col(new_col(sz-ii))
721 s%a (j1+ii) = old_a(new_col(sz-ii))
734 type(
csr) ,
intent(in) :: LLT
735 real(RP),
dimension(:),
intent(in) :: b
736 real(RP),
dimension(:),
intent(out) :: x
738 real(RP),
dimension(:),
allocatable :: z
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')
752 call invlt(x, llt, z)
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
765 integer :: ii,jj, jDeb,jFin
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')
777 x(lt%nl) = b(lt%nl) / lt%a(lt%nnz)
781 jfin = lt%row(ii+1)-1
785 x(ii) = x(ii) - lt%a(jj) * x(lt%col(jj))
788 x(ii) = x(ii) / lt%a(jdeb)
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
800 integer :: ii,jj, jDeb,jFin
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')
819 x(ii) = x(ii) - l%a(jj) * x(l%col(jj))
822 x(ii) = x(ii) / l%a(jfin)
835 type(
csr),
intent(inout) :: A
837 integer :: i, k, kDeb, kFin
841 &
'csr_mod: getDiag: 1')
846 if (
allocated( a%diag) )
then 847 if (
size(a%diag,1) == a%nl )
return 858 if (a%col(k)==i)
then 870 type(
csr),
intent(inout) :: g2
871 type(
csr),
intent(in) :: g1
873 integer,
dimension(:),
allocatable :: nnz
874 integer :: nc, jj, ii, j1, j2
878 &
'csr_mod: csr_transpose: 1')
890 nnz(jj) = nnz(jj) + 1
901 call addentry(g2, g1%col(jj), ii, g1%a(jj))
913 type(
csr) ,
intent(inout) :: m
914 type(
csr) ,
intent(in) :: m2
916 integer :: j1, j2, ii, jj
920 &
'csr_mod: add_subMatrix1: 1')
922 &
'csr_mod: add_subMatrix1: 2')
931 call addentry(m, ii, m2%col(jj), m2%a(jj))
943 type(
csr) ,
intent(inout) :: m
944 type(
csr) ,
intent(in) :: m2
945 real(RP) ,
intent(in) :: a
947 integer :: j1, j2, ii, jj
951 &
'csr_mod: add_subMatrix1: 1')
953 &
'csr_mod: add_subMatrix1: 2')
962 call addentry(m, ii, m2%col(jj), a*m2%a(jj))
974 type(
csr),
intent(in) :: g
976 write(*,*)
'csr_mod : print' 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
subroutine csr_add_2(m, m1, a, m2, b)
m <= a*m1 + b*m2 m1, m2 with identical patterns
subroutine csr_write(m, fileName, fileFormat)
Write matrix m to file 'fileName' with format 'fileFormat'.
is the csr matrix well defined ?
subroutine csr_addentry(m, rw, cl, val)
Filling matrix m entries.
subroutine csr_sortrows(s)
Sort column indexing per line (increasingly)
deallocate memory for real(RP) arrays
logical function csr_equal(g1, g2)
test equality
type(csr) function csr_create(nnz)
Constructor for the type csr
write the matrix to a file for a given format
subroutine csr_setrow(m, val, col, ii)
Set line ii.
subroutine csr_invllt(x, LLT, b)
x solution of LLT x = b, L lower diagonal LT = transpose(L)
x solution of LLT x = b, L lower diagonal LT = transpose(L)
subroutine add_submatrix_2(m, m2, a)
m <= m + a*m2 m2(i,j) /=0 ==> m(i,j)/=0
addition of csr matrices (with the same pattern)
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
transpose matrix of a sparse csr matrix
subroutine csr_clear(m)
Destructor for csr type.
print a short description
subroutine invlt(x, LT, b)
x solution of LT x = b, L lower diagonal LT = transpose(L)
subroutine csr_reorder(s, perm, permInv)
reorder s according to the permutation perm, permInv = perm^{-1}
test csr matrices equality
subroutine csr_transpose(g2, g1)
Transpose csr.
subroutine csr_scale(m1, a)
m1 <= a*m1
Derived type for sparse matrices in CSR format.
ALGEBRAIC OPERATIONS ON SETS
subroutine invl(x, L, b)
x solution of L x = b
subroutine matlab_write(m, fileName)
Write matrix m to file 'fileName' matlab ASCII format.
allocate memory for real(RP) arrays
subroutine, public getdiag(A)
Computes the array Adiag of the indices in Acol for the diagonal coefficients.
reorder rows and cols of a matrix following a permutation of the indexes
add a CSR matrix M1 with pattern P1 to a CSR matrix M with pattern P such that P1 is included in P ...
subroutine csr_add_1(m1, a, m2, b)
m1 = a*m1 + b*m2 m1, m2 with identical patterns
subroutine csr_print(g)
Print a short description.
subroutine csr_getrow(sz, v, p, n, g, ii)
Extract row ii of the csr matrix g.
subroutine csr_copy(s2, s1)
s2 = s1, copy : s2 is nullified, reallocated and set to s1
subroutine csr_matvecprod(v, m, u)
Matrix Vector Product v <= m*u.
subroutine, public shellsort_dec(t, new_i, n)
Sort integer array : shell sort.
DERIVED TYPE csr for sparse matrices
subroutine, public addentry_2(m, rw, cl, val)
Filling matrix m entries.
logical function csr_valid(m)
Check if structure contents is correct.
subroutine add_submatrix_1(m, m2)
m <= m + m2 m2(i,j) /=0 ==> m(i,j)/=0