← Back to branch summary

~airpollution/fluidity/fluidity_airpollution

« back to all changes in this revision

Viewing changes to femtools/tests/test_matmul_t_sparse.F90

  • Committer: ziyouzhj
  • Date: 2013-12-09 16:51:29 UTC
  • Revision ID: ziyouzhj@gmail.com-20131209165129-ucoetc3u0atyy05c
airpolution

Show diffs side-by-side

added added

removed removed

Lines of Context:
femtools/tests/test_matmul_t_sparse.F90
 
1
subroutine test_matmul_t_sparse
 
2
 
 
3
  use vector_tools
 
4
  use sparse_tools
 
5
  use unittest_tools
 
6
  implicit none
 
7
 
 
8
  integer, parameter :: size_mat = 50
 
9
  integer, parameter :: n_samples = 100
 
10
  real, dimension(size_mat,size_mat) :: A, B, C, D
 
11
  type(dynamic_csr_matrix) :: A_d, B_d, C_d
 
12
  type(csr_matrix) :: A_c, B_c, C_dc
 
13
  integer :: i, j, k, n
 
14
  logical :: fail,fail1,fail2, fail3
 
15
  character(len=size_mat) :: buf
 
16
  real, dimension(4) :: rand0
 
17
  
 
18
  do k=1,5
 
19
 
 
20
     call allocate(A_d, size_mat,size_mat)
 
21
     call allocate(B_d, size_mat,size_mat)
 
22
 
 
23
    write(buf,'(i0)') k
 
24
    fail = .false.
 
25
    fail1 = .false.
 
26
    fail2 = .false.
 
27
    fail3 = .false.
 
28
 
 
29
    A = 0.
 
30
    B = 0.
 
31
    call zero(A_d)
 
32
    call zero(B_d)
 
33
 
 
34
    do n = 1, n_samples
 
35
       call random_number(rand0)
 
36
       i = ceiling(rand0(1)*size_mat)
 
37
       j = ceiling(rand0(2)*size_mat)
 
38
       A(i,j) = A(i,j) +1.0
 
39
       B(i,j) = B(i,j) + 1.0
 
40
       !call addto(A_d,i,j,rand0(3))
 
41
       !call addto(B_d,i,j,rand0(4))
 
42
       call addto(A_d,i,j,1.0)
 
43
       call addto(B_d,i,j,1.0)
 
44
    end do
 
45
 
 
46
    if (any(A /= dense(A_d))) then
 
47
       write(0,*) 'A assembly bungled'
 
48
       write(0,*) '---'
 
49
       write(0,*) A
 
50
       write(0,*) '---'
 
51
       write(0,*) dense(A_d)
 
52
       stop
 
53
    end if
 
54
    if (any(B /= dense(B_d))) then
 
55
       write(0,*) 'B assembly bungled'
 
56
       stop
 
57
    end if
 
58
 
 
59
    A_c = dcsr2csr(A_d)
 
60
    B_c = dcsr2csr(B_d)
 
61
 
 
62
    C_d = matmul_T(A_d,B_d)
 
63
    !C_c = matmul_T(A_c,B_c)
 
64
 
 
65
    C_dc = dcsr2csr(C_d)
 
66
 
 
67
    C=dense(C_d)
 
68
    !E = dense(C_c)
 
69
 
 
70
    if (any(abs(C-dense(C_dc))>1.0e-14)) fail3 = .true.
 
71
 
 
72
 
 
73
    do i = 1, size(B_d%colm)
 
74
       if(any(C_d%val(i)%ptr<0.5)) fail2=.true.
 
75
    end do
 
76
 
 
77
    call deallocate(C_d)
 
78
    call deallocate(A_c)
 
79
    call deallocate(B_c)
 
80
    !call deallocate(C_c)
 
81
    call deallocate(C_dc)
 
82
 
 
83
    D = matmul(A, transpose(B))
 
84
 
 
85
    if (any(abs(C-D)>1.0e-14)) fail = .true.
 
86
    !if (any(abs(E-D)>1.0e-14)) fail2 = .true.
 
87
    do i = 1, size(A_d%colm)
 
88
       do j = 2, size(A_d%colm(i)%ptr)
 
89
          if(A_d%colm(i)%ptr(j).le.A_d%colm(i)%ptr(j-1)) then
 
90
             fail1 = .true.
 
91
          end if
 
92
       end do
 
93
    end do
 
94
 
 
95
    do i = 1, size(B_d%colm)
 
96
       do j = 2, size(B_d%colm(i)%ptr)
 
97
          if(B_d%colm(i)%ptr(j).le.B_d%colm(i)%ptr(j-1)) then
 
98
             fail1 = .true.
 
99
          end if
 
100
       end do
 
101
    end do
 
102
 
 
103
    call report_test("[matmul_t_sparse dcsr " // trim(buf) // "]", fail, .false., "The output of matmul_t and matmul should be identical.")
 
104
    call report_test("[matmul_t_sparse dcsr dense " // trim(buf) // "]", fail3, .false., "The output of dense(dscr_matrix) and dense(dcsr2csr(dcsr_matrix)) should be identical.")
 
105
    call report_test("[matmul_t_sparse zeros " // trim(buf) // "]", fail2, .false., "The sparsity pattern should not contain zeros")
 
106
    call report_test("[matmul_t_sparse " // trim(buf) // " ordering]", fail1, .false., "We expect the rows to be incrementally ordered in dcsr matrices")
 
107
 
 
108
     call deallocate(A_d)
 
109
     call deallocate(B_d)
 
110
 
 
111
  end do
 
112
 
 
113
end subroutine test_matmul_t_sparse