因为taichi 中只有3*3的特征值分解，我需要做更大的矩阵的特征值分解。就自己写了个jacob特征值算法。结果遇到了一个奇怪的问题，
import taichi as ti
ti.init(ti.cpu,dynamic_index=True)
@ti.func
def maxElem(a): # Find largest off-diag. element a[k,l]
aMax = -1.0
k=0
l=0
for i in range(a.n-1):
for j in range(i+1,a.n):
if ti.abs(a[i,j]) >= aMax:
aMax = ti.abs(a[i,j])
k = i
l = j
return aMax,k,l
@ti.func
def jacob_eigen_test(a:ti.template()):
maxRot = 5*(a.n**2) # Set limit on number of rotations
p = ti.math.eye(a.n)
tol = 1.0e-7
sig = ti.Vector.zero(ti.f32,a.n)
aMax = 1.0
t = 0.0
N = 0
# for _ in range(1):
print('p1',p[0,0])
while aMax > tol:
# while True:
print('p2',p[0,0])
aMax,k,l = maxElem(a)
# print(aMax,k,l)
aDiff = a[l,l] - a[k,k]
# print(aDiff)
if ti.abs(a[k,l]) < ti.abs(aDiff)*1.0e-12:
t = a[k,l]/aDiff
else:
phi = aDiff/(2.0*a[k,l])
# print(phi)
t = 1.0/(ti.abs(phi) + ti.sqrt(phi**2 + 1.0))
if phi < 0.0:
t = -t
# print('pppp')
c = 1.0/ti.sqrt(t**2 + 1.0); s = t*c
tau = s/(1.0 + c)
temp = a[k,l]
a[k,l] = 0.0
a[k,k] = a[k,k] - t*temp
a[l,l] = a[l,l] + t*temp
for i in range(k): # Case of i < k
temp1 = a[i,k]
a[i,k] = temp1 - s*(a[i,l] + tau*temp1)
a[i,l] = a[i,l] + s*(temp1 - tau*a[i,l])
for i in range(k+1,l): # Case of k < i < l
temp2 = a[k,i]
a[k,i] = temp2 - s*(a[i,l] + tau*a[k,i])
a[i,l] = a[i,l] + s*(temp2 - tau*a[i,l])
for i in range(l+1,a.n): # Case of i > l
temp3 = a[k,i]
a[k,i] = temp3 - s*(a[l,i] + tau*temp3)
a[l,i] = a[l,i] + s*(temp3 - tau*a[l,i])
for i in range(a.n): # Update transformation matrix
temp4 = p[i,k]
p[i,k] = temp4 - s*(p[i,l] + tau*p[i,k])
p[i,l] = p[i,l] + s*(temp4 - tau*p[i,l])
N+=1
if N>maxRot:
break
# print(p[0,0],p[0,1],p[0,2],p[0,3])
for ii in range(a.n):
sig[ii] = a[ii,ii]
return sig, p
@ti.kernel
def test():
test_S =ti.Matrix([[1.0,-1.0,2.0,3.0],[-1.0,1.0,-1.5,-0.5],[2.0,-1.5,1.0,1.7],[3.0,-0.5,1.7,1.0]],ti.f32)
Sig, P = jacob_eigen_test(test_S)
# print(Sig,P)
@ti.kernel
def test2():
test_S =ti.Matrix([[1.0,-1.0,2.0,3.0],[-1.0,1.0,-1.5,-0.5],[2.0,-1.5,1.0,1.7],[3.0,-0.5,1.7,1.0]],ti.f32)
Sig, P = jacob_eigen_test(test_S)
print(Sig,P)
test()
print('test2')
test2()
test函数的显示
p1 1.000000
p2 1.000000
test2函数显示
p1 1.000000
p2 0.000000
完全同样的代码，就加了一个print，为什么中间的变量值会改变。
另外，我如果使用f64的就不会出现这个问题。将while换成for循环也不会出现这个问题。
卡了我一下午，实在是无法理解。