The format statement for outputting the mass, stiffness and damping matrices results (subroutines psmass.f and pptang.f) in the numbers having data exponent of type "d". In loading these files in Matlab or Octave using "load TANG_001" for example results in truncating the exponent from the values (this is of course wrong). I found that replacing "d" with "e" in the format statement in these subroutines resolves this issue and makes their loading in Matlab correct.