Skip to content

Commit

Permalink
Merge pull request #21 from NLESC-JCER/free
Browse files Browse the repository at this point in the history
moved miscellaneous to utils
  • Loading branch information
felipeZ authored Feb 22, 2019
2 parents a9eb0e3 + 5f623f6 commit 853e460
Show file tree
Hide file tree
Showing 13 changed files with 552 additions and 711 deletions.
4 changes: 2 additions & 2 deletions src/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@

# create binary
set(UTILS numeric_kinds.f90 lapack_wrapper.f90)
set(SOURCES main.f90 davidson.f90 benchmark.f90 ${UTILS})
set(UTILS array_utils.f90 numeric_kinds.f90 lapack_wrapper.f90)
set(SOURCES main.f90 davidson.f90 ${UTILS})
message (STATUS "SOURCES: " ${SOURCES})
add_executable(main ${SOURCES})

Expand Down
135 changes: 135 additions & 0 deletions src/array_utils.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,135 @@
module array_utils

use numeric_kinds, only: dp
use lapack_wrapper, only: lapack_generalized_eigensolver, lapack_matmul, lapack_matrix_vector, &
lapack_qr, lapack_solver
implicit none

!> \private
private
!> \public
public :: concatenate, diagonal,eye, generate_diagonal_dominant, norm

contains

pure function eye(m, n, alpha)
!> Create a matrix with ones in the diagonal and zero everywhere else
!> \param m: number of rows
!> \param n: number of colums
!> \param alpha: optional diagonal value
!> \return matrix of size n x m
integer, intent(in) :: n, m
real(dp), dimension(m, n) :: eye
real(dp), intent(in), optional :: alpha

!local variable
integer :: i, j
real(dp) :: x

! check optional values
x = 1.d0
if (present(alpha)) x = alpha

do i=1, m
do j=1, n
if (i /= j) then
eye(i, j) = 0
else
eye(i, i) = x
end if
end do
end do

end function eye

pure function norm(vector)
!> compute the norm-2 of a vector
real(dp), dimension(:), intent(in) :: vector
real(dp) :: norm

norm = sqrt(sum(vector ** 2.d0))

end function norm

subroutine concatenate(arr, brr)

!> Concatenate two matrices
!> \param arr: first array
!> \param brr: second array
!> \return arr concatenate brr (overwrites arr)

real(dp), dimension(:, :), intent(inout), allocatable :: arr
real(dp), dimension(:, :), intent(in) :: brr
real(dp), dimension(:, :), allocatable :: tmp_array
integer :: new_dim, dim_cols, dim_rows

! dimension
dim_rows = size(arr, 1)
dim_cols = size(arr, 2)

! Number of columns of the new matrix
new_dim = dim_cols + size(brr, 2)

! move to temporal array
allocate(tmp_array(dim_rows, new_dim))
tmp_array(:, :dim_cols) = arr

! Move to new expanded matrix
deallocate(arr)
call move_alloc(tmp_array, arr)

arr(:, dim_cols + 1:) = brr

end subroutine concatenate

function generate_diagonal_dominant(m, sparsity, diag_val) result(arr)
!> Generate a diagonal dominant square matrix of dimension m
!> \param m dimension of the matrix
!> \param sparsity magnitude order of the off-diagonal values

integer, intent(in) :: m ! size of the square matrix
real(dp), optional :: diag_val
integer :: i, j
real(dp) :: sparsity
real(dp), dimension(m, m) :: arr
call random_number(arr)

arr = arr * sparsity
do j=1, m
do i=1, m
if (i > j) then
arr(i, j) = arr(j, i)
else if(i == j) then
if (present(diag_val))then
arr(i,i) = diag_val
else
arr(i, i) = i
end if
end if
end do
end do

end function generate_diagonal_dominant

function diagonal(matrix)
!> return the diagonal of a matrix
real(dp), dimension(:, :), intent(in) :: matrix
real(dp), dimension(size(matrix, 1)) :: diagonal

! local variables
integer :: i, j, m

! dimension of the matrix
m = size(matrix, 1)

do i=1,m
do j=1,m
if (i == j) then
diagonal(i) = matrix(i, j)
end if
end do
end do

end function diagonal

end module array_utils
94 changes: 0 additions & 94 deletions src/benchmark.f90

This file was deleted.

3 changes: 2 additions & 1 deletion src/benchmark_free.f90
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,8 @@ end module matrix_free

program main
use numeric_kinds, only: dp
use davidson, only: generalized_eigensolver, generate_diagonal_dominant, norm
use davidson, only: generalized_eigensolver
use array_utils, only: generate_diagonal_dominant, norm
use matrix_free, only: compute_matrix_on_the_fly, compute_stx_on_the_fly

implicit none
Expand Down
Loading

0 comments on commit 853e460

Please sign in to comment.