贴上一段代码,分解对角矩阵用的
注释里面是和你一样想用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