fixed Fix() function, glaring error before that had different phases for different processors. Now, however, -np 3 gives inconsistent results
didn't really change anything, but currently getting inconsistent results fornp >= 3, and sometimes the Newton doesn't converge after the NullSolve. Interestingly, the JtJ is exactly the same always, even if the results of the solve are different
Added a iteration limit to Newton.

Some discoveries:

- The slowness of my code compared to past commits was not real. This was really due to me using Nc = 3 instead of Nc = 1 before.

- The non-convergence of the NewtonSolve after NullSolve in passive is due to the approximate null-vector being bad, i.e. imaginary part too small or something. This is probably a limitation of the NullSolve method.

- The inconsistency of results, both in my NewtonSolve and in the EPSSolve in NullSolve, might be a real issue. Also, it does happen, albeit much more rarely, for -np 2. This is probably a memory leak in my code, or much less likely, an actual feature of one of the solvers I'm using. Maybe Pastix uses Monte Carlo or something (doubtful).
all three problems not solved, but accounted for. Feel free to continue
improved makefile, changed some names around
* operators for functions
fixed MAJOR bug with the setlast2 in FormJf for psifix not getting the correct psi_fix.imag() due to the parallel functions.
added * operations to Jacobian
completely cleaned up Jacobian, removed ugly temp functions
added "GetValue" for parallel Vec Get Values
changed some file names around
fixed important typo, the Bcast for c in Mode's read constructor used MPI_INT instead of MPI_DOUBLE
changed some minor stuff around, now ready to code above-threshold
wrote FieldSquare function, seems to work
changed const Medium to Medium in a bunch of places, since to compute H we modify cv
wrote code for above threshold derivative terms and SetBlock Tensor, probably doesn't work yet
nope, as usual, doesn't work
fixed some TODO stuff relating to w.imag() versus c, now |f| stays at around 1e+1 instead of blowing up quickly to 1e+8
fixed minor typo, forgot c in denominator of H. Still doesn't work though
holy ****ing ****! It finally converges. It doesn't do so cubically, and takes bout 6-7 iterations, but at least now the c doesn't blow up. It turns out the problem was due to the MODE READ CONSTRUCTOR of all things, not doing the MPI_BCAST of Dout correctly, I was putting it in **Dout, not *Dout.

Now all I have to do is fix up the other derivatives, dfdk and dfdpsi to get it to converge cubically!