Skip to content

Commit

Permalink
mtx: eliminate all quad precision code
Browse files Browse the repository at this point in the history 
None of this code is used in the other parts of MESA.
  • Loading branch information
@VincentVanlaer
VincentVanlaer committed Sep 30, 2024
1 parent c065f89 commit 97e531b
Show file tree
Hide file tree
Showing 9 changed files with 17 additions and 1,040 deletions.
296 changes: 0 additions & 296 deletions mtx/private/mtx_support.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -558,39 +558,6 @@ subroutine do_dense_to_col_sparse_0_based_qp( &
end subroutine do_dense_to_col_sparse_0_based_qp


subroutine do_quad_dense_to_col_sparse_0_based( &
n,ndim,a,nzmax,nz,colptr,rowind,values,diags,ierr)
integer, intent(in) :: n,ndim,nzmax
real(qp), intent(in) :: a(:,:) ! (ndim,n)
integer, intent(inout) :: colptr(:) ! (n+1)
integer, intent(inout) :: rowind(:) ! (nzmax)
real(qp), intent(out) :: values(:) ! (nzmax)
logical, intent(in) :: diags
integer, intent(out) :: nz,ierr
integer :: i,j
ierr = 0
nz = 0
do j=1,n
colptr(j) = nz ! index in values of first entry in this column
do i=1,n
if (a(i,j) == 0) then
if (i /= j) cycle ! not a diagonal
if (.not. diags) cycle
! else include diagonals even if are == 0
end if
nz = nz+1
if (nz > nzmax) then
ierr = j
return
end if
values(nz) = a(i,j)
rowind(nz) = i-1
end do
end do
colptr(n+1) = nz
end subroutine do_quad_dense_to_col_sparse_0_based


subroutine do_column_sparse_to_dense(n,ndim,a,nz,colptr,rowind,values,ierr)
integer, intent(in) :: n,ndim,nz
real(dp), intent(inout) :: a(ndim,n)
Expand All @@ -617,28 +584,6 @@ subroutine do_column_sparse_to_dense(n,ndim,a,nz,colptr,rowind,values,ierr)
end subroutine do_column_sparse_to_dense


subroutine do_quad_column_sparse_to_dense(n,ndim,a,nz,colptr,rowind,values,ierr)
integer, intent(in) :: n,ndim,nz
real(qp), intent(out) :: a(ndim,n)
integer, intent(in) :: colptr(n+1),rowind(nz)
real(qp), intent(in) :: values(nz)
integer, intent(out) :: ierr
integer :: i,j,k
ierr = 0
a = 0
do j=1,n
do k=colptr(j),colptr(j+1)-1
i = rowind(k)
if (i > n) then
ierr = j
return
endif
a(i,j) = values(k)
end do
end do
end subroutine do_quad_column_sparse_to_dense


subroutine do_col_sparse_0_based_to_dense(n,ndim,a,nz,colptr,rowind,values,ierr)
integer, intent(in) :: n,ndim,nz
real(dp), intent(inout) :: a(ndim,n)
Expand Down Expand Up @@ -685,56 +630,6 @@ subroutine do_block_dble_mv(nvar, nz, lblk, dblk, ublk, b, prod)
end subroutine do_block_dble_mv


subroutine do_block_mv_quad(lblk, dblk, ublk, b, prod)
real(qp), pointer, dimension(:,:,:), intent(in) :: lblk, dblk, ublk ! (nvar,nvar,nz)
real(qp), pointer, dimension(:,:), intent(in) :: b ! (nvar,nz)
real(qp), pointer, dimension(:,:), intent(inout) :: prod ! (nvar,nz)
integer :: nvar, nz, k
include 'formats'
nvar = size(b,dim=1)
nz = size(b,dim=2)
!$OMP PARALLEL DO PRIVATE(k)
do k = 1, nz
prod(:,k) = 0
call qdgemv(nvar,nvar,dblk(:,:,k),nvar,b(:,k),prod(:,k))
if (k > 1) then
call qdgemv(nvar,nvar,lblk(:,:,k),nvar,b(:,k-1),prod(:,k))
end if
if (k < nz) then
call qdgemv(nvar,nvar,ublk(:,:,k),nvar,b(:,k+1),prod(:,k))
end if
end do
!$OMP END PARALLEL DO

contains

subroutine qdgemv(m,n,a,lda,x,y)
! y := alpha*a*x + beta*y
use const_def, only: dp

integer lda,m,n
real(qp) a(lda,*),x(*),y(*)
real(qp) :: tmp
! trans = 'n'
! alpha = 1
! beta = 1
! incx = 1
! incy = 1
integer :: j, i
do j = 1,n
tmp = x(j)
if (tmp.ne.0d0) then
do i = 1,m
y(i) = y(i) + tmp*a(i,j)
end do
end if
end do
end subroutine qdgemv


end subroutine do_block_mv_quad


subroutine do_LU_factored_block_dble_mv(lblk, dblk, ublk, b, ipiv, prod)
real(dp), pointer, dimension(:,:,:), intent(in) :: lblk, dblk, ublk ! (nvar,nvar,nz)
real(dp), pointer, dimension(:,:), intent(in) :: b ! (nvar,nz)
Expand Down Expand Up @@ -836,21 +731,6 @@ subroutine do_multiply_xa(n, A1, x, b)
end subroutine do_multiply_xa


subroutine do_quad_multiply_xa(n, A1, x, b)
! calculates b = x*A
integer, intent(in) :: n
real(qp), pointer, intent(in) :: A1(:) ! =(n, n)
real(qp), pointer, intent(in) :: x(:) ! (n)
real(qp), pointer, intent(inout) :: b(:) ! (n)
integer :: j
real(qp), pointer :: A(:,:) ! (n, n)
A(1:n,1:n) => A1(1:n*n)
do j = 1, n
b(j) = dot_product(x(1:n),A(1:n,j))
end do
end subroutine do_quad_multiply_xa


subroutine do_multiply_xa_plus_c(n, A1, x, c, b)
! calculates b = x*A + c
integer, intent(in) :: n
Expand All @@ -867,22 +747,6 @@ subroutine do_multiply_xa_plus_c(n, A1, x, c, b)
end subroutine do_multiply_xa_plus_c


subroutine do_quad_multiply_xa_plus_c(n, A1, x, c, b)
! calculates b = x*A + c
integer, intent(in) :: n
real(qp), pointer, intent(in) :: A1(:) ! =(n, n)
real(qp), pointer, intent(in) :: x(:) ! (n)
real(qp), pointer, intent(in) :: c(:) ! (n)
real(qp), pointer, intent(inout) :: b(:) ! (n)
integer :: j
real(qp), pointer :: A(:,:) ! (n, n)
A(1:n,1:n) => A1(1:n*n)
do j = 1, n
b(j) = dot_product(x(1:n),A(1:n,j)) + c(j)
end do
end subroutine do_quad_multiply_xa_plus_c


subroutine do_block_multiply_xa(nvar, nz, lblk1, dblk1, ublk1, x1, b1)
! calculates b = x*A
integer, intent(in) :: nvar, nz
Expand Down Expand Up @@ -919,42 +783,6 @@ subroutine do_block_multiply_xa(nvar, nz, lblk1, dblk1, ublk1, x1, b1)
end subroutine do_block_multiply_xa


subroutine do_quad_block_multiply_xa(nvar, nz, lblk1, dblk1, ublk1, x1, b1)
! calculates b = x*A
integer, intent(in) :: nvar, nz
real(qp), dimension(:), intent(in), pointer :: lblk1, dblk1, ublk1 ! =(nvar,nvar,nz)
real(qp), intent(in), pointer :: x1(:) ! =(nvar,nz)
real(qp), intent(inout), pointer :: b1(:) ! =(nvar,nz)
integer :: k, shift, shift2, nvar2
real(qp), pointer, dimension(:) :: p1, p2, p3, p4
nvar2 = nvar*nvar
!$OMP PARALLEL DO PRIVATE(k,shift,shift2,p1,p2,p3,p4)
do k = 1, nz
shift = nvar*(k-1)
shift2 = nvar2*(k-1)
p1(1:nvar2) => dblk1(shift2+1:shift2+nvar2)
p2(1:nvar) => x1(shift+1:shift+nvar)
p3(1:nvar) => b1(shift+1:shift+nvar)
call do_quad_multiply_xa(nvar,p1,p2,p3)
if (k > 1) then
p1(1:nvar2) => ublk1(shift2+1:shift2+nvar2)
p2(1:nvar) => x1(shift+1:shift+nvar)
p3(1:nvar) => b1(shift+1:shift+nvar)
p4(1:nvar) => b1(shift+1:shift+nvar)
call do_quad_multiply_xa_plus_c(nvar,p1,p2,p3,p4)
end if
if (k < nz) then
p1(1:nvar2) => lblk1(shift2+1:shift2+nvar2)
p2(1:nvar) => x1(shift+1+nvar:shift+2*nvar)
p3(1:nvar) => b1(shift+1:shift+nvar)
p4(1:nvar) => b1(shift+1:shift+nvar)
call do_quad_multiply_xa_plus_c(nvar,p1,p2,p3,p4)
end if
end do
!$OMP END PARALLEL DO
end subroutine do_quad_block_multiply_xa


subroutine do_band_multiply_xa(n, kl, ku, ab1, ldab, x, b)
! calculates b = x*a = transpose(a)*x
integer, intent(in) :: n
Expand Down Expand Up @@ -988,38 +816,6 @@ subroutine do_band_multiply_xa(n, kl, ku, ab1, ldab, x, b)
end subroutine do_band_multiply_xa


subroutine do_quad_band_multiply_xa(n, kl, ku, ab, ldab, x, b)
! calculates b = x*a = transpose(a)*x
integer, intent(in) :: n
! the number of linear equations, i.e., the order of the
! matrix a. n >= 0.
integer, intent(in) :: kl
! the number of subdiagonals within the band of a. kl >= 0.
integer, intent(in) :: ku
! the number of superdiagonals within the band of a. ku >= 0.
integer, intent(in) :: ldab
! the leading dimension of the array ab. ldab >= kl+ku+1.
real(qp), intent(in) :: ab(:,:) ! (ldab, n)
! the matrix a in band storage, in rows 1 to kl+ku+1;
! the j-th column of a is stored in the j-th column of the
! array ab as follows:
! ab(ku+1+i-j, j) = a(i, j) for max(1, j-ku)<=i<=min(n, j+kl)
real(qp), intent(in) :: x(:) ! (n)
! the input vector to be multiplied by the matrix.
real(qp), intent(inout) :: b(:) ! (n)
! on exit, set to matrix product of x*a = b
integer :: i, j, k
do j = 1, n
k = ku+1-j
b(j) = 0
do i = max(1, j-ku), min(n, j+kl)
b(j) = b(j) + x(i)*ab(k+i, j)
end do
end do
end subroutine do_quad_band_multiply_xa



subroutine do_clip_blocks( &
mblk, clip_limit, lmat, dmat, umat, dmat_nnz, total_nnz)
integer, intent(in) :: mblk
Expand Down Expand Up @@ -1106,52 +902,6 @@ subroutine read_hbcode1(iounit, nrow, ncol, nnzero, values, rowind, colptr, ierr
end subroutine read_hbcode1


subroutine read_hbcode1_quad(iounit, nrow, ncol, nnzero, values, rowind, colptr, ierr)

CHARACTER TITLE*72 , KEY*8 , MXTYPE*3 , &
PTRFMT*16, INDFMT*16, VALFMT*20, RHSFMT*20

INTEGER TOTCRD, PTRCRD, INDCRD, VALCRD, RHSCRD, &
iounit, NROW , NCOL , NNZERO, NELTVL

INTEGER, pointer :: COLPTR (:), ROWIND (:)

REAL(qp), pointer :: VALUES (:)
integer, intent(out) :: ierr

integer :: i
ierr = 0
READ (iounit, 1000, iostat=ierr ) TITLE , KEY , &
TOTCRD, PTRCRD, INDCRD, VALCRD, RHSCRD, &
MXTYPE, NROW , NCOL , NNZERO, NELTVL, &
PTRFMT, INDFMT, VALFMT, RHSFMT
if (ierr /= 0) return
1000 FORMAT ( A72, A8 / 5I14 / A3, 11X, 4I14 / 2A16, 2A20 )

allocate(VALUES(NNZERO), ROWIND(NNZERO), COLPTR(NCOL+1), stat=ierr)
if (ierr /= 0) return

READ (iounit, PTRFMT, iostat=ierr ) ( COLPTR (I), I = 1, NCOL+1 )
if (ierr /= 0) return

READ (iounit, INDFMT, iostat=ierr ) ( ROWIND (I), I = 1, NNZERO )
if (ierr /= 0) return

IF ( VALCRD .GT. 0 ) THEN

! ----------------------
! ... READ MATRIX VALUES
! ----------------------

READ (iounit, VALFMT, iostat=ierr ) ( VALUES (I), I = 1, NNZERO )
if (ierr /= 0) return

ENDIF


end subroutine read_hbcode1_quad


subroutine write_hbcode1(iounit, nrow, ncol, nnzero, values, rowind, colptr, ierr)

CHARACTER TITLE*72 , KEY*8 , MXTYPE*3 , &
Expand Down Expand Up @@ -1304,37 +1054,6 @@ subroutine read_block_tridiagonal(iounit,nvar,nblk,lblk1,dblk1,ublk1,ierr)
end subroutine read_block_tridiagonal


subroutine read_quad_block_tridiagonal(iounit,nvar,nblk,lblk1,dblk1,ublk1,ierr)
integer, intent(in) :: iounit
integer, intent(out) :: nvar, nblk
real(qp), pointer, dimension(:) :: lblk1,dblk1,ublk1 ! will be allocated
integer, intent(out) :: ierr
integer :: k
real(qp), pointer, dimension(:,:,:) :: lblk,dblk,ublk
ierr = 0
read(iounit,*,iostat=ierr) nvar, nblk
if (ierr /= 0) return
allocate(lblk1(nvar*nvar*nblk), dblk1(nvar*nvar*nblk), ublk1(nvar*nvar*nblk), stat=ierr)
if (ierr /= 0) return
lblk(1:nvar,1:nvar,1:nblk) => lblk1(1:nvar*nvar*nblk)
dblk(1:nvar,1:nvar,1:nblk) => dblk1(1:nvar*nvar*nblk)
ublk(1:nvar,1:nvar,1:nblk) => ublk1(1:nvar*nvar*nblk)
do k=1,nblk
if (k > 1) then
call read1_sparse_block_quad(iounit, nvar, lblk(:,:,k), ierr)
if (ierr /= 0) return
end if
call read1_sparse_block_quad(iounit, nvar, dblk(:,:,k), ierr)
if (ierr /= 0) return
if (k < nblk) then
call read1_sparse_block_quad(iounit, nvar, ublk(:,:,k), ierr)
if (ierr /= 0) return
end if
end do

end subroutine read_quad_block_tridiagonal


subroutine read1_sparse_block(iounit, nvar, blk, ierr)
integer, intent(in) :: iounit, nvar
real(dp) :: blk(:,:) ! (nvar,nvar)
Expand All @@ -1350,21 +1069,6 @@ subroutine read1_sparse_block(iounit, nvar, blk, ierr)
end subroutine read1_sparse_block


subroutine read1_sparse_block_quad(iounit, nvar, blk, ierr)
integer, intent(in) :: iounit, nvar
real(qp) :: blk(:,:) ! (nvar,nvar)
integer, intent(out) :: ierr
integer :: nnz, nrow, ncol
integer, pointer :: rowind(:), colptr(:)
real(qp), pointer :: values(:)
ierr = 0
call read_hbcode1_quad(iounit, nrow, ncol, nnz, values, rowind, colptr,ierr)
if (ierr /= 0 .or. nrow /= nvar .or. nrow /= ncol) return
call do_quad_column_sparse_to_dense(nrow,ncol,blk,nnz,colptr,rowind,values,ierr)
deallocate(colptr,rowind,values)
end subroutine read1_sparse_block_quad


subroutine write_block_tridiagonal(iounit,nvar,nblk,lblk,dblk,ublk,ierr)
integer, intent(in) :: iounit, nvar, nblk
real(dp), intent(in), dimension(:,:,:) :: lblk,dblk,ublk
Expand Down
11 changes: 0 additions & 11 deletions mtx/public/mtx_decsolblk_quad.dek
View file

This file was deleted.

1 change: 0 additions & 1 deletion mtx/public/mtx_def.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,6 @@ module mtx_def
! matrix solver options
integer, parameter :: lapack = 1
integer, parameter :: block_thomas_dble = 2
integer, parameter :: block_thomas_quad = 3
integer, parameter :: block_thomas_refine = 4
integer, parameter :: bcyclic_dble = 5

Expand Down
Loading

0 comments on commit 97e531b

Please sign in to comment.