Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! copyright (c) 2012 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 my_lapack95
21 :
22 : use const_def, only: dp
23 : use utils_lib, only: mesa_error
24 : implicit none
25 : integer, parameter :: fltp = dp
26 :
27 : contains
28 :
29 7 : subroutine my_gemv(m,n,a,lda,x,y) ! y = y - a*x
30 : integer :: lda,m,n
31 : real(fltp) :: a(:,:) ! (lda,*)
32 : real(fltp) :: x(:), y(:)
33 7 : real(fltp) :: tmp
34 : real(fltp), parameter :: zero=0
35 : ! trans = 'n'
36 : ! alpha = -1
37 : ! beta = 1
38 : ! incx = 1
39 : ! incy = 1
40 : integer :: j, i
41 182 : do j = 1,n
42 175 : tmp = x(j)
43 182 : if (tmp/=zero) then
44 4550 : do i = 1,m
45 4550 : y(i) = y(i) - tmp*a(i,j)
46 : end do
47 : end if
48 : end do
49 7 : end subroutine my_gemv
50 :
51 :
52 14 : subroutine my_gemv_mv(m,n,a,x,b,z,y) ! y = y - a*x - b*z
53 : integer :: m, n
54 : real(fltp) :: a(:,:), b(:,:)
55 : real(fltp) :: x(:), z(:), y(:)
56 14 : real(fltp) :: tmp_x, tmp_z
57 : real(fltp), parameter :: zero=0
58 : integer :: j, i
59 364 : do j = 1,n
60 350 : tmp_x = x(j)
61 350 : tmp_z = z(j)
62 364 : if (tmp_x/=zero) then
63 350 : if (tmp_z /= zero) then
64 9100 : do i = 1,m
65 9100 : y(i) = y(i) - tmp_x*a(i,j) - tmp_z*b(i,j)
66 : end do
67 : else
68 0 : do i = 1,m
69 0 : y(i) = y(i) - tmp_x*a(i,j)
70 : end do
71 : end if
72 0 : else if (tmp_z /= zero) then
73 0 : do i = 1,m
74 0 : y(i) = y(i) - tmp_z*b(i,j)
75 : end do
76 : end if
77 : end do
78 14 : end subroutine my_gemv_mv
79 :
80 :
81 61 : subroutine my_gemv_p1(m,n,a,lda,x,y) ! y = y + a*x
82 : integer :: lda,m,n
83 : real(fltp) :: a(:,:) ! (lda,*)
84 : real(fltp) :: x(:), y(:)
85 61 : real(fltp) :: tmp
86 : real(fltp), parameter :: zero=0
87 : ! trans = 'n'
88 : ! alpha = -1
89 : ! beta = 1
90 : ! incx = 1
91 : ! incy = 1
92 : integer :: j, i
93 1586 : do j = 1,n
94 1525 : tmp = x(j)
95 1586 : if (tmp/=zero) then
96 39650 : do i = 1,m
97 39650 : y(i) = y(i) + tmp*a(i,j)
98 : end do
99 : end if
100 : end do
101 61 : end subroutine my_gemv_p1
102 :
103 :
104 16 : subroutine my_gemv_p_mv(m,n,a,x,b,z,y) ! y = y + a*x + b*z
105 : integer :: m, n
106 : real(fltp) :: a(:,:), b(:,:)
107 : real(fltp) :: x(:), z(:), y(:)
108 16 : real(fltp) :: tmp_x, tmp_z
109 : real(fltp), parameter :: zero=0
110 : integer :: j, i
111 416 : do j = 1,n
112 400 : tmp_x = x(j)
113 400 : tmp_z = z(j)
114 416 : if (tmp_x/=zero) then
115 400 : if (tmp_z /= zero) then
116 10400 : do i = 1,m
117 10400 : y(i) = y(i) + tmp_x*a(i,j) + tmp_z*b(i,j)
118 : end do
119 : else
120 0 : do i = 1,m
121 0 : y(i) = y(i) + tmp_x*a(i,j)
122 : end do
123 : end if
124 0 : else if (tmp_z /= zero) then
125 0 : do i = 1,m
126 0 : y(i) = y(i) + tmp_z*b(i,j)
127 : end do
128 : end if
129 : end do
130 16 : end subroutine my_gemv_p_mv
131 :
132 :
133 0 : subroutine my_gemm(m,n,k,a,lda,b,ldb,c,ldc) ! c := c - a*b
134 : integer, intent(in) :: k,lda,ldb,ldc,m,n
135 : real(fltp), dimension(:,:) :: a, b, c ! a(lda,*),b(ldb,*),c(ldc,*)
136 0 : real(fltp) :: tmp
137 : real(fltp), parameter :: zero=0
138 : integer :: j, i, l
139 : ! transa = 'n'
140 : ! transb = 'n'
141 : ! alpha = -1
142 : ! beta = 1
143 : ! assumes other args are valid
144 0 : do j = 1,n
145 0 : do l = 1,k
146 0 : tmp = b(l,j)
147 0 : if (tmp /= zero) then
148 0 : do i = 1,m
149 0 : c(i,j) = c(i,j) - tmp*a(i,l)
150 : end do
151 : end if
152 : end do
153 : end do
154 0 : end subroutine my_gemm
155 :
156 :
157 42 : subroutine my_gemm_p1(m,n,k,a,lda,b,ldb,c,ldc) ! c := c + a*b
158 : integer, intent(in) :: k,lda,ldb,ldc,m,n
159 : real(fltp), dimension(:,:) :: a, b, c ! a(lda,*),b(ldb,*),c(ldc,*)
160 42 : real(fltp) :: tmp
161 : real(fltp), parameter :: zero=0
162 : integer :: j, i, l
163 : ! transa = 'n'
164 : ! transb = 'n'
165 : ! alpha = 1
166 : ! beta = 1
167 : ! assumes other args are valid
168 1092 : do j = 1,n
169 27342 : do l = 1,k
170 26250 : tmp = b(l,j)
171 27300 : if (tmp /= zero) then
172 57876 : do i = 1,m
173 57876 : c(i,j) = c(i,j) + tmp*a(i,l)
174 : end do
175 : end if
176 : end do
177 : end do
178 42 : end subroutine my_gemm_p1
179 :
180 :
181 14 : subroutine my_gemm_plus_mm(m,n,k,a,b,d,e,c) ! c := c + a*b + d*e
182 : integer, intent(in) :: k,m,n
183 : real(fltp), dimension(:,:) :: a, b, c, d, e
184 14 : real(fltp) :: tmp_b, tmp_e
185 : real(fltp), parameter :: zero=0
186 : integer :: j, i, l
187 364 : do j = 1,n
188 9114 : do l = 1,k
189 8750 : tmp_b = b(l,j)
190 8750 : tmp_e = e(l,j)
191 9100 : if (tmp_b /= zero) then
192 742 : if (tmp_e /= zero) then
193 14196 : do i = 1,m
194 14196 : c(i,j) = c(i,j) + tmp_b*a(i,l) + tmp_e*d(i,l)
195 : end do
196 : else
197 5096 : do i = 1,m
198 5096 : c(i,j) = c(i,j) + tmp_b*a(i,l)
199 : end do
200 : end if
201 8008 : else if (tmp_e /= zero) then
202 5096 : do i = 1,m
203 5096 : c(i,j) = c(i,j) + tmp_e*d(i,l)
204 : end do
205 : end if
206 : end do
207 : end do
208 14 : end subroutine my_gemm_plus_mm
209 :
210 :
211 0 : subroutine my_gemm0(m,n,k,a,lda,b,ldb,c,ldc)
212 : ! c := -a*b
213 : integer, intent(in) :: k,lda,ldb,ldc,m,n
214 : real(fltp), dimension(:,:) :: a, b, c ! a(lda,*),b(ldb,*),c(ldc,*)
215 : integer :: j, i
216 : real(fltp), parameter :: zero=0
217 : include 'formats'
218 : ! transa = 'n'
219 : ! transb = 'n'
220 : ! alpha = -1
221 : ! beta = 0
222 : ! assumes other args are valid
223 0 : do j=1,n
224 0 : do i=1,m
225 0 : c(i,j) = zero
226 : end do
227 : end do
228 0 : call my_gemm(m,n,k,a,lda,b,ldb,c,ldc)
229 0 : end subroutine my_gemm0
230 :
231 :
232 35 : subroutine my_gemm0_p1(m,n,k,a,lda,b,ldb,c,ldc)
233 : ! c := -a*b
234 : integer, intent(in) :: k,lda,ldb,ldc,m,n
235 : real(fltp), dimension(:,:) :: a, b, c ! a(lda,*),b(ldb,*),c(ldc,*)
236 : integer :: j, i
237 : real(fltp), parameter :: zero=0
238 : include 'formats'
239 : ! transa = 'n'
240 : ! transb = 'n'
241 : ! alpha = -1
242 : ! beta = 0
243 : ! assumes other args are valid
244 910 : do j=1,n
245 22785 : do i=1,m
246 22750 : c(i,j) = zero
247 : end do
248 : end do
249 35 : call my_gemm_p1(m,n,k,a,lda,b,ldb,c,ldc)
250 35 : end subroutine my_gemm0_p1
251 :
252 :
253 20 : subroutine my_getf2(m, a, lda, ipiv, sfmin, info)
254 : integer :: info, lda, m
255 : integer :: ipiv(:)
256 : real(fltp) :: a(:,:)
257 : real(fltp) :: sfmin
258 : real(fltp), parameter :: one=1, zero=0
259 : integer :: i, j, jp, ii, jj, n, mm
260 20 : real(fltp) :: tmp, da
261 20 : if (m == 4) then
262 0 : call my_getf2_4_by_4(a, lda, ipiv, sfmin, info)
263 0 : return
264 20 : else if (m == 5) then
265 0 : call my_getf2_5_by_5(a, lda, ipiv, sfmin, info)
266 0 : return
267 : end if
268 520 : do j = 1, m
269 500 : info = 0
270 7000 : jp = j - 1 + maxloc(abs(a(j:lda,j)),dim=1)
271 500 : ipiv( j ) = jp
272 500 : if( a( jp, j )/=zero ) then
273 500 : if( jp/=j ) then ! swap a(j,:) and a(jp,:)
274 0 : do i=1,m
275 0 : tmp = a(j,i)
276 0 : a(j,i) = a(jp,i)
277 0 : a(jp,i) = tmp
278 : end do
279 : end if
280 500 : if( j<m ) then
281 480 : if( abs(a( j, j )) >= sfmin ) then
282 480 : da = one / a( j, j )
283 480 : n = m-j
284 480 : mm = mod(n,5)
285 480 : if (mm /= 0) then
286 1400 : do i = 1,mm
287 1400 : a(j+i,j) = da*a(j+i,j)
288 : end do
289 : end if
290 480 : if (n >= 5) then
291 400 : do i = mm + 1,n,5
292 1000 : a(j+i,j) = da*a(j+i,j)
293 1000 : a(j+i+1,j) = da*a(j+i+1,j)
294 1000 : a(j+i+2,j) = da*a(j+i+2,j)
295 1000 : a(j+i+3,j) = da*a(j+i+3,j)
296 1000 : a(j+i+4,j) = da*a(j+i+4,j)
297 : end do
298 : end if
299 : else ! no scale
300 0 : do i = 1, m-j
301 0 : a( j+i, j ) = a( j+i, j ) / a( j, j )
302 : end do
303 : end if
304 : end if
305 : else if( info==0 ) then
306 0 : info = j
307 : end if
308 520 : if( j<m ) then
309 : !call dger( m-j, m-j, -one, a( j+1, j ), 1, a( j, j+1 ), lda, a( j+1, j+1 ), lda )
310 6480 : do jj = j+1, m
311 104480 : do ii = j+1, m
312 104000 : a(ii,jj) = a(ii,jj) - a(ii,j)*a(j,jj)
313 : end do
314 : end do
315 : end if
316 : end do
317 : end subroutine my_getf2
318 :
319 :
320 0 : subroutine my_getf2_4_by_4(a, lda, ipiv, sfmin, info)
321 : integer :: info, lda ! m=4
322 : integer :: ipiv(:)
323 : real(fltp) :: a(:,:)
324 : real(fltp) :: sfmin
325 : real(fltp), parameter :: one=1, zero=0
326 : integer :: jp
327 0 : real(fltp) :: tmp, da
328 0 : info = 0
329 :
330 0 : jp = maxloc(abs(a(1:lda,1)),dim=1)
331 0 : ipiv( 1 ) = jp
332 0 : if( a( jp, 1 )/=zero ) then
333 0 : if( jp/=1 ) then ! swap a(1,:) and a(jp,:)
334 0 : tmp = a(1,1)
335 0 : a(1,1) = a(jp,1)
336 0 : a(jp,1) = tmp
337 0 : tmp = a(1,2)
338 0 : a(1,2) = a(jp,2)
339 0 : a(jp,2) = tmp
340 0 : tmp = a(1,3)
341 0 : a(1,3) = a(jp,3)
342 0 : a(jp,3) = tmp
343 0 : tmp = a(1,4)
344 0 : a(1,4) = a(jp,4)
345 0 : a(jp,4) = tmp
346 : end if
347 0 : if( abs(a( 1, 1 )) >= sfmin ) then
348 0 : da = one / a( 1, 1 )
349 0 : a(2,1) = da*a(2,1)
350 0 : a(3,1) = da*a(3,1)
351 0 : a(4,1) = da*a(4,1)
352 : else ! no scale
353 0 : a( 2, 1 ) = a( 2, 1 ) / a( 1, 1 )
354 0 : a( 3, 1 ) = a( 3, 1 ) / a( 1, 1 )
355 0 : a( 4, 1 ) = a( 4, 1 ) / a( 1, 1 )
356 : end if
357 : else if( info==0 ) then
358 0 : info = 1
359 : end if
360 0 : a(2,2) = a(2,2) - a(2,1)*a(1,2)
361 0 : a(3,2) = a(3,2) - a(3,1)*a(1,2)
362 0 : a(4,2) = a(4,2) - a(4,1)*a(1,2)
363 0 : a(2,3) = a(2,3) - a(2,1)*a(1,3)
364 0 : a(3,3) = a(3,3) - a(3,1)*a(1,3)
365 0 : a(4,3) = a(4,3) - a(4,1)*a(1,3)
366 0 : a(2,4) = a(2,4) - a(2,1)*a(1,4)
367 0 : a(3,4) = a(3,4) - a(3,1)*a(1,4)
368 0 : a(4,4) = a(4,4) - a(4,1)*a(1,4)
369 :
370 0 : jp = 1 + maxloc(abs(a(2:lda,2)),dim=1)
371 0 : ipiv( 2 ) = jp
372 0 : if( a( jp, 2 )/=zero ) then
373 0 : if( jp/=2 ) then ! swap a(2,:) and a(jp,:)
374 0 : tmp = a(2,1)
375 0 : a(2,1) = a(jp,1)
376 0 : a(jp,1) = tmp
377 0 : tmp = a(2,2)
378 0 : a(2,2) = a(jp,2)
379 0 : a(jp,2) = tmp
380 0 : tmp = a(2,3)
381 0 : a(2,3) = a(jp,3)
382 0 : a(jp,3) = tmp
383 0 : tmp = a(2,4)
384 0 : a(2,4) = a(jp,4)
385 0 : a(jp,4) = tmp
386 : end if
387 0 : if( abs(a( 2, 2 )) >= sfmin ) then
388 0 : da = one / a( 2, 2 )
389 0 : a(3,2) = da*a(3,2)
390 0 : a(4,2) = da*a(4,2)
391 : else ! no scale
392 0 : a( 3, 2 ) = a( 3, 2 ) / a( 2, 2 )
393 0 : a( 4, 2 ) = a( 4, 2 ) / a( 2, 2 )
394 : end if
395 0 : else if( info==0 ) then
396 0 : info = 2
397 : end if
398 0 : a(3,3) = a(3,3) - a(3,2)*a(2,3)
399 0 : a(4,3) = a(4,3) - a(4,2)*a(2,3)
400 0 : a(3,4) = a(3,4) - a(3,2)*a(2,4)
401 0 : a(4,4) = a(4,4) - a(4,2)*a(2,4)
402 :
403 0 : jp = 2 + maxloc(abs(a(3:lda,3)),dim=1)
404 0 : ipiv( 3 ) = jp
405 0 : if( a( jp, 3 )/=zero ) then
406 0 : if( jp/=3 ) then ! swap a(3,:) and a(jp,:)
407 0 : tmp = a(3,1)
408 0 : a(3,1) = a(jp,1)
409 0 : a(jp,1) = tmp
410 0 : tmp = a(3,2)
411 0 : a(3,2) = a(jp,2)
412 0 : a(jp,2) = tmp
413 0 : tmp = a(3,3)
414 0 : a(3,3) = a(jp,3)
415 0 : a(jp,3) = tmp
416 0 : tmp = a(3,4)
417 0 : a(3,4) = a(jp,4)
418 0 : a(jp,4) = tmp
419 : end if
420 0 : if( abs(a( 3, 3 )) >= sfmin ) then
421 0 : da = one / a( 3, 3 )
422 0 : a(4,3) = da*a(4,3)
423 : else ! no scale
424 0 : a( 4, 3 ) = a( 4, 3 ) / a( 3, 3 )
425 : end if
426 0 : else if( info==0 ) then
427 0 : info = 3
428 : end if
429 0 : a(4,4) = a(4,4) - a(4,3)*a(3,4)
430 :
431 0 : jp = 3 + maxloc(abs(a(4:lda,4)),dim=1)
432 0 : ipiv( 4 ) = jp
433 0 : if( a( jp, 4 )/=zero ) then
434 0 : if( jp/=4 ) then ! swap a(4,:) and a(jp,:)
435 0 : tmp = a(4,1)
436 0 : a(4,1) = a(jp,1)
437 0 : a(jp,1) = tmp
438 0 : tmp = a(4,2)
439 0 : a(4,2) = a(jp,2)
440 0 : a(jp,2) = tmp
441 0 : tmp = a(4,3)
442 0 : a(4,3) = a(jp,3)
443 0 : a(jp,3) = tmp
444 0 : tmp = a(4,4)
445 0 : a(4,4) = a(jp,4)
446 0 : a(jp,4) = tmp
447 : end if
448 0 : else if( info==0 ) then
449 0 : info = 4
450 : end if
451 :
452 0 : end subroutine my_getf2_4_by_4
453 :
454 :
455 0 : subroutine my_getf2_5_by_5(a, lda, ipiv, sfmin, info)
456 : integer :: info, lda ! m=5
457 : integer :: ipiv(:)
458 : real(fltp) :: a(:,:)
459 : real(fltp) :: sfmin
460 : real(fltp), parameter :: one=1, zero=0
461 : integer :: jp
462 0 : real(fltp) :: tmp, da
463 0 : info = 0
464 0 : jp = maxloc(abs(a(1:lda,1)),dim=1)
465 0 : ipiv( 1 ) = jp
466 0 : if( a( jp, 1 )/=zero ) then
467 0 : if( jp/=1 ) then ! swap a(1,:) and a(jp,:)
468 0 : tmp = a(1,1)
469 0 : a(1,1) = a(jp,1)
470 0 : a(jp,1) = tmp
471 0 : tmp = a(1,2)
472 0 : a(1,2) = a(jp,2)
473 0 : a(jp,2) = tmp
474 0 : tmp = a(1,3)
475 0 : a(1,3) = a(jp,3)
476 0 : a(jp,3) = tmp
477 0 : tmp = a(1,4)
478 0 : a(1,4) = a(jp,4)
479 0 : a(jp,4) = tmp
480 0 : tmp = a(1,5)
481 0 : a(1,5) = a(jp,5)
482 0 : a(jp,5) = tmp
483 : end if
484 0 : if( abs(a( 1, 1 )) >= sfmin ) then
485 0 : da = one / a( 1, 1 )
486 0 : a(2,1) = da*a(2,1)
487 0 : a(3,1) = da*a(3,1)
488 0 : a(4,1) = da*a(4,1)
489 0 : a(5,1) = da*a(5,1)
490 : else ! no scale
491 0 : a( 2, 1 ) = a( 2, 1 ) / a( 1, 1 )
492 0 : a( 3, 1 ) = a( 3, 1 ) / a( 1, 1 )
493 0 : a( 4, 1 ) = a( 4, 1 ) / a( 1, 1 )
494 0 : a( 5, 1 ) = a( 5, 1 ) / a( 1, 1 )
495 : end if
496 : else if( info==0 ) then
497 0 : info = 1
498 : end if
499 0 : a(2,2) = a(2,2) - a(2,1)*a(1,2)
500 0 : a(3,2) = a(3,2) - a(3,1)*a(1,2)
501 0 : a(4,2) = a(4,2) - a(4,1)*a(1,2)
502 0 : a(5,2) = a(5,2) - a(5,1)*a(1,2)
503 0 : a(2,3) = a(2,3) - a(2,1)*a(1,3)
504 0 : a(3,3) = a(3,3) - a(3,1)*a(1,3)
505 0 : a(4,3) = a(4,3) - a(4,1)*a(1,3)
506 0 : a(5,3) = a(5,3) - a(5,1)*a(1,3)
507 0 : a(2,4) = a(2,4) - a(2,1)*a(1,4)
508 0 : a(3,4) = a(3,4) - a(3,1)*a(1,4)
509 0 : a(4,4) = a(4,4) - a(4,1)*a(1,4)
510 0 : a(5,4) = a(5,4) - a(5,1)*a(1,4)
511 0 : a(2,5) = a(2,5) - a(2,1)*a(1,5)
512 0 : a(3,5) = a(3,5) - a(3,1)*a(1,5)
513 0 : a(4,5) = a(4,5) - a(4,1)*a(1,5)
514 0 : a(5,5) = a(5,5) - a(5,1)*a(1,5)
515 :
516 0 : jp = 1 + maxloc(abs(a(2:lda,2)),dim=1)
517 0 : ipiv( 2 ) = jp
518 0 : if( a( jp, 2 )/=zero ) then
519 0 : if( jp/=2 ) then ! swap a(2,:) and a(jp,:)
520 0 : tmp = a(2,1)
521 0 : a(2,1) = a(jp,1)
522 0 : a(jp,1) = tmp
523 0 : tmp = a(2,2)
524 0 : a(2,2) = a(jp,2)
525 0 : a(jp,2) = tmp
526 0 : tmp = a(2,3)
527 0 : a(2,3) = a(jp,3)
528 0 : a(jp,3) = tmp
529 0 : tmp = a(2,4)
530 0 : a(2,4) = a(jp,4)
531 0 : a(jp,4) = tmp
532 0 : tmp = a(2,5)
533 0 : a(2,5) = a(jp,5)
534 0 : a(jp,5) = tmp
535 : end if
536 0 : if( abs(a( 2, 2 )) >= sfmin ) then
537 0 : da = one / a( 2, 2 )
538 0 : a(3,2) = da*a(3,2)
539 0 : a(4,2) = da*a(4,2)
540 0 : a(5,2) = da*a(5,2)
541 : else ! no scale
542 0 : a( 3, 2 ) = a( 3, 2 ) / a( 2, 2 )
543 0 : a( 4, 2 ) = a( 4, 2 ) / a( 2, 2 )
544 0 : a( 5, 2 ) = a( 5, 2 ) / a( 2, 2 )
545 : end if
546 0 : else if( info==0 ) then
547 0 : info = 2
548 : end if
549 0 : a(3,3) = a(3,3) - a(3,2)*a(2,3)
550 0 : a(4,3) = a(4,3) - a(4,2)*a(2,3)
551 0 : a(5,3) = a(5,3) - a(5,2)*a(2,3)
552 0 : a(3,4) = a(3,4) - a(3,2)*a(2,4)
553 0 : a(4,4) = a(4,4) - a(4,2)*a(2,4)
554 0 : a(5,4) = a(5,4) - a(5,2)*a(2,4)
555 0 : a(3,5) = a(3,5) - a(3,2)*a(2,5)
556 0 : a(4,5) = a(4,5) - a(4,2)*a(2,5)
557 0 : a(5,5) = a(5,5) - a(5,2)*a(2,5)
558 :
559 0 : jp = 2 + maxloc(abs(a(3:lda,3)),dim=1)
560 0 : ipiv( 3 ) = jp
561 0 : if( a( jp, 3 )/=zero ) then
562 0 : if( jp/=3 ) then ! swap a(3,:) and a(jp,:)
563 0 : tmp = a(3,1)
564 0 : a(3,1) = a(jp,1)
565 0 : a(jp,1) = tmp
566 0 : tmp = a(3,2)
567 0 : a(3,2) = a(jp,2)
568 0 : a(jp,2) = tmp
569 0 : tmp = a(3,3)
570 0 : a(3,3) = a(jp,3)
571 0 : a(jp,3) = tmp
572 0 : tmp = a(3,4)
573 0 : a(3,4) = a(jp,4)
574 0 : a(jp,4) = tmp
575 0 : tmp = a(3,5)
576 0 : a(3,5) = a(jp,5)
577 0 : a(jp,5) = tmp
578 : end if
579 0 : if( abs(a( 3, 3 )) >= sfmin ) then
580 0 : da = one / a( 3, 3 )
581 0 : a(4,3) = da*a(4,3)
582 0 : a(5,3) = da*a(5,3)
583 : else ! no scale
584 0 : a( 4, 3 ) = a( 4, 3 ) / a( 3, 3 )
585 0 : a( 4, 3 ) = a( 4, 3 ) / a( 3, 3 )
586 : end if
587 0 : else if( info==0 ) then
588 0 : info = 3
589 : end if
590 0 : a(4,4) = a(4,4) - a(4,3)*a(3,4)
591 0 : a(5,4) = a(5,4) - a(5,3)*a(3,4)
592 0 : a(4,5) = a(4,5) - a(4,3)*a(3,5)
593 0 : a(5,5) = a(5,5) - a(5,3)*a(3,5)
594 :
595 0 : jp = 3 + maxloc(abs(a(4:lda,4)),dim=1)
596 0 : ipiv( 4 ) = jp
597 0 : if( a( jp, 4 )/=zero ) then
598 0 : if( jp/=4 ) then ! swap a(4,:) and a(jp,:)
599 0 : tmp = a(4,1)
600 0 : a(4,1) = a(jp,1)
601 0 : a(jp,1) = tmp
602 0 : tmp = a(4,2)
603 0 : a(4,2) = a(jp,2)
604 0 : a(jp,2) = tmp
605 0 : tmp = a(4,3)
606 0 : a(4,3) = a(jp,3)
607 0 : a(jp,3) = tmp
608 0 : tmp = a(4,4)
609 0 : a(4,4) = a(jp,4)
610 0 : a(jp,4) = tmp
611 0 : tmp = a(4,5)
612 0 : a(4,5) = a(jp,5)
613 0 : a(jp,5) = tmp
614 : end if
615 0 : if( abs(a( 4, 4 )) >= sfmin ) then
616 0 : da = one / a( 4, 4 )
617 0 : a(5,4) = da*a(5,4)
618 : else ! no scale
619 0 : a( 5, 4 ) = a( 5, 4 ) / a( 4, 4 )
620 : end if
621 0 : else if( info==0 ) then
622 0 : info = 4
623 : end if
624 0 : a(5,5) = a(5,5) - a(5,4)*a(4,5)
625 :
626 0 : jp = 4 + maxloc(abs(a(5:lda,5)),dim=1)
627 0 : ipiv( 5 ) = jp
628 0 : if( a( jp, 5 )/=zero ) then
629 0 : if( jp/=5 ) then ! swap a(5,:) and a(jp,:)
630 0 : tmp = a(5,1)
631 0 : a(5,1) = a(jp,1)
632 0 : a(jp,1) = tmp
633 0 : tmp = a(5,2)
634 0 : a(5,2) = a(jp,2)
635 0 : a(jp,2) = tmp
636 0 : tmp = a(5,3)
637 0 : a(5,3) = a(jp,3)
638 0 : a(jp,3) = tmp
639 0 : tmp = a(5,4)
640 0 : a(5,4) = a(jp,4)
641 0 : a(jp,4) = tmp
642 0 : tmp = a(5,5)
643 0 : a(5,5) = a(jp,5)
644 0 : a(jp,5) = tmp
645 : end if
646 0 : else if( info==0 ) then
647 0 : info = 5
648 : end if
649 :
650 0 : end subroutine my_getf2_5_by_5
651 :
652 :
653 58 : subroutine my_laswp( n, a, lda, k1, k2, ipiv, incx )
654 : integer :: incx, k1, k2, lda, n
655 : integer :: ipiv(:)
656 : real(fltp) :: a(:,:) ! a( lda, * )
657 : integer :: i, i1, i2, inc, ip, ix, ix0, j, k, n32
658 58 : real(fltp) :: temp
659 : ! interchange row i with row ipiv(i) for each of rows k1 through k2.
660 58 : if( incx>0 ) then
661 58 : ix0 = k1
662 58 : i1 = k1
663 58 : i2 = k2
664 58 : inc = 1
665 0 : else if( incx<0 ) then
666 0 : ix0 = 1 + ( 1-k2 )*incx
667 0 : i1 = k2
668 0 : i2 = k1
669 0 : inc = -1
670 : else
671 : return
672 : end if
673 58 : n32 = ( n / 32 )*32
674 58 : if( n32/=0 ) then
675 0 : do j = 1, n32, 32
676 0 : ix = ix0
677 0 : do i = i1, i2, inc
678 0 : ip = ipiv( ix )
679 0 : if( ip/=i ) then
680 0 : do k = j, j + 31
681 0 : temp = a( i, k )
682 0 : a( i, k ) = a( ip, k )
683 0 : a( ip, k ) = temp
684 : end do
685 : end if
686 0 : ix = ix + incx
687 : end do
688 : end do
689 : end if
690 58 : if( n32/=n ) then
691 58 : n32 = n32 + 1
692 58 : ix = ix0
693 58 : do i = i1, i2, inc
694 1450 : ip = ipiv( ix )
695 1450 : if( ip/=i ) then
696 0 : do k = n32, n
697 0 : temp = a( i, k )
698 0 : a( i, k ) = a( ip, k )
699 0 : a( ip, k ) = temp
700 : end do
701 : end if
702 1450 : ix = ix + incx
703 : end do
704 : end if
705 : end subroutine my_laswp
706 :
707 :
708 0 : subroutine my_laswp_4_by_1( a, lda, ipiv )
709 : ! n == 1, k1 == 1, k2 == 4, incx == 1
710 : integer :: lda
711 : integer :: ipiv(:)
712 : real(fltp) :: a(:,:) ! a( lda, * )
713 : integer :: ip
714 0 : real(fltp) :: temp
715 : ! interchange row i with row ipiv(i) for each of rows k1 through k2.
716 0 : ip = ipiv( 1 )
717 0 : if( ip/=1 ) then
718 0 : temp = a( 1, 1 )
719 0 : a( 1, 1 ) = a( ip, 1 )
720 0 : a( ip, 1 ) = temp
721 : end if
722 0 : ip = ipiv( 2 )
723 0 : if( ip/=2 ) then
724 0 : temp = a( 2, 1 )
725 0 : a( 2, 1 ) = a( ip, 1 )
726 0 : a( ip, 1 ) = temp
727 : end if
728 0 : ip = ipiv( 3 )
729 0 : if( ip/=3 ) then
730 0 : temp = a( 3, 1 )
731 0 : a( 3, 1 ) = a( ip, 1 )
732 0 : a( ip, 1 ) = temp
733 : end if
734 0 : ip = ipiv( 4 )
735 0 : if( ip/=4 ) then
736 0 : temp = a( 4, 1 )
737 0 : a( 4, 1 ) = a( ip, 1 )
738 0 : a( ip, 1 ) = temp
739 : end if
740 0 : end subroutine my_laswp_4_by_1
741 :
742 :
743 0 : subroutine my_laswp_5_by_1( a, lda, ipiv )
744 : ! n == 1, k1 == 1, k2 == 5, incx == 1
745 : integer :: lda
746 : integer :: ipiv(:)
747 : real(fltp) :: a(:,:) ! a( lda, * )
748 : integer :: ip
749 0 : real(fltp) :: temp
750 : ! interchange row i with row ipiv(i) for each of rows k1 through k2.
751 0 : ip = ipiv( 1 )
752 0 : if( ip/=1 ) then
753 0 : temp = a( 1, 1 )
754 0 : a( 1, 1 ) = a( ip, 1 )
755 0 : a( ip, 1 ) = temp
756 : end if
757 0 : ip = ipiv( 2 )
758 0 : if( ip/=2 ) then
759 0 : temp = a( 2, 1 )
760 0 : a( 2, 1 ) = a( ip, 1 )
761 0 : a( ip, 1 ) = temp
762 : end if
763 0 : ip = ipiv( 3 )
764 0 : if( ip/=3 ) then
765 0 : temp = a( 3, 1 )
766 0 : a( 3, 1 ) = a( ip, 1 )
767 0 : a( ip, 1 ) = temp
768 : end if
769 0 : ip = ipiv( 4 )
770 0 : if( ip/=4 ) then
771 0 : temp = a( 4, 1 )
772 0 : a( 4, 1 ) = a( ip, 1 )
773 0 : a( ip, 1 ) = temp
774 : end if
775 0 : ip = ipiv( 5 )
776 0 : if( ip/=5 ) then
777 0 : temp = a( 5, 1 )
778 0 : a( 5, 1 ) = a( ip, 1 )
779 0 : a( ip, 1 ) = temp
780 : end if
781 0 : end subroutine my_laswp_5_by_1
782 :
783 :
784 0 : subroutine my_laswp_4_by_4( a, lda, ipiv )
785 : integer :: lda
786 : integer :: ipiv(:)
787 : real(fltp) :: a(:,:) ! a( lda, * )
788 0 : real(fltp) :: temp
789 : integer :: ip
790 : ! interchange row i with row ipiv(i) for each of rows 1 through 4.
791 0 : ip = ipiv( 1 )
792 0 : if( ip/=1 ) then
793 0 : temp = a( 1, 1 )
794 0 : a( 1, 1 ) = a( ip, 1 )
795 0 : a( ip, 1 ) = temp
796 0 : temp = a( 1, 2 )
797 0 : a( 1, 2 ) = a( ip, 2 )
798 0 : a( ip, 2 ) = temp
799 0 : temp = a( 1, 3 )
800 0 : a( 1, 3 ) = a( ip, 3 )
801 0 : a( ip, 3 ) = temp
802 0 : temp = a( 1, 4 )
803 0 : a( 1, 4 ) = a( ip, 4 )
804 0 : a( ip, 4 ) = temp
805 : end if
806 0 : ip = ipiv( 2 )
807 0 : if( ip/=2 ) then
808 0 : temp = a( 2, 1 )
809 0 : a( 2, 1 ) = a( ip, 1 )
810 0 : a( ip, 1 ) = temp
811 0 : temp = a( 2, 2 )
812 0 : a( 2, 2 ) = a( ip, 2 )
813 0 : a( ip, 2 ) = temp
814 0 : temp = a( 2, 3 )
815 0 : a( 2, 3 ) = a( ip, 3 )
816 0 : a( ip, 3 ) = temp
817 0 : temp = a( 2, 4 )
818 0 : a( 2, 4 ) = a( ip, 4 )
819 0 : a( ip, 4 ) = temp
820 : end if
821 0 : ip = ipiv( 3 )
822 0 : if( ip/=3 ) then
823 0 : temp = a( 3, 1 )
824 0 : a( 3, 1 ) = a( ip, 1 )
825 0 : a( ip, 1 ) = temp
826 0 : temp = a( 3, 2 )
827 0 : a( 3, 2 ) = a( ip, 2 )
828 0 : a( ip, 2 ) = temp
829 0 : temp = a( 3, 3 )
830 0 : a( 3, 3 ) = a( ip, 3 )
831 0 : a( ip, 3 ) = temp
832 0 : temp = a( 3, 4 )
833 0 : a( 3, 4 ) = a( ip, 4 )
834 0 : a( ip, 4 ) = temp
835 : end if
836 0 : ip = ipiv( 4 )
837 0 : if( ip/=4 ) then
838 0 : temp = a( 4, 1 )
839 0 : a( 4, 1 ) = a( ip, 1 )
840 0 : a( ip, 1 ) = temp
841 0 : temp = a( 4, 2 )
842 0 : a( 4, 2 ) = a( ip, 2 )
843 0 : a( ip, 2 ) = temp
844 0 : temp = a( 4, 3 )
845 0 : a( 4, 3 ) = a( ip, 3 )
846 0 : a( ip, 3 ) = temp
847 0 : temp = a( 4, 4 )
848 0 : a( 4, 4 ) = a( ip, 4 )
849 0 : a( ip, 4 ) = temp
850 : end if
851 0 : end subroutine my_laswp_4_by_4
852 :
853 :
854 0 : subroutine my_laswp_5_by_5( a, lda, ipiv )
855 : integer :: lda
856 : integer :: ipiv(:)
857 : real(fltp) :: a(:,:) ! a( lda, * )
858 0 : real(fltp) :: temp
859 : integer :: ip
860 : ! interchange row i with row ipiv(i) for each of rows 1 through 5.
861 0 : ip = ipiv( 1 )
862 0 : if( ip/=1 ) then
863 0 : temp = a( 1, 1 )
864 0 : a( 1, 1 ) = a( ip, 1 )
865 0 : a( ip, 1 ) = temp
866 0 : temp = a( 1, 2 )
867 0 : a( 1, 2 ) = a( ip, 2 )
868 0 : a( ip, 2 ) = temp
869 0 : temp = a( 1, 3 )
870 0 : a( 1, 3 ) = a( ip, 3 )
871 0 : a( ip, 3 ) = temp
872 0 : temp = a( 1, 4 )
873 0 : a( 1, 4 ) = a( ip, 4 )
874 0 : a( ip, 4 ) = temp
875 0 : temp = a( 1, 5 )
876 0 : a( 1, 5 ) = a( ip, 5 )
877 0 : a( ip, 5 ) = temp
878 : end if
879 0 : ip = ipiv( 2 )
880 0 : if( ip/=2 ) then
881 0 : temp = a( 2, 1 )
882 0 : a( 2, 1 ) = a( ip, 1 )
883 0 : a( ip, 1 ) = temp
884 0 : temp = a( 2, 2 )
885 0 : a( 2, 2 ) = a( ip, 2 )
886 0 : a( ip, 2 ) = temp
887 0 : temp = a( 2, 3 )
888 0 : a( 2, 3 ) = a( ip, 3 )
889 0 : a( ip, 3 ) = temp
890 0 : temp = a( 2, 4 )
891 0 : a( 2, 4 ) = a( ip, 4 )
892 0 : a( ip, 4 ) = temp
893 0 : temp = a( 2, 5 )
894 0 : a( 2, 5 ) = a( ip, 5 )
895 0 : a( ip, 5 ) = temp
896 : end if
897 0 : ip = ipiv( 3 )
898 0 : if( ip/=3 ) then
899 0 : temp = a( 3, 1 )
900 0 : a( 3, 1 ) = a( ip, 1 )
901 0 : a( ip, 1 ) = temp
902 0 : temp = a( 3, 2 )
903 0 : a( 3, 2 ) = a( ip, 2 )
904 0 : a( ip, 2 ) = temp
905 0 : temp = a( 3, 3 )
906 0 : a( 3, 3 ) = a( ip, 3 )
907 0 : a( ip, 3 ) = temp
908 0 : temp = a( 3, 4 )
909 0 : a( 3, 4 ) = a( ip, 4 )
910 0 : a( ip, 4 ) = temp
911 0 : temp = a( 3, 5 )
912 0 : a( 3, 5 ) = a( ip, 5 )
913 0 : a( ip, 5 ) = temp
914 : end if
915 0 : ip = ipiv( 4 )
916 0 : if( ip/=4 ) then
917 0 : temp = a( 4, 1 )
918 0 : a( 4, 1 ) = a( ip, 1 )
919 0 : a( ip, 1 ) = temp
920 0 : temp = a( 4, 2 )
921 0 : a( 4, 2 ) = a( ip, 2 )
922 0 : a( ip, 2 ) = temp
923 0 : temp = a( 4, 3 )
924 0 : a( 4, 3 ) = a( ip, 3 )
925 0 : a( ip, 3 ) = temp
926 0 : temp = a( 4, 4 )
927 0 : a( 4, 4 ) = a( ip, 4 )
928 0 : a( ip, 4 ) = temp
929 0 : temp = a( 4, 5 )
930 0 : a( 4, 5 ) = a( ip, 5 )
931 0 : a( ip, 5 ) = temp
932 : end if
933 0 : ip = ipiv( 5 )
934 0 : if( ip/=5 ) then
935 0 : temp = a( 5, 1 )
936 0 : a( 5, 1 ) = a( ip, 1 )
937 0 : a( ip, 1 ) = temp
938 0 : temp = a( 5, 2 )
939 0 : a( 5, 2 ) = a( ip, 2 )
940 0 : a( ip, 2 ) = temp
941 0 : temp = a( 5, 3 )
942 0 : a( 5, 3 ) = a( ip, 3 )
943 0 : a( ip, 3 ) = temp
944 0 : temp = a( 5, 4 )
945 0 : a( 5, 4 ) = a( ip, 4 )
946 0 : a( ip, 4 ) = temp
947 0 : temp = a( 5, 5 )
948 0 : a( 5, 5 ) = a( ip, 5 )
949 0 : a( ip, 5 ) = temp
950 : end if
951 0 : end subroutine my_laswp_5_by_5
952 :
953 :
954 58 : subroutine my_getrs( n, nrhs, a, lda, ipiv, b, ldb, info )
955 : integer :: info, lda, ldb, n, nrhs
956 : integer, pointer :: ipiv(:)
957 : real(fltp), pointer :: a(:,:), b(:,:) ! a( lda, * ), b( ldb, * )
958 : real(fltp), parameter :: one=1, zero=0
959 : integer :: i, j, k
960 58 : info = 0
961 :
962 58 : if (nrhs == 1) then
963 20 : if (n == 4) then
964 0 : call my_getrs_4_by_1( a, lda, ipiv, b, ldb, info )
965 0 : return
966 : end if
967 20 : if (n == 5) then
968 0 : call my_getrs_5_by_1( a, lda, ipiv, b, ldb, info )
969 0 : return
970 : end if
971 38 : else if (nrhs == 4 .and. n == 4) then
972 0 : call my_getrs_4_by_4( a, lda, ipiv, b, ldb, info )
973 0 : return
974 38 : else if (nrhs == 5 .and. n == 5) then
975 0 : call my_getrs_5_by_5( a, lda, ipiv, b, ldb, info )
976 0 : return
977 : end if
978 :
979 58 : call my_laswp(nrhs, b, ldb, 1, n, ipiv, 1 )
980 : !call dtrsm( 'left', 'lower', 'no transpose', 'unit', n, nrhs, one, a, lda, b, ldb )
981 1028 : do j = 1,nrhs
982 25278 : do k = 1,n
983 25220 : if (b(k,j)/=zero) then
984 31567 : do i = k + 1,n
985 31567 : b(i,j) = b(i,j) - b(k,j)*a(i,k)
986 : end do
987 : end if
988 : end do
989 : end do
990 : !call dtrsm( 'left', 'upper', 'no transpose', 'non-unit', n, nrhs, one, a, lda, b, ldb )
991 1028 : do j = 1,nrhs
992 25278 : do k = n,1,-1
993 25220 : if (b(k,j)/=zero) then
994 2514 : b(k,j) = b(k,j)/a(k,k)
995 23562 : do i = 1,k - 1
996 23562 : b(i,j) = b(i,j) - b(k,j)*a(i,k)
997 : end do
998 : end if
999 : end do
1000 : end do
1001 :
1002 : end subroutine my_getrs
1003 :
1004 :
1005 0 : subroutine my_getrs_5_by_5( a, lda, ipiv, b, ldb, info )
1006 : integer :: info, lda, ldb ! , n=5, nrhs=5
1007 : integer, pointer :: ipiv(:)
1008 : real(fltp), pointer :: a(:,:), b(:,:) ! a( lda, * ), b( ldb, * )
1009 : real(fltp), parameter :: zero=0
1010 :
1011 0 : info = 0
1012 :
1013 : !call my_laswp(5, b, ldb, 1, 5, ipiv, 1 )
1014 0 : call my_laswp_5_by_5( b, ldb, ipiv )
1015 :
1016 0 : b(2,1) = b(2,1) - b(1,1)*a(2,1)
1017 0 : b(3,1) = b(3,1) - b(1,1)*a(3,1)
1018 0 : b(4,1) = b(4,1) - b(1,1)*a(4,1)
1019 0 : b(5,1) = b(5,1) - b(1,1)*a(5,1)
1020 0 : b(3,1) = b(3,1) - b(2,1)*a(3,2)
1021 0 : b(4,1) = b(4,1) - b(2,1)*a(4,2)
1022 0 : b(5,1) = b(5,1) - b(2,1)*a(5,2)
1023 0 : b(4,1) = b(4,1) - b(3,1)*a(4,3)
1024 0 : b(5,1) = b(5,1) - b(3,1)*a(5,3)
1025 0 : b(5,1) = b(5,1) - b(4,1)*a(5,4)
1026 :
1027 0 : b(2,2) = b(2,2) - b(1,2)*a(2,1)
1028 0 : b(3,2) = b(3,2) - b(1,2)*a(3,1)
1029 0 : b(4,2) = b(4,2) - b(1,2)*a(4,1)
1030 0 : b(5,2) = b(5,2) - b(1,2)*a(5,1)
1031 0 : b(3,2) = b(3,2) - b(2,2)*a(3,2)
1032 0 : b(4,2) = b(4,2) - b(2,2)*a(4,2)
1033 0 : b(5,2) = b(5,2) - b(2,2)*a(5,2)
1034 0 : b(4,2) = b(4,2) - b(3,2)*a(4,3)
1035 0 : b(5,2) = b(5,2) - b(3,2)*a(5,3)
1036 0 : b(5,2) = b(5,2) - b(4,2)*a(5,4)
1037 :
1038 0 : b(2,3) = b(2,3) - b(1,3)*a(2,1)
1039 0 : b(3,3) = b(3,3) - b(1,3)*a(3,1)
1040 0 : b(4,3) = b(4,3) - b(1,3)*a(4,1)
1041 0 : b(5,3) = b(5,3) - b(1,3)*a(5,1)
1042 0 : b(3,3) = b(3,3) - b(2,3)*a(3,2)
1043 0 : b(4,3) = b(4,3) - b(2,3)*a(4,2)
1044 0 : b(5,3) = b(5,3) - b(2,3)*a(5,2)
1045 0 : b(4,3) = b(4,3) - b(3,3)*a(4,3)
1046 0 : b(5,3) = b(5,3) - b(3,3)*a(5,3)
1047 0 : b(5,3) = b(5,3) - b(4,3)*a(5,4)
1048 :
1049 0 : b(2,4) = b(2,4) - b(1,4)*a(2,1)
1050 0 : b(3,4) = b(3,4) - b(1,4)*a(3,1)
1051 0 : b(4,4) = b(4,4) - b(1,4)*a(4,1)
1052 0 : b(5,4) = b(5,4) - b(1,4)*a(5,1)
1053 0 : b(3,4) = b(3,4) - b(2,4)*a(3,2)
1054 0 : b(4,4) = b(4,4) - b(2,4)*a(4,2)
1055 0 : b(5,4) = b(5,4) - b(2,4)*a(5,2)
1056 0 : b(4,4) = b(4,4) - b(3,4)*a(4,3)
1057 0 : b(5,4) = b(5,4) - b(3,4)*a(5,3)
1058 0 : b(5,4) = b(5,4) - b(4,4)*a(5,4)
1059 :
1060 0 : b(2,5) = b(2,5) - b(1,5)*a(2,1)
1061 0 : b(3,5) = b(3,5) - b(1,5)*a(3,1)
1062 0 : b(4,5) = b(4,5) - b(1,5)*a(4,1)
1063 0 : b(5,5) = b(5,5) - b(1,5)*a(5,1)
1064 0 : b(3,5) = b(3,5) - b(2,5)*a(3,2)
1065 0 : b(4,5) = b(4,5) - b(2,5)*a(4,2)
1066 0 : b(5,5) = b(5,5) - b(2,5)*a(5,2)
1067 0 : b(4,5) = b(4,5) - b(3,5)*a(4,3)
1068 0 : b(5,5) = b(5,5) - b(3,5)*a(5,3)
1069 0 : b(5,5) = b(5,5) - b(4,5)*a(5,4)
1070 :
1071 : !call dtrsm( 'left', 'upper', 'no transpose', 'non-unit', n, nrhs, one, a, lda, b, ldb )
1072 0 : b(5,1) = b(5,1)/a(5,5)
1073 0 : b(1,1) = b(1,1) - b(5,1)*a(1,5)
1074 0 : b(2,1) = b(2,1) - b(5,1)*a(2,5)
1075 0 : b(3,1) = b(3,1) - b(5,1)*a(3,5)
1076 0 : b(4,1) = b(4,1) - b(5,1)*a(4,5)
1077 0 : b(4,1) = b(4,1)/a(4,4)
1078 0 : b(1,1) = b(1,1) - b(4,1)*a(1,4)
1079 0 : b(2,1) = b(2,1) - b(4,1)*a(2,4)
1080 0 : b(3,1) = b(3,1) - b(4,1)*a(3,4)
1081 0 : b(3,1) = b(3,1)/a(3,3)
1082 0 : b(1,1) = b(1,1) - b(3,1)*a(1,3)
1083 0 : b(2,1) = b(2,1) - b(3,1)*a(2,3)
1084 0 : b(2,1) = b(2,1)/a(2,2)
1085 0 : b(1,1) = b(1,1) - b(2,1)*a(1,2)
1086 0 : b(1,1) = b(1,1)/a(1,1)
1087 :
1088 0 : b(5,2) = b(5,2)/a(5,5)
1089 0 : b(1,2) = b(1,2) - b(5,2)*a(1,5)
1090 0 : b(2,2) = b(2,2) - b(5,2)*a(2,5)
1091 0 : b(3,2) = b(3,2) - b(5,2)*a(3,5)
1092 0 : b(4,2) = b(4,2) - b(5,2)*a(4,5)
1093 0 : b(4,2) = b(4,2)/a(4,4)
1094 0 : b(1,2) = b(1,2) - b(4,2)*a(1,4)
1095 0 : b(2,2) = b(2,2) - b(4,2)*a(2,4)
1096 0 : b(3,2) = b(3,2) - b(4,2)*a(3,4)
1097 0 : b(3,2) = b(3,2)/a(3,3)
1098 0 : b(1,2) = b(1,2) - b(3,2)*a(1,3)
1099 0 : b(2,2) = b(2,2) - b(3,2)*a(2,3)
1100 0 : b(2,2) = b(2,2)/a(2,2)
1101 0 : b(1,2) = b(1,2) - b(2,2)*a(1,2)
1102 0 : b(1,2) = b(1,2)/a(1,1)
1103 :
1104 0 : b(5,3) = b(5,3)/a(5,5)
1105 0 : b(1,3) = b(1,3) - b(5,3)*a(1,5)
1106 0 : b(2,3) = b(2,3) - b(5,3)*a(2,5)
1107 0 : b(3,3) = b(3,3) - b(5,3)*a(3,5)
1108 0 : b(4,3) = b(4,3) - b(5,3)*a(4,5)
1109 0 : b(4,3) = b(4,3)/a(4,4)
1110 0 : b(1,3) = b(1,3) - b(4,3)*a(1,4)
1111 0 : b(2,3) = b(2,3) - b(4,3)*a(2,4)
1112 0 : b(3,3) = b(3,3) - b(4,3)*a(3,4)
1113 0 : b(3,3) = b(3,3)/a(3,3)
1114 0 : b(1,3) = b(1,3) - b(3,3)*a(1,3)
1115 0 : b(2,3) = b(2,3) - b(3,3)*a(2,3)
1116 0 : b(2,3) = b(2,3)/a(2,2)
1117 0 : b(1,3) = b(1,3) - b(2,3)*a(1,2)
1118 0 : b(1,3) = b(1,3)/a(1,1)
1119 :
1120 0 : b(5,4) = b(5,4)/a(5,5)
1121 0 : b(1,4) = b(1,4) - b(5,4)*a(1,5)
1122 0 : b(2,4) = b(2,4) - b(5,4)*a(2,5)
1123 0 : b(3,4) = b(3,4) - b(5,4)*a(3,5)
1124 0 : b(4,4) = b(4,4) - b(5,4)*a(4,5)
1125 0 : b(4,4) = b(4,4)/a(4,4)
1126 0 : b(1,4) = b(1,4) - b(4,4)*a(1,4)
1127 0 : b(2,4) = b(2,4) - b(4,4)*a(2,4)
1128 0 : b(3,4) = b(3,4) - b(4,4)*a(3,4)
1129 0 : b(3,4) = b(3,4)/a(3,3)
1130 0 : b(1,4) = b(1,4) - b(3,4)*a(1,3)
1131 0 : b(2,4) = b(2,4) - b(3,4)*a(2,3)
1132 0 : b(2,4) = b(2,4)/a(2,2)
1133 0 : b(1,4) = b(1,4) - b(2,4)*a(1,2)
1134 0 : b(1,4) = b(1,4)/a(1,1)
1135 :
1136 0 : b(5,5) = b(5,5)/a(5,5)
1137 0 : b(1,5) = b(1,5) - b(5,5)*a(1,5)
1138 0 : b(2,5) = b(2,5) - b(5,5)*a(2,5)
1139 0 : b(3,5) = b(3,5) - b(5,5)*a(3,5)
1140 0 : b(4,5) = b(4,5) - b(5,5)*a(4,5)
1141 0 : b(4,5) = b(4,5)/a(4,4)
1142 0 : b(1,5) = b(1,5) - b(4,5)*a(1,4)
1143 0 : b(2,5) = b(2,5) - b(4,5)*a(2,4)
1144 0 : b(3,5) = b(3,5) - b(4,5)*a(3,4)
1145 0 : b(3,5) = b(3,5)/a(3,3)
1146 0 : b(1,5) = b(1,5) - b(3,5)*a(1,3)
1147 0 : b(2,5) = b(2,5) - b(3,5)*a(2,3)
1148 0 : b(2,5) = b(2,5)/a(2,2)
1149 0 : b(1,5) = b(1,5) - b(2,5)*a(1,2)
1150 0 : b(1,5) = b(1,5)/a(1,1)
1151 :
1152 0 : end subroutine my_getrs_5_by_5
1153 :
1154 :
1155 0 : subroutine my_getrs_5_by_1( a, lda, ipiv, b, ldb, info )
1156 : integer :: info, lda, ldb ! , n=5, nrhs=1
1157 : integer, pointer :: ipiv(:)
1158 : real(fltp), pointer :: a(:,:), b(:,:) ! a( lda, * ), b( ldb, * )
1159 : real(fltp), parameter :: zero=0
1160 0 : info = 0
1161 0 : call my_laswp_5_by_1( b, ldb, ipiv )
1162 : !call dtrsm( 'left', 'lower', 'no transpose', 'unit', n, nrhs, one, a, lda, b, ldb )
1163 0 : b(2,1) = b(2,1) - b(1,1)*a(2,1)
1164 0 : b(3,1) = b(3,1) - b(1,1)*a(3,1)
1165 0 : b(4,1) = b(4,1) - b(1,1)*a(4,1)
1166 0 : b(5,1) = b(5,1) - b(1,1)*a(5,1)
1167 0 : b(3,1) = b(3,1) - b(2,1)*a(3,2)
1168 0 : b(4,1) = b(4,1) - b(2,1)*a(4,2)
1169 0 : b(5,1) = b(5,1) - b(2,1)*a(5,2)
1170 0 : b(4,1) = b(4,1) - b(3,1)*a(4,3)
1171 0 : b(5,1) = b(5,1) - b(3,1)*a(5,3)
1172 0 : b(5,1) = b(5,1) - b(4,1)*a(5,4)
1173 :
1174 : !call dtrsm( 'left', 'upper', 'no transpose', 'non-unit', n, nrhs, one, a, lda, b, ldb )
1175 0 : b(5,1) = b(5,1)/a(5,5)
1176 0 : b(1,1) = b(1,1) - b(5,1)*a(1,5)
1177 0 : b(2,1) = b(2,1) - b(5,1)*a(2,5)
1178 0 : b(3,1) = b(3,1) - b(5,1)*a(3,5)
1179 0 : b(4,1) = b(4,1) - b(5,1)*a(4,5)
1180 0 : b(4,1) = b(4,1)/a(4,4)
1181 0 : b(1,1) = b(1,1) - b(4,1)*a(1,4)
1182 0 : b(2,1) = b(2,1) - b(4,1)*a(2,4)
1183 0 : b(3,1) = b(3,1) - b(4,1)*a(3,4)
1184 0 : b(3,1) = b(3,1)/a(3,3)
1185 0 : b(1,1) = b(1,1) - b(3,1)*a(1,3)
1186 0 : b(2,1) = b(2,1) - b(3,1)*a(2,3)
1187 0 : b(2,1) = b(2,1)/a(2,2)
1188 0 : b(1,1) = b(1,1) - b(2,1)*a(1,2)
1189 0 : b(1,1) = b(1,1)/a(1,1)
1190 :
1191 0 : end subroutine my_getrs_5_by_1
1192 :
1193 :
1194 0 : subroutine my_getrs_4_by_4( a, lda, ipiv, b, ldb, info )
1195 : integer :: info, lda, ldb ! , n=4, nrhs=4
1196 : integer, pointer :: ipiv(:)
1197 : real(fltp), pointer :: a(:,:), b(:,:) ! a( lda, * ), b( ldb, * )
1198 : real(fltp), parameter :: zero=0
1199 0 : info = 0
1200 :
1201 0 : call my_laswp_4_by_4( b, ldb, ipiv )
1202 :
1203 : !call dtrsm( 'left', 'lower', 'no transpose', 'unit', n, nrhs, one, a, lda, b, ldb )
1204 0 : b(2,1) = b(2,1) - b(1,1)*a(2,1)
1205 0 : b(3,1) = b(3,1) - b(1,1)*a(3,1)
1206 0 : b(4,1) = b(4,1) - b(1,1)*a(4,1)
1207 0 : b(3,1) = b(3,1) - b(2,1)*a(3,2)
1208 0 : b(4,1) = b(4,1) - b(2,1)*a(4,2)
1209 0 : b(4,1) = b(4,1) - b(3,1)*a(4,3)
1210 :
1211 0 : b(2,2) = b(2,2) - b(1,2)*a(2,1)
1212 0 : b(3,2) = b(3,2) - b(1,2)*a(3,1)
1213 0 : b(4,2) = b(4,2) - b(1,2)*a(4,1)
1214 0 : b(3,2) = b(3,2) - b(2,2)*a(3,2)
1215 0 : b(4,2) = b(4,2) - b(2,2)*a(4,2)
1216 0 : b(4,2) = b(4,2) - b(3,2)*a(4,3)
1217 :
1218 0 : b(2,3) = b(2,3) - b(1,3)*a(2,1)
1219 0 : b(3,3) = b(3,3) - b(1,3)*a(3,1)
1220 0 : b(4,3) = b(4,3) - b(1,3)*a(4,1)
1221 0 : b(3,3) = b(3,3) - b(2,3)*a(3,2)
1222 0 : b(4,3) = b(4,3) - b(2,3)*a(4,2)
1223 0 : b(4,3) = b(4,3) - b(3,3)*a(4,3)
1224 :
1225 0 : b(2,4) = b(2,4) - b(1,4)*a(2,1)
1226 0 : b(3,4) = b(3,4) - b(1,4)*a(3,1)
1227 0 : b(4,4) = b(4,4) - b(1,4)*a(4,1)
1228 0 : b(3,4) = b(3,4) - b(2,4)*a(3,2)
1229 0 : b(4,4) = b(4,4) - b(2,4)*a(4,2)
1230 0 : b(4,4) = b(4,4) - b(3,4)*a(4,3)
1231 :
1232 : !call dtrsm( 'left', 'upper', 'no transpose', 'non-unit', n, nrhs, one, a, lda, b, ldb )
1233 0 : b(4,1) = b(4,1)/a(4,4)
1234 0 : b(1,1) = b(1,1) - b(4,1)*a(1,4)
1235 0 : b(2,1) = b(2,1) - b(4,1)*a(2,4)
1236 0 : b(3,1) = b(3,1) - b(4,1)*a(3,4)
1237 0 : b(3,1) = b(3,1)/a(3,3)
1238 0 : b(1,1) = b(1,1) - b(3,1)*a(1,3)
1239 0 : b(2,1) = b(2,1) - b(3,1)*a(2,3)
1240 0 : b(2,1) = b(2,1)/a(2,2)
1241 0 : b(1,1) = b(1,1) - b(2,1)*a(1,2)
1242 0 : b(1,1) = b(1,1)/a(1,1)
1243 :
1244 0 : b(4,2) = b(4,2)/a(4,4)
1245 0 : b(1,2) = b(1,2) - b(4,2)*a(1,4)
1246 0 : b(2,2) = b(2,2) - b(4,2)*a(2,4)
1247 0 : b(3,2) = b(3,2) - b(4,2)*a(3,4)
1248 0 : b(3,2) = b(3,2)/a(3,3)
1249 0 : b(1,2) = b(1,2) - b(3,2)*a(1,3)
1250 0 : b(2,2) = b(2,2) - b(3,2)*a(2,3)
1251 0 : b(2,2) = b(2,2)/a(2,2)
1252 0 : b(1,2) = b(1,2) - b(2,2)*a(1,2)
1253 0 : b(1,2) = b(1,2)/a(1,1)
1254 :
1255 0 : b(4,3) = b(4,3)/a(4,4)
1256 0 : b(1,3) = b(1,3) - b(4,3)*a(1,4)
1257 0 : b(2,3) = b(2,3) - b(4,3)*a(2,4)
1258 0 : b(3,3) = b(3,3) - b(4,3)*a(3,4)
1259 0 : b(3,3) = b(3,3)/a(3,3)
1260 0 : b(1,3) = b(1,3) - b(3,3)*a(1,3)
1261 0 : b(2,3) = b(2,3) - b(3,3)*a(2,3)
1262 0 : b(2,3) = b(2,3)/a(2,2)
1263 0 : b(1,3) = b(1,3) - b(2,3)*a(1,2)
1264 0 : b(1,3) = b(1,3)/a(1,1)
1265 :
1266 0 : b(4,4) = b(4,4)/a(4,4)
1267 0 : b(1,4) = b(1,4) - b(4,4)*a(1,4)
1268 0 : b(2,4) = b(2,4) - b(4,4)*a(2,4)
1269 0 : b(3,4) = b(3,4) - b(4,4)*a(3,4)
1270 0 : b(3,4) = b(3,4)/a(3,3)
1271 0 : b(1,4) = b(1,4) - b(3,4)*a(1,3)
1272 0 : b(2,4) = b(2,4) - b(3,4)*a(2,3)
1273 0 : b(2,4) = b(2,4)/a(2,2)
1274 0 : b(1,4) = b(1,4) - b(2,4)*a(1,2)
1275 0 : b(1,4) = b(1,4)/a(1,1)
1276 :
1277 0 : end subroutine my_getrs_4_by_4
1278 :
1279 :
1280 0 : subroutine my_getrs_4_by_1( a, lda, ipiv, b, ldb, info )
1281 : integer :: info, lda, ldb ! , n=4, nrhs=1
1282 : integer, pointer :: ipiv(:)
1283 : real(fltp), pointer :: a(:,:), b(:,:) ! a( lda, * ), b( ldb, * )
1284 : real(fltp), parameter :: zero=0
1285 :
1286 0 : info = 0
1287 0 : call my_laswp_4_by_1( b, ldb, ipiv )
1288 :
1289 : !call dtrsm( 'left', 'lower', 'no transpose', 'unit', n, nrhs, one, a, lda, b, ldb )
1290 0 : b(2,1) = b(2,1) - b(1,1)*a(2,1)
1291 0 : b(3,1) = b(3,1) - b(1,1)*a(3,1)
1292 0 : b(4,1) = b(4,1) - b(1,1)*a(4,1)
1293 0 : b(3,1) = b(3,1) - b(2,1)*a(3,2)
1294 0 : b(4,1) = b(4,1) - b(2,1)*a(4,2)
1295 0 : b(4,1) = b(4,1) - b(3,1)*a(4,3)
1296 :
1297 : !call dtrsm( 'left', 'upper', 'no transpose', 'non-unit', n, nrhs, one, a, lda, b, ldb )
1298 0 : b(4,1) = b(4,1)/a(4,4)
1299 0 : b(1,1) = b(1,1) - b(4,1)*a(1,4)
1300 0 : b(2,1) = b(2,1) - b(4,1)*a(2,4)
1301 0 : b(3,1) = b(3,1) - b(4,1)*a(3,4)
1302 0 : b(3,1) = b(3,1)/a(3,3)
1303 0 : b(1,1) = b(1,1) - b(3,1)*a(1,3)
1304 0 : b(2,1) = b(2,1) - b(3,1)*a(2,3)
1305 0 : b(2,1) = b(2,1)/a(2,2)
1306 0 : b(1,1) = b(1,1) - b(2,1)*a(1,2)
1307 0 : b(1,1) = b(1,1)/a(1,1)
1308 :
1309 0 : end subroutine my_getrs_4_by_1
1310 :
1311 :
1312 0 : subroutine my_getrs_dbg( n, nrhs, a, lda, ipiv, b, ldb, info )
1313 : integer :: info, lda, ldb, n, nrhs
1314 : integer, pointer :: ipiv(:)
1315 : real(fltp), pointer :: a(:,:), b(:,:) ! a( lda, * ), b( ldb, * )
1316 : real(fltp), parameter :: one=1, zero=0
1317 : integer :: i, j, k
1318 0 : info = 0
1319 0 : call my_laswp_dbg(nrhs, b, ldb, 1, n, ipiv, 1 )
1320 : !call dtrsm( 'left', 'lower', 'no transpose', 'unit', n, nrhs, one, a, lda, b, ldb )
1321 0 : do j = 1,nrhs
1322 0 : do k = 1,n
1323 0 : if (b(k,j)/=zero) then
1324 0 : do i = k + 1,n
1325 0 : b(i,j) = b(i,j) - b(k,j)*a(i,k)
1326 : end do
1327 : end if
1328 : end do
1329 : end do
1330 : !call dtrsm( 'left', 'upper', 'no transpose', 'non-unit', n, nrhs, one, a, lda, b, ldb )
1331 0 : do j = 1,nrhs
1332 0 : do k = n,1,-1
1333 0 : if (b(k,j)/=zero) then
1334 0 : b(k,j) = b(k,j)/a(k,k)
1335 0 : do i = 1,k - 1
1336 0 : b(i,j) = b(i,j) - b(k,j)*a(i,k)
1337 : end do
1338 : end if
1339 : end do
1340 : end do
1341 0 : end subroutine my_getrs_dbg
1342 :
1343 :
1344 0 : subroutine my_laswp_dbg( n, a, lda, k1, k2, ipiv, incx )
1345 : integer :: incx, k1, k2, lda, n
1346 : integer :: ipiv(:)
1347 : real(fltp) :: a(:,:) ! a( lda, * )
1348 : integer :: i, i1, i2, inc, ip, ix, ix0, j, k, n32
1349 0 : real(fltp) :: temp
1350 : include 'formats'
1351 0 : write(*,2) 'n', n
1352 0 : write(*,2) 'incx', incx
1353 0 : write(*,2) 'k1', k1
1354 0 : write(*,2) 'k2', k2
1355 0 : do j = 1, n
1356 0 : write(*,3) 'ipiv(j)', j, ipiv(j)
1357 : end do
1358 : ! interchange row i with row ipiv(i) for each of rows k1 through k2.
1359 0 : if( incx>0 ) then
1360 0 : ix0 = k1
1361 0 : i1 = k1
1362 0 : i2 = k2
1363 0 : inc = 1
1364 0 : else if( incx<0 ) then
1365 0 : ix0 = 1 + ( 1-k2 )*incx
1366 0 : i1 = k2
1367 0 : i2 = k1
1368 0 : inc = -1
1369 : else
1370 0 : return
1371 : end if
1372 0 : n32 = ( n / 32 )*32
1373 0 : if( n32/=0 ) then
1374 0 : do j = 1, n32, 32
1375 0 : ix = ix0
1376 0 : do i = i1, i2, inc
1377 0 : ip = ipiv( ix )
1378 0 : if( ip/=i ) then
1379 0 : do k = j, j + 31
1380 0 : temp = a( i, k )
1381 0 : a( i, k ) = a( ip, k )
1382 0 : a( ip, k ) = temp
1383 : end do
1384 : end if
1385 0 : ix = ix + incx
1386 : end do
1387 : end do
1388 : end if
1389 0 : if( n32/=n ) then
1390 0 : n32 = n32 + 1
1391 0 : ix = ix0
1392 0 : do i = i1, i2, inc
1393 0 : ip = ipiv( ix )
1394 :
1395 :
1396 0 : if (ip == 0) then
1397 :
1398 0 : stop 'my_lapack95 ip == 0'
1399 :
1400 : end if
1401 :
1402 0 : if( ip/=i ) then
1403 0 : do k = n32, n
1404 0 : temp = a( i, k )
1405 0 : a( i, k ) = a( ip, k )
1406 0 : a( ip, k ) = temp
1407 : end do
1408 : end if
1409 0 : ix = ix + incx
1410 : end do
1411 : end if
1412 : end subroutine my_laswp_dbg
1413 :
1414 : end module my_lapack95
|
