!***************************************** !* Rasmus Handberg !* Department of Physics and Astronomy !* University of Aarhus !***************************************** module matrix interface function backsub(T, c) implicit none real, intent(in) :: T(:,:) real, intent(in) :: c(:) real :: backsub(size(c)) end function backsub subroutine qrdec(matrix, Q, R) implicit none real, dimension(:,:) :: matrix real, dimension(size(matrix,1),size(matrix,2)) :: Q real, dimension(size(matrix,2),size(matrix,2)) :: R end subroutine qrdec real function determinant(R) implicit none real, dimension(:,:) :: R end function determinant function lineqnsystem(Q, R, b) implicit none real, dimension(:,:) :: Q real, dimension(:,:) :: R real, dimension(:) :: b real, dimension(size(b)) :: lineqnsystem end function lineqnsystem function inverse(Q, R) implicit none real, dimension(:,:) :: Q real, dimension(:,:) :: R real, dimension(size(Q,1),size(Q,2)) :: inverse end function inverse function inv(A) implicit none real, dimension(:,:) :: A real, dimension(size(A,1),size(A,2)) :: inv end function inv function unit_matrix(n) implicit none integer, intent(in) :: n real :: unit_matrix(n,n) end function unit_matrix function upper_triangle(n,m,incdiag) implicit none integer, intent(in) :: n, m integer, optional, intent(in) :: incdiag logical, dimension(n,m) :: upper_triangle end function upper_triangle real function norm(v) implicit none real, intent(in) :: v(:) end function norm end interface ! Generic function print_matrix: interface print_matrix module procedure MatOut module procedure VecOut end interface contains subroutine MatOut(streng, A) implicit none character(len=*) :: streng real, dimension(:,:) :: A integer :: i, n, m n = size(A,1) m = size(A,2) print '(1x,a)', streng(1:len_trim(streng)) do i = 1,n print '(10f10.3)', A(i,1:m) enddo end subroutine MatOut subroutine VecOut(streng, A) implicit none character(len=*) :: streng real, dimension(:) :: A integer :: i, n n = size(A) print '(1x,a)', streng(1:len_trim(streng)) do i = 1,n print '(f10.3)', A(i) enddo end subroutine VecOut end module matrix ! Back-substitution. function backsub(T, c) implicit none real, intent(in) :: T(:,:) real, intent(in) :: c(:) real :: backsub(size(c)) integer :: i, n n = size(c) do i = n,1,-1 backsub(i) = (c(i) - dot_product(T(i,i+1:n),backsub(i+1:n)))/T(i,i) enddo end function backsub ! QR-decomposition. ! Input: matrix - n*m-matrix to be decomposed. ! Q - n*m-matrix ! R - m*m-matrix subroutine qrdec(matrix, Q, R) implicit none real, intent(in) :: matrix(:,:) real, intent(out) :: Q(size(matrix,1),size(matrix,2)) real, intent(out) :: R(size(matrix,2),size(matrix,2)) real :: A(size(matrix,1),size(matrix,2)) integer :: i, j, n, m n = size(matrix,1) m = size(matrix,2) A = matrix Q = 0.0 R = 0.0 do i = 1,m R(i,i) = sqrt(dot_product(A(:,i),A(:,i))) Q(:,i) = A(:,i)/R(i,i) do j = i+1,m R(i,j) = dot_product(Q(:,i),A(:,j)) A(:,j) = A(:,j) - Q(:,i)*R(i,j) enddo enddo end subroutine qrdec ! Calcultaes the determinant, given the R-matrix ! from the QR-decomposition. real function determinant(R) implicit none real, intent(in) :: R(:,:) integer :: i determinant = 1.0 do i = 1,size(R,1) determinant = determinant*R(i,i) enddo end function determinant ! Solve linear equation system Ax=b, given the ! Q and R-matrices from the QR-decomposition and ! the vector b. function lineqnsystem(Q, R, b) use matrix, only : backsub implicit none real, intent(in) :: Q(:,:) real, intent(in) :: R(:,:) real, intent(in) :: b(:) real :: lineqnsystem(size(Q,2)) ! size(R,2) real :: c(size(Q,2)) c = matmul(transpose(Q), b) lineqnsystem = backsub(R, c) end function lineqnsystem ! Calculate the inverse: function inverse(Q, R) use matrix, only : lineqnsystem implicit none real, intent(in) :: Q(:,:) real, intent(in) :: R(:,:) real :: inverse(size(Q,1),size(Q,2)) real :: z(size(Q,1)) integer :: i inverse = 0.0 do i = 1,size(Q,1) z = 0.0 z(i) = 1.0 inverse(:,i) = lineqnsystem(Q, R, z) enddo end function inverse function inv(A) use matrix, only : qrdec, lineqnsystem implicit none real, intent(in) :: A(:,:) real :: Q(size(A,1),size(A,2)) real :: R(size(A,2),size(A,2)) real :: inv(size(A,1),size(A,2)) real :: z(size(A,1)) integer :: i call qrdec(A,Q,R) inv = 0.0 do i = 1,size(Q,1) z = 0.0 z(i) = 1.0 inv(:,i) = lineqnsystem(Q, R, z) enddo end function inv ! Create the n*n unit-matrix. function unit_matrix(n) implicit none integer, intent(in) :: n real :: unit_matrix(n,n) integer :: i unit_matrix = 0.0 do i = 1,n unit_matrix(i,i) = 1.0 enddo end function unit_matrix ! Returns an upper-triangle matrix-mask of size n*m. function upper_triangle(n, m, incdiag) implicit none integer, intent(in) :: n, m integer, optional, intent(in) :: incdiag logical, dimension(n,m) :: upper_triangle integer :: i, j, d d = 0 if (present(incdiag)) d = incdiag do i = 1,n do j = 1,m upper_triangle(i,j) = (i-j < d) enddo enddo end function upper_triangle ! Calcultae the norm of a vector: real function norm(v) implicit none real, intent(in) :: v(:) norm = sqrt(dot_product(v,v)) end function norm