在kernel中print变量为什么会改变func中的变量的值

因为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循环也不会出现这个问题。
卡了我一下午,实在是无法理解。

1 Like

看起来是个 bug,我这里复现了一下,打印的时候写 print(Sig) 或者 print(P) 都是正常的,两个写在一起 print(Sig, P) 会有这个问题。

不过 32 位精度下面由于浮点误差你的 tol 设置成 1e-7 然后比较已经失去意义了。

1 Like

代码最后是用了64位精度的,tol忘了改过去了。