namespace ngfem {

template <int M>
void FastMat (int n, double * pa, double * pb, double * pc)
{
  double * hpa = pa;

  for (int i = 0; i < n; i++)
    {
      double * hpb = pb;
      double * hpc = pc+n*i;

      for (int j = 0; j < i; j++)
	{
	  double sum = *hpc;
	  
	  for (int k = 0; k < M; k++)
	    sum += hpa[k] * hpb[k];
	  
	  *hpc = sum;
	  pc[i+n*j] = sum;

	  hpb += M;
	  hpc++;
	}

      double sum = *hpc;
      
      for (int k = 0; k < M; k++)
	sum += hpa[k] * hpb[k];
      
      *hpc = sum;

      hpa += M;
    }
}

template void FastMat<1> (int n, double * pa, double * pb, double * pc);
template void FastMat<2> (int n, double * pa, double * pb, double * pc);
template void FastMat<3> (int n, double * pa, double * pb, double * pc);
template void FastMat<4> (int n, double * pa, double * pb, double * pc);
template void FastMat<5> (int n, double * pa, double * pb, double * pc);
template void FastMat<6> (int n, double * pa, double * pb, double * pc);
template void FastMat<9> (int n, double * pa, double * pb, double * pc);



template void FastMat<25> (int n, double * pa, double * pb, double * pc);
template void FastMat<26> (int n, double * pa, double * pb, double * pc);
template void FastMat<27> (int n, double * pa, double * pb, double * pc);
template void FastMat<28> (int n, double * pa, double * pb, double * pc);
template void FastMat<29> (int n, double * pa, double * pb, double * pc);
template void FastMat<30> (int n, double * pa, double * pb, double * pc);


template void FastMat<12> (int n, double * pa, double * pb, double * pc);
template void FastMat<18> (int n, double * pa, double * pb, double * pc);
template void FastMat<24> (int n, double * pa, double * pb, double * pc);
template void FastMat<36> (int n, double * pa, double * pb, double * pc);
template void FastMat<48> (int n, double * pa, double * pb, double * pc);

template void FastMat<100> (int n, double * pa, double * pb, double * pc);






void FastMatN (int n, int M, double * pa, double * pb, double * pc)
{
  double * hpa = pa;

  for (int i = 0; i < n; i++)
    {
      double * hpb = pb;
      double * hpc = pc+n*i;

      for (int j = 0; j < i; j++)
	{
	  double sum = *hpc;
	  
	  for (int k = 0; k < M; k++)
	    sum += hpa[k] * hpb[k];
	  
	  *hpc = sum;
	  pc[i+n*j] = sum;

	  hpb += M;
	  hpc++;
	}

      double sum = *hpc;
      
      for (int k = 0; k < M; k++)
	sum += hpa[k] * hpb[k];
      
      *hpc = sum;

      hpa += M;
    }
}


}