CurrentValue[$FrontEnd,{"PrintAction"}]={"PrintToNotebook"}
CurrentValue[$FrontEnd,{"MessageOptions","KernelMessageAction"}]={"PrintToNotebook","Beep"}

(*Change history *)
(*bugs to nma@12000.org*)
(*may 10, 2023. Initial Version *)
(*may 10, 2023. Added check for geometric multiplicity >1 for now until supported*)
(*may 11, 2023. Added support for geometric multiplicity >1, improved formatting*)
(*may 11, 2023. Print final result at end*)

findV[mat_,V_]:=Module[{data,keys,result,header,p},
   {data,keys}=generateInitialTable[mat];
   result=processTable[mat,data,keys,V];
   
   
   Print[Style["The following is summary of the result found",14,Bold]];
   header={"λ","P matrix of eigenvectors"};
   p=SortBy[result,#[[1]]&];
   Print[Grid[Join[{header},
               MapAt[MatrixForm[Transpose[#]]&,p,{All,2}]
               ],
               Frame->All]];
   Print["Therefore the generalized P matrix of all the eigenvectors is"];
   p=Transpose@Flatten[p[[All,2]],1];
   Print[MatrixForm@p];
         
   result
]
 
(*****************************)
printTable[data_]:=Module[{header,new},

    header={"λ","algebraic multiplicity","geometric multiplicity","normal eigenvectors","Defective λ?"};
    Print@Text[Style[
    "We first use Mathematica to find the eigenvalues and the norma\n"<>
    "eigenvectors and examine the result to determine which eigenvalues\n"<>
    "are defective and which are not. This results in this table ",TextAlignment->Left]];
   
    Print[Grid[Join[{header},
               MapAt[MatrixForm[Transpose[#]]&,data,{All,4}]
               ],
               Frame->All]]
]
(******************************)
generateInitialTable[mat_]:=Module[{lambda,v,p,keys,k,n,data,idx,theV,header,x},
   {lambda,v}=Eigensystem[mat];
   Print["The Matrix is ",MatrixForm@mat];
   Print["its eigenvalues are ",MatrixForm@lambda];
   Print["its normal Eigenvectors are ",MatrixForm@Transpose[v]];

   p=PositionIndex[lambda];
   keys=Keys[p];

   (*this table will hold summary of Eigensystem finding. One row per each lambda*)
   (*1=lambda,2=algebraic multiplicity,3=geometric multiplicity,4=normal V,5=defective?*)

   data=Table[{0,0,0,0,0},Length@keys];

   (*Now build the table *)
   Do[
     k=keys[[n]];
     data[[n,1]]=k;
     data[[n,2]]=Length@p[k];(*algebraic multiplicity*)
     idx=p[k];
     theV=v[[idx[[1]];;idx[[-1]]]];
     theV=Cases[theV,x_/;Not@AllTrue[x,#==0&]];
     data[[n,3]]=Length@theV;(*geometric multiplicity*)
     data[[n,4]]=theV;
     data[[n,5]]=If[data[[n,2]]>data[[n,3]],True,False]
     ,
     {n,Length@keys}
   ];
   
    printTable[data];
   {data,keys}
]
(******************************)
(*geometric multiplicity =1, easy case *)
handleCaseOne[data_,mat_,V_,n_]:=Module[{howMuchDefective,currentV,g,eq,m,size=Length@mat,sol,setOfV},
  howMuchDefective=data[[n,2]]-data[[n,3]];
   Print["Since eigenvalue ",data[[n,1]]," is defective, and the degree of defect is ",
        howMuchDefective," then we need to find ",howMuchDefective,
        " additional generalized eigenvectors for this eigenvalue."];
 currentV=data[[n,4]]; (*normal V found*)
 setOfV=currentV;
 
 Do[
   g=Table[V[i],{i,1,size}];
   eq=mat.g==currentV;
   Print["Now we solve for generalized eigevector ",MatrixForm[g],
          " Using (A-λI).g=v1. This gives"];
   Print[MatrixForm[mat],"-(",data[[n,1]],")",
         MatrixForm[IdentityMatrix[size]],".",MatrixForm[g],"=",MatrixForm[Transpose@currentV]];

   Print["The above equation becomes\n",MatrixForm[mat-data[[n,1]]*IdentityMatrix[size]],
         ".",MatrixForm[g],"=",MatrixForm[Transpose@currentV]];

   sol=LinearSolve[(mat-data[[n,1]]*IdentityMatrix[size]),currentV[[1]]];
   Print["Using LinearSolve gives ",MatrixForm[g],"=",Style[MatrixForm@sol,Bold]];
   currentV={sol};
   setOfV=Join[setOfV, currentV]
   ,
   {m,1,howMuchDefective}
   ];
   
   setOfV
]
  
(******************************)
(* when geometric multiplicity >1, hard case *)   
handleCaseTwo[data_,mat_,V_,n_]:=Module[{howMuchDefective,currentV,g,eq,m,size=Length@mat,sol,newV,i},
   Print["This case of geometric multiplicity >1."];
   
   howMuchDefective=data[[n,2]]-data[[n,3]];
   Print["Since eigenvalue ",data[[n,1]]," is defective, and the\n",
         "degree of defect is ",howMuchDefective," then we need to find ",howMuchDefective,"\n",
         "additional generalized eigenvectors for this eigenvalue."];
  currentV=data[[n,4]]; (*normal V's found. More than one*)
  g=Table[V[i],{i,1,size}];
  
  Do[
     Print["Now we solve need to find a generlized eigenvector ",MatrixForm@g,"\n",
          "which is linearly independent to the set of eigenvectors already found ",
          MatrixForm[Transpose[currentV]],"\n",
          "Which we do by looking at the NullSpace of the transpose of the above as a matrix\n",
          "Which gives ",MatrixForm[NullSpace[currentV]],"\n",
          "And we need to just pick one of these rows as our new eigenvector.\n"
     ];
     newV=NullSpace[currentV];
     newV=newV[[-1,All]];
     currentV=Join[currentV,{newV}];
     Print["Taking the last row, and transposing it gives the generalized eigenvector as ",MatrixForm@g," = ",
            Style[MatrixForm[Transpose[{newV}]],Bold],"\nwhich we now add to the set of known eigenvectors found giving\n",
            MatrixForm[Transpose@currentV]
            ];
    
    
   ,
   {m,1,howMuchDefective}
   ];
   
   currentV  
]
(******************************)
processTable[mat_,data_,keys_,V_]:=Module[{n,setOfV,result},

   (* this table will hold each eigenvalue with its correspoding P matrix*)
   (* The P matrix will contains the Liearly indepenent eigenvectors generated*)
   (* for each lambda. This way we can make Jordan form at end *)
   
   result=Table[{0,0},{Length@keys}];
   
   Do[
      result[[n,1]]=data[[n,1]];
      If[data[[n,5]]==True,(*i.e. defective eigenvalue*)
          (*check if geometric multiplicity=1 or not. These cases are handled differently*)
          
          If[data[[n,3]]==1,(*standard case, where one normal vector is found*)
          
              (*geometric multiplicity =1 case*)
              setOfV=handleCaseOne[data,mat,V,n];
              Print["Hence the set of eigenvectors associated with eigenvalue ",data[[n,1]]," is ",
                    Style[MatrixForm[Transpose@setOfV],Bold]
              ]
          ,
              (*geometric multiplicity >1 case*)
              setOfV=handleCaseTwo[data,mat,V,n];
              Print["Hence the set of eigenvectors associated with eigenvalue ",data[[n,1]]," is ",
                    Style[MatrixForm[Transpose@setOfV],Bold]
              ]
          ];    
          result[[n,2]]=setOfV
     ,
        Print["Since eigenvalue ",data[[n,1]]," is not defective, then we already have its normal "];
        Print["eigenvectors from the table above,  which was found by Eigensystem. Nothing to do."];
        result[[n,2]]=data[[n,4]]
    ]
    ,
    {n,Length@keys}
  ];
  
  result
]

mat={{1,2,0},{1,1,2},{0,-1,1}};
mat={{39,8,-16},{-36,-5,16},{72,16,-29}};

mat={{5,1,-2,4},{0,5,2,2},{0,0,5,3},{0,0,0,4}};
mat={{0,1,-1,1},{0,0,0,0},{0,-1,2,2},{0,1,-2,2}};
mat={{3,1,0,0},{0,3,1,0},{0,0,3,0},{0,0,0,3}};
mat={{0,1,-1,1},{0,0,0,0},{0,-1,2,2},{0,1,-2,2}};
mat={{3,1,0,0},{0,3,1,0},{0,0,3,0},{0,0,0,3}};
mat={{0,1,-1,1},{0,0,0,0},{0,-1,2,2},{0,1,-2,2}};
mat={{1,0,0,0,0},{3,1,0,0,0},{6,3,2,0,0},{10,6,3,2,0},{15,10,6,3,2}};
mat={{3,1,0,0},{0,3,1,0},{0,0,3,0},{0,0,0,3}};
(*4*)
mat={{4,1,1},{-2,1,-2},{1,1,4}};
mat={{5,3,2},{-1,8,1},{-1,2,8}};
mat={{3,1,0,0},{0,3,1,0},{0,0,3,0},{0,0,0,3}};
mat={{1,0,0,0,0},{1,-1,0,0,-1},{1,-1,0,0,-1},{0,0,0,0,-1},{-1,1,0,0,1}};

mat={{3,1,0,0},{0,3,1,0},{0,0,3,0},{0,0,0,3}};(*correct*)
mat={{1,0,0,0,0},{1,-1,0,0,-1},{1,-1,0,0,-1},{0,0,0,0,-1},{-1,1,0,0,1}};(*wrong*)
mat={{4,1,1},{-2,1,-2},{1,1,4}};
result=findV[mat,V];