//  (c) Drobotenko    http://drobotenko.com
namespace La
{    inline static  void  inversCovarianceMtrix(DMat m)
     {
       int N= dimA(m);
       int *rn=(int*)alloca((sizeof(int))*N);
       double *str= (double*)alloca((sizeof(double))*N)
             ,*strm=(double*)alloca((sizeof(double))*N);

       int j,k;

       for(j=N;j--;)
         rn[j]=j;
       gaus_minved=1e99;
       gaus_deter=1;

       for(gaus_ostatok=N;gaus_ostatok;gaus_ostatok--)
       {  int jved;
          double aved=-1;
          for(j=N;j--;)
          if(~rn[j])
          if(aved<fabs(m[j][j]))
                aved=fabs(m[j][j]),jved=j;
          rn[jved]=-1;
          double mved=m[jved][jved];
          if(gaus_minved>aved)
             gaus_minved=aved;
          gaus_deter*=mved;

          if(aved<1e-99)
          {  for(j=N;j--;)
             {  if(~rn[j]) for(k=N;k--;)
                  m[j][k]=m[k][j]=_NaN();
             }
             break;
          }

          for(j=0;j<jved;j++)
            strm[j]=(str[j]=m[j][jved])/mved;
          for(j=jved;j<N;j++)
            strm[j]=(str[j]=m[jved][j])/mved;

          #if 0
          for(j=0;j<N;j++)
          for(k=j;k<N;k++)
            m[j][k]-= strm[j]*str[k];
          #else
          for(j=0;j<N;j++)
          {  double *mj= &m[j][0];
             for(k=j;k<N;k++)
               mj[k]-= strm[j]*str[k];
               //  оптимизация  m[j][k] заменили на mj[k]

          }
          #endif

          for(j=0;j<jved;j++)
            m[j][jved]= strm[j];
          for(j=jved;j<N;j++)
            m[jved][j]= strm[j];
          m[jved][jved]=-1./mved;

       }
       for(j=0;j<N;m[j][j]=-m[j][j],j++)
       for(k=j+1;k<N;k++)
           m[k][j]= m[j][k]= -m[j][k];

     }
}