output_mat_tri.f90

Go to the documentation of this file.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Time-stamp: <2016-10-18 15:35:01 pbrowne>
!!!
!!!    Subroutine to output triangular matrix
!!!    Copyright (C) 2015  Philip A. Browne
!!!
!!!    This program is free software: you can redistribute it and/or modify
!!!    it under the terms of the GNU General Public License as published by
!!!    the Free Software Foundation, either version 3 of the License, or
!!!    (at your option) any later version.
!!!
!!!    This program 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 General Public License
!!!    along with this program.  If not, see <http://www.gnu.org/licenses/>.
!!!
!!!    Email: p.browne @ reading.ac.uk
!!!    Mail:  School of Mathematical and Physical Sciences,
!!!           University of Reading,
!!!       Reading, UK
!!!       RG6 6BB
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine output_mat_tri(n,A,filename,output_type)
  use output_empire, only : unit_mat_tri,emp_e
  implicit none
  integer, parameter :: rk = kind(1.0d0)
  integer, intent(in) :: n
  real(kind=rk), intent(in), dimension(n*(n+1)/2) :: A
  character(40), intent(in) :: filename
  integer, intent(in) :: output_type


  real(kind=rk), dimension(n*(n+1)/2) :: Aout
  integer :: outunit=unit_mat_tri
  logical :: opend
  character(11) :: fm
  integer :: err

  if(output_type .eq. 0) then
     write(emp_e,*) 'Error in output_mat_tri. output_type = 0&
          & unsupported. Stopping.'
     stop '-4'
  elseif(output_type .gt. 0) then
     fm='unformatted'
  else
     fm='formatted'
  end if

  


  !ensure unit is not opened then open file=filename
  do
     inquire(unit=outunit,opened=opend)
     if(.not. opend) then
        open(outunit,file=filename,action='write',form=fm)
        exit
     else
        outunit=outunit-1
     end if
  end do

  select case(abs(output_type))
  case(1) ! standard packed format (TP)
     call dtfttp('N','U',n,a,aout,err)
  case(2) ! rectangular full packed format (TF)
     aout = a
  case default
     write(emp_e,*) 'Error in output_mat_tri, unsupported output_type'
     write(emp_e,*) 'output_type=',output_type,'. Stopping'
     stop '-4'
  end select


  if(output_type .gt. 0) then
     write(outunit) aout
  else
     write(outunit,*) aout
  end if

  close(outunit)


end subroutine output_mat_tri