// pde_abc_eq.java  see WEB page for development, CMSC 455 supplemental  
//        The PDE to be solved is:                             
//        a1(x,y)*uxx(x,y) + b1(x,y)*uxy(x,y) + c1(x,y)*uyy(x,y) +    
//        d1(x,y)*ux(x,y) + e1(x,y)*uy(x,y) + f1(x,y)*u(x,y) = c(x,y) 
//                                                                    
// The functions and other user information is in the files:          
//    abc.java                                              
//    compile with abc.java and nderiv.java and simeq.java rderiv.java

 
// uses simeq.java
// uses deriv.java
import java.text.*;

class pde_abc_eq
{
  abc ABC = new abc();
  final int nd = ABC.nd; // degree of discretization
  final int ifcheck = ABC.ifcheck; // non zero for checking 
  final int nx = ABC.nx; // problem parameters
  final int ny = ABC.ny;
  final double xmin = ABC.xmin;
  final double xmax = ABC.xmax;
  final double ymin = ABC.ymin;
  final double ymax = ABC.ymax;
  double us[] = new double[(nx-2)*(ny-2)];  // solution being computed 
  double uu[] = new double[(nx-2)*(ny-2)];  // solution if checking 
  double ut[][] = new double[(nx-2)*(ny-2)][(nx-2)*(ny-2)+1];
                         // equations to be solved  
                         // nx-2 by ny-2  internal, non boundary cells 
  // all elements 0..nx-1 by 0..ny-1 includes boundary elements 
  // xmin, xmax, ymin, ymax are defined in  abc.java
  
  double hx; // x step = (xmax-xmin)/(nx-1) 
  double hy; // y step = (ymax-ymin)/(ny-1) 
  double hx2, hy2, hxy; // hx^2, hy^2, hx*hy 
  int cs;    // index of constant vector 
  double pdxx[][] = new double[nd][nd];
  double pdxy[][] = new double[nd][nd];
  double pdyy[][] = new double[nd][nd];
  double pdx[][] = new double[nd][nd];
  double pdy[][] = new double[nd][nd];
  double pd[][] = new double[nd][nd];
  int md;  // mid point correction  nd/2 
  double error;         // sum of absolute exact-computed  
  double avgerror;      // average error                   
  double maxerror;      // maximum cell error              
  DecimalFormat f6 = new DecimalFormat("0.000");
  DecimalFormat f7 = new DecimalFormat(" 0.000");
  DecimalFormat f8 = new DecimalFormat("  0.000");
  DecimalFormat fe = new DecimalFormat("0.000E00");
  // sample  System.out.println("Ua("+f6.format(xmin)+","+
  nderiv N = new nderiv();

  pde_abc_eq()
  {
    int i, j, ii, jj;
  
    System.out.println("pde_abc_eq.java running with user file  abc.java");
    md = nd/2; // mid point correction 
    hx = (xmax-xmin)/(nx-1); // step size in x 
    hy = (ymax-ymin)/(ny-1); // step size in y 
    System.out.println("xmin="+f6.format(xmin)+", xmax="+f6.format(xmax)+
		       ", nx="+nx+", hx="+f6.format(hx));
    System.out.println("ymin="+f6.format(ymin)+", ymax="+f6.format(ymax)+
		       ", ny="+ny+", hy="+f6.format(hy));
    hx2 = hx*hx;
    hy2 = hy*hy;
    hxy = hx*hy;
    cs = (nx-2)*(ny-2); // constant column subscript 

    clear();
    init_matrix();
    System.out.println("matrix initialized");
    if(ifcheck>0) print_mat();

    new simeq((nx-2)*(ny-2), ut, us);

    printus();
    if(ifcheck>0)
    {
      printexact(); // us vector is the solution 
      printdiff();
      check();
      avgerror = error/((nx-2)*(ny-2));
      System.out.println("total error="+error);
      System.out.println("average error="+avgerror);
      System.out.println("max error="+maxerror);
    }
    else
    {
      // addedderivative computation on solution and check against PDE 
	check_soln();
    }
    check_soln(); // do it always for this demo
  } // end pde_abc_eq  constructor

  void clear()
  {
    int ii;

    for(int i=0; i<nx-2; i++) // solution coordinates 
    {
      for(int j=0; j<ny-2; j++)
      {
        ii = i+(nx-2)*j;
        us[ii] = 0.0;
        if(ifcheck>0) uu[ii] = ABC.u(xmin+(i+1)*hx,ymin+(j+1)*hy);
      }
    }
  } // end clear 

  void printus() // print the computed solution 
  {
    System.out.println(" computed solution");
    System.out.println("us[x,y] ");
    for(int i=1; i<nx-1; i++)
    {
	System.out.println("row i="+i+", at x="+f6.format(xmin+i*hx)+
			   ", first y="+f6.format(ymin+hy));
      for(int k=1; k<ny-1; k=k+8)
      {
        for(int j=k; j<Math.min(ny-1,k+8); j++)
        {
	  System.out.print(f8.format(us[(i-1)+(nx-2)*(j-1)])+" ");
        }
        System.out.println(" ");
      }
    }
    System.out.println(" ");
  } // end printus 

  void check() // typically not known, instrumentation 
  {
    double uexact, err;
  
    error = 0.0;
    maxerror = 0.0;
    for(int i=1; i<ny-1; i++)
    {
      for(int j=1; j<nx-1; j++)
      {
        uexact = ABC.u(xmin+i*hx, ymin+j*hy);
        err = Math.abs(us[(i-1)+(nx-2)*(j-1)]-uexact);
        error = error + err;
        if(err>maxerror) maxerror=err;
      }
    }
  } // end check 

  void printexact() // typically not known 
  {
    System.out.println(" exact solution");
    System.out.println("u(x,y) ");
    for(int i=1; i<nx-1; i++)
    {
      for(int j=1; j<Math.min(ny-1,9); j++)
      {
	System.out.print(f8.format(ABC.u(xmin+i*hx, ymin+j*hy))+" ");
      }
      System.out.println(" ");
    }
    if(Math.min(ny,9)<ny) System.out.println("Not all columns printed ");
    System.out.println(" ");
  } // end printexact 

  void printdiff() // only when checking 
  {
    System.out.println(" difference exact-computed solution");
    System.out.println("");
    System.out.println("u(x,y)-us[x,y] ");
    for(int i=1; i<nx-1; i++)
    {
      for(int j=1; j<Math.min(ny-1,8); j++)
      {
        System.out.print(f8.format(ABC.u(xmin+i*hx, ymin+j*hy)-
			   us[(i-1)+(nx-2)*(j-1)])+" ");
      }
      System.out.println(" ");
    }
    System.out.println(" ");
  } // end printdiff 

  void print_coef()
  {
    System.out.println("pdxx     j=0     j=1     j=2     j=3     j=4  y");
    for(int i=0; i<nd; i++)
    {
      System.out.print("i="+i+"   ");
      for(int j=0; j<nd; j++)
      {
	  System.out.print(f7.format(pdxx[i][j])+"  ");
      }
      System.out.println(" ");
    }
    System.out.println("x");

    System.out.println("pdxy     j=0     j=1     j=2     j=3     j=4  y");
    for(int i=0; i<nd; i++)
    {
      System.out.print("i="+i+"   ");
      for(int j=0; j<nd; j++)
      {
	System.out.print(f7.format(pdxy[i][j])+"  ");
      }
      System.out.println(" ");
    }
    System.out.println(" ");

    System.out.println("pdyy     j=0     j=1     j=2     j=3     j=4  y");
    for(int i=0; i<nd; i++)
    {
      System.out.print("i="+i+"   ");
      for(int j=0; j<nd; j++)
      {
	System.out.print(f7.format(pdyy[i][j])+"  ");
      }
      System.out.println(" ");
    }
    System.out.println(" ");

    System.out.println("pdx      j=0     j=1     j=2     j=3     j=4  y");
    for(int i=0; i<nd; i++)
    {
      System.out.print("i="+i+"   ");
      for(int j=0; j<nd; j++)
      {
	System.out.print(f7.format(pdx[i][j])+"  ");
      }
      System.out.println(" ");
    }
    System.out.println(" ");

    System.out.println("pdy      j=0     j=1     j=2     j=3     j=4  y");
    for(int i=0; i<nd; i++)
    {
      System.out.print("i="+i+"   ");
      for(int j=0; j<nd; j++)
      {
	System.out.print(f7.format(pdy[i][j])+"  ");
      }
      System.out.println(" ");
    }
    System.out.println(" ");
  } // end print_coef 

  void build_coef(int ix, int iy, int ic)
  {
    int a[] = new int[25];
    int bb[] = new int[1];

    if((nd&1)==0 || nd<3)
    { 
      System.out.println("nd must be odd and >1 , bye.");
      System.exit(1);
    }
    if(nd!=5)
    {
      System.out.println("init_matrix needs enhancement, bye.");
      System.exit(1);
    }
    // build discretization matrices 
    if(ifcheck>0)
    {
      System.out.println("build_coef for point ix="+ix+", iy="+iy);
    }
    for(int i=0; i<nd; i++)
    {
      for(int j=0; j<nd; j++)
      {
        pdxx[i][j]=0.0; // coefficient of a1(x,y) 
        pdxy[i][j]=0.0; // coefficient of b1(x,y) 
        pdyy[i][j]=0.0; // coefficient of c1(x,y) 
        pdx [i][j]=0.0; // coefficient of d1(x,y) 
        pdy [i][j]=0.0; // coefficient of e1(x,y) 
        pd  [i][j]=0.0; // coefficient of f1(x,y) 
      }
    }
    pd[ix][iy] = 1.0;
    N.deriv(2, nd, ix, a, bb); // for pdxx 
    if(ifcheck>5)
    {
      System.out.println("deriv pdxx returns bb="+bb[0]+", aa="+
		         a[0]+" "+a[1]+" "+a[2]+" "+a[3]+" "+a[4]);
    }
    for(int i=0; i<nd; i++)
      pdxx[i][iy] = (double)a[i]/((double)bb[0]*hx2);

    N.deriv(2, nd, iy, a, bb); // for pdyy 
    if(ifcheck>5)
    {
      System.out.println("deriv pdyy returns bb="+bb[0]+", aa="+
		         a[0]+" "+a[1]+" "+a[2]+" "+a[3]+" "+a[4]);
    }
    for(int j=0; j<nd; j++)
      pdyy[ix][j] = (double)a[j]/((double)bb[0]*hy2);

    N.deriv(1, nd, ix, a, bb); // for pdx 
    if(ifcheck>5)
    {
      System.out.println("deriv pdx returns bb="+bb[0]+", aa="+
		         a[0]+" "+a[1]+" "+a[2]+" "+a[3]+" "+a[4]);
    }
    for(int i=0; i<nd; i++)
      pdx[i][iy] = (double)a[i]/((double)bb[0]*hx);

    N.deriv(1, nd, iy, a, bb); // for pdy 
    if(ifcheck>5)
    {
      System.out.println("deriv pdy returns bb="+bb[0]+", aa="+
		         a[0]+" "+a[1]+" "+a[2]+" "+a[3]+" "+a[4]);
    }
    for(int j=0; j<nd; j++)
      pdy[ix][j] = (double)a[j]/((double)bb[0]*hy);

    for(int i=0; i<nd; i++)
      for(int j=0; j<nd; j++)
        pdxy[i][j] = pdx[i][iy]*pdy[ix][j];

    if(ifcheck>2) print_coef();
  } // end build_coef 

  void check_row(int row)
  {
    int ii;
    double sum = 0.0;

    for(int i=0; i<nx-2; i++)
    {
      for(int j=0; j<nx-2; j++)
      {
        ii = i+(nx-2)*j;
        sum += ut[row][ii]*uu[ii];
      }
    }
    sum -= ut[row][cs];
    if(Math.abs(sum)>0.001)
       System.out.println("row="+row+" is bad err="+sum);
  } // end check_row 

  void build_row(int i, int j, int bi, int bj)
  {
    // i,j are in element coordinates, ir,jr are all element coord 
    int ii, jj, ij, ir, jr, ia, ja;
    double  p, a1c, b1c, c1c, d1c, e1c, f1c, cc;

    ir = i;
    jr = j;
    ii = (i-1)+(nx-2)*(j-1); // row of matrix 
    if(ifcheck>2)
      System.out.println("working row "+ii+", i="+i+", j="+j+", bi="+bi+
                         ", bj="+bj);
      a1c = ABC.a1(xmin+ir*hx,ymin+jr*hy);
      b1c = ABC.b1(xmin+ir*hx,ymin+jr*hy);
      c1c = ABC.c1(xmin+ir*hx,ymin+jr*hy);
      d1c = ABC.d1(xmin+ir*hx,ymin+jr*hy);
      e1c = ABC.e1(xmin+ir*hx,ymin+jr*hy);
      f1c = ABC.f1(xmin+ir*hx,ymin+jr*hy);
      cc  =  ABC.c(xmin+ir*hx,ymin+jr*hy);
      if(ifcheck>5)
      {
        System.out.println("a1="+a1c+", b1="+b1c+", c1="+c1c);
        System.out.println("d1="+d1c+", e1="+e1c);
        System.out.println("f1="+f1c+", c="+cc);
      }
      for(int ki=0; ki<nd; ki++)
      {
        for(int kj=0; kj<nd; kj++)
	{
          ia = ir-bi+ki-md; // real i,j in all element coordinates 
          ja = jr-bj+kj-md; // may be boundary 
          ij = (ia-1)+(nx-2)*(ja-1); // in ut coordinates 
          p =  a1c*pdxx[ki][kj];
          p += b1c*pdxy[ki][kj];
          p += c1c*pdyy[ki][kj];
          p += d1c* pdx[ki][kj];
          p += e1c* pdy[ki][kj];
          p += f1c*  pd[ki][kj];
          if(ia<0)    System.out.println("***-error in ia="+ia);
          if(ia>nx-1) System.out.println("***+error in ia="+ia);
          if(ja<0)    System.out.println("***-error in ja="+ja);
          if(ja>ny-1) System.out.println("***+error in ja="+ja);
          if(ia==0)  // boundary value tests 
            ut[ii][cs] -= p * ABC.u(xmin+(ia)*hx, ymin+(ja)*hy);
          else if(ia==nx-1)
            ut[ii][cs] -= p * ABC.u(xmin+(ia)*hx, ymin+(ja)*hy);
          else if(ja==0)
            ut[ii][cs] -= p * ABC.u(xmin+(ia)*hx, ymin+(ja)*hy);
          else if(ja==ny-1)
            ut[ii][cs] -= p * ABC.u(xmin+(ia)*hx, ymin+(ja)*hy);
          else // non boundary case 
	  {
            ut[ii][ij] = p;
            if(ifcheck>5)
		System.out.println("ij="+ij+", p=u[ii][ij]="+ut[ii][ij]);
	  } 
	} 
      }
      ut[ii][cs] += cc; 
      System.out.println("ii="+ii+", ut[ii][ii]="+f8.format(ut[ii][ii])+
			 ", ut[ii][cs]="+f8.format(ut[ii][cs]));
      if(ifcheck>0) check_row(ii);
  } // end build_row 

  void init_matrix()
  {
      int i, j, ii, jj;

    // zero internal cells, the matrix does not include boundary cells 
    for(i=1; i<nx-1; i++)
      for(j=1; j<ny-1; j++)
        for(ii=1; ii<nx-1; ii++)
          for(jj=1; jj<ny-1; jj++)
            ut[(i-1)+(nx-2)*(j-1)][(ii-1)+(nx-2)*(jj-1)] = 0.0;

    System.out.println("internal cells zeroed ");          
    build_coef(md, md, md);
    for(i=2; i<nx-2; i++) // compute standard central cells 
    {
      for(j=2; j<ny-2; j++)
      {
        build_row(i, j, 0, 0);
      } // end j loop 
    } // end i loop 
    System.out.println("standard central cells generated ");

    // do four sides, one element in. special cases 
    build_coef(md-1, md, md);
    for(j=2; j<ny-2; j++) // do i=1 
    {
      i = 1;
      if(ifcheck>2) System.out.println("working row edge i="+i+", j="+j);
      build_row(i, j, -1, 0);
    } // end j loop on edge 

    build_coef(md+1, md, md);
    for(j=2; j<ny-2; j++) // do  i=nx-2 
    {
      i = nx-2;
      if(ifcheck>2) System.out.println("working row edge i="+i+", j="+j);
      build_row(i, j, 1, 0);
    } // end j loop on edge 

    build_coef(md, md-1, md);
    for(i=2; i<nx-2; i++) // do j=1 
    {
      j = 1;
      if(ifcheck>2) System.out.println("working row edge i="+i+", j="+j);
      build_row(i, j, 0 , -1);
    } // end i loop on edge 

    build_coef(md, md+1, md);
    for(i=2; i<nx-2; i++) // do  j=ny-2 
    {
      j = ny-2;
      if(ifcheck>2) System.out.println("working row edge i="+i+", j="+j);
      build_row(i, j, 0 , 1);
    } // end i loop on edge 
    System.out.println("four sides, one element in, special case, generated ");
    // now do four corners 
    i=1;
    j=1;
    build_coef(md-1, md-1, md);
    build_row(i, j, -1, -1);
    i=1;
    j=ny-2;
    build_coef(md-1, md+1, md);
    build_row(i, j, -1, 1);
    i=nx-2;
    j=1;
    build_coef(md+1, md-1, md);
    build_row(i, j, 1, -1);
    i=nx-2;
    j=ny-2;
    build_coef(md+1, md+1, md);
    build_row(i, j, 1, 1);
  } // end init_matrix 

  void print_mat()
  {
    System.out.println(" initial matrix ");
    for(int i=0; i<(nx-2)*(ny-2); i++)
    {
      System.out.println("i="+i+", j=0 .. "+((nx-2)*(ny-2)));
      for(int k=0; k<=(nx-2)*(ny-2); k=k+7)
      {
        for(int j=0; j<7 && j+k<=(nx-2)*(ny-2); j++)
        {
	  System.out.print(f8.format(ut[i][j+k])+" ");
        }
        System.out.println(" ");
      }
      System.out.println(" ");
    } 
    System.out.println(" ");
  } // end print_mat 

  void check_soln() // against PDE, added later, used rderiv and xg,yg 
  {
               // for all interior points 
    double x, y, err, maxerr, avgerr, rmserr;
    double sumerr, sumsqerr;
    double cx[][] = new double[nx][nx];
    double cxx[][] = new double[nx][nx];
    double cy[][] = new double[ny][ny];
    double cyy[][] = new double[ny][ny];
    double tcxy[] = new double[ny];
    double dxx, dxy, dyy, dx, dy; // a1, b1, c1, d1, e1 coefficients  
    double xg[] = new double[nx]; // may be available 
    double yg[] = new double[ny]; 
    double ug[][] = new double[nx][ny]; // full, including boundaries, copied from us[] 
    int cnt = 0;
    int ii;

    if(ifcheck>5) System.out.println("check_soln setup");
    for(int i=0; i<nx; i++) xg[i] = xmin + (double)i*hx;
    for(int j=0; j<ny; j++) yg[j] = ymin + (double)j*hy;
    // transform us[] into ug, add boundaries, make code simpler 
    for(int i=1; i<nx-1; i++) // get from us[] 
    {
      for(int j=1; j<ny-1; j++)
      {
        ii = (i-1)+(nx-2)*(j-1); // ii for  us[ii] computed solution at i,j 
        ug[i][j] = us[ii];
        if(ifcheck>5) System.out.println("ug["+i+"]["+j+"]="+ug[i][j]);
      }
    }
    for(int i=0; i<nx; i++) // fill boundaries, loops below easier 
    {
      ug[i][0]    = ABC.u(xg[i],yg[0]);
      ug[i][ny-1] = ABC.u(xg[i],yg[ny-1]);
    }
    for(int j=0; j<ny; j++) // fill boundaries 
    {
      ug[0][j]    = ABC.u(xg[0],   yg[j]);
      ug[nx-1][j] = ABC.u(xg[nx-1],yg[j]);
    }

    for(int i=0; i<nx; i++)
    {
      new rderiv(1,nx,i,hx,cx[i]); // fills row of derivative matrix 
      new rderiv(2,nx,i,hx,cxx[i]);
    }
    for(int j=0; j<ny; j++)
    {
      new rderiv(1,ny,j,hy,cy[j]);
      new rderiv(2,ny,j,hy,cyy[j]);
    }
    if(ifcheck>5)
    {
      System.out.println("cx=");
      for(int i=0; i<nx; i++)
      {
        for(int k=0; k<nx; k++)
        {
	  System.out.print(f8.format(cx[i][k])+" ");
        }
        System.out.println(" ");
      }
      System.out.println(" ");
      System.out.println("cxx=");
      for(int i=0; i<nx; i++)
      {
        for(int k=0; k<nx; k++)
        {
	  System.out.print(f8.format(cxx[i][k])+" ");
        }
        System.out.println(" ");
      }
      System.out.println(" ");
      System.out.println("cy=");
      for(int j=0; j<ny; j++)
      {
        for(int k=0; k<ny; k++)
        {
	  System.out.print(f8.format(cy[j][k])+" ");
        }
        System.out.println(" ");
      }
      System.out.println(" ");
      System.out.println("cyy=");
      for(int j=0; j<ny; j++)
      {
        for(int k=0; k<ny; k++)
        {
	    System.out.print(f8.format(cyy[j][k])+" ");
        }
        System.out.println(" ");
      }
      System.out.println(" ");
    } // end ifcheck 

    sumerr = 0.0;
    sumsqerr = 0.0;
    maxerr = 0.0;
    for(int i=1; i<nx-1; i++) // check full equation, at every point 
    {
      x = xg[i];
      for(int j=1; j<ny-1; j++)
      {
        y = yg[j];
        // compute derivatives, plug into PDE, subtract RHS, should be zero 
        dxx = 0.0;
        dx  = 0.0;
        for(int k=0; k<nx; k++)
        {
          // PDE  a1(x,y)*uxx(x,y) term 
	  dxx += cxx[i][k]*ug[k][j];
          // PDE  d1(x,y)*ux(x,y) term 
	  dx += cx[i][k]*ug[k][j];
        } 
        dyy = 0.0;
        dy  = 0.0;
        for(int k=0; k<ny; k++)
        {
          // PDE c1(x,y)*uyy(x,y) term 
	  dyy += cyy[j][k]*ug[i][k];
          // PDE e1(x,y)*uy(x,y) term 
	  dy += cy[j][k]*ug[i][k];
        }
        // PDE b1(x,y)*uxy(x,y) term 
        dxy = 0.0;
        for(int jj=0; jj<ny; jj++) // walk jj in y, computing uxx 
        {
	  tcxy[jj] = 0.0; // numeric  ux part of uxy 
          for(int k=0; k<nx; k++)
          {
	    tcxy[jj] += cx[i][k]*ug[k][jj];
	  }
        }
        for(int k=0; k<ny; k++)  // uy part of uxy 
        {
          dxy += cy[j][k]*tcxy[k];  // PDE uxy term 
        }
        // check 
        if(ifcheck>5)
        {
	  System.out.println("numerical dxx="+f8.format(dxx)+
                             ", dxy="+f8.format(dxy)+
                             ", dyy="+f8.format(dyy));
	  System.out.println("analytic  dxx="+f8.format(6.0*x+6.0*y)+
                             ", dxy="+f8.format(6.0*x+8.0*y+5.0)+
                             ", dyy="+f8.format(12.0*y+8.0*x));
	  System.out.println("numerical dx= "+f8.format(dx)+
                             ", dy= "+f8.format(dy)+
                             ", u = "+f8.format(ug[i][j]));
	  System.out.println("analytic  dx= "+f8.format(3.0*x*x+6.0*x*y+4.0*y*y+5.0*y+6.0)+
                             ", dy= "+f8.format(6.0*y*y+3.0*x*x+8.0*x*y+5.0*x+7.0)+
                             ", u = "+f8.format(ABC.u(x,y)));
        }
        err = ABC.a1(xg[i],yg[j])*dxx +
	      ABC.b1(xg[i],yg[j])*dxy +
              ABC.c1(xg[i],yg[j])*dyy +
              ABC.d1(xg[i],yg[j])*dx  +
              ABC.e1(xg[i],yg[j])*dy  +
              ABC.f1(xg[i],yg[j])*ug[i][j] - // PDE u term at i,j
	      ABC.c(xg[i],yg[j]);           // RHS 
        err = Math.abs(err);
        sumerr += err;
        sumsqerr += err*err;
        maxerr = Math.max(err,maxerr);
        cnt++;  
      } // end j 
    } // end i 
    rmserr = Math.sqrt(sumsqerr/(double)cnt);
    avgerr = sumerr/(double)cnt;
    System.out.println("check_soln values against PDE");
    System.out.println("maxerr="+maxerr+", rmserr="+rmserr+
                       ", avgerr="+avgerr);
  } // end check_soln  

  public static void main(String[] args)
  {
    new pde_abc_eq();
  }
} // end class pde_abc_eq