!>
!> @file tmatrix_gb.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
PROGRAM main
!
!   Test some routines of module matrix
!
  USE matrix
  IMPLICIT NONE
  TYPE(gbmat) :: mata
  INTEGER, PARAMETER :: n=5, ku=3, kl=3
  DOUBLE PRECISION :: arr(n), fulla(n,n), base
  DOUBLE PRECISION, DIMENSION(:,:), POINTER :: p
  INTEGER :: i, j, iaway, pow
  CHARACTER(len=32) :: str
!
  CALL init(ku, ku, n, 0, mata)
  CALL getvalp(mata, p)
  PRINT*, 'shape of A: ', SHAPE(p)
!
!   Test updtmat
  p = 0.0d0
  DO j=1,n
     DO i=1,n
        arr(i) = 10*i + j
        IF( ABS(j-i) .LT. ku+1 ) CALL updtmat(mata, i, j, arr(i))
     END DO
  END DO
  CALL prntmat('Test of UPDTMAT', p)
!
!   Test PUTCOL
  p = 0.0d0
  DO j=1,n
     DO i=1,n
        arr(i) = 10*i + j
     END DO
     CALL putcol(mata, j, arr)
  END DO
  CALL prntmat('Test of PUTCOL', p)
!
!
!   Test PUTROW
  p = 0.0d0
  DO i=1,n
     DO j=1,n
        arr(j) = 10*i + j
     END DO
     CALL putrow(mata, i, arr)
  END DO
  CALL prntmat('Test of PUTROW', p)
!
!   Test GETCOL
  fulla = 0.0
  DO j=1,n
     CALL getcol(mata, j, fulla(:,j))
  END DO
  CALL prntmat('Test of GETCOL', fulla)
!
  iaway=4
  arr = 0.0d0
  arr(iaway) =1.0
  CALL putrow(mata, iaway, arr)
  CALL putcol(mata, iaway, arr)
  WRITE(str,'(a,i3)') 'Away on i = j =',iaway
  CALL prntmat(TRIM(str), p)
!
!   Test GETCOL
  fulla = 0.0
  DO j=1,n
     CALL getcol(mata, j, fulla(:,j))
  END DO
  CALL prntmat('Matrix full', fulla)
!
!   Test of determinant
  CALL determinant(mata, base, pow)
  CALL prntmat('Factored A (gb)', p)
  PRINT*, 'Prod. of factored A diagnonals', PRODUCT(p(kl+ku+1,:))
  WRITE(*,'(a,1pe15.6,i6)') 'Determinant(bas,power) =', base, pow
  PRINT*, 'Pivots ', mata%piv
!
  call destroy(mata)
CONTAINS
  SUBROUTINE prntmat(str, a)
    DOUBLE PRECISION, DIMENSION(:,:) :: a
    CHARACTER(len=*) :: str
    INTEGER :: i
    WRITE(*,'(a)') TRIM(str)
    DO i=1,SIZE(a,1)
       WRITE(*,'(10f8.1)') a(i,:)
    END DO
  END SUBROUTINE prntmat
END PROGRAM main