@@ -25,7 +25,8 @@ contains
25
25
new_unittest('sellc', test_sellc), &
26
26
new_unittest('symmetries', test_symmetries), &
27
27
new_unittest('diagonal', test_diagonal), &
28
- new_unittest('add_get_values', test_add_get_values) &
28
+ new_unittest('add_get_values', test_add_get_values), &
29
+ new_unittest('tridiagonal', test_tridiagonal) &
29
30
]
30
31
end subroutine
31
32
@@ -383,6 +384,43 @@ contains
383
384
#:endfor
384
385
end subroutine
385
386
387
+ subroutine test_tridiagonal(error)
388
+ !> Error handling
389
+ type(error_type), allocatable, intent(out) :: error
390
+ #:for k1, t1, s1 in (KINDS_TYPES)
391
+ #:if k1 != "qp" and k1 != "xdp"
392
+ block
393
+ integer, parameter :: wp = ${k1}$
394
+ integer, parameter :: n = 5
395
+ type(Tridiagonal_${s1}$_type) :: A
396
+ ${t1}$, allocatable :: Amat(:,:), dl(:), dv(:), du(:)
397
+ ${t1}$, allocatable :: x(:)
398
+ ${t1}$, allocatable :: y1(:), y2(:)
399
+
400
+ ! Initialize matrix.
401
+ allocate(dl(n-1), dv(n), du(n-1))
402
+ call random_number(dl) ; call random_number(dv) ; call random_number(du)
403
+ A = Tridiagonal(dl, dv, du) ; Amat = dense(A)
404
+
405
+ ! Random vectors.
406
+ allocate(x(n), source = 0.0_wp) ; call random_number(x)
407
+ allocate(y1(n), source = 0.0_wp) ; allocate(y2(n), source=0.0_wp)
408
+
409
+ ! Test y = A @ x
410
+ y1 = matmul(Amat, x) ; call spmv(A, x, y2)
411
+ call check(error, all(y1 == y2))
412
+ if (allocated(error)) return
413
+
414
+ ! Test y = A.T @ x
415
+ y1 = 0.0_wp ; y2 = 0.0_wp
416
+ y1 = matmul(transpose(Amat), x) ; call spmv(A, x, y2, op="T")
417
+ call check(error, all(y1 == y2))
418
+ if (allocated(error)) return
419
+ end block
420
+ #:endif
421
+ #:endfor
422
+ end subroutine
423
+
386
424
end module
387
425
388
426
0 commit comments