Skip to content

Commit f6e6fb4

Browse files
cgcgcgcgcgcg
committed
Fix for transpose operators
1 parent cdee592 commit f6e6fb4

File tree

3 files changed

+90
-15
lines changed

3 files changed

+90
-15
lines changed

‎base/PyNucleus_base/LinearOperator_{SCALAR}.pxi

+6-5
Original file line numberDiff line numberDiff line change
@@ -651,14 +651,12 @@ cdef class {SCALAR_label}Transpose_Linear_Operator({SCALAR_label}LinearOperator)
651651
cdef INDEX_t matvec(self,
652652
{SCALAR}_t[::1] x,
653653
{SCALAR}_t[::1] y) except -1:
654-
self.A.matvecTrans(x, y)
655-
return 0
654+
return self.A.matvecTrans(x, y)
656655

657656
cdef INDEX_t matvec_no_overwrite(self,
658657
{SCALAR}_t[::1] x,
659658
{SCALAR}_t[::1] y) except -1:
660-
self.A.matvecTrans_no_overwrite(x, y)
661-
return 0
659+
return self.A.matvecTrans_no_overwrite(x, y)
662660

663661
def isSparse(self):
664662
return self.A.isSparse()
@@ -677,7 +675,10 @@ cdef class {SCALAR_label}Transpose_Linear_Operator({SCALAR_label}LinearOperator)
677675
return Bcsr
678676

679677
def toarray(self):
680-
return self.A.transpose().toarray()
678+
try:
679+
return self.A.transpose().toarray()
680+
except AttributeError:
681+
return np.ascontiguousarray(self.A.toarray().T)
681682

682683
def get_diagonal(self):
683684
return np.array(self.A.diagonal, copy=False)

‎base/PyNucleus_base/linear_operators.pyx

+83-9
Original file line numberDiff line numberDiff line change
@@ -1213,26 +1213,57 @@ cdef class sumMultiplyOperator(LinearOperator):
12131213
cdef:
12141214
INDEX_t i
12151215
LinearOperator op
1216+
int ret
12161217
op = self.ops[0]
1217-
op.matvec(x, y)
1218+
ret = op.matvec(x, y)
12181219
scaleScalar(y, self.coeffs[0])
12191220
for i in range(1, self.coeffs.shape[0]):
12201221
op = self.ops[i]
1221-
op.matvec(x, self.z)
1222+
ret = min(ret, op.matvec(x, self.z))
12221223
assign3(y, y, 1.0, self.z, self.coeffs[i])
1223-
return 0
1224+
return ret
12241225

12251226
cdef INDEX_t matvec_no_overwrite(self,
12261227
REAL_t[::1] x,
12271228
REAL_t[::1] y) except -1:
1229+
cdef:
1230+
INDEX_t i
1231+
LinearOperator op = 0
1232+
int ret = 0
1233+
for i in range(self.coeffs.shape[0]):
1234+
op = self.ops[i]
1235+
ret = min(op.matvec(x, self.z), ret)
1236+
assign3(y, y, 1.0, self.z, self.coeffs[i])
1237+
return ret
1238+
1239+
cdef INDEX_t matvecTrans(self,
1240+
REAL_t[::1] x,
1241+
REAL_t[::1] y) except -1:
12281242
cdef:
12291243
INDEX_t i
12301244
LinearOperator op
1245+
int ret
1246+
op = self.ops[0]
1247+
ret = op.matvecTrans(x, y)
1248+
scaleScalar(y, self.coeffs[0])
1249+
for i in range(1, self.coeffs.shape[0]):
1250+
op = self.ops[i]
1251+
ret = min(ret, op.matvecTrans(x, self.z))
1252+
assign3(y, y, 1.0, self.z, self.coeffs[i])
1253+
return ret
1254+
1255+
cdef INDEX_t matvecTrans_no_overwrite(self,
1256+
REAL_t[::1] x,
1257+
REAL_t[::1] y) except -1:
1258+
cdef:
1259+
INDEX_t i
1260+
LinearOperator op
1261+
int ret = 0
12311262
for i in range(self.coeffs.shape[0]):
12321263
op = self.ops[i]
1233-
op.matvec(x, self.z)
1264+
ret = min(ret, op.matvecTrans(x, self.z))
12341265
assign3(y, y, 1.0, self.z, self.coeffs[i])
1235-
return 0
1266+
return ret
12361267

12371268
def toarray(self):
12381269
return sum([c*op.toarray() for c, op in zip(self.coeffs, self.ops)])
@@ -1436,8 +1467,34 @@ cdef class multiIntervalInterpolationOperator(LinearOperator):
14361467
interpolationOperator op
14371468
assert self.selected != -1
14381469
op = self.ops[self.selected]
1439-
op.matvec(x, y)
1440-
return 0
1470+
return op.matvec(x, y)
1471+
1472+
cdef INDEX_t matvec_no_overwrite(self,
1473+
REAL_t[::1] x,
1474+
REAL_t[::1] y) except -1:
1475+
cdef:
1476+
interpolationOperator op
1477+
assert self.selected != -1
1478+
op = self.ops[self.selected]
1479+
return op.matvec_no_overwrite(x, y)
1480+
1481+
cdef INDEX_t matvecTrans(self,
1482+
REAL_t[::1] x,
1483+
REAL_t[::1] y) except -1:
1484+
cdef:
1485+
interpolationOperator op
1486+
assert self.selected != -1
1487+
op = self.ops[self.selected]
1488+
return op.matvecTrans(x, y)
1489+
1490+
cdef INDEX_t matvecTrans_no_overwrite(self,
1491+
REAL_t[::1] x,
1492+
REAL_t[::1] y) except -1:
1493+
cdef:
1494+
interpolationOperator op
1495+
assert self.selected != -1
1496+
op = self.ops[self.selected]
1497+
return op.matvecTrans_no_overwrite(x, y)
14411498

14421499
def toarray(self):
14431500
assert self.selected != -1
@@ -1521,8 +1578,25 @@ cdef class delayedConstructionOperator(LinearOperator):
15211578
REAL_t[::1] x,
15221579
REAL_t[::1] y) except -1:
15231580
self.assure_constructed()
1524-
self.A.matvec(x, y)
1525-
return 0
1581+
return self.A.matvec(x, y)
1582+
1583+
cdef INDEX_t matvec_no_overwrite(self,
1584+
REAL_t[::1] x,
1585+
REAL_t[::1] y) except -1:
1586+
self.assure_constructed()
1587+
return self.A.matvec_no_overwrite(x, y)
1588+
1589+
cdef INDEX_t matvecTrans(self,
1590+
REAL_t[::1] x,
1591+
REAL_t[::1] y) except -1:
1592+
self.assure_constructed()
1593+
return self.A.matvecTrans(x, y)
1594+
1595+
cdef INDEX_t matvecTrans_no_overwrite(self,
1596+
REAL_t[::1] x,
1597+
REAL_t[::1] y) except -1:
1598+
self.assure_constructed()
1599+
return self.A.matvecTrans_no_overwrite(x, y)
15261600

15271601
def toarray(self):
15281602
self.assure_constructed()

‎base/PyNucleus_base/utilsFem.py

+1-1
Original file line numberDiff line numberDiff line change
@@ -646,7 +646,7 @@ def diff(self, d):
646646
result[p.label] = (p.value, d[p.label])
647647
elif isinstance(p.value, (int, INDEX, REAL, float)):
648648
if not np.allclose(p.value, d[p.label],
649-
rtol=rTol, atol=aTol):
649+
rtol=rTol, atol=aTol) and not (np.isnan(p.value) and np.isnan(d[p.label])):
650650
print(p.label, p.value, d[p.label], rTol, aTol, p.rTol, p.aTol)
651651
result[p.label] = (p.value, d[p.label])
652652
else:

0 commit comments

Comments
 (0)