Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010 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 :
21 : module mtx_support
22 :
23 : use const_def, only: dp, qp
24 : use utils_lib, only: mesa_error
25 :
26 : implicit none
27 :
28 : integer, parameter :: num_chunks = 4
29 :
30 : contains
31 :
32 :
33 :
34 1 : subroutine do_dense_to_band(n,ndim,a,ml,mu,ab,ldab,ierr)
35 : integer, intent(in) :: n,ndim,ml,mu,ldab
36 : real(dp), intent(in) :: a(:,:) ! (ndim,n)
37 : real(dp), intent(inout) :: ab(:,:) ! (ldab,n)
38 : integer, intent(out) :: ierr
39 : integer :: i, j
40 1 : if (ml+mu+1 > n) then
41 0 : ierr = -1
42 0 : return
43 : end if
44 1 : ierr = 0
45 37 : ab = 0
46 7 : do j=1,n
47 27 : do i=max(1,j-mu),min(n,j+ml)
48 26 : ab(ldab-ml+i-j,j) = a(i,j)
49 : end do
50 : end do
51 : end subroutine do_dense_to_band
52 :
53 :
54 1 : subroutine do_band_to_dense(n,ml,mu,ab,ldab,ndim,a,ierr)
55 : integer, intent(in) :: n,ndim,ml,mu,ldab
56 : real(dp), intent(in) :: ab(:,:) ! (ldab,n)
57 : real(dp), intent(inout) :: a(:,:) ! (ndim,n)
58 : integer, intent(out) :: ierr
59 : integer :: i, j
60 1 : if (ml+mu+1 > n) then
61 0 : ierr = -1
62 0 : return
63 : end if
64 1 : ierr = 0
65 43 : a = 0
66 7 : do j=1,n
67 27 : do i=max(1,j-mu),min(n,j+ml)
68 26 : a(i,j) = ab(ldab-ml+i-j,j)
69 : end do
70 : end do
71 : end subroutine do_band_to_dense
72 :
73 :
74 1 : subroutine do_band_to_column_sparse(n,ml,mu,ab,ldab,nzmax,nz,colptr,rowind,values,diags,ierr)
75 : integer, intent(in) :: n,ml,mu,nzmax,ldab
76 : real(dp), intent(in) :: ab(ldab,n)
77 : integer, intent(inout) :: colptr(n+1),rowind(nzmax)
78 : real(dp), intent(inout) :: values(nzmax)
79 : ! real(dp), intent(in) :: ab(:,:) ! (ldab,n)
80 : ! integer, intent(inout) :: colptr(:) ! (n+1)
81 : ! integer, intent(inout) :: rowind(:) ! (nzmax)
82 : ! real(dp), intent(inout) :: values(:) ! (nzmax)
83 : logical, intent(in) :: diags
84 : integer, intent(out) :: nz,ierr
85 : integer :: i, j
86 1 : if (ml+mu+1 > n) then
87 0 : ierr = -1
88 0 : return
89 : end if
90 1 : ierr = 0
91 1 : nz = 0
92 7 : do j=1,n
93 6 : colptr(j) = nz+1 ! index in values of first entry in this column
94 27 : do i=max(1,j-mu),min(n,j+ml)
95 20 : if (ab(ldab-ml+i-j,j) == 0) then
96 7 : if (i /= j) cycle ! not a diagonal
97 0 : if (.not. diags) cycle
98 : ! else include diagonals even if are == 0
99 : end if
100 13 : nz = nz+1
101 13 : if (nz > nzmax) then
102 0 : ierr = j
103 0 : return
104 : end if
105 13 : values(nz) = ab(ldab-ml+i-j,j)
106 26 : rowind(nz) = i
107 : end do
108 : end do
109 1 : colptr(n+1) = nz+1
110 : end subroutine do_band_to_column_sparse
111 :
112 :
113 1 : subroutine do_column_sparse_to_band(n,ml,mu,ab,ldab,nz,colptr,rowind,values,ierr)
114 : integer, intent(in) :: n,ml,mu,nz,ldab
115 :
116 : real(dp), intent(inout) :: ab(ldab,n)
117 : integer, intent(in) :: colptr(n+1),rowind(nz)
118 : real(dp), intent(in) :: values(nz)
119 : ! real(dp), intent(inout) :: ab(:,:) ! (ldab,n)
120 : ! integer, intent(inout) :: colptr(:) ! (n+1)
121 : ! integer, intent(inout) :: rowind(:) ! (nzmax)
122 : ! real(dp), intent(in) :: values(:) ! (nz)
123 : integer, intent(out) :: ierr
124 : integer :: i,j,k
125 1 : ierr = 0
126 37 : ab = 0
127 7 : do j=1,n
128 20 : do k=colptr(j),colptr(j+1)-1
129 13 : i = rowind(k)
130 13 : if (i > j+ml .or. i < j-mu) then
131 0 : ierr = j
132 0 : return
133 : end if
134 19 : ab(ldab-ml+i-j,j) = values(k)
135 : end do
136 : end do
137 : end subroutine do_column_sparse_to_band
138 :
139 :
140 1 : subroutine do_band_to_row_sparse(n,ml,mu,ab,ldab,nzmax,nz,rowptr,colind,values,diags,ierr)
141 : integer, intent(in) :: n,ml,mu,nzmax,ldab
142 :
143 : real(dp), intent(in) :: ab(ldab,n)
144 : integer, intent(inout) :: rowptr(n+1),colind(nzmax)
145 : real(dp), intent(inout) :: values(nzmax)
146 : ! real(dp), intent(in) :: ab(:,:) ! (ldab,n)
147 : ! integer, intent(out) :: rowptr(:) ! (n+1)
148 : ! integer, intent(out) :: colind(:) ! (nzmax)
149 : ! real(dp), intent(inout) :: values(:) ! (nzmax)
150 : logical, intent(in) :: diags
151 : integer, intent(out) :: ierr, nz
152 : integer :: idiag, op_err, j1, j2, k, i, nz1
153 : integer, dimension(num_chunks) :: nz_per_chunk, nz_start, nz_max, i_lo, i_hi
154 :
155 : logical, parameter :: dbg = .false.
156 :
157 : include 'formats'
158 :
159 : if (dbg) write(*,*) 'enter do_band_to_row_sparse'
160 :
161 1 : if (ml+mu+1 > n) then
162 0 : ierr = -1
163 : if (dbg) then
164 : write(*,*) 'do_band_to_row_sparse'
165 : write(*,*) 'ml+mu+1', ml+mu+1
166 : write(*,*) 'n', n
167 : end if
168 0 : return
169 : end if
170 :
171 1 : ierr = 0
172 1 : nz = 0
173 1 : idiag = ldab - ml
174 :
175 1 : nz_start(1) = 1
176 1 : i_lo(1) = 1
177 4 : do k = 2, num_chunks
178 3 : nz_start(k) = ((k-1)*nzmax)/num_chunks
179 3 : nz_max(k-1) = nz_start(k) - 1
180 3 : i_lo(k) = ((k-1)*n)/num_chunks
181 4 : i_hi(k-1) = i_lo(k) - 1
182 : end do
183 1 : nz_max(num_chunks) = nzmax
184 1 : i_hi(num_chunks) = n
185 :
186 : if (dbg) write(*,*) 'do_band_to_row_sparse - do chunks'
187 1 : op_err = 0
188 1 : !$OMP PARALLEL DO PRIVATE(k,op_err)
189 : do k = 1, num_chunks
190 : call do_chunk_band_to_row_sparse( &
191 : k, ldab, n, nzmax, nz_per_chunk, nz_start, nz_max, ml, mu, idiag, diags, &
192 : i_lo, i_hi, ab, rowptr, colind, values, op_err)
193 : if (op_err /= 0) ierr = op_err
194 : end do
195 : !$OMP END PARALLEL DO
196 : if (dbg) write(*,*) 'do_band_to_row_sparse - done chunks'
197 :
198 1 : if (ierr /= 0) then
199 :
200 0 : write(*,*) 'do_band_to_row_sparse: failed to fit in chunks'
201 0 : write(*,*) 'please increase the max fill factor for your sparse matrix'
202 0 : call mesa_error(__FILE__,__LINE__)
203 :
204 : else
205 :
206 : if (dbg) then
207 : do k=1,num_chunks
208 : write(*,2) 'k', k
209 : write(*,2) 'nz_per_chunk(k)', nz_per_chunk(k)
210 : write(*,2) 'nz_start(k)', nz_start(k)
211 : write(*,2) 'nz_max(k)', nz_max(k)
212 : write(*,2) 'i_lo(k)', i_lo(k)
213 : write(*,2) 'i_hi(k)', i_hi(k)
214 : write(*,'(A)')
215 : end do
216 : end if
217 :
218 : ! reposition the chunk results
219 : if (dbg) write(*,*) 'reposition the chunk results'
220 1 : i = nz_per_chunk(1)
221 4 : do k = 2, num_chunks
222 3 : nz1 = nz_per_chunk(k)
223 : if (dbg) write(*,2) 'i', i
224 : if (dbg) write(*,2) 'nz1', nz1
225 : if (dbg) write(*,2) 'k', k
226 : if (dbg) write(*,2) 'nz_start(k)', nz_start(k)
227 3 : j2 = nz_start(k)
228 16 : do j1 = i+1, i+nz1 ! avoid ifort segfault
229 13 : values(j1) = values(j2)
230 13 : colind(j1) = colind(j2)
231 16 : j2 = j2+1
232 : end do
233 : if (dbg) write(*,2) 'i_lo(k)', i_lo(k)
234 : if (dbg) write(*,2) 'i_hi(k)', i_hi(k)
235 : if (dbg) write(*,2) 'nz_start(k)-i-1', nz_start(k)-i-1
236 3 : j2 = (nz_start(k)-i-1)
237 9 : do j1 = i_lo(k), i_hi(k)
238 9 : rowptr(j1) = rowptr(j1) - j2
239 : end do
240 4 : i = i+nz1
241 : end do
242 1 : nz = i
243 : end if
244 :
245 :
246 1 : rowptr(n+1) = nz+1
247 : !write(*,*) 'done do_band_to_row_sparse - fill fraction', dble(nz)/dble(nzmax)
248 :
249 : end subroutine do_band_to_row_sparse
250 :
251 :
252 4 : subroutine do_chunk_band_to_row_sparse( &
253 : num, ldab, n, nzmax, nz_per_chunk, nz_start, nz_max, ml, mu, idiag, diags, &
254 4 : i_lo, i_hi, ab, rowptr, colind, values, ierr)
255 : integer, intent(in) :: num, ldab, n, nzmax, ml, mu, idiag
256 : logical, intent(in) :: diags
257 :
258 : real(dp), intent(in) :: ab(ldab,n)
259 : integer, intent(inout) :: rowptr(n+1), colind(nzmax)
260 : real(dp), intent(inout) :: values(nzmax)
261 : ! real(dp), intent(in) :: ab(:,:) ! (ldab,n)
262 : ! integer, intent(inout) :: rowptr(:) ! (n+1)
263 : ! integer, intent(inout) :: colind(:) ! (nzmax)
264 : ! real(dp), intent(inout) :: values(:) ! (nzmax)
265 : integer, dimension(num_chunks) :: nz_per_chunk, nz_start, nz_max, i_lo, i_hi
266 : integer, intent(out) :: ierr
267 :
268 : integer :: i, j, nz
269 4 : real(dp) :: val
270 :
271 : logical, parameter :: dbg = .false.
272 : include 'formats'
273 :
274 4 : ierr = 0
275 4 : nz = nz_start(num) - 1
276 10 : do i = i_lo(num), i_hi(num)
277 6 : rowptr(i) = nz+1 ! index in values of first entry in this row
278 : !write(*,*) 'set rowptr(i)', i, rowptr(i)
279 30 : do j = max(1,i-ml), min(n,i+mu)
280 20 : val = ab(idiag+i-j,j)
281 20 : if (val == 0) then
282 7 : if (i /= j) cycle ! not a diagonal, so skip if 0
283 0 : if (.not. diags) cycle
284 : ! if (diags) then include diagonals even if are == 0
285 : end if
286 13 : nz = nz+1
287 13 : if (nz > nz_max(num)) then
288 0 : ierr = i
289 : if (dbg) then
290 : write(*,*) 'nz > nz_max(num)', nz, nz_max(num), num
291 : end if
292 : return
293 : end if
294 13 : values(nz) = val
295 26 : colind(nz) = j
296 : end do
297 : end do
298 4 : nz_per_chunk(num) = nz - nz_start(num) + 1 ! number of non-zero values for this chunk
299 : end subroutine do_chunk_band_to_row_sparse
300 :
301 :
302 :
303 0 : subroutine do_row_sparse_to_band(n,ml,mu,ab,ldab,nz,rowptr,colind,values,ierr)
304 : integer, intent(in) :: n,ml,mu,nz,ldab
305 : real(dp), intent(inout) :: ab(ldab,n)
306 : integer, intent(in) :: rowptr(n+1),colind(nz)
307 : real(dp), intent(in) :: values(nz)
308 : ! real(dp), intent(inout) :: ab(:,:) ! (ldab,n)
309 : ! integer, intent(inout) :: rowptr(:) ! (n+1)
310 : ! integer, intent(inout) :: colind(:) ! (nz)
311 : ! real(dp), intent(in) :: values(:) ! (nz)
312 : integer, intent(out) :: ierr
313 : integer :: i,j,k
314 0 : ierr = 0
315 0 : ab = 0
316 0 : do i=1,n
317 0 : do k=rowptr(i),rowptr(i+1)-1
318 0 : j = colind(k)
319 0 : if (i > j+ml .or. i < j-mu) then
320 0 : ierr = j
321 0 : return
322 : end if
323 0 : ab(ldab-ml+i-j,j) = values(k)
324 : end do
325 : end do
326 : end subroutine do_row_sparse_to_band
327 :
328 :
329 : ! sparse conversion based on similar routines from sparskit_src/formats.f
330 :
331 1 : subroutine do_dense_to_row_sparse(n,ndim,a,nzmax,nz,rowptr,colind,values,diags,ierr)
332 : integer, intent(in) :: n,ndim,nzmax
333 : !real(dp), intent(in) :: a(ndim,n)
334 : !integer, intent(inout) :: rowptr(n+1),colind(nzmax)
335 : !real(dp), intent(inout) :: values(nzmax)
336 : real(dp), intent(in) :: a(:,:) ! (ndim,n)
337 : integer, intent(inout) :: rowptr(:) ! (n+1)
338 : integer, intent(inout) :: colind(:) ! (nzmax)
339 : real(dp), intent(inout) :: values(:) ! (nzmax)
340 : logical, intent(in) :: diags
341 : integer, intent(out) :: nz,ierr
342 : integer :: i,j
343 1 : ierr = 0
344 1 : nz = 0
345 7 : do i=1,n
346 6 : rowptr(i) = nz+1 ! index in values of first entry in this row
347 43 : do j=1,n
348 36 : if (a(i,j) == 0) then
349 21 : if (i /= j) cycle ! not a diagonal
350 0 : if (.not. diags) cycle
351 : ! else include diagonals even if are == 0
352 : end if
353 15 : nz = nz+1
354 15 : if (nz > nzmax) then
355 0 : ierr = i
356 0 : return
357 : end if
358 15 : values(nz) = a(i,j)
359 42 : colind(nz) = j
360 : end do
361 : end do
362 1 : rowptr(n+1) = nz+1
363 : end subroutine do_dense_to_row_sparse
364 :
365 :
366 0 : subroutine do_dense_to_row_sparse_0_based( &
367 0 : n,ndim,a,nzmax,nz,rowptr,colind,values,diags,ierr)
368 : integer, intent(in) :: n,ndim,nzmax
369 : !real(dp), intent(in) :: a(ndim,n)
370 : !integer, intent(inout) :: rowptr(n+1),colind(nzmax)
371 : !real(dp), intent(inout) :: values(nzmax)
372 : real(dp), intent(in) :: a(:,:) ! (ndim,n)
373 : integer, intent(inout) :: rowptr(:) ! (n+1)
374 : integer, intent(inout) :: colind(:) ! (nzmax)
375 : real(dp), intent(inout) :: values(:) ! (nzmax)
376 : logical, intent(in) :: diags
377 : integer, intent(out) :: nz,ierr
378 : integer :: i,j
379 0 : ierr = 0
380 0 : nz = 0
381 0 : do i=1,n
382 0 : rowptr(i) = nz ! index in values of first entry in this row
383 0 : do j=1,n
384 0 : if (a(i,j) == 0) then
385 0 : if (i /= j) cycle ! not a diagonal
386 0 : if (.not. diags) cycle
387 : ! else include diagonals even if are == 0
388 : end if
389 0 : nz = nz+1
390 0 : if (nz > nzmax) then
391 0 : ierr = i
392 0 : return
393 : end if
394 0 : values(nz) = a(i,j)
395 0 : colind(nz) = j-1
396 : end do
397 : end do
398 0 : rowptr(n+1) = nz
399 : end subroutine do_dense_to_row_sparse_0_based
400 :
401 :
402 0 : subroutine do_row_sparse_to_dense(n,ndim,a,nz,rowptr,colind,values,ierr)
403 : integer, intent(in) :: n,ndim,nz
404 : real(dp), intent(inout) :: a(ndim,n)
405 : integer, intent(in) :: rowptr(n+1),colind(nz)
406 : real(dp), intent(in) :: values(nz)
407 : ! real(dp), intent(inout) :: a(:,:) ! (ndim,n)
408 : ! integer, intent(inout) :: rowptr(:) ! (n+1)
409 : ! integer, intent(inout) :: colind(:) ! (nz)
410 : ! real(dp), intent(inout) :: values(:) ! (nz)
411 : integer, intent(out) :: ierr
412 : integer :: i,j,k
413 0 : ierr = 0
414 0 : a = 0
415 0 : do i=1,n
416 0 : do k=rowptr(i),rowptr(i+1)-1
417 0 : j = colind(k)
418 0 : if (j > n) then
419 0 : ierr = i
420 0 : return
421 : end if
422 0 : a(i,j) = values(k)
423 : end do
424 : end do
425 : end subroutine do_row_sparse_to_dense
426 :
427 :
428 0 : subroutine do_row_sparse_0_based_to_dense(n,ndim,a,nz,rowptr,colind,values,ierr)
429 : integer, intent(in) :: n,ndim,nz
430 : real(dp), intent(inout) :: a(ndim,n)
431 : integer, intent(in) :: rowptr(0:n),colind(0:nz-1)
432 : real(dp), intent(in) :: values(nz)
433 : integer, intent(out) :: ierr
434 : integer :: i,j,k
435 0 : ierr = 0
436 0 : a = 0
437 0 : do i=1,n
438 0 : do k=rowptr(i),rowptr(i+1)-1
439 0 : j = colind(k)
440 0 : if (j > n) then
441 0 : ierr = i
442 0 : return
443 : end if
444 0 : a(i,j) = values(k)
445 : end do
446 : end do
447 : end subroutine do_row_sparse_0_based_to_dense
448 :
449 :
450 0 : subroutine do_dense_to_column_sparse(n,ndim,a,nzmax,nz,colptr,rowind,values,diags,ierr)
451 : integer, intent(in) :: n,ndim,nzmax
452 : !real(dp), intent(in) :: a(ndim,n)
453 : !integer, intent(inout) :: colptr(n+1),rowind(nzmax)
454 : !real(dp), intent(inout) :: values(nzmax)
455 : real(dp), intent(in) :: a(:,:) ! (ndim,n)
456 : integer, intent(inout) :: colptr(:) ! (n+1)
457 : integer, intent(inout) :: rowind(:) ! (nzmax)
458 : real(dp), intent(inout) :: values(:) ! (nzmax)
459 : logical, intent(in) :: diags
460 : integer, intent(out) :: nz,ierr
461 : integer :: i,j
462 0 : ierr = 0
463 0 : nz = 0
464 0 : do j=1,n
465 0 : colptr(j) = nz+1 ! index in values of first entry in this column
466 0 : do i=1,n
467 0 : if (a(i,j) == 0) then
468 0 : if (i /= j) cycle ! not a diagonal
469 0 : if (.not. diags) cycle
470 : ! else include diagonals even if are == 0
471 : end if
472 0 : nz = nz+1
473 0 : if (nz > nzmax) then
474 0 : ierr = j
475 0 : return
476 : end if
477 0 : values(nz) = a(i,j)
478 0 : rowind(nz) = i
479 : end do
480 : end do
481 0 : colptr(n+1) = nz+1
482 : end subroutine do_dense_to_column_sparse
483 :
484 :
485 0 : subroutine do_dense_to_col_sparse_0_based( &
486 0 : n,ndim,a,nzmax,nz,colptr,rowind,values,diags,ierr)
487 : integer, intent(in) :: n,ndim,nzmax
488 : !real(dp), intent(in) :: a(ndim,n)
489 : !integer, intent(inout) :: colptr(n+1),rowind(nzmax)
490 : !real(dp), intent(inout) :: values(nzmax)
491 : real(dp), intent(in) :: a(:,:) ! (ndim,n)
492 : integer, intent(inout) :: colptr(:) ! (n+1)
493 : integer, intent(inout) :: rowind(:) ! (nzmax)
494 : real(dp), intent(inout) :: values(:) ! (nzmax)
495 : logical, intent(in) :: diags
496 : integer, intent(out) :: nz,ierr
497 : integer :: i,j
498 0 : ierr = 0
499 0 : nz = 0
500 0 : do j=1,n
501 0 : colptr(j) = nz ! index in values of first entry in this column
502 0 : do i=1,n
503 0 : if (a(i,j) == 0) then
504 0 : if (i /= j) cycle ! not a diagonal
505 0 : if (.not. diags) cycle
506 : ! else include diagonals even if are == 0
507 : end if
508 0 : nz = nz+1
509 0 : if (nz > nzmax) then
510 0 : ierr = j
511 0 : return
512 : end if
513 0 : values(nz) = a(i,j)
514 0 : rowind(nz) = i-1
515 : end do
516 : end do
517 0 : colptr(n+1) = nz
518 : end subroutine do_dense_to_col_sparse_0_based
519 :
520 :
521 0 : subroutine do_dense_to_col_sparse_0_based_qp( &
522 0 : n,ndim,a,nzmax,nz,colptr,rowind,values,diags,ierr)
523 : integer, intent(in) :: n,ndim,nzmax
524 : !real(dp), intent(in) :: a(ndim,n)
525 : !integer, intent(inout) :: colptr(n+1),rowind(nzmax)
526 : !real(dp), intent(inout) :: values(nzmax)
527 : real(dp), intent(in) :: a(:,:) ! (ndim,n)
528 : integer, intent(inout) :: colptr(:) ! (n+1)
529 : integer, intent(inout) :: rowind(:) ! (nzmax)
530 : real(qp), intent(out) :: values(:) ! (nzmax)
531 : logical, intent(in) :: diags
532 : integer, intent(out) :: nz,ierr
533 : integer :: i,j
534 0 : ierr = 0
535 0 : nz = 0
536 0 : do j=1,n
537 0 : colptr(j) = nz ! index in values of first entry in this column
538 0 : do i=1,n
539 0 : if (a(i,j) == 0) then
540 0 : if (i /= j) cycle ! not a diagonal
541 0 : if (.not. diags) cycle
542 : ! else include diagonals even if are == 0
543 : end if
544 0 : nz = nz+1
545 0 : if (nz > nzmax) then
546 0 : ierr = j
547 0 : return
548 : end if
549 0 : values(nz) = a(i,j)
550 0 : rowind(nz) = i-1
551 : end do
552 : end do
553 0 : colptr(n+1) = nz
554 : end subroutine do_dense_to_col_sparse_0_based_qp
555 :
556 :
557 2935 : subroutine do_column_sparse_to_dense(n,ndim,a,nz,colptr,rowind,values,ierr)
558 : integer, intent(in) :: n,ndim,nz
559 : real(dp), intent(inout) :: a(ndim,n)
560 : integer, intent(in) :: colptr(n+1),rowind(nz)
561 : real(dp), intent(in) :: values(nz)
562 : ! real(dp), intent(inout) :: a(:,:) ! (ndim,n)
563 : ! integer, intent(in) :: colptr(:) ! (n+1)
564 : ! integer, intent(in) :: rowind(:) ! (nz)
565 : ! real(dp), intent(in) :: values(:) ! (nz)
566 : integer, intent(out) :: ierr
567 : integer :: i,j,k
568 2935 : ierr = 0
569 1910685 : a = 0
570 76310 : do j=1,n
571 200137 : do k=colptr(j),colptr(j+1)-1
572 123827 : i = rowind(k)
573 123827 : if (i > n) then
574 0 : ierr = j
575 0 : return
576 : end if
577 197202 : a(i,j) = values(k)
578 : end do
579 : end do
580 : end subroutine do_column_sparse_to_dense
581 :
582 :
583 0 : subroutine do_col_sparse_0_based_to_dense(n,ndim,a,nz,colptr,rowind,values,ierr)
584 : integer, intent(in) :: n,ndim,nz
585 : real(dp), intent(inout) :: a(ndim,n)
586 : integer, intent(in) :: colptr(n+1),rowind(nz)
587 : real(dp), intent(in) :: values(nz)
588 : ! real(dp), intent(inout) :: a(:,:) ! (ndim,n)
589 : ! integer, intent(in) :: colptr(:) ! (n+1)
590 : ! integer, intent(in) :: rowind(:) ! (nz)
591 : ! real(dp), intent(in) :: values(:) ! (nz)
592 : integer, intent(out) :: ierr
593 : integer :: i,j,k
594 0 : ierr = 0
595 0 : a = 0
596 0 : do j=1,n
597 0 : do k=colptr(j)+1,colptr(j+1)
598 0 : i = rowind(k)+1
599 0 : if (i > n) then
600 0 : ierr = j
601 0 : return
602 : end if
603 0 : a(i,j) = values(k)
604 : end do
605 : end do
606 : end subroutine do_col_sparse_0_based_to_dense
607 :
608 :
609 1 : subroutine do_block_dble_mv(nvar, nz, lblk, dblk, ublk, b, prod)
610 : use my_lapack95, only: my_gemv_p1
611 : integer, intent(in) :: nvar, nz
612 : real(dp), pointer, dimension(:,:,:), intent(in) :: lblk, dblk, ublk ! (nvar,nvar,nz)
613 : real(dp), pointer, dimension(:,:), intent(in) :: b ! (nvar,nz)
614 : real(dp), pointer, dimension(:,:), intent(inout) :: prod ! (nvar,nz)
615 : integer :: k
616 21 : do k = 1, nz
617 520 : prod(1:nvar,k) = 0
618 20 : call my_gemv_p1(nvar,nvar,dblk(1:nvar,1:nvar,k),nvar,b(1:nvar,k),prod(1:nvar,k))
619 20 : if (k > 1) then
620 19 : call my_gemv_p1(nvar,nvar,lblk(1:nvar,1:nvar,k),nvar,b(1:nvar,k-1),prod(1:nvar,k))
621 : end if
622 21 : if (k < nz) then
623 19 : call my_gemv_p1(nvar,nvar,ublk(1:nvar,1:nvar,k),nvar,b(1:nvar,k+1),prod(1:nvar,k))
624 : end if
625 : end do
626 1 : end subroutine do_block_dble_mv
627 :
628 :
629 0 : subroutine do_LU_factored_block_dble_mv(lblk, dblk, ublk, b, ipiv, prod)
630 : real(dp), pointer, dimension(:,:,:), intent(in) :: lblk, dblk, ublk ! (nvar,nvar,nz)
631 : real(dp), pointer, dimension(:,:), intent(in) :: b ! (nvar,nz)
632 : integer, intent(in) :: ipiv(:,:) ! (nvar,nz)
633 : real(dp), pointer, dimension(:,:), intent(inout) :: prod ! (nvar,nz)
634 : integer :: nvar, nz, k, incx, incy
635 0 : nvar = size(b,dim=1)
636 0 : nz = size(b,dim=2)
637 0 : incx = 1
638 0 : incy = 1
639 0 : !$OMP PARALLEL DO PRIVATE(k)
640 : do k = 1, nz
641 : call do_LU_factored_square_mv(nvar,dblk(:,:,k),b(:,k),ipiv(:,k),prod(:,k))
642 : if (k > 1) then
643 : call dgemv('N',nvar,nvar,1d0,lblk(:,:,k),nvar,b(:,k-1),incx,1d0,prod(:,k),incy)
644 : end if
645 : if (k < nz) then
646 : call dgemv('N',nvar,nvar,1d0,ublk(:,:,k),nvar,b(:,k+1),incx,1d0,prod(:,k),incy)
647 : end if
648 : end do
649 : !$OMP END PARALLEL DO
650 0 : end subroutine do_LU_factored_block_dble_mv
651 :
652 :
653 0 : subroutine do_LU_factored_square_mv(m,a,b,ipiv,x) ! set x = A*b
654 : ! A factored in LU manner = P*L*U.
655 : integer, intent(in) :: m
656 : real(dp), intent(in) :: a(:,:) ! (lda,m), lda >= m
657 : real(dp), intent(in) :: b(:) ! (m)
658 : integer, intent(in) :: ipiv(:) ! (m)
659 : real(dp), intent(inout) :: x(:) ! (m)
660 : integer :: i, j
661 0 : real(dp), dimension(m) :: y
662 : include 'formats'
663 : ! y = U*b
664 0 : do i=1,m
665 0 : y(i) = 0
666 0 : do j=i,m
667 0 : y(i) = y(i) + a(i,j)*b(j)
668 : end do
669 : end do
670 : ! x = L*y
671 0 : do i=1,m
672 0 : x(i) = y(i)
673 0 : do j=1,i-1
674 0 : x(i) = x(i) + a(i,j)*y(j)
675 : end do
676 : end do
677 : ! x = P*x
678 0 : call dlaswp(1, x, m, 1, m, ipiv, -1)
679 0 : end subroutine do_LU_factored_square_mv
680 :
681 :
682 0 : subroutine do_LU_factored_square_mm(m,A,B,ipiv,C) ! set C = A*B
683 : ! A factored in LU manner = P*L*U.
684 : integer, intent(in) :: m
685 : real(dp), intent(in) :: A(:,:) ! (m,m)
686 : real(dp), intent(in) :: B(:,:) ! (m,m)
687 : integer, intent(in) :: ipiv(:) ! (m)
688 : real(dp), intent(inout) :: C(:,:) ! (m,m)
689 : integer :: i, j, k
690 0 : real(dp), dimension(m,m) :: Y
691 : include 'formats'
692 : ! Y = U*B
693 0 : do i=1,m
694 0 : do j=1,m
695 0 : Y(i,j) = 0
696 0 : do k=i,m
697 0 : Y(i,j) = Y(i,j) + A(i,k)*B(k,j)
698 : end do
699 : end do
700 : end do
701 : ! C = L*Y
702 0 : do i=1,m
703 0 : do j=1,m
704 0 : C(i,j) = Y(i,j)
705 0 : do k=1,i-1
706 0 : C(i,j) = C(i,j) + A(i,k)*Y(k,j)
707 : end do
708 : end do
709 : end do
710 : ! C = P*C
711 0 : call dlaswp(m, C, m, 1, m, ipiv, -1)
712 0 : end subroutine do_LU_factored_square_mm
713 :
714 :
715 0 : subroutine do_multiply_xa(n, A1, x, b)
716 : ! calculates b = x*A
717 : integer, intent(in) :: n
718 : real(dp), pointer, intent(in) :: A1(:) ! =(n, n)
719 : real(dp), pointer, intent(in) :: x(:) ! (n)
720 : real(dp), pointer, intent(inout) :: b(:) ! (n)
721 : integer :: j
722 0 : real(dp), pointer :: A(:,:) ! (n, n)
723 0 : A(1:n,1:n) => A1(1:n*n)
724 0 : do j = 1, n
725 0 : b(j) = dot_product(x(1:n),A(1:n,j))
726 : end do
727 0 : end subroutine do_multiply_xa
728 :
729 :
730 0 : subroutine do_multiply_xa_plus_c(n, A1, x, c, b)
731 : ! calculates b = x*A + c
732 : integer, intent(in) :: n
733 : real(dp), pointer, intent(in) :: A1(:) ! =(n,n)
734 : real(dp), pointer, intent(in) :: x(:) ! (n)
735 : real(dp), pointer, intent(in) :: c(:) ! (n)
736 : real(dp), pointer, intent(inout) :: b(:) ! (n)
737 : integer :: j
738 0 : real(dp), pointer :: A(:,:) ! (n,n)
739 0 : A(1:n,1:n) => A1(1:n*n)
740 0 : do j = 1, n
741 0 : b(j) = dot_product(x(1:n),A(1:n,j)) + c(j)
742 : end do
743 0 : end subroutine do_multiply_xa_plus_c
744 :
745 :
746 0 : subroutine do_block_multiply_xa(nvar, nz, lblk1, dblk1, ublk1, x1, b1)
747 : ! calculates b = x*A
748 : integer, intent(in) :: nvar, nz
749 : real(dp), dimension(:), intent(in), pointer :: lblk1, dblk1, ublk1 ! =(nvar,nvar,nz)
750 : real(dp), intent(in), pointer :: x1(:) ! =(nvar,nz)
751 : real(dp), intent(inout), pointer :: b1(:) ! =(nvar,nz)
752 : integer :: k, shift, shift2, nvar2
753 0 : real(dp), pointer, dimension(:) :: p1, p2, p3, p4
754 0 : nvar2 = nvar*nvar
755 0 : !$OMP PARALLEL DO PRIVATE(k,shift,shift2,p1,p2,p3,p4)
756 : do k = 1, nz
757 : shift = nvar*(k-1)
758 : shift2 = nvar2*(k-1)
759 : p1(1:nvar2) => dblk1(shift2+1:shift2+nvar2)
760 : p2(1:nvar) => x1(shift+1:shift+nvar)
761 : p3(1:nvar) => b1(shift+1:shift+nvar)
762 : call do_multiply_xa(nvar,p1,p2,p3)
763 : if (k > 1) then
764 : p1(1:nvar2) => ublk1(shift2+1:shift2+nvar2)
765 : p2(1:nvar) => x1(shift+1:shift+nvar)
766 : p3(1:nvar) => b1(shift+1:shift+nvar)
767 : p4(1:nvar) => b1(shift+1:shift+nvar)
768 : call do_multiply_xa_plus_c(nvar,p1,p2,p3,p4)
769 : end if
770 : if (k < nz) then
771 : p1(1:nvar2) => lblk1(shift2+1:shift2+nvar2)
772 : p2(1:nvar) => x1(shift+1+nvar:shift+2*nvar)
773 : p3(1:nvar) => b1(shift+1:shift+nvar)
774 : p4(1:nvar) => b1(shift+1:shift+nvar)
775 : call do_multiply_xa_plus_c(nvar,p1,p2,p3,p4)
776 : end if
777 : end do
778 : !$OMP END PARALLEL DO
779 0 : end subroutine do_block_multiply_xa
780 :
781 :
782 20 : subroutine do_band_multiply_xa(n, kl, ku, ab1, ldab, x, b)
783 : ! calculates b = x*a = transpose(a)*x
784 : integer, intent(in) :: n
785 : ! the number of linear equations, i.e., the order of the
786 : ! matrix a. n >= 0.
787 : integer, intent(in) :: kl
788 : ! the number of subdiagonals within the band of a. kl >= 0.
789 : integer, intent(in) :: ku
790 : ! the number of superdiagonals within the band of a. ku >= 0.
791 : integer, intent(in) :: ldab
792 : ! the leading dimension of the array ab. ldab >= kl+ku+1.
793 : real(dp), intent(in), pointer :: ab1(:) ! =(ldab, n)
794 : ! the matrix a in band storage, in rows 1 to kl+ku+1;
795 : ! the j-th column of a is stored in the j-th column of the
796 : ! array ab as follows:
797 : ! ab(ku+1+i-j, j) = a(i, j) for max(1, j-ku)<=i<=min(n, j+kl)
798 : real(dp), intent(in), pointer :: x(:) ! (n)
799 : ! the input vector to be multiplied by the matrix.
800 : real(dp), intent(inout), pointer :: b(:) ! (n)
801 : ! on exit, set to matrix product of x*a = b
802 : integer :: i, j, k
803 20 : real(dp), pointer :: ab(:,:)
804 20 : ab(1:ldab,1:n) => ab1(1:ldab*n)
805 40060 : do j = 1, n
806 40040 : k = ku+1-j
807 40040 : b(j) = 0
808 320100 : do i = max(1,j-ku), min(n,j+kl)
809 320080 : b(j) = b(j) + x(i)*ab(k+i,j)
810 : end do
811 : end do
812 20 : end subroutine do_band_multiply_xa
813 :
814 :
815 0 : subroutine do_clip_blocks( &
816 0 : mblk, clip_limit, lmat, dmat, umat, dmat_nnz, total_nnz)
817 : integer, intent(in) :: mblk
818 : real(dp), intent(in) :: clip_limit
819 : real(dp), intent(inout) :: lmat(:,:), dmat(:,:), umat(:,:)
820 : integer, intent(inout) :: dmat_nnz, total_nnz
821 : integer :: i, j
822 0 : dmat_nnz = 0; total_nnz = 0
823 0 : do j=1,mblk
824 0 : do i=1,mblk
825 0 : if (i /= j .and. abs(lmat(i,j)) < clip_limit) lmat(i,j) = 0d0
826 0 : if (lmat(i,j) /= 0) total_nnz = total_nnz + 1
827 0 : if (i /= j .and. abs(dmat(i,j)) < clip_limit) dmat(i,j) = 0d0
828 0 : if (dmat(i,j) /= 0) then
829 0 : total_nnz = total_nnz + 1
830 0 : dmat_nnz = dmat_nnz + 1
831 : end if
832 0 : if (i /= j .and. abs(umat(i,j)) < clip_limit) umat(i,j) = 0d0
833 0 : if (umat(i,j) /= 0) total_nnz = total_nnz + 1
834 : end do
835 : end do
836 0 : end subroutine do_clip_blocks
837 :
838 :
839 0 : subroutine do_clip_block(mblk, clip_limit, dmat, dmat_nnz)
840 : integer, intent(in) :: mblk
841 : real(dp), intent(in) :: clip_limit
842 : real(dp), intent(inout) :: dmat(:,:)
843 : integer, intent(inout) :: dmat_nnz
844 : integer :: i, j
845 0 : dmat_nnz = 0
846 0 : do j=1,mblk
847 0 : do i=1,mblk
848 0 : if (i /= j .and. abs(dmat(i,j)) < clip_limit) dmat(i,j) = 0d0
849 0 : if (dmat(i,j) /= 0) dmat_nnz = dmat_nnz + 1
850 : end do
851 : end do
852 0 : end subroutine do_clip_block
853 :
854 :
855 2935 : subroutine read_hbcode1(iounit, nrow, ncol, nnzero, values, rowind, colptr, ierr)
856 :
857 : character :: TITLE*72, KEY*8, MXTYPE*3, PTRFMT*16, INDFMT*16, VALFMT*20, RHSFMT*20
858 :
859 : integer :: TOTCRD, PTRCRD, INDCRD, VALCRD, RHSCRD, iounit, NROW , NCOL , NNZERO, NELTVL
860 :
861 : integer, pointer :: COLPTR (:), ROWIND (:)
862 :
863 : real(dp), pointer :: VALUES (:)
864 : integer, intent(out) :: ierr
865 :
866 : integer :: i
867 2935 : ierr = 0
868 2935 : READ (iounit, 1000, iostat=ierr ) TITLE , KEY , &
869 2935 : TOTCRD, PTRCRD, INDCRD, VALCRD, RHSCRD, &
870 2935 : MXTYPE, NROW , NCOL , NNZERO, NELTVL, &
871 5870 : PTRFMT, INDFMT, VALFMT, RHSFMT
872 2935 : if (ierr /= 0) return
873 : 1000 FORMAT ( A72, A8 / 5I14 / A3, 11X, 4I14 / 2A16, 2A20 )
874 :
875 2935 : allocate(VALUES(NNZERO), ROWIND(NNZERO), COLPTR(NCOL+1), stat=ierr)
876 2935 : if (ierr /= 0) return
877 :
878 2935 : READ (iounit, PTRFMT, iostat=ierr ) ( COLPTR (I), I = 1, NCOL+1 )
879 2935 : if (ierr /= 0) return
880 :
881 2935 : READ (iounit, INDFMT, iostat=ierr ) ( ROWIND (I), I = 1, NNZERO )
882 2935 : if (ierr /= 0) return
883 :
884 2935 : IF ( VALCRD > 0 ) THEN
885 :
886 : ! ----------------------
887 : ! ... READ MATRIX VALUES
888 : ! ----------------------
889 :
890 2935 : READ (iounit, VALFMT, iostat=ierr ) ( VALUES (I), I = 1, NNZERO )
891 2935 : if (ierr /= 0) return
892 :
893 : end if
894 :
895 :
896 : end subroutine read_hbcode1
897 :
898 :
899 0 : subroutine write_hbcode1(iounit, nrow, ncol, nnzero, values, rowind, colptr, ierr)
900 :
901 : character :: TITLE*72 , KEY*8, MXTYPE*3, PTRFMT*16, INDFMT*16, use_VALFMT*20, VALFMT*20, RHSFMT*20
902 :
903 : integer :: TOTCRD, PTRCRD, INDCRD, VALCRD, RHSCRD, iounit, NROW , NCOL , NNZERO, NELTVL
904 :
905 : integer :: COLPTR (*), ROWIND (*), ierr
906 :
907 : real(dp) :: VALUES (*)
908 :
909 : integer :: i
910 :
911 0 : ierr = 0
912 :
913 : ! ------------------------
914 : ! ... WRITE HEADER BLOCK
915 : ! ------------------------
916 :
917 : ! Line 1 (A72, A8)
918 : ! Col. 1 - 72 Title (TITLE)
919 : ! Col. 73 - 80 Matrix name / identifier (MTRXID)
920 : !
921 : ! Line 2 (I14, 3(1X, I13))
922 : ! Col. 1 - 14 Total number of lines excluding header (TOTCRD)
923 : ! Col. 16 - 28 Number of lines for pointers (PTRCRD)
924 : ! Col. 30 - 42 Number of lines for row (or variable) indices (INDCRD)
925 : ! Col. 44 - 56 Number of lines for numerical values (VALCRD)
926 : !
927 : ! Line 3 (A3, 11X, 4(1X, I13))
928 : ! Col. 1 - 3 Matrix type (see below) (MXTYPE)
929 : ! Col. 15 - 28 Compressed Column: Number of rows (NROW)
930 : ! Elemental: Largest integer used to index variable (MVAR)
931 : ! Col. 30 - 42 Compressed Column: Number of columns (NCOL)
932 : ! Elemental: Number of element matrices (NELT)
933 : ! Col. 44 - 56 Compressed Column: Number of entries (NNZERO)
934 : ! Elemental: Number of variable indices (NVARIX)
935 : ! Col. 58 - 70 Compressed Column: Unused, explicitly zero
936 : ! Elemental: Number of elemental matrix entries (NELTVL)
937 : !
938 : ! Line 4 (2A16, A20)
939 : ! Col. 1 - 16 Fortran format for pointers (PTRFMT)
940 : ! Col. 17 - 32 Fortran format for row (or variable) indices (INDFMT)
941 : ! Col. 33 - 52 Fortran format for numerical values of coefficient matrix
942 : ! (VALFMT)
943 : ! (blank in the case of matrix patterns)
944 : !
945 : ! The three character type field on line 3 describes the matrix type.
946 : ! The following table lists the permitted values for each of the three
947 : ! characters. As an example of the type field, RSA denotes that the matrix
948 : ! is real, symmetric, and assembled.
949 : !
950 : ! First Character:
951 : ! R Real matrix
952 : ! C Complex matrix
953 : ! I integer matrix
954 : ! P Pattern only (no numerical values supplied)
955 : ! Q Pattern only (numerical values supplied in associated auxiliary value
956 : ! file)
957 : !
958 : ! Second Character:
959 : ! S Symmetric
960 : ! U Unsymmetric
961 : ! H Hermitian
962 : ! Z Skew symmetric
963 : ! R Rectangular
964 : !
965 : ! Third Character:
966 : ! A Compressed column form
967 : ! E Elemental form
968 : !
969 :
970 :
971 :
972 0 : TITLE = ''
973 0 : KEY = ''
974 :
975 0 : PTRFMT = '(10I8)'
976 0 : INDFMT = '(12I6)'
977 0 : use_VALFMT = '(5(1pE27.16))'
978 0 : VALFMT = '(5E27.16)'
979 0 : RHSFMT = ''
980 :
981 0 : PTRCRD = (NCOL+1)/10 + 1 ! number of lines for COLPTR
982 0 : INDCRD = NNZERO/12 + 1 ! number of lines for ROWIND
983 0 : VALCRD = NNZERO/5 + 1 ! number of lines for VALUES
984 0 : RHSCRD = 0
985 0 : TOTCRD = 3 + PTRCRD + INDCRD + VALCRD + RHSCRD
986 :
987 0 : MXTYPE = 'RUA'
988 0 : NELTVL = 0
989 :
990 0 : WRITE (iounit, 1000 ) TITLE , KEY , &
991 0 : TOTCRD, PTRCRD, INDCRD, VALCRD, RHSCRD, &
992 0 : MXTYPE, NROW , NCOL , NNZERO, NELTVL, &
993 0 : PTRFMT, INDFMT, VALFMT, RHSFMT
994 : 1000 FORMAT ( A72, A8 / 5I14 / A3, 11X, 4I14 / 2A16, 2A20 )
995 :
996 : ! -------------------------
997 : ! ... WRITE MATRIX STRUCTURE
998 : ! -------------------------
999 :
1000 0 : WRITE (iounit, PTRFMT ) ( COLPTR (I), I = 1, NCOL+1 )
1001 :
1002 0 : WRITE (iounit, INDFMT ) ( ROWIND (I), I = 1, NNZERO )
1003 :
1004 0 : IF ( VALCRD > 0 ) THEN
1005 :
1006 : ! ----------------------
1007 : ! ... WRITE MATRIX VALUES
1008 : ! ----------------------
1009 :
1010 0 : WRITE (iounit, use_VALFMT ) ( VALUES (I), I = 1, NNZERO )
1011 :
1012 : end if
1013 :
1014 0 : return
1015 : end subroutine write_hbcode1
1016 :
1017 :
1018 1 : subroutine read_block_tridiagonal(iounit,nvar,nblk,lblk1,dblk1,ublk1,ierr)
1019 : integer, intent(in) :: iounit
1020 : integer, intent(out) :: nvar, nblk
1021 : real(dp), pointer, dimension(:) :: lblk1,dblk1,ublk1 ! will be allocated
1022 : integer, intent(out) :: ierr
1023 : integer :: k
1024 1 : real(dp), pointer, dimension(:,:,:) :: lblk,dblk,ublk
1025 1 : ierr = 0
1026 1 : read(iounit,*,iostat=ierr) nvar, nblk
1027 1 : if (ierr /= 0) return
1028 1 : allocate(lblk1(nvar*nvar*nblk), dblk1(nvar*nvar*nblk), ublk1(nvar*nvar*nblk), stat=ierr)
1029 1 : if (ierr /= 0) return
1030 1 : lblk(1:nvar,1:nvar,1:nblk) => lblk1(1:nvar*nvar*nblk)
1031 1 : dblk(1:nvar,1:nvar,1:nblk) => dblk1(1:nvar*nvar*nblk)
1032 1 : ublk(1:nvar,1:nvar,1:nblk) => ublk1(1:nvar*nvar*nblk)
1033 980 : do k=1,nblk
1034 979 : if (k > 1) then
1035 978 : call read1_sparse_block(iounit, nvar, lblk(:,:,k), ierr)
1036 978 : if (ierr /= 0) return
1037 : end if
1038 979 : call read1_sparse_block(iounit, nvar, dblk(:,:,k), ierr)
1039 979 : if (ierr /= 0) return
1040 980 : if (k < nblk) then
1041 978 : call read1_sparse_block(iounit, nvar, ublk(:,:,k), ierr)
1042 978 : if (ierr /= 0) return
1043 : end if
1044 : end do
1045 :
1046 1 : end subroutine read_block_tridiagonal
1047 :
1048 :
1049 2935 : subroutine read1_sparse_block(iounit, nvar, blk, ierr)
1050 : integer, intent(in) :: iounit, nvar
1051 : real(dp) :: blk(:,:) ! (nvar,nvar)
1052 : integer, intent(out) :: ierr
1053 : integer :: nnz, nrow, ncol
1054 2935 : integer, pointer :: rowind(:), colptr(:)
1055 2935 : real(dp), pointer :: values(:)
1056 : ierr = 0
1057 2935 : call read_hbcode1(iounit, nrow, ncol, nnz, values, rowind, colptr,ierr)
1058 2935 : if (ierr /= 0 .or. nrow /= nvar .or. nrow /= ncol) return
1059 2935 : call do_column_sparse_to_dense(nrow,ncol,blk,nnz,colptr,rowind,values,ierr)
1060 2935 : deallocate(colptr,rowind,values)
1061 2935 : end subroutine read1_sparse_block
1062 :
1063 :
1064 0 : subroutine write_block_tridiagonal(iounit,nvar,nblk,lblk,dblk,ublk,ierr)
1065 : integer, intent(in) :: iounit, nvar, nblk
1066 : real(dp), intent(in), dimension(:,:,:) :: lblk,dblk,ublk
1067 : integer, intent(out) :: ierr
1068 : integer :: k
1069 0 : ierr = 0
1070 0 : write(iounit,*) nvar, nblk
1071 0 : do k=1,nblk
1072 0 : if (k > 1) then
1073 0 : call write1_sparse_block(iounit, nvar, lblk(:,:,k), ierr)
1074 0 : if (ierr /= 0) return
1075 : end if
1076 0 : call write1_sparse_block(iounit, nvar, dblk(:,:,k), ierr)
1077 0 : if (ierr /= 0) return
1078 0 : if (k < nblk) then
1079 0 : call write1_sparse_block(iounit, nvar, ublk(:,:,k), ierr)
1080 0 : if (ierr /= 0) return
1081 : end if
1082 : end do
1083 : end subroutine write_block_tridiagonal
1084 :
1085 :
1086 0 : subroutine write1_sparse_block(iounit, nvar, blk, ierr)
1087 : integer, intent(in) :: iounit, nvar
1088 : real(dp), intent(in) :: blk(:,:) ! (nvar,nvar)
1089 : integer, intent(out) :: ierr
1090 0 : integer :: nnz, rowind(nvar*nvar), colptr(nvar+1)
1091 0 : real(dp) :: values(nvar*nvar)
1092 : call do_dense_to_column_sparse( &
1093 0 : nvar, nvar, blk, nvar*nvar, nnz, colptr, rowind, values, .true., ierr)
1094 0 : if (ierr /= 0) return
1095 0 : call write_hbcode1(iounit, nvar, nvar, nnz, values, rowind, colptr, ierr)
1096 : end subroutine write1_sparse_block
1097 :
1098 :
1099 : end module mtx_support
|
