@tomleslie 

1. The "extraneous crap" inverts the isomorphism G2-->G1. This would be necessary when IsIsomorphic(G1,G2,f) gives an error but  G1 and G2 are still isomorphic and f is wanted.
2. I said: seems to work. Anyway it works for the original question.
3. Do you have an example of a workaround for a Maple command which is guaranteed to work in each and every situation? I'd say that IsIso it's a typical workaround; do you have a counterexample for it?

@tomleslie 

It seems that my workaround is invisible :-)

IsIso(G1,G2);
      false

@mehdibaghaee 

It is not easy to construct good examples. This one is not good enough. How do you know that your eigenvalues are correct? (Keep in mind that the determinant is very complicated to compute and would need exact entries which you don't have).
An ideal example would have rational (or even integer) entries and all the eigenvalues be known (not necessarily rational, but relatively simple).
Maple is able to compute Eigenvalues(K,M)  but I simply don't know how far are these eigenvalues from the exact ones for such huge matrices. So, good examples are mandatory.

Note that for standard eigenvalues Maple works very well and probably it is fine for generalized ones too; for your example they are not close but I suspect that your eigenvalues are not correct.

restart;
with(LinearAlgebra):
n:=700:
L:=[seq(3.333*k,k=1..n)]:
A:=DiagonalMatrix(L,datatype=float[8]):
Q:=RandomMatrix(n,n, generator=rand(-1. .. 1.)):
Q1:=Q^(-1):
B:=Q1.A.Q:
LL:=sort(convert(Eigenvalues(B),list)):
max(abs~(L-LL));
#                  HFloat(4.88817022414878e-8)

Have you tried the above workaround? (Just use IsIso instead of IsIsomorphic). E.g.

restart;
IsIso := proc(G1::GRAPHLN, G2::GRAPHLN, phi::name)
local r;
try  GraphTheory:-IsIsomorphic(args); 
  catch:
  if nargs=2 then GraphTheory:-IsIsomorphic(G2,G1);
  else r:=GraphTheory:-IsIsomorphic(G2,G1,phi);
       phi:=sort(map(u -> rhs(u)=lhs(u), eval(phi)));
  r;
  fi;  
end try;
end proc:
interface(rtablesize=infinity):
n:= 3: N:= n^2:
<ListTools:-Categorize(
   (M1,M2)-> IsIso(GraphTheory:-Digraph~([M1,M2], n)[]),
   map(L-> Matrix(n,n,L), combinat:-permute([0$N, 1$N], N)   )
)>; 

@tomasmedici 

This is not possible without examples with known eigenvalues, because the accuracy of the command (working with hardware floats) is not documented for such large matrices. But you refuse to provide them.

@mehdibaghaee 

The workaround above seems to work now in all situations (i.e. including the correct isomorphism if requested).

@dharr 

@tomleslie 

An ad-hoc workaround:

IsIso := proc(G1::GRAPHLN, G2::GRAPHLN, phi::name)
local r;
try  GraphTheory:-IsIsomorphic(args); 
  catch:
  if nargs=2 then GraphTheory:-IsIsomorphic(G2,G1);
  else r:=GraphTheory:-IsIsomorphic(G2,G1,phi);
       phi:=sort(map(u -> rhs(u)=lhs(u), eval(phi)));
  r;
  fi;  
end try;
end proc:

Edit:  seems to work in all situations and returns the correct isomorphism.

@mehdibaghaee 

OK, I wish you good luck.

@mehdibaghaee 
It is your turn to say something about the points 1-2 above. After that I will be able to comment. (How do you know that the method works without examples with known eigenvalues?)

@mehdibaghaee 

The problem is changed now.
You should test the problem with small matrices. Nobody will want to see (and work with) your 700x700 matrices.
Note that Eigenvalues (unlike Determinant) works for big matrices (due to special algorithms) so probably the change of variables is not needed.

Download lista1.mw

@sand15 

assuming positive  means that all the variables are assumed to be >0.
On my system (Maple 2016.2, 64 bit) I get FAIL  for     is( r^2 - r*s + s^2> 0) assuming r >0, s>0;

positive means real and >0. So,  is(I, positive)  is of course false (not being real).

r^2+s^2 > 0   is simple for Maple (sum of two positive numbers); s^2 - s*r + r^2  seems to be too complicated for the moment.

Edit (concerning "simple conditions"):

is( (r - s)^2 + r^2  + s^2  > 0)  assuming positive;
                              FAIL
is( (r - s/2)^2 + 3/4*s^2  > 0) assuming positive;
                              FAIL
is( (r - s/2)^2 + 3/4*s^2  > 0) assuming positive, r>s/2;
                              true

@mehdibaghaee 

My question was about the type of the problem. Where/how  those 10^6 variables are they used?
For instance, do you have a system of equations with 10^6 unknowns?

@mehdibaghaee 

It would be interesting to know the type of the problem which needs 10^6 variables and can be solved.

@Markiyan Hirnyk 

But you have spent a rather long time with this problem as I can see. So, selectively busy.