| 1 | ! (C) Copyright 2014- ECMWF. |
|---|
| 2 | ! |
|---|
| 3 | ! This software is licensed under the terms of the Apache Licence Version 2.0 |
|---|
| 4 | ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. |
|---|
| 5 | ! |
|---|
| 6 | ! In applying this licence, ECMWF does not waive the privileges and immunities |
|---|
| 7 | ! granted to it by virtue of its status as an intergovernmental organisation |
|---|
| 8 | ! nor does it submit to any jurisdiction. |
|---|
| 9 | |
|---|
| 10 | ! This simple program tests the functioning of the matrix operations |
|---|
| 11 | ! in the radiation_matrix module |
|---|
| 12 | |
|---|
| 13 | program test_spartacus_math |
|---|
| 14 | |
|---|
| 15 | use parkind1 |
|---|
| 16 | use radiation_matrix |
|---|
| 17 | |
|---|
| 18 | implicit none |
|---|
| 19 | |
|---|
| 20 | integer, parameter :: n = 10, m = 3 |
|---|
| 21 | |
|---|
| 22 | real(jprb), dimension(n,m,m) :: A, B, C |
|---|
| 23 | real(jprb), dimension(n,m) :: v, w |
|---|
| 24 | |
|---|
| 25 | ! Terms in exchange matrix |
|---|
| 26 | real(jprb), dimension(n) :: aa, bb, cc, dd |
|---|
| 27 | |
|---|
| 28 | integer :: j |
|---|
| 29 | |
|---|
| 30 | A = -1.0_jprb |
|---|
| 31 | A(:,1,1) = 10.0_jprb |
|---|
| 32 | A(:,2,2) = 5.0_jprb |
|---|
| 33 | A(:,2,1) = -0.5_jprb |
|---|
| 34 | |
|---|
| 35 | B = 2.0_jprb |
|---|
| 36 | B(:,1,1) = -1.4_jprb |
|---|
| 37 | B(:,2,1) = 1.0e4_jprb |
|---|
| 38 | |
|---|
| 39 | v(:,1) = 2.0_jprb |
|---|
| 40 | v(:,2) = 3.0_jprb |
|---|
| 41 | if (m == 3) then |
|---|
| 42 | v(:,m) = -1.5_jprb |
|---|
| 43 | A(:,m,m) = 4.5 |
|---|
| 44 | end if |
|---|
| 45 | |
|---|
| 46 | write(*,*) 'A =' |
|---|
| 47 | do j = 1,m |
|---|
| 48 | write(*,*) A(1,j,:) |
|---|
| 49 | end do |
|---|
| 50 | |
|---|
| 51 | write(*,*) 'B =' |
|---|
| 52 | do j = 1,m |
|---|
| 53 | write(*,*) B(1,j,:) |
|---|
| 54 | end do |
|---|
| 55 | |
|---|
| 56 | write(*,*) 'v = ', v(1,:) |
|---|
| 57 | |
|---|
| 58 | C = mat_x_mat(n,n,m,A,B) |
|---|
| 59 | w = mat_x_vec(n,n,m,A,v) |
|---|
| 60 | |
|---|
| 61 | write(*,*) 'C = A*B =' |
|---|
| 62 | do j = 1,m |
|---|
| 63 | write(*,*) C(1,j,:) |
|---|
| 64 | end do |
|---|
| 65 | write(*,*) 'w = A*v = ', w(1,:) |
|---|
| 66 | |
|---|
| 67 | w = solve_vec(n,n,m,A,v) |
|---|
| 68 | C = solve_mat(n,n,m,A,B) |
|---|
| 69 | |
|---|
| 70 | |
|---|
| 71 | write(*,*) 'C = A\B =' |
|---|
| 72 | do j = 1,m |
|---|
| 73 | write(*,*) C(1,j,:) |
|---|
| 74 | end do |
|---|
| 75 | write(*,*) 'w = A\v = ', w(1,:) |
|---|
| 76 | |
|---|
| 77 | call expm(n,n,m,A,IMatrixPatternDense) |
|---|
| 78 | |
|---|
| 79 | write(*,*) 'expm(A) =' |
|---|
| 80 | do j = 1,m |
|---|
| 81 | write(*,*) A(1,j,:) |
|---|
| 82 | end do |
|---|
| 83 | |
|---|
| 84 | ! Test fast_expm_exchange_3 |
|---|
| 85 | aa = 2.0_jprb |
|---|
| 86 | bb = 3.0_jprb |
|---|
| 87 | cc = 5.0_jprb |
|---|
| 88 | dd = 7.0_jprb |
|---|
| 89 | |
|---|
| 90 | if (m == 2) then |
|---|
| 91 | |
|---|
| 92 | call fast_expm_exchange(n,n,aa,bb,A) |
|---|
| 93 | |
|---|
| 94 | write(*,*) 'fast_expm(A) = ' |
|---|
| 95 | do j = 1,m |
|---|
| 96 | write(*,*) A(1,j,:) |
|---|
| 97 | end do |
|---|
| 98 | |
|---|
| 99 | A = 0.0_jprb |
|---|
| 100 | A(:,1,1) = -aa |
|---|
| 101 | A(:,2,1) = aa |
|---|
| 102 | A(:,1,2) = bb |
|---|
| 103 | A(:,2,2) = -bb |
|---|
| 104 | |
|---|
| 105 | call expm(n,n,m,A,IMatrixPatternDense) |
|---|
| 106 | |
|---|
| 107 | write(*,*) 'expm(A) = ' |
|---|
| 108 | do j = 1,m |
|---|
| 109 | write(*,*) A(1,j,:) |
|---|
| 110 | end do |
|---|
| 111 | |
|---|
| 112 | ! Test zeros lead to identity matrix |
|---|
| 113 | aa = 0.0_jprb |
|---|
| 114 | call fast_expm_exchange(n,n,aa,aa,A) |
|---|
| 115 | |
|---|
| 116 | write(*,*) 'expm(zeros) = ' |
|---|
| 117 | do j = 1,m |
|---|
| 118 | write(*,*) A(1,j,:) |
|---|
| 119 | end do |
|---|
| 120 | |
|---|
| 121 | else |
|---|
| 122 | |
|---|
| 123 | call fast_expm_exchange(n,n,aa,bb,cc,dd,A) |
|---|
| 124 | |
|---|
| 125 | write(*,*) 'fast_expm(A) = ' |
|---|
| 126 | do j = 1,m |
|---|
| 127 | write(*,*) A(1,j,:) |
|---|
| 128 | end do |
|---|
| 129 | |
|---|
| 130 | A = 0.0_jprb |
|---|
| 131 | A(:,1,1) = -aa |
|---|
| 132 | A(:,2,1) = aa |
|---|
| 133 | A(:,1,2) = bb |
|---|
| 134 | A(:,2,2) = -bb-cc |
|---|
| 135 | A(:,3,2) = cc |
|---|
| 136 | A(:,2,3) = dd |
|---|
| 137 | A(:,3,3) = -dd |
|---|
| 138 | |
|---|
| 139 | call expm(n,n,m,A,IMatrixPatternDense) |
|---|
| 140 | |
|---|
| 141 | write(*,*) 'expm(A) = ' |
|---|
| 142 | do j = 1,m |
|---|
| 143 | write(*,*) A(1,j,:) |
|---|
| 144 | end do |
|---|
| 145 | |
|---|
| 146 | ! Test zeros lead to identity matrix |
|---|
| 147 | aa = 0.0_jprb |
|---|
| 148 | call fast_expm_exchange(n,n,aa,aa,aa,aa,A) |
|---|
| 149 | |
|---|
| 150 | write(*,*) 'expm(zeros) = ' |
|---|
| 151 | do j = 1,m |
|---|
| 152 | write(*,*) A(1,j,:) |
|---|
| 153 | end do |
|---|
| 154 | |
|---|
| 155 | end if |
|---|
| 156 | |
|---|
| 157 | end program test_spartacus_math |
|---|
