我可能遇到了一个bug

今天打算学习一下利用牛顿法解线性方程组,于是考虑着先用高斯消元法实现一下,在一步一步弄成牛顿的。

这段代码产生了这样一个错误
UnboundLocalError: local variable ‘j’ referenced before assignment

我想这我的 j 也定义了呀,找了半天发现只要把我15行的 range(i + 1,3)改成range(i,3)就可以运行了。
我想倘若是语法不允许这样做那报错信息也不应该跟 j 有关,所以我把我这个“奇丑无比”的代码丢上来请大佬看一下…
kernel里的if(true)是希望它单线程运行,因为不是很熟悉并行计算,就先这样。


这个文档编辑器我好像不太会用 :confounded:
麻烦各位 …将就着…看一下
看来我确实得花点时间学markdown了

  • 据我对并行计算的了解
  • taichi使用static时,就相当于在hlsl里面使用unroll关键字
  • [https://docs.microsoft.com/zh-tw/previous-versions/windows/xna/bb313989(v=xnagamestudio.10)]
  • 如果你对C了解的话,这个展开就是在编译阶段,拿替代,相当于把你的代码拷贝3份,i被0 1 2替代即可
  • 第二层中i不存在了,变成了0 1 2,但是必须保证range里面是静态变量
  • i+1操作,相当于要有一个缓存,这里是j=i+1
  • 在编译器层面引入了一个local变量,那么下面就展不开了
  • 而j=i,那么j也是静态变量,就能展开了

···
个人理解,仅供参考
我之前也遇到过这个问题
编译器展开不了,就自己手动展开一下咯,人肉拷贝三份,先把外层展开了
如果有懂编译器的朋友也可以发表下见解

贴上一段代码,分解对角矩阵用的
注释里面是和你一样想用static展开的过程,但是引入了局部变量就不行了
所以人肉展开了


@ti.func
def LUPDecompose(A) :
    Tol = 1e-15
    ret = 1
    P   = ti.Vector([0,1,2,3])

    '''
    for i in ti.static(range(3)):
        maxA = 0.0
        imax = i

        for k in ti.static(range(i,3)):
            absA = ti.abs(A[k,i])
            if absA > maxA:
                maxA = absA
                imax = k

        if (maxA < Tol): 
            ret = 0
            break

        if (imax != i) :
            P = swap_V(P, i, imax)
            A = swap_M(A, i, imax)
            P[3] += 1

        for j in ti.static(range(i+1,3)):
            A[j,i] /= A[i,i]
            for k in ti.static(range(i+1,3)):
                A[j,k] -= A[j,i] * A[i,k]
    '''

    ########### i =0 ##############################
    maxA = 0.0
    imax = 0

    for k in ti.static(range(3)):
        absA = ti.abs(A[k,0])
        if absA > maxA:
            maxA = absA
            imax = k

    if (maxA < Tol): 
        ret = 0

    if (imax != 0) :
        if imax == 1:
            P = ti.Vector([P[1],P[0],P[2],P[3]])
            A = ti.Matrix([[A[1,0], A[1,1], A[1,2]], [A[0,0], A[0,1], A[0,2]], [A[2,0], A[2,1], A[2,2]]]) 
        else:#IMAX = 2
            P = ti.Vector([P[2],P[1],P[0],P[3]])
            A = ti.Matrix([[A[2,0], A[2,1], A[2,2]],[A[1,0], A[1,1], A[1,2]],[A[0,0], A[0,1], A[0,2]]]) 
        P[3] += 1

    # j = 1
    A[1,0] /= A[0,0]
    A[1,1] -= A[1,0] * A[0,1]
    A[1,2] -= A[1,0] * A[0,2]

    # j = 2
    A[2,0] /= A[0,0]
    A[2,1] -= A[2,0] * A[0,1]
    A[2,2] -= A[2,0] * A[0,2]

    
    ########### i =1 ##############################
    maxA = 0.0
    imax = 1

    for k in ti.static(range(1,3)):
        absA = ti.abs(A[k,0])
        if absA > maxA:
            maxA = absA
            imax = k

    if (maxA < Tol): 
        ret = 0

    if (imax != 1) :#IMAX = 2
        P = ti.Vector([P[0],P[2],P[1],P[3]])
        A = ti.Matrix([[A[0,0], A[0,1], A[0,2]],[A[2,0], A[2,1], A[2,2]],[A[1,0], A[1,1], A[1,2]]]) 
        P[3] += 1

    # j = 2
    A[2,1] /= A[1,1]
    A[2,2] -= A[2,1] * A[1,2]

    #print(A)

    ########### i =2 ##############################
    if (ti.abs(A[2,2]) < Tol): 
        ret = 0

    
    B = A

    '''
    for k in ti.static(range(3)):
        if P[k] == 0:
            B[0, 0] = A[k,0] 
            B[0, 1] = A[k,1] 
            B[0, 2] = A[k,2] 
        elif P[k] == 1:
            B[1, 0] = A[k,0] 
            B[1, 1] = A[k,1] 
            B[1, 2] = A[k,2] 
        else:
            B[2, 0] = A[k,0] 
            B[2, 1] = A[k,1] 
            B[2, 2] = A[k,2] 
    '''
    return ret,B,P

感谢大佬回复,虽然还是有点小懵逼,不过已经基本了解问题所在了,果然我还有好多东西需要学习呀 :laughing:

1 个赞