From 3e60500db16a815a48a06952901133b573c5e581 Mon Sep 17 00:00:00 2001
From: ShilohSanity <51039238+ShilohSanity@users.noreply.github.com>
Date: Tue, 5 May 2020 15:22:40 +0800
Subject: [PATCH] Drazin-inverse & sympy

This is the code to caculate the Drazin inverse of matrix A. Be careful, k is the Drazin index of matrix A.
---
 Drazin-inverse & sympy | 27 +++++++++++++++++++++++++++
 1 file changed, 27 insertions(+)
 create mode 100644 Drazin-inverse & sympy

diff --git a/Drazin-inverse & sympy b/Drazin-inverse & sympy
new file mode 100644
index 0000000..3738979
--- /dev/null
+++ b/Drazin-inverse & sympy	
@@ -0,0 +1,27 @@
+import sympy as sp
+def Shiloh_GI1(A,k):
+    n=A.shape[1]
+    p=1                                      
+    F=sp.eye(n)  
+    Y=sp.eye(n)
+    while(k!=0):
+        Y=Y*A
+        k=k-1
+    YA=Y*A                                   
+    j=1
+    while p!=0:                              
+        YAF=YA*F
+        p=sp.trace(YAF)/j                    
+        A_pinv=F*Y/p                          
+        F=YAF-p*sp.eye(n)                    
+        j+=1
+        p=sp.trace(YA*F)/j                   
+    return A_pinv
+def Shiloh_GI(A,k):
+    (m,n)=A.shape
+    if m>=n:
+        A_pinv,DIM= Shiloh_GI1(A,k)
+    else:
+        AH_pinv,DIM= Shiloh_GI1(A.H,k)
+        A_pinv=AH_pinv.H
+    return A_pinv,DIM