Attributes[buildmatrix] = {HoldRest}
buildmatrix::usage = "buildmatrix[\"f\",m] builds a matrix m from the file
f \n   with the structure:\n
dimension \n
row index column index value \n
..."
buildmatrix/: buildmatrix[file_, matrix_] := 
   Block[{data, f, i, j, l, n}, 
    f = OpenRead[file]; n = Read[f, Number]; matrix = IdentityMatrix[n]; 
     data = ReadList[f, {Number, Number, Number}]; 
     Do[i = data[[l,1]]; j = data[[l,2]]; matrix[[i,j]] = data[[l,3]], 
      {l, Length[data]}]; Close[f]]

Attributes[buildpat] = {HoldRest}
buildpat::usage = "buildpat[\"f\",m] builds a matrix pattern m
from the file f \n   with the structure:\n
dimension \n
row index column index value \n
..."
buildpat/: buildpat[file_, matrix_] := 
   Block[{data, f, i, j, l, n}, 
    f = OpenRead[file]; n = Read[f, Number]; matrix = IdentityMatrix[n]; 
     data = ReadList[f, {Number, Number, Number}]; 
     Do[i = data[[l,1]]; j = data[[l,2]]; matrix[[i,j]] = 1.0, 
      {l, Length[data]}]; Close[f];
      matrix = Reverse[ Map[ If[# != 0, 0, 1]&, matrix, {2}] ]; ]




lcp/: lcp[m_, q_, n_] := Block[{l, indices, t, x}, indices = 
  Flatten[Apply[Outer, Prepend[Table[{l, l + n}, {l, n}], List]], n-1]; 
  t = Table[Join[IdentityMatrix[n][[l]],-m[[l]] ], {l,n}];
  Do[x = Inverse[Transpose[Transpose[t][[indices[[l]] ]] ] ] . q; 
  Print[indices[[l]], " ", x]; Null, {l, 2^n}] ]
lcp::usage = "lcp prints all the basic sequences and basic variables
for the n x n nondegenerate LCP (m/q)."




(* Formulas on 21 point stencil. *)
v[x_,y_] = Normal[Series[u[x,y],{x,0,5},{y,0,5}]];
sumx= Expand[c1 v[-1,2] +c2 v[0,2] +c3 v[1,2] +
c4 v[-2,1] + c5 v[-1,1] + c6 v[0,1] + c7 v[1,1] + c8 v[2,1] +
c9 v[-2,0] + c10 v[-1,0] + c11 v[0,0] + c12 v[1,0] + c13 v[2,0] +
c14 v[-2,-1] + c15 v[-1,-1] + c16 v[0,-1] + c17 v[1,-1] + c18 v[2,-1] +
c19 v[-1,-2] + c20 v[0,-2] + c21 v[1,-2] ] ;
eq[ 1]=Coefficient[sumx,Derivative[0,0][u][0,0]] == 0;
eq[ 2]=Coefficient[sumx,Derivative[1,0][u][0,0]] == 0;
eq[ 3]=Coefficient[sumx,Derivative[0,1][u][0,0]] == 0;
eq[ 4]=Coefficient[sumx,Derivative[2,0][u][0,0]] == 0;
eq[ 5]=Coefficient[sumx,Derivative[1,1][u][0,0]] == 0;
eq[ 6]=Coefficient[sumx,Derivative[0,2][u][0,0]] == 0;
eq[ 7]=Coefficient[sumx,Derivative[3,0][u][0,0]] == 0;
eq[ 8]=Coefficient[sumx,Derivative[2,1][u][0,0]] == 0;
eq[ 9]=Coefficient[sumx,Derivative[1,2][u][0,0]] == 0;
eq[10]=Coefficient[sumx,Derivative[0,3][u][0,0]] == 0;
eq[11]=Coefficient[sumx,Derivative[4,0][u][0,0]] == 0;
eq[12]=Coefficient[sumx,Derivative[3,1][u][0,0]] == 0;
eq[13]=Coefficient[sumx,Derivative[2,2][u][0,0]] == 0;
eq[14]=Coefficient[sumx,Derivative[1,3][u][0,0]] == 1;
eq[15]=Coefficient[sumx,Derivative[0,4][u][0,0]] == 0;
eq[16]=Coefficient[sumx,Derivative[5,0][u][0,0]] == 0;
eq[17]=Coefficient[sumx,Derivative[4,1][u][0,0]] == 0;
eq[18]=Coefficient[sumx,Derivative[3,2][u][0,0]] == 0;
eq[19]=Coefficient[sumx,Derivative[2,3][u][0,0]] == 0;
eq[20]=Coefficient[sumx,Derivative[1,4][u][0,0]] == 0;
eq[21]=Coefficient[sumx,Derivative[0,5][u][0,0]] == 0;

Solve[Prepend[Table[eq[i],{i,21}],c3==0]] (*Set some c's to eliminate
					some of the free parameters.*)
 
(* Differential operator formulas for equation transformed from
equilateral (x-y coordinates) triangle to right triangle (s-t
coordinates). *)
s/: s[x_, y_] := x - (1*y)/Sqrt[3]
 
t/: t[x_, y_] := (2*y)/Sqrt[3]
 
u/: u[x_, y_] := v[s[x, y], t[x, y]]
 
ux/: ux[x_, y_] = Derivative[1, 0][v][x - y/3^(1/2), (2*y)/3^(1/2)]
 
uy/: uy[x_, y_] = 
   (2*Derivative[0, 1][v][x - y/3^(1/2), (2*y)/3^(1/2)])/3^(1/2) - 
    Derivative[1, 0][v][x - y/3^(1/2), (2*y)/3^(1/2)]/3^(1/2)
 
xytost = {x -> s + t/2, y -> (3^(1/2)*t)/2, v[s, t] -> vV[0, 0], 
   Derivative[0, 1][v][s, t] -> vV[0, 1], 
   Derivative[0, 2][v][s, t] -> vV[0, 2], 
   Derivative[0, 3][v][s, t] -> vV[0, 3], 
   Derivative[0, 4][v][s, t] -> vV[0, 4], 
   Derivative[1, 0][v][s, t] -> vV[1, 0], 
   Derivative[1, 1][v][s, t] -> vV[1, 1], 
   Derivative[1, 2][v][s, t] -> vV[1, 2], 
   Derivative[1, 3][v][s, t] -> vV[1, 3], 
   Derivative[1, 4][v][s, t] -> vV[1, 4], 
   Derivative[2, 0][v][s, t] -> vV[2, 0], 
   Derivative[2, 1][v][s, t] -> vV[2, 1], 
   Derivative[2, 2][v][s, t] -> vV[2, 2], 
   Derivative[2, 3][v][s, t] -> vV[2, 3], 
   Derivative[2, 4][v][s, t] -> vV[2, 4], 
   Derivative[3, 0][v][s, t] -> vV[3, 0], 
   Derivative[3, 1][v][s, t] -> vV[3, 1], 
   Derivative[3, 2][v][s, t] -> vV[3, 2], 
   Derivative[3, 3][v][s, t] -> vV[3, 3], 
   Derivative[3, 4][v][s, t] -> vV[3, 4], 
   Derivative[4, 0][v][s, t] -> vV[4, 0], 
   Derivative[4, 1][v][s, t] -> vV[4, 1], 
   Derivative[4, 2][v][s, t] -> vV[4, 2], 
   Derivative[4, 3][v][s, t] -> vV[4, 3], 
   Derivative[4, 4][v][s, t] -> vV[4, 4]}
del2/: del2[f_] := Through[((D[#1, {x, 2}] & ) + (D[#1, {y, 2}] & ))[f], Plus]
 	(* Laplacian operator with respect to x, y, applied to function
	f. *)

(* Differential operator formulas for equation transformed from
trapezoid (x-y coordinates) to unit square (s-t coordinates). *)
s/: s[x_, y_] := (d y + 3 x - 3 d)/(2(d y + 3 Sqrt[3] - 3 d))

t/: t[x_, y_] := y/3

u/: u[x_, y_] := v[s[x, y], t[x, y]]

ux/: ux[x_, y_] = (Derivative[1, 0][v][
  (d y + 3 x - 3 d)/(2(d y + 3 Sqrt[3] - 3 d)), y/3] D[s[x,y], x])

uy/: uy[x_, y_] = (Derivative[1, 0][v][
  (d y + 3 x - 3 d)/(2(d y + 3 Sqrt[3] - 3 d)), y/3] D[s[x,y], y]) +
  (1/3) (Derivative[0, 1][v][
  (d y + 3 x - 3 d)/(2(d y + 3 Sqrt[3] - 3 d)), y/3] )

xytost = Flatten[{ (d y + 3 x - 3 d)/(2(d y + 3 Sqrt[3] - 3 d)) -> ss,
  x -> d + 2 Sqrt[3] ss - 2 d ss - d tt + 2 d ss tt, y -> 3 tt,
 Table[Derivative[i, j][v][ss, tt] -> vV[i, j], {i, 0, 4}, {j, 0, 4}]
 }]
del2/: del2[f_] := Through[((D[#1, {x, 2}] & ) + (D[#1, {y, 2}] & ))[f], Plus]
 	(* Laplacian operator with respect to x, y, applied to function
	f. *)

(* Compute and save coefficients of partials in full differential
	operator.
LU = del2[del2[u[x,y]]] + aA del2[ux[x,y]] + bB del2[uy[x,y]] +
	cC ux[x,y] + dD uy[x,y] //. xytost ;
Do[Do[ cvV[i,j] = Simplify[Coefficient[Collect[ LU, vV[i,j]], vV[i,j]]];
  OutputForm[SequenceForm["cvV[",i,",",j,"]=",
  FortranForm[cvV[i,j]]]] >>> coef.f , {j, 0, 4-i, 1}], {i, 0, 4, 1}];
*)
