flush to zero denormalized floating point numbers

[veselic]

Dear Martin

thanks indeed. But my program still does not flush.
Maybe you just kindly look at the whole stuff --
it is not much. The otput is still

denormalized =      1.175002775E-40

This will help me indeed, thanks again, Kresimir Veselic.

program axb;
const

  MM_MaskDAZMode   = %0000000001000000;
  MM_MaskInvalidOp = %0000000010000000;
  MM_MaskDenorm    = %0000000100000000;
  MM_MaskDivZero   = %0000001000000000;
  MM_MaskOverflow  = %0000010000000000;
  MM_MaskUnderflow = %0000100000000000;
  MM_MaskPrecision = %0001000000000000;
  MM_RoundingCtrl  = %0110000000000000;
  MM_MaskFTZ       = %1000000000000000;

  mxcsr : dword = MM_MaskUnderflow
               or MM_MaskPrecision
               or MM_MaskDenorm
//               or MM_MaskFTZ
               or MM_MaskDAZMode;
Procedure SysResetFPU;
var
  localmxcsr: dword;
begin
  if has_sse_support then
    begin
      localmxcsr:=mxcsr;
      asm
        ldmxcsr localmxcsr
      end;
    end;
end;
var a,c: single;
begin
SysResetFPU;
a := 1.175E-38;
c := a/100;
writeln('denormalized = ',c:20);
end.

[Nitorami]

Veselic

The compiler per default uses the x87 FPU which does not support FTZ or DAZ.
You need to tell it to use SSE. Luckily, this can be done by a compiler switch

{$FPUTYPE SSE2}

Insert this at the top of your code.

Regards
Martin

[veselic]

Dear Martin,

excellent. Thanks again. Maybe this trick could be of
some interest eleswhere.


Regards, K. Veselic

[engkin]

According to the "Guidelines for Writing to the MXCSR Register" (IA32 Vol. 1) setting a reserved bit in MXCSR causes general-protection exception. We are supposed to use MXCSR_MASK located at offset 28 from the 512 bytes fxsave gives (the buffer should be aligned to 16byte boundary), applying these guidelines on Nitorami's code:

Code: [Select]

program project1;

const

  MM_MaskDAZMode   = %0000000001000000;
  MM_MaskInvalidOp = %0000000010000000;
  MM_MaskDenorm    = %0000000100000000;
  MM_MaskDivZero   = %0000001000000000;
  MM_MaskOverflow  = %0000010000000000;
  MM_MaskUnderflow = %0000100000000000;
  MM_MaskPrecision = %0001000000000000;
  MM_RoundingCtrl  = %0110000000000000;
  MM_MaskFTZ       = %1000000000000000;

  mxcsr : dword = MM_MaskUnderflow
               or MM_MaskPrecision
               or MM_MaskDenorm
               or MM_MaskFTZ
               or MM_MaskDAZMode;


{$asmmode intel}

procedure fxsave(P: pointer);assembler;
asm
  fxsave [p];
end;

function get_mxcsr_mask: dword;
var
  buf: array[0..511+15] of byte;
  pAlignedBuf: PByte;
  pMXCSR_MASK: pdword;
begin
  pAlignedBuf := pointer((PtrUInt(@Buf[0])+15) and not $F);
  FillByte(pAlignedBuf^, 512, 0);
  fxsave(pAlignedBuf);
  pMXCSR_MASK := @pAlignedBuf[28];
  Result := pMXCSR_MASK^;
end;

function has_daz_support: boolean;
begin
  Result := (get_mxcsr_mask and MM_MaskDAZMode)<>0;
end;

function safe_mxcsr(vMXCSR: dword): dword;
var
  lMask: dword;
begin
  lMask := get_mxcsr_mask;
  if lMask = 0 then
   lMask := $FFBF; { value came from IA32 Vol. 1 "Guidelines for Writing to the MXCSR Register"}
  Result := vMXCSR and lMask;
end;

Procedure SysResetFPU;
var
  localmxcsr: dword;
begin
  if has_sse_support then
    begin
      localmxcsr := safe_mxcsr(mxcsr);
      {$asmmode intel}
      asm
        ldmxcsr localmxcsr
      end;
    end;
end;

{$FPUTYPE sse2}

var a,c: single;
begin
  SysResetFPU;
  a := 1.175E-38;
  c := a/100;
  writeln('denormalized = ',c:20);
  readln;
end.


[Nitorami]

I see.
Well, it could be guessed that writing to the FPU register willy-nilly may be risky, but it worked, so I did not bother.

BTW - can't we achieve the purpose of this awkward construct

Code: [Select]

  pAlignedBuf := pointer((PtrUInt(@Buf[0])+15) and not $F);

simply by

Code: [Select]

{$CODEALIGN LOCALMIN=16}

?

[engkin]

I see.
Well, it could be guessed that writing to the FPU register willy-nilly may be risky, but it worked, so I did not bother.

I hope you did not read my post as a criticism, because it is not. You went straight to the core of the code, which IMHO is impressive and should work fine on all new CPUs. I just tried to following the little mentioned steps to include older CPUs.

BTW - can't we achieve the purpose of this awkward construct

Code: [Select]

  pAlignedBuf := pointer((PtrUInt(@Buf[0])+15) and not $F);

simply by

Code: [Select]

{$CODEALIGN LOCALMIN=16}

?

I tried it first, but it did not work on my 32bit system and old 2.7.1 FPC.

[Nitorami]

I hope you did not read my post as a criticism,

No. If it came across like this, it is because english is not my native language.

Codealign:

I tried it first, but it did not work on my 32bit system and old 2.7.1 FPC.

Odd, it works here, FPC 2.6.4 on my AMD Atlon 64 under Win XP. With codealign, the starting address of buf is a multiple of 16; without, it is (mostly) not.

[veselic]

Itried to compile the modified program by engkin but got
errors:

flush_to_zero2.pas(39,10) Error: Identifier not found "Result"
flush_to_zero2.pas(44,10) Error: Identifier not found "Result"
flush_to_zero2.pas(54,10) Error: Identifier not found "Result"
flush_to_zero2.pas(81) Fatal: There were 3 errors compiling module, stopping

What was wrong? Thanks in advance, K. Veselic.

[Laksen]

You are not compiling in ObjFpc mode

[veselic]

Laksen:

According to your comment and the man-page instruction I tried to
compile as

 fpc -Mmode objfpc flush_to_zero2.pas

but got

Error: Illegal parameter: -Mmode

I apparently  made error in command line syntax...
What is wrong?

Thanks, K. Veselic.

[rvk]

I tried to compile the modified program by engkin but got errors:

You could just add

Code: [Select]

{$mode objfpc}
at the top of your program (just below your "program flush_to_zero2;" line.

Otherwise you need to do:

Code: [Select]

fpc -Mobjfpc flush_to_zero2.pas



See also the help in FPC:

Code: [Select]

  -M<x>  Set language mode to <x>
      -Mfpc      Free Pascal dialect (default)
      -Mobjfpc   FPC mode with Object Pascal support
      -Mdelphi   Delphi 7 compatibility mode
      -Mtp       TP/BP 7.0 compatibility mode
      -Mmacpas   Macintosh Pascal dialects compatibility mode


[veselic]

Laksen,

Yes. Thanks to all of you, also for using your time to help
ignorant matehmaticians. Just to explain myself:
the LAPACK routine (s,d)gesvj, which I coauthored
is based of my own extensive previous experimental work
in (Turbo)pascal.  With (s,d)gesvj flush-to-zero
mode was essential, and we had quite a nasty time to deliver
a (numerical) workaround for processors without this mode
on their disposal.  I am now resuming
the resaerch in this direction and find pascal quite
comfortable for this purpose.
Although in this phase the efficiancy is not an issue, it is good to
be able to experiment with flush-to-zero in order to observe
possible side effiects e.g. in accuracy.

So, maybe you could explain in two words, is there essential
difference between the two codes above: martin's and engkin's?

K. Veselic.

[Nitorami]

Not in mathematical terms.

Engkin's code follows the textbook, read the the FPU registers first, then modify what is allowed, then write back. While my code is quick and dirty. It might cause an exception on some FPUs. Not sure.

If in doubt, use Engkin's code. But the mathematical result will certainly be the same.

[veselic]

Yes, thanks once more that seems to settle it. K. Veselic.