double **init_a(int n) {
      double **a=new double*[n];
      for (int i=0;i<n;i++)
          a[i]=new double[n];
      for (int i=0;i<n;i++) {
          for (int j=0;j<n;j++) {
              a[i][j]=(i==j)?1:0;
          }
       }
       return a;
  }            
