T_RealArr = array of real;
T_Mul_s_s = class   
  private    
  fSrc1: T_RealArr;          
  fSrc2: T_RealArr;     
  fRes: T_RealArr;     
  procedure exec_SSE(const i0,i1: integer);   
end;
…
// IMPLEMENTATION
const FloatSize = 8;
procedure T_Mul_s_s.exec_SSE(const i0,i1: integer); 
{$IFDEF CPUX64}
asm   
// Self = RCX   
  .PUSHNV R13
  .PUSHNV R14
  .PUSHNV R15
  XOR R15,R15
  XOR R14,R14
  XOR R13,R13
  mov R13, Self.fRes      // R13 = @Res[0]
  mov R14, Self.fSrc1     // R14 = @Src1[0]
  mov R15, self.fSrc2     // R15 = @Src2[0]
  // расчет смещения Elem[i0] относительно начала массива
  XOR RAX, RAX            // очистили RAX
  mov EAX, i0             // i0- int32, потому кладем ее в EAX, для RAX она маловата
  imul RAX, FloatSize     // а умножать можно весь RAX
  // начальные адреса блоков данных, которые будем обрабатывать.
  add R13, RAX           // R13 = @Res[i0]
  add R14, RAX           // R14 = @Src1[i0]
  add R15, RAX           // R15 = @Src2[i0]
  // а теперь закинем куда-нибудь число оставшихся элементов:
  mov EAX, i1
  sub EAX, i0
  inc EAX // EAX = 1+i1-i0
// to be continued…

// continue…
 @new_pair:
    movapd XMM0, oWORD[R14]    // XMM0 = src1[i]   | src1[i+1]
    movapd XMM1, oWORD[R15]    // XMM1 = src2[i]   | src2[i+1]
    mulpd  XMM0, XMM1          // XMM0 = src1[i+0]src2[i+0] | src1[i+1]src2[i+1]
    movapd oWORD[R13], XMM0    // save to Res[i], Res[i+1]
    add R15, FloatSize*2  // next src2 pair
    add R14, FloatSize*2  // next src1 pair
    add R13, FloatSize*2  // next res pair
    sub RAX, 2
    cmp RAX, 1
  JA @new_pair
  cmp RAX, 0
  JE @finish
  @unaligned_end:
      movlpd XMM0, qWORD[R14]    // XMM0 = src1[i]   | garbage0
      movlpd XMM1, qWORD[R15]    // XMM1 = src2[i]   | garbage1
      mulsd  XMM0, XMM1          // XMM0 = src1[i]* src2[i]  | garbage0
      movlpd qWORD[R13], XMM0
      add R15, FloatSize         // move to the next element
      add R14, FloatSize
      add R13, FloatSize
      dec RAX
      cmp RAX,0
  JA @unaligned_end // if RAX > 0- jump to next elem
@finish:
{$ELSE CPUX64} begin {$ENDIF CPUX64}
end;

T_Vect = record    
  x: real;     
  y: real;     
  z: real; 
  {$IFDEF ALIGNED_VECTORS_ON}
  t: real;
  {$ENDIF }
end;
T_VectArr = array of T_Vect;
{$IFDEF ALIGNED_VECTORS_ON}
const VectSize = FloatSize * 4;
{$ELSE}
const VectSize = FloatSize * 3;
{$ENDIF }
 T_Mul_v_p = class
  private
    fCoeff: mreal;
    fSrc: T_VectArr;
    procedure exec_SSE(const i0,i1: integer);
  end;
  procedure T_Mul_v_p.exec_SSE_fast(const i0,i1: integer);
{$IFDEF CPUX64}
  asm
  // Self = RCX
  /// Result = Source!
     .NOFRAME
     .PUSHNV R14
      XOR RAX, RAX
      mov EAX, i0
      imul RAX, VectSize
      XOR R14, R14
      mov R14, Self.fSrc      // R14 = @Src[0]
      add R14, RAX            // R14 = @Src[i0]
      mov EAX, i1
      sub EAX, i0
      movlpd  XMM2, Self.fCoeff // XMM2 = fCoeff | x
      movlhps XMM2, XMM2        // продублировали в верхнюю половину
  @new_vect:
        movapd XMM0, oWORD[R14]                // XMM0 = v[i].x   | v[i].y
        movapd XMM1, oWORD[R14+2 * FloatSize]  // XMM0 = v[i].z   | v[i].t
        mulpd  XMM0, XMM2
        mulsd  XMM1, XMM2                      // только нижняя половина!
        movapd oWORD[R14              ], XMM0
        {$IFDEF ALIGNED_VECTORS_ON}
        movapd oWORD[R14+2 * FloatSize], XMM1
        {$ELSE}
        movlpd qWORD[R14+2 * FloatSize], XMM1
        {$ENDIF }
        add R14, VectSize  // next vector
        dec RAX
        cmp RAX, 0
  JGE @new_vect // Jump if RAX >=0
{$ENDIF CPUX64}
end;

T_Tens = packed record
   x: T_Vect;
   y: T_Vect;
   z: T_Vect;
end;

class procedure T_SSE.Add(var a_i: T_Vect; const B_ij: T_Tens; const c_j: T_Vect)
; static; register;
{$IFDEF CPUX64} 
asm 
{   a_i.x := a_i.x+b_ij.x.x * C_j.x+b_ij.x.y * C_j.y+b_ij.x.z * C_j.z ;   
    a_i.y := a_i.y+b_ij.y.x * C_j.x+b_ij.y.y * C_j.y+b_ij.y.z * C_j.z ;   
    a_i.z := a_i.z+b_ij.z.x * C_j.x+b_ij.z.y * C_j.y+b_ij.z.z * C_j.z ;}   
    // C_j: RAX,   
    // b_ij:  RDX, r8   (? compiler dependent! do not use)   
    // a_i: RCX
.NOFRAME
.SAVENV XMM5
.SAVENV XMM6
.SAVENV XMM7
movupd  XMM6, oWORD[C_j+0 * FloatSize]      // XMM6 = c0 | c1
movsd   XMM7, qWORD[C_j+2 * FloatSize]      // XMM7 = c2 | w
movupd  XMM4, oWORD[a_i+0 * FloatSize]      // XMM4 = a0 | a1
movsd   XMM5, qWORD[a_i+2 * FloatSize]      // XMM5 = a2 | w = waste
movupd  XMM0, oWORD[ b_ij+0 * FloatSize ]   //  XMM0 = m00 | m01
movsd   XMM1, qWORD[ b_ij+2 * FloatSize ]   //  XMM1 = m02 | w
movupd  XMM2, oWORD[ b_ij+4 * FloatSize ]   //  XMM2 = m10 | m11
movsd   XMM3, qWORD[ b_ij+6 * FloatSize ]   //  XMM3 = m12 | w
mulpd   XMM0, XMM6            // xmm0 = m00c0 | m01c1
mulsd   XMM1, XMM7            // xmm1 = m02c2 | 0
haddpd  XMM0, XMM1            // xmm0 = m00c0 + m01c1 | m02c2
mulpd   XMM2, XMM6            // xmm2 = m10c0 | m11c1
mulsd   XMM3, XMM7            // xmm3 = m22c2 | 0
haddpd  XMM2, XMM3            // xmm2 = m10c0 + m11c1 | m12c2
haddpd  XMM0, XMM2            // xmm0 = m00c0 + m01c1 + m02c2 | m10c0 + m11c1 + m12c2
addpd   XMM4, XMM0            // xmm4 = a0 + m0jcj | a1 + m1jcj
movupd  XMM2, oWORD[ b_ij+8 * FloatSize  ]   //  XMM2 = m20 | m21
movsd   XMM3, qWORD[ b_ij+10 * FloatSize  ]  //  XMM3 = m22 | w
mulpd   XMM2, XMM6            // XMM2 = m20c0 | m21c1
mulsd   XMM3, XMM7            // XMM3 = m22c2 | 0
movlhps XMM5, XMM3            // XMM5 = a2, m22c2
haddpd  XMM5, XMM2            // xmm3 = m22c2+a2 | m20c0 + m21c1
haddpd  XMM5, XMM5            // xmm5 = m22c2+a2 + m20c0 + m21c1 | m22c2+a2 + m20c0 + m21c1
movupd  oWORD[ a_i ], XMM4
movsd   qWORD[ a_i+2 * FloatSize ], XMM5
{$ELSE} 
begin   
a_i.x := a_i.x+b_ij.x.x * C_j.x+b_ij.x.y * C_j.y+b_ij.x.z * C_j.z;
a_i.y := a_i.y+b_ij.y.x * C_j.x+b_ij.y.y * C_j.y+b_ij.y.z * C_j.z;
a_i.z := a_i.z+b_ij.z.x * C_j.x+b_ij.z.y * C_j.y+b_ij.z.z * C_j.z; 
{$IFEND CPUX64}
end;

class procedure T_SSE.Add(var a_j: T_Vect; const c_i: T_Vect; const B_ij:
    T_Tens); static; register;
 {$IFDEF CPUX64}
asm
{  a_i.x := a_i.x+b_ij.x.xC_j.x+b_ij.x.yC_j.y+b_ij.x.zC_j.z ;
  a_i.y := a_i.y+b_ij.y.xC_j.x+b_ij.y.yC_j.y+b_ij.y.zC_j.z ;
  a_i.z := a_i.z+b_ij.z.xC_j.x+b_ij.z.yC_j.y+b_ij.z.zC_j.z ;}
  // C_i: RDX,
  // b_ij:  r8   (? compiler dependent! do not use)
  // a_j: RCX
    .NOFRAME
    .SAVENV XMM5
    .SAVENV XMM6
    .SAVENV XMM7
    movupd  XMM4, oWORD[ a_j+0 * FloatSize  ]   // XMM4 = a0 | a1
    movsd   XMM5, qWORD[ a_j+2 * FloatSize  ]   // XMM5 = a2 | w
   movapd  XMM0, oWORD[ b_ij  ]     //  XMM0 = m00 | m01
    ADD b_ij, 2 * FloatSize            
    movsd   XMM1, qWORD[ b_ij  ]     //  XMM1 = m02 | w
    ADD b_ij, 2 * FloatSize            
    movsd   XMM7, qWORD[ C_i  ]      //  XMM7 = c0 | 0
    ADD     C_i, FloatSize           
    movlhps XMM7, XMM7               //  XMM7 = c0 | c0
    mulpd   XMM0, XMM7               //  xmm0 = m00c0 | m01c0
    mulsd   XMM1, XMM7               //  xmm1 = m02c0 | 0
    addpd   XMM4, XMM0               //  a0+c0m00 | a1 + c0m01
    addsd   XMM5, XMM1               //  a2+c0m02 | 0
    movapd  XMM0, oWORD[ b_ij  ]     //  XMM0 = m10 | m11
    ADD b_ij, 2 * FloatSize            
    movsd   XMM1, qWORD[ b_ij  ]     //  XMM1 = m12 | w
    ADD b_ij, 2 * FloatSize            
    movsd   XMM7, qWORD[ C_i    ]    // XMM7 = c1 | 0
    ADD     C_i, FloatSize           //  RDX = C_i + 2 elem
    movlhps XMM7, XMM7               // XMM7 = c1 | c1
    mulpd   XMM0, XMM7               // xmm0 = m10c1 | m11c1
    mulsd   XMM1, XMM7               // xmm1 = m12c1 | 0
    addpd   XMM4, XMM0               //  a0+c0m00+m10c1 | a1 + c0m01+m11c1
    addsd   XMM5, XMM1               //  a2+c0m02+m12c1 | 0
    movapd  XMM0, oWORD[ b_ij  ]     //  XMM0 = m20 | m21
    ADD b_ij, 2 * FloatSize            
    movsd   XMM1, qWORD[ b_ij  ]     //  XMM1 = m22 | w
    movsd   XMM7, qWORD[ C_i   ]     // XMM7 = c2 | 0
    movlhps XMM7, XMM7               // XMM7 = c2 | c2
    mulpd   XMM0, XMM7               // xmm0 = m20c2 | m21c2
    mulsd   XMM1, XMM7               // xmm1 = m22c2 | 0
    addpd   XMM4, XMM0               //  a0+c0m00+m10c1+m20c2 | a1 + c0m01+m11c1+m21c2
    addsd   XMM5, XMM1               //  a2+c0m02+m12c1+m22c2 | 0
    movupd  oWORD[a_j   ], XMM4
    movsd   qWORD[a_j+2 * FloatSize], XMM5
{$ELSE}
begin
  a_j.x := a_j.x+b_ij.x.x * C_i.x+b_ij.y.x * C_i.y+b_ij.z.x * C_i.z ;
  a_j.y := a_j.y+b_ij.x.y * C_i.x+b_ij.y.y * C_i.y+b_ij.z.y * C_i.z ;
  a_j.z := a_j.z+b_ij.x.z * C_i.x+b_ij.y.z * C_i.y+b_ij.z.z * C_i.z ;
{$ENDIF CPUX64}
end;

class procedure T_SSE.Invert(const S: T_Tens; const i0, i1: integer);
{$IFDEF CPUX64}
  const One: double = 1.0;
        Zero: double = 0.0;
        eps: double = 1.0e-10;
        sign: Uint64 = $7FFFFFFFFFFFFFFF;
  asm
    .NOFRAME
    .SAVENV XMM5   // S.z.z
    .SAVENV XMM6   // Res.x.x | Res.x.y
    .SAVENV XMM7   // Res.x.z | w
    .SAVENV XMM8   // Res.y.x | Res.y.y
    .SAVENV XMM9   // Res.y.z | w
    .SAVENV XMM10  // Res.z.x | Res.z.y
    .SAVENV XMM11  // Res.z.z | w
    .SAVENV XMM12  // tmp1
    .SAVENV XMM13  // tmp2
    XOR RAX, RAX
    mov EAX, i0
    imul RAX, 3*VectorSize // Size of Tensor
    ADD RCX, RAX // RCX = @S[i0]
    XOR RAX, RAX
    mov EAX, i1
    sub EAX, i0
    inc RAX
    PXOR XMM6, XMM6
    PXOR XMM7, XMM7
    PXOR XMM8, XMM8
    PXOR XMM9, XMM9
    PXOR XMM10, XMM10
    PXOR XMM11, XMM11
@begin:
    prefetch [RCX+ 12*FloatSize  ]     // предзагрузка следующего тензора в кэш
    prefetch [RCX+ 16*FloatSize  ]
    prefetch [RCX+ 20*FloatSize  ]
    movupd XMM0, [RCX              ]   //   S.x.x | S.x.y
    movupd XMM1, [RCX+  2*FloatSize]   //   S.x.z | w
    movupd XMM2, [RCX+  4*FloatSize]   //   S.y.x | S.y.y
    movupd XMM3, [RCX+  6*FloatSize]   //   S.y.z | w
    movupd XMM4, [RCX+  8*FloatSize]   //   S.z.x | S.z.y
    movupd XMM5, [RCX+ 10*FloatSize]   //   S.z.z | w
  // Рассчет алгебраических дополнений (с учетом знака перестановки)
  {   Res.x.x :=    ( + S.y.y*S.z.z - S.y.z*S.z.y  );
      Res.x.y :=    ( + S.x.z*S.z.y - S.x.y*S.z.z  );          }
    movhlps XMM6, XMM2     // XMM6.L = S.y.y
    movlhps XMM6, XMM3     // XMM6.h = S.y.z
    blendpd XMM5, XMM4, 2  // XMM5= S.z.z | S.z.y
    mulpd XMM6, XMM5       // XMM6 = S.y.yS.z.z |  S.y.zS.z.y
    movhlps XMM12, XMM4    // XMM12.L = S.z.y
    movlhps XMM12, XMM5    // XMM12.h = S.z.z
    blendpd XMM1,  XMM0, 2 // XMM1 = S.x.z | S.x.y
    mulpd XMM12, XMM1      // XMM13 = S.x.z*S.z.y | S.x.y*S.z.z
    hsubpd XMM6, XMM12     // XMM6 = Res.x.x | Res.x.y     }
  {  Res.y.x :=    ( + S.z.x*S.y.z - S.y.x*S.z.z  );
     Res.y.y :=    ( + S.x.x*S.z.z - S.x.z*S.z.x  );   }
    movsd   XMM8, XMM4     // XMM8.L = S.z.x
    movlhps XMM8, XMM2     // XMM8.h = S.y.x
    movlhps XMM3 , XMM5    // XMM3 = S.y.z | S.z.z
    mulpd   XMM8, XMM3     // XMM8 = S.z.x*S.y.z |  S.y.x*S.z.z
    movsd   XMM13, XMM0    // XMM12.L = S.x.x
    movlhps XMM13, XMM1    // XMM13 = S.x.x | S.x.z
    movlhps XMM5, XMM4     // XMM5 = S.z.z  | S.z.x
    mulpd XMM13, XMM5      // XMM13 = S.x.x*S.z.z | S.x.z*S.z.x
    hsubpd XMM8, XMM13     // XMM8 = Res.y.x | Res.y.y
  {   Res.z.x :=    ( + S.y.x*S.z.y - S.y.y*S.z.x  );
      Res.z.y :=    ( + S.x.y*S.z.x - S.x.x*S.z.y  );   }
    movhlps XMM10, XMM4
    movlhps XMM10, XMM4    // XMM10 = S.z.y | S.z.x
    pxor XMM12, XMM12
    subpd  XMM12, XMM10    // XMM10 = -S.z.y | -S.z.x
    mulpd  XMM12, XMM0     // XMM12 = S.z.y*S.y.x | S.z.x*S.y.y
    mulpd  XMM10, XMM2     // XMM10 = -S.x.x*S.z.y | -S.x.y*S.z.x
    hsubpd XMM10, XMM12    // XMM10 = Res.z.x | Res.z.y
  {  Res.x.z :=    ( + S.x.y*S.y.z - S.x.z*S.y.y  ); //  XMM7.L
     Res.y.z :=    ( + S.y.x*S.x.z - S.x.x*S.y.z  ); // XMM9.L
    // Reorder:
    Res.x.z :=    ( + S.x.y*S.y.z - S.y.y*S.x.z  ); //  XMM7.L
    Res.y.z :=   -(   S.x.x*S.y.z - S.y.x*S.x.z  ); // XMM9.L       }
    movlhps XMM3, XMM3     // XMM3 = S.y.z | S.y.z
    mulpd   XMM3, XMM0     // XMM3 = S.x.x*S.y.z | S.x.y*S.y.z
    movapd  XMM12, XMM2    // XMM12 = S.y.x | S.y.y
    movlhps XMM1, XMM1     // XMM1  = S.x.z | S.x.z
    mulpd   XMM12, XMM1    // XMM12 = S.y.x*S.x.z | S.y.y*S.x.z
    subpd   XMM3, XMM12    // XMM3 = S.x.x*S.y.z-S.y.x*S.x.z | S.x.y*S.y.z - S.y.y*S.x.z
                           //      = -Res.y.z                | Res.x.z
    pxor  XMM9, XMM9       // XMM9 = 0 | 0
    subsd XMM9, XMM3       // XMM9 = Res.y.z
    movhlps XMM7, XMM3     // XMM11 = Res.x.z
  {  Res.z.z :=    ( + S.x.x*S.y.y - S.x.y*S.y.x  );     }
    movhlps XMM11, XMM2    // XMM10 = S.y.y | w
    movlhps XMM11, XMM2    // XMM10 = S.y.y | S.y.x
    mulpd   XMM11, XMM0    // XMM11 = S.x.x*S.y.y | S.x.y*S.y.x
    hsubpd  XMM11, XMM11   // XMM11 = S.x.x*S.y.y - S.x.y*S.y.x
  {   D := S.x.x*Res.x.x +   S.x.y*Res.y.x +   S.x.z*Res.z.x ==
           Res.x.x*S.x.x +   Res.x.y*S.y.x +   Res.x.z*S.z.x         }
    movlhps XMM0, XMM2
    dppd    XMM0, XMM6,  49 // dppd- dot product packed double
                            // imm8 = bit0 | ~bit1 | bit4 | bit5
                            // bit0- save res to low half
                            // bit1- save res to high half
                            // bit4- * low parts
                            // bit5- * high parts
                            // XXM0 = (s.xx | s.yx)*(Res.xx | Res.xy)
    mulsd XMM4, XMM7        // XMM4 =  S.zx * Res.xz
    addsd XMM4, XMM0        // now XMM4 = D = Determinant(S)!
    // Check if D< eps:
    movsd XMM0, XMM4
    movsd XMM1, sign
    pand  XMM0, XMM1        // убрали знак D- получили XMM0 = Abs(D)
    movsd XMM1, eps
    comisd XMM1, XMM0       // if ABS (D) > eps then PROCEED else EXIT
    jnbe @exit
  //  Res now is the Matrix of algebraic complements,
  //  Out := mult(Res,1.0/D);
    movsd   XMM3, One
    divsd   XMM3, XMM4
    movlhps XMM3, XMM3                       // XMM3 = 1.0/D | 1.0/D
    mulpd   XMM6, XMM3
    movupd  oWORD[RCX              ], XMM6   //   S.x.x | S.x.y
    mulpd   XMM7, XMM3
    movupd  oWORD[RCX+  2*FloatSize], XMM7   //   S.x.z |
    mulpd   XMM8, XMM3
    movupd  oWORD[RCX+  4*FloatSize], XMM8   //   S.y.x | S.y.y
    mulpd   XMM9, XMM3
    movupd  oWORD[RCX+  6*FloatSize], XMM9   //   S.y.z |
    mulpd   XMM10, XMM3
    movupd  oWORD[RCX+  8*FloatSize], XMM10  //   S.z.x | S.z.y
    mulpd   XMM11, XMM3
    movupd  oWORD[RCX+ 10*FloatSize], XMM11  //   S.z.z |
  //  sfence
    @exit:
    Add RCX, 3*VectorSize  // next matrix
    dec RAX
    cmp RAX, 0
    jg @begin              // if RAX > 0- invert next matrix
{$ELSE}
var i: integer;
begin
  for i := i0 to i1 do
    Invert(  T_TensorArr(@S)[i] );
{$IFEND}
end;

Выбираем строку с номером i=1.
сохраняем D.ii с главной диагонали этой строки, например- “xx”, 
на его место записываем 1.0. 
вычисляем 1/D, 
всю строку с этим элементом умножаем на 1/D. 
  выбираем j-ю строку j=i+1
  выбираем в ней i-й элемент (‘yx’), копируем его.
  вместо него записываем 0. 
  вычисляем 1/’yx’
  умножаем на это значение всю j-ю строку, 
  вычитаем из строки j строку i. 
  переходим к j+1
  если строки кончились- увеличиваем i=i+1, и повторяем. 
когда обработаем так все строки- аналогичным образом выполняем обратный ход метода Гаусса

type
  T_M4 = array [0..3,0..3] of real;
class procedure T_SSE.Invert_gauss(var S: T_M4);
{$IFDEF CPUX64}
const FloatSize = 8;
      One: double = 1.0;
      Zero: double = 0.0;
asm
  .NOFRAME
  .SAVENV XMM5
  .SAVENV XMM6
  .SAVENV XMM7
  .SAVENV XMM8
  .SAVENV XMM9
  .SAVENV XMM10
  movupd XMM0, [RCX              ]   //   S.x.x | S.x.y
  movupd XMM1, [RCX+  2*FloatSize]   //   S.x.z | S.x.t
  prefetch [RCX + 4*FloatSize]       // slightly increase overall performance (~1%)
  PXOR XMM10,XMM10
//  FORWARD MOVE
// STEP  #1
// STEP #1.prep:  D = 1.0/ S.x.x    S.x.x = 1; S.x := S.xD;
  movsd XMM8, One             // XMM8.L = 1.0;
  movsd XMM9, XMM0            // EXTRACT DIVIDER!
  movsd XMM0, XMM8            //XMM0 = 1 | S.x.y
  divsd XMM8, XMM9            // D := 1.0 / S.x.x;
  movlhps XMM8, XMM8          // XMM8 = D  | D
  mulpd XMM0, XMM8            // S.x = S.xD
  mulpd XMM1, XMM8            // XMM0 XMM1= a00d | a01d   ||  a02d | a03d
  movupd XMM2, [RCX+  4*FloatSize]   //   S.y.x | S.y.y
  movupd XMM3, [RCX+  6*FloatSize]   //   S.y.z | S.y.t
  prefetch [RCX + 8*FloatSize]
// STEP  #1.y:   D := S.y.x;  s.y.x :=  0; S.y = S.y- S.xD
  movsd   XMM8, XMM2   // XMM8 = (D= S.y.x) | 0
  movlhps XMM8, XMM8   // XMM8 =  D          |   D
  movsd XMM2, XMM10    // S.y.x := 0;  XMM2 = 0 | S.y.y
  movapd XMM9, XMM0    // copy S.x.x, S.x.y
  mulpd XMM9, XMM8     // mul
  subpd XMM2, XMM9     // XMM2 = -S.x.xD         | S.y.y - S.x.yD
  movapd XMM9, XMM1    // copy S.x.z, S.x.t
  mulpd XMM9, XMM8     // mul
  subpd XMM3, XMM9     // XMM3 = S.y.z - S.x.zd |  S.y.t - S.x.td        //}
// STEP  #1.z:   D := S.z.x;  s.z.x :=  0; S.z = S.z- S.xD
  movupd XMM4, [RCX+  8*FloatSize]   //   S.z.x | S.z.y
  movupd XMM5, [RCX+ 10*FloatSize]   //   S.z.z | S.z.t
  prefetch [RCX + 12*FloatSize]
  movsd   XMM8, XMM4   // XMM8 = (D= S.y.x) | 0
  movlhps XMM8, XMM8   // XMM8 =  D          |   D
  movsd XMM4, XMM10    // S.z.x:= 0;  XMM4 = 0 | S.z.y
  movapd XMM9, XMM0    // Copy S.x.x, S.x.y
  mulpd XMM9, XMM8     // mul
  subpd XMM4, XMM9     // XMM4 = -S.x.xD         | S.y.y - S.x.yD
  movapd XMM9, XMM1    // copy S.x.z , S.x.t
  mulpd XMM9, XMM8     // mul
  subpd XMM5, XMM9     // XMM5 = S.z.z - S.z.zd |  S.z.t - S.z.td        //}
// STEP  #1.t:   D := S.t.x;  s.t.x :=  0; S.t = S.t- S.xD
  movupd XMM6, [RCX+ 12*FloatSize]   //   S.t.x | S.t.y
  movupd XMM7, [RCX+ 14*FloatSize]   //   S.t.z | S.t.t
  movsd   XMM8, XMM6   // XMM8 = (D= S.y.x) | 0
  movlhps XMM8, XMM8   // XMM8 =  D          |   D
  movsd XMM6, XMM10    // S.t.x:= 0;  XMM6 = 0 | S.z.y
  movapd XMM9, XMM0    // Copy S.y.x, S.y.y
  mulpd XMM9, XMM8     // mul
  subpd XMM6, XMM9     // XMM6 = -S.x.xD         | S.y.y - S.x.yD
  movapd XMM9, XMM1    // copy S.z.z , S.z.t
  mulpd XMM9, XMM8     // mul
  subpd XMM7, XMM9     // XMM7 = S.z.z - S.z.zd |  S.z.t - S.z.td        //}
// STEP  #4.prep  D := 1.0/S.y.y;    S.y.y = 1;  S.y := S.yD;
  movsd   XMM8, One
  movhlps XMM9, XMM2    // XMM9 = S.y.y    |    w
  movlhps XMM2, XMM8    // XMM2 = S.y.x    |    S.y.y=1
  divsd   XMM8, XMM9    // XMM8 = 1.0/ S.y.y   = D
  movlhps XMM8, XMM8
  mulpd   XMM2, XMM8    // XMM2 = S.y.xD  |    D= 1/S.y.y
  mulpd   XMM3, XMM8    // XMM3 = S.y.zD  |    S.y.tD
// STEP  #4.z  D := S.z.y; S.z.y := 0; S.z := S.z-S.yD
  movhlps XMM8, XMM4   // XMM9= D=S.z.y   |   w
  movlhps XMM8, XMM8   // XMM8= D         |   D
  movlhps XMM4, XMM10  // XMM4= S.z.x     | S.z.y=0
  movapd  XMM9, XMM2   // XMM9= S.y.x      | S.y.y
  mulpd   XMM9, XMM8   // XMM9= S.y.xD    | S.y.yD
  subpd   XMM4, XMM9   // XMM4= S.z.x - S.y.xD | 0 -S.y.yD
  movapd  XMM9, XMM3   // XMM9= S.y.z      | S.y.t
  mulpd   XMM9, XMM8   // XMM9= S.y.zD    | S.y.tD
  subpd   XMM5, XMM9   // XMM4= S.z.z - S.y.zD | S.z.t - S.y.tD
// STEP  #4.t  D := S.t.y; S.t.y := 0; S.t := S.t-S.yD
  movhlps XMM8, XMM6   // XMM9= D=S.t.y   |   w
  movlhps XMM8, XMM8   // XMM8= D         |   D
  movlhps XMM6, XMM10  // XMM4= S.t.x     | S.t.y=0
  movapd  XMM9, XMM2   // XMM9= S.y.x      | S.y.y
  mulpd   XMM9, XMM8   // XMM9= S.y.xD    | S.y.yD
  subpd   XMM6, XMM9   // XMM4= S.t.x - S.y.xD | 0 -S.y.yD
  movapd  XMM9, XMM3   // XMM9= S.y.z      | S.y.t
  mulpd   XMM9, XMM8   // XMM9= S.y.zD    | S.y.tD
  subpd   XMM7, XMM9   // XMM4= S.t.z - S.y.zD | S.t.t - S.y.tD
// STEP  #5.prep  D := 1.0/S.z.z;    S.z.z = 1;  S.z := S.zD;
  movsd   XMM8, One
  movsd   XMM9, XMM5    // XMM9 = S.z.z    |    w
  movsd   XMM5, XMM8    // XMM5 = S.z.z=1  |    S.z.t
  divsd   XMM8, XMM9    // XMM8 = 1.0/ S.z.z   = D
  movlhps XMM8, XMM8
  mulpd   XMM4, XMM8    // XMM4 = S.z.xD  |    D= 1/S.z.z
  mulpd   XMM5, XMM8    // XMM5 = S.z.zD  |    S.z.tD
// STEP  #5.t  D := S.t.z; S.t.z := 0; S.t := S.t-S.zD
  movsd   XMM8, XMM7   // XMM8= D=S.t.z   |   w
  movlhps XMM8, XMM8   // XMM8= D         |   D
  movsd   XMM7, XMM10  // XMM7= S.t.z=0   | S.t.t
  movapd  XMM9, XMM4   // XMM9= S.z.x      | S.z.y
  mulpd   XMM9, XMM8   // XMM9= S.z.xD    | S.z.yD
  subpd   XMM6, XMM9   // XMM4= S.t.x - S.z.xD | S.t.y -S.z.yD
  movapd  XMM9, XMM5   // XMM9= S.z.z      | S.z.t
  mulpd   XMM9, XMM8   // XMM9= S.z.zD    | S.z.tD
  subpd   XMM7, XMM9   // XMM7= 0 - S.z.zD | S.t.t - S.z.tD
// STEP  #6  D := 1/S.t.t; S.t.t := 1.0; S.t := S.tD
  movsd   XMM8, One
  movhlps XMM9, XMM7    // XMM9 = S.t.t    |    w
  movlhps XMM7, XMM8    // XMM7 = S.t.z    |    S.t.t=1
  divsd   XMM8, XMM9    // XMM8 = 1.0/ S.z.z   = D
  movlhps XMM8, XMM8
  mulpd XMM6, XMM8
  mulpd XMM7, XMM8
  // save S.t
  movupd oWORD[RCX+ 12*FloatSize], XMM6   //   S.z.x | S.z.y
  movupd oWORD[RCX+ 14*FloatSize], XMM7   //   S.z.z |
// BACKWARD MOVEMENT!
// STEP  #7.x-y  D := S.x.y; S.x.y := 0; S.x -= S.yD
  movhlps XMM8, XMM0     // XMM8= D=S.x.y | w
  movlhps XMM8, XMM8     // XMM8 = D | D
  movlhps XMM0, XMM10    // S.x.y = 0
  movapd XMM9, XMM2      // XMM9= S.y.x      | S.y.y
  mulpd XMM9, XMM8       // XMM9= S.y.xD    | S.y.yD
  subpd XMM0, XMM9       // XMM0= S.x.x+S.y.xD | S.y.yD
  movapd XMM9, XMM3       // XMM9= S.y.z      | t
  mulpd XMM9, XMM8       // XMM9= S.y.zD    | t
  subpd XMM1, XMM9       // XMM0= S.x.z+S.y.zD | t
// STEP  #7.y-z  D := S.y.z; S.y.z := 0; S.y -= S.zD
  movsd   XMM8, XMM3     // XMM8= D=S.y.z | w
  movsd   XMM3, XMM10    // S.y.z := 0
  movlhps XMM8, XMM8     // XMM8 = D | D
  movapd XMM9, XMM4      // XMM9= S.z.x      | S.z.y
  mulpd XMM9, XMM8       // XMM9= S.z.xD    | S.z.yD
  subpd XMM2, XMM9       // XMM0= S.y.x-S.z.xD | S.y.y-S.z.yD
  movapd XMM9, XMM5       // XMM9= S.y.z      | t
  mulpd XMM9, XMM8       // XMM9= S.y.zD    | t
  subpd XMM3, XMM9       // XMM0= S.x.z+S.y.zD | t
// STEP  #7.x-z  D := S.x.z; S.x.z := 0; S.x -= S.zD
  movsd   XMM8, XMM1     // XMM8= D=S.x.z | w
  movsd   XMM1, XMM10    // S.x.z := 0
  movlhps XMM8, XMM8     // XMM8 = D | D
  movapd XMM9, XMM4      // XMM9= S.t.x      | S.t.y
  mulpd XMM9, XMM8       // XMM9= S.t.xD    | S.t.yD
  subpd XMM0, XMM9       // XMM0= S.z.x-S.t.xD | S.z.y-S.t.yD
  movapd XMM9, XMM5      // XMM9= S.t.z      | S.t.t
  mulpd XMM9, XMM8       // XMM9= S.t.zD    | S.t.tD
  subpd XMM1, XMM9       // XMM0= S.z.z-S.t.zD | S.z.t-S.t.tD
// STEP  #7.z-t  D := S.z.t; S.z.t := 0; S.z -= S.tD
  movhlps XMM8, XMM5     // XMM8= D=S.z.t | w
  movlhps XMM5, XMM10    // S.z.t := 0
  movlhps XMM8, XMM8     // XMM8 = D | D
  movapd XMM9, XMM6      // XMM9= S.t.x      | S.t.y
  mulpd XMM9, XMM8       // XMM9= S.t.xD    | S.t.yD
  subpd XMM4, XMM9       // XMM0= S.z.x-S.t.xD | S.z.y-S.t.yD
  movapd XMM9, XMM7      // XMM9= S.t.z      | S.t.t
  mulpd XMM9, XMM8       // XMM9= S.t.zD    | S.t.tD
  subpd XMM5, XMM9       // XMM0= S.z.z-S.t.zD | S.z.t-S.t.tD
  // SAVE S.z
  movupd oWORD[RCX+  8*FloatSize], XMM4   //   S.z.x | S.z.y
  movupd oWORD[RCX+ 10*FloatSize], XMM5   //   S.z.z |
// STEP  #8.y-t  D := S.y.t; S.y.t := 0; S.y -= S.tD
  movhlps XMM8, XMM3     // XMM8= D=S.y.t | w
  movlhps XMM3, XMM10    // S.y.t = 0
  movlhps XMM8, XMM8     // XMM8 = D | D
  movapd XMM9, XMM6      // XMM9= S.y.x      | S.y.y
  mulpd XMM9, XMM8       // XMM9= S.y.xD    | S.y.yD
  subpd XMM2, XMM9       // XMM0= S.x.x+S.y.xD | S.y.yD
  movapd XMM9, XMM7      // XMM9= S.y.z      | t
  mulpd XMM9, XMM8       // XMM9= S.y.zD    | t
  subpd XMM3, XMM9       // XMM0= S.x.z+S.y.zD | t
//   save S.y
  movupd oWORD[RCX+  4*FloatSize], XMM2   //   S.y.x | S.y.y
  movupd oWORD[RCX+  6*FloatSize], XMM3   //   S.y.z |
// STEP  #8.x-t  D := S.x.t; S.x.t := 0; S.x -= S.tD
  movhlps XMM8, XMM1     // XMM8= D=S.x.t | w
  movlhps XMM1, XMM10    // S.x.t := 0
  movlhps XMM8, XMM8     // XMM8 = D | D
  movapd XMM9, XMM6      // XMM9= S.t.x      | S.t.y
  mulpd XMM9, XMM8       // XMM9= S.t.xD    | S.t.yD
  subpd XMM0, XMM9       // XMM0= S.x.x-S.t.xD | S.x.y-S.t.yD
  movapd XMM9, XMM7      // XMM9= S.t.z         | S.t.t
  mulpd XMM9, XMM8       // XMM9= S.t.zD       | S.t.tD
  subpd XMM1, XMM9       // XMM0= S.x.z-S.t.zD | S.x.t-S.t.tD
  movupd oWORD[RCX              ], XMM0   //   S.x.x | S.x.y
  movupd oWORD[RCX+  2*FloatSize], XMM1   //   S.x.z |
{$ELSE}
begin
  Invert_m4( S );
{$ENDIF}
end;

Dot product: compiler :  5240,8 MB/s
Dot product: SSE      :  7717,2 MB/s
Dot product: SSE, arr :  9704,8 MB/s

V3 += V3M33 compiler : 3599,5 MB/s
V3 += V3M33 SSE      : 6810,8 MB/s
V3 += M33V3 compiler : 4304,6 MB/s
V3 += M33V3 SSE      : 8644,2 MB/s

M3х3: Invert compiler    : 3384,8 MB/s
M3х3: T_SSE.Invert_gauss : 3386,8 MB/s
M3х3: T_SSE.Invert       : 4990,0 MB/s

M4х4: Invert compiler    :  333,2 MB/s
M4х4: Invert_SSE_gauss   :  2958,2 MB/s