! calculate generalised inverse of real symmetric, halfstored matrix
! input file: non-zero elements of matrix, one line per element
!             -> row no., column no., element  (space separated)
! program will ask for file name & no. of rows/columns
!                                                        karin meyer

program ginverse

implicit none
real(8), dimension(:), allocatable :: avec
real(8)                            :: xx, det, zero=1.d-9
character(len=30)                  :: fname
logical                            :: lex
integer                            :: neqns, i, j, n, ii, nrank, ij
integer, external                  :: ihmssf

! get file name
write(*,*)'Give name of input file (max 30 char.s): '
read(*,'(a)',iostat=ii) fname
if(ii/=0) then
   write(*,*)'error trying to read file name'; stop
end if

! open file
inquire(file=trim(fname),exist=lex)
if( .not. lex )then
   write(*,*)'file "',trim(fname),'" does not seem to exist'; stop
end if
open(1,file=trim(fname),status='old',form='formatted',action='read')
write(*,*)'"',trim(fname),'" opened '

write(*,*)'Give size (no. of rows/columns) of matrix: '
read(*,*,iostat=ii) neqns
if(ii/=0) then
   write(*,*)'error trying to read number'; stop
end if
write(*,*)'No. of rows/columns specified  = ',neqns

! allocate matrix
allocate(avec(neqns*(neqns+1)/2), stat=ii)
if( ii/= 0 ) then
    write(*,*)'error trying to allocate matrix'; stop
end if

! read matrix
n=0; avec=0.d0
do
  read(1,*,iostat=ii) i, j, xx;  if( ii/=0 )exit
  if(i < 1 .or. i > neqns .or. j<1 .or. j> neqns ) then
     write(*,*)'invalid row or column no. encountered', i,j,xx; stop
  end if
  avec(ihmssf(i,j,neqns)) = xx;   n=n+1
end do
write(*,*)'No. of elements in matrix found =',n
close(1)

! invert
call dkmxhf( avec, det, zero, nrank, neqns, .true.)
write(*,*)'Inversion done'
write(*,*)'Rank of matrix                  =',nrank
write(*,*)'Log determinant                 =', det
n = count(mask=(avec/=0))
write(*,*)'No. of non-zero elements        =',n

! write out matrix: 1st line has determinant & no. of elements!
open(2,file='ginverse.out')
write(2,*)det,n
! ... following lines have non-zero elements
do i=1,neqns
   do j=1,i
      ij = ihmssf(i,j,neqns); if( avec(ij) == 0 ) cycle
      write(2,'(2i10,g18.8)')j,i,avec(ij)
   end do
end do
write(*,*)'Output file is: "ginverse.out" '
stop
end program ginverse

! ===============================
  INTEGER FUNCTION IHMSSF(I,J,N)
! ===============================

  implicit none
  integer, intent(in) :: i, j, n
  integer             :: i1, j1

  IF(I.LE.J)THEN
     I1=I-1; IHMSSF=N*I1-I*I1/2+J
  ELSE
     J1=J-1; IHMSSF=N*J1-J*J1/2+I
  END IF

  RETURN
  END function ihmssf