Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2011 Bill Paxton & The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module test_square
21 :
22 : use mtx_lib
23 : use mtx_def
24 : use utils_lib, only: mesa_error
25 :
26 : implicit none
27 :
28 : contains
29 :
30 1 : subroutine do_test_square
31 1 : call test_square1
32 1 : call test_square2
33 1 : call test_square_inv
34 1 : end subroutine do_test_square
35 :
36 1 : subroutine test_square_inv
37 :
38 : integer, parameter :: n = 3, nrhs = 1
39 : integer :: info, ipiv(n)
40 33 : double precision :: A1(n, n), B1(n), A2(n, n), B2(n)
41 33 : double precision :: A1_init(n, n), A2_init(n, n), X(n), prod(n)
42 :
43 : include 'formats'
44 :
45 1 : write (*, *) 'test_square_inv'
46 1 : write (*, *)
47 :
48 4 : A1(1, 1:n) = [3.14d0, 7.5d0, 0.00d0]
49 4 : A1(2, 1:n) = [4.1d0, 3.2d0, 0.3d0]
50 4 : A1(3, 1:n) = [0.00d0, 1d0, 4.1d0]
51 1 : A1_init = A1
52 :
53 4 : A2(1, 1:n) = [0d0, 3.1d0, 0d0]
54 4 : A2(2, 1:n) = [4.7d0, 6.2d0, 0d0]
55 4 : A2(3, 1:n) = [3.2d0, 0d0, 0.31d0]
56 1 : A2_init = A2
57 :
58 1 : B1(1:n) = [1.0d0, 2.0d0, 3.0d0]
59 1 : B2(1:n) = [1.1d0, 2.1d0, 3.1d0]
60 :
61 1 : call DGETRF(n, n, A1, n, ipiv, info)
62 1 : if (info /= 0) then
63 0 : write (*, *) 'singular matrix?', info
64 0 : call mesa_error(__FILE__, __LINE__)
65 : end if
66 1 : X = B1
67 : ! solve A1*X = B1
68 1 : call DGETRS('N', n, nrhs, A1, n, ipiv, X, n, info)
69 1 : if (info /= 0) call mesa_error(__FILE__, __LINE__)
70 :
71 1 : write (*, 1) 'B1', B1(1:n)
72 : ! prod = A1_init*X; should get prod == B1
73 1 : call dgemv('N', n, n, 1d0, A1_init, n, X, 1, 0d0, prod, 1)
74 1 : write (*, 1) 'A1_init*X', prod(1:n)
75 : ! prod = A1*X; should get prod == B1
76 :
77 1 : call DGETRF(n, n, A2, n, ipiv, info)
78 1 : if (info /= 0) then
79 0 : write (*, *) 'singular matrix?', info
80 0 : call mesa_error(__FILE__, __LINE__)
81 : end if
82 1 : X = B2
83 : ! solve A2*X = B2
84 1 : call DGETRS('N', n, nrhs, A2, n, ipiv, X, n, info)
85 1 : if (info /= 0) call mesa_error(__FILE__, __LINE__)
86 :
87 1 : write (*, 1) 'B2', B2(1:n)
88 : ! prod = A2_init*X; should get prod == B2
89 1 : call dgemv('N', n, n, 1d0, A2_init, n, X, 1, 0d0, prod, 1)
90 1 : write (*, 1) 'A2_init*X', prod(1:n)
91 : ! prod = A2*X; should get prod == B2
92 :
93 1 : end subroutine test_square_inv
94 :
95 1 : subroutine test_square2
96 :
97 : integer, parameter :: n = 3, nrhs = 1
98 : integer :: info, ipiv(n)
99 33 : double precision :: A1(n, n), B1(n, nrhs), A2(n, n), B2(n, nrhs)
100 :
101 : include 'formats'
102 :
103 1 : write (*, *) 'test_square2'
104 1 : write (*, *)
105 :
106 4 : A1(1, 1:n) = [3.14d0, 7.5d0, 0.00d0]
107 4 : A1(2, 1:n) = [4.1d0, 3.2d0, 0.3d0]
108 4 : A1(3, 1:n) = [0.00d0, 1d0, 4.1d0]
109 :
110 4 : A2(1, 1:n) = [4.7d0, 6.2d0, 0d0]
111 4 : A2(2, 1:n) = [3.2d0, 0d0, 0.31d0]
112 4 : A2(3, 1:n) = [0d0, 3.1d0, 0d0]
113 :
114 1 : B1(1:n, 1) = [1.0d0, 2.0d0, 3.0d0]
115 1 : B2(1:n, 1) = [1.1d0, 2.1d0, 3.1d0]
116 :
117 1 : call DGETRF(n, n, A1, n, ipiv, info)
118 1 : if (info /= 0) then
119 0 : write (*, *) 'singular matrix?', info
120 0 : call mesa_error(__FILE__, __LINE__)
121 : end if
122 1 : call DGETRS('N', n, nrhs, A1, n, ipiv, B1, n, info)
123 1 : if (info /= 0) call mesa_error(__FILE__, __LINE__)
124 :
125 1 : write (*, 1) 'B1', B1(1:n, 1)
126 :
127 1 : call DGETRF(n, n, A2, n, ipiv, info)
128 1 : if (info /= 0) then
129 0 : write (*, *) 'singular matrix?', info
130 0 : call mesa_error(__FILE__, __LINE__)
131 : end if
132 1 : call DGETRS('N', n, nrhs, A2, n, ipiv, B2, n, info)
133 1 : if (info /= 0) call mesa_error(__FILE__, __LINE__)
134 :
135 1 : write (*, 1) 'B2', B2(1:n, 1)
136 1 : write (*, *)
137 :
138 1 : end subroutine test_square2
139 :
140 1 : subroutine test_square1
141 :
142 : integer, parameter :: n = 4, nrhs = 1
143 : integer :: info, ipiv(n)
144 46 : double precision :: A(n, n), B(n, nrhs), A2(n, n)
145 :
146 : include 'formats'
147 :
148 5 : A(1, 1:n) = [1.80d0, 2.88d0, 2.05d0, 0.00d0]
149 5 : A(2, 1:n) = [5.25d0, -2.95d0, -0.95d0, -3.80d0]
150 5 : A(3, 1:n) = [0.00d0, 0.00d0, -2.90d0, -1.04d0]
151 5 : A(4, 1:n) = [-1.11d0, 0.00d0, -0.59d0, 0.80d0]
152 1 : B(1:n, 1) = [4.35d0, 5.05d0, 3.04d0, -2.05d0]
153 :
154 1 : A2 = A
155 :
156 1 : write (*, *) ' test_square1'
157 1 : write (*, *)
158 :
159 1 : call DGETRF(n, n, A, n, ipiv, info)
160 1 : if (info /= 0) then
161 0 : write (*, *) 'singular matrix?', info
162 0 : call mesa_error(__FILE__, __LINE__)
163 : end if
164 :
165 1 : call DGETRS('N', n, nrhs, A, n, ipiv, B, n, info)
166 1 : if (info /= 0) call mesa_error(__FILE__, __LINE__)
167 :
168 1 : write (*, 1) 'B', B(1:n, 1)
169 1 : write (*, *)
170 :
171 1 : end subroutine test_square1
172 :
173 : end module test_square
|
