Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- Eigenvectors@Transpose[A]
- SetAttributes[makeVecs, HoldAll];
- makeVecs[vals_, vecs_, s_, rQ_] := With[{n = Length[vals]},
- If[TrueQ[rQ],
- Do[Switch[Sign[Im[vals[[k]]]], 0, Null,
- 1, vecs[[k]] += s I vecs[[k + 1]],
- -1, vecs[[k]] = Conjugate[vecs[[k - 1]]]], {k, n}]]]
- Options[getEigensystem] = {Mode -> Automatic};
- getEigensystem[mat_?SquareMatrixQ, opts : OptionsPattern[]] :=
- Module[{m = mat, chk, ei, ev, lm, lv, n, rm, rQ, rv}, n = Length[mat];
- Switch[OptionValue[Mode],
- Right | Automatic, {lm, rm} = {"N", "V"},
- Left, {lm, rm} = {"V", "N"},
- All, {lm, rm} = {"V", "V"},
- _, {lm, rm} = {"N", "V"}];
- LinearAlgebra`LAPACK`GEEV[lm, rm, m, ev, ei, lv, rv, chk, rQ];
- If[! TrueQ[chk], Message[getEigensystem::eivec0]; Return[$Failed, Module]];
- If[rQ, ev += I ei];
- Switch[OptionValue[Mode],
- Right | Automatic, rv = ConjugateTranspose[rv];
- makeVecs[ev, rv, 1, rQ]; {ev, rv},
- Left, lv = ConjugateTranspose[lv];
- makeVecs[ev, lv, -1, rQ]; {ev, lv},
- All, {lv, rv} = ConjugateTranspose /@ {lv, rv};
- makeVecs[ev, rv, 1, rQ]; makeVecs[ev, lv, -1, rQ];
- {ev, rv, lv},
- _, rv = ConjugateTranspose[rv];
- makeVecs[ev, rv, 1, rQ]; {ev, rv}]]
- m = {{1., -2., 1., 9., 6.},
- {8., 5., 3., 7., 7.},
- {-9., -3., 9., 5., -2.},
- {-5., -4., 6., 8., -6.},
- {-3., 9., -2., 6., -3.}};
- {vals, rvecs, lvecs} = getEigensystem[m, Mode -> All];
- Norm[m.Transpose[rvecs] - Transpose[rvecs].DiagonalMatrix[vals]]
- 4.16257*10^-14
- Norm[lvecs.m - DiagonalMatrix[vals].lvecs]
- 2.73754*10^-14
Add Comment
Please, Sign In to add comment