Guest User

Untitled

a guest
Mar 22nd, 2018
99
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 1.79 KB | None | 0 0
  1. Eigenvectors@Transpose[A]
  2.  
  3. SetAttributes[makeVecs, HoldAll];
  4. makeVecs[vals_, vecs_, s_, rQ_] := With[{n = Length[vals]},
  5. If[TrueQ[rQ],
  6. Do[Switch[Sign[Im[vals[[k]]]], 0, Null,
  7. 1, vecs[[k]] += s I vecs[[k + 1]],
  8. -1, vecs[[k]] = Conjugate[vecs[[k - 1]]]], {k, n}]]]
  9.  
  10. Options[getEigensystem] = {Mode -> Automatic};
  11.  
  12. getEigensystem[mat_?SquareMatrixQ, opts : OptionsPattern[]] :=
  13. Module[{m = mat, chk, ei, ev, lm, lv, n, rm, rQ, rv}, n = Length[mat];
  14. Switch[OptionValue[Mode],
  15. Right | Automatic, {lm, rm} = {"N", "V"},
  16. Left, {lm, rm} = {"V", "N"},
  17. All, {lm, rm} = {"V", "V"},
  18. _, {lm, rm} = {"N", "V"}];
  19. LinearAlgebra`LAPACK`GEEV[lm, rm, m, ev, ei, lv, rv, chk, rQ];
  20. If[! TrueQ[chk], Message[getEigensystem::eivec0]; Return[$Failed, Module]];
  21. If[rQ, ev += I ei];
  22. Switch[OptionValue[Mode],
  23. Right | Automatic, rv = ConjugateTranspose[rv];
  24. makeVecs[ev, rv, 1, rQ]; {ev, rv},
  25. Left, lv = ConjugateTranspose[lv];
  26. makeVecs[ev, lv, -1, rQ]; {ev, lv},
  27. All, {lv, rv} = ConjugateTranspose /@ {lv, rv};
  28. makeVecs[ev, rv, 1, rQ]; makeVecs[ev, lv, -1, rQ];
  29. {ev, rv, lv},
  30. _, rv = ConjugateTranspose[rv];
  31. makeVecs[ev, rv, 1, rQ]; {ev, rv}]]
  32.  
  33. m = {{1., -2., 1., 9., 6.},
  34. {8., 5., 3., 7., 7.},
  35. {-9., -3., 9., 5., -2.},
  36. {-5., -4., 6., 8., -6.},
  37. {-3., 9., -2., 6., -3.}};
  38.  
  39. {vals, rvecs, lvecs} = getEigensystem[m, Mode -> All];
  40.  
  41. Norm[m.Transpose[rvecs] - Transpose[rvecs].DiagonalMatrix[vals]]
  42. 4.16257*10^-14
  43.  
  44. Norm[lvecs.m - DiagonalMatrix[vals].lvecs]
  45. 2.73754*10^-14
Add Comment
Please, Sign In to add comment