double uquaf(int n, double *x, double *a, int lda)
#define A(I,J) a[(I) + (J)* ( lda)]
{
   int i, j;
   float aa=0.0, bb=0.0;
   
   for (i=0; i<=n-1; i++) 
      for (j=0; j<=n-1; j++) {
         bb=x[i]*A(j,i)*x[j];
         aa=aa+bb;
      }
   return(aa);
}