Skip to content

Commit

Permalink
replace from_data by as_source/range_array or make_array where approp…
Browse files Browse the repository at this point in the history 
…riate
  • Loading branch information
@sdrave
sdrave committed Dec 2, 2016
1 parent 0a16f69 commit c331d14
Show file tree
Hide file tree
Showing 2 changed files with 21 additions and 33 deletions.
46 changes: 17 additions & 29 deletions src/pymor/discretizations/iosys.py
View file
Expand Up @@ -421,16 +421,14 @@ def eval_tf(self, s):

sEmA = LincombOperator((E, A), (s, -1))
if self.m <= self.p:
I_m = B.source.from_data(sp.eye(self.m))
tfs = C.apply(sEmA.apply_inverse(B.apply(I_m))).data.T
tfs = C.apply(sEmA.apply_inverse(B.as_range_array())).data.T
else:
I_p = C.range.from_data(sp.eye(self.p))
tfs = B.apply_adjoint(sEmA.apply_inverse_adjoint(C.apply_adjoint(I_p))).data.conj()
tfs = B.apply_adjoint(sEmA.apply_inverse_adjoint(C.as_source_array())).data.conj()
if not isinstance(D, ZeroOperator):
if self.m <= self.p:
tfs += D.apply(I_m).data.T
tfs += D.as_range_array().data.T
else:
tfs += D.apply_adjoint(I_p).data
tfs += D.as_source_array().data
return tfs

def eval_dtf(self, s):
Expand Down Expand Up @@ -463,12 +461,10 @@ def eval_dtf(self, s):

sEmA = LincombOperator((E, A), (s, -1))
if self.m <= self.p:
I_m = B.source.from_data(sp.eye(self.m))
dtfs = -C.apply(sEmA.apply_inverse(E.apply(sEmA.apply_inverse(B.apply(I_m))))).data.T
dtfs = -C.apply(sEmA.apply_inverse(E.apply(sEmA.apply_inverse(B.as_range_array())))).data.T
else:
I_p = C.range.from_data(sp.eye(self.p))
dtfs = B.apply_adjoint(sEmA.apply_inverse_adjoint(E.apply_adjoint(sEmA.apply_inverse_adjoint(
C.apply_adjoint(I_p))))).data.conj()
C.as_source_array())))).data.conj()
return dtfs

@cached
Expand Down Expand Up @@ -760,13 +756,11 @@ def eval_tf(self, s):

s2MpsDpK = LincombOperator((M, D, K), (s ** 2, s, 1))
if self.m <= self.p:
I_m = B.source.from_data(sp.eye(self.m))
CppsCv = LincombOperator((Cp, Cv), (1, s))
tfs = CppsCv.apply(s2MpsDpK.apply_inverse(B.apply(I_m))).data.T
tfs = CppsCv.apply(s2MpsDpK.apply_inverse(B.as_range_array())).data.T
else:
I_p = Cp.range.from_data(sp.eye(self.p))
tfs = B.apply_adjoint(s2MpsDpK.apply_inverse_adjoint(Cp.apply_adjoint(I_p) +
Cv.apply_adjoint(I_p) * s.conj())).data.conj()
tfs = B.apply_adjoint(s2MpsDpK.apply_inverse_adjoint(Cp.as_source_array() +
Cv.as_source_array() * s.conj())).data.conj()
return tfs

def eval_dtf(self, s):
Expand Down Expand Up @@ -804,15 +798,13 @@ def eval_dtf(self, s):
s2MpsDpK = LincombOperator((M, D, K), (s ** 2, s, 1))
sM2pD = LincombOperator((M, D), (2 * s, 1))
if self.m <= self.p:
I_m = B.source.from_data(sp.eye(self.m))
dtfs = Cv.apply(s2MpsDpK.apply_inverse(B.apply(I_m))).data.T * s
dtfs = Cv.apply(s2MpsDpK.apply_inverse(B.as_range_array())).data.T * s
CppsCv = LincombOperator((Cp, Cv), (1, s))
dtfs -= CppsCv.apply(s2MpsDpK.apply_inverse(sM2pD.apply(s2MpsDpK.apply_inverse(B.apply(I_m))))).data.T
dtfs -= CppsCv.apply(s2MpsDpK.apply_inverse(sM2pD.apply(s2MpsDpK.apply_inverse(B.as_range_array())))).data.T
else:
I_p = Cp.range.from_data(sp.eye(self.p))
dtfs = B.apply_adjoint(s2MpsDpK.apply_inverse_adjoint(Cv.apply_adjoint(I_m))).data.conj() * s
dtfs = B.apply_adjoint(s2MpsDpK.apply_inverse_adjoint(Cv.as_source_array())).data.conj() * s
dtfs -= B.apply_adjoint(s2MpsDpK.apply_inverse_adjoint(sM2pD.apply_adjoint(s2MpsDpK.apply_inverse_adjoint(
Cp.apply_adjoint(I_p) + Cv.apply_adjoint(I_p) * s.conj())))).data.conj()
Cp.as_source_array() + Cv.as_source_array() * s.conj())))).data.conj()
return dtfs


Expand Down Expand Up @@ -928,11 +920,9 @@ def eval_tf(self, s):

middle = LincombOperator((E, A) + Ad, (s, -1) + tuple(-np.exp(-taui * s) for taui in self.tau))
if self.m <= self.p:
I_m = B.source.from_data(sp.eye(self.m))
tfs = C.apply(middle.apply_inverse(B.apply(I_m))).data.T
tfs = C.apply(middle.apply_inverse(B.as_range_array())).data.T
else:
I_p = C.range.from_data(sp.eye(self.p))
tfs = B.apply_adjoint(middle.apply_inverse_adjoint(C.apply_adjoint(I_p))).data.conj()
tfs = B.apply_adjoint(middle.apply_inverse_adjoint(C.as_source_array())).data.conj()
return tfs

def eval_dtf(self, s):
Expand Down Expand Up @@ -969,13 +959,11 @@ def eval_dtf(self, s):
left_and_right = LincombOperator((E, A) + Ad, (s, -1) + tuple(-np.exp(-taui * s) for taui in self.tau))
middle = LincombOperator((E,) + Ad, (s,) + tuple(taui * np.exp(-taui * s) for taui in self.tau))
if self.m <= self.p:
I_m = B.source.from_data(sp.eye(self.m))
dtfs = C.apply(left_and_right.apply_inverse(middle.apply(left_and_right.apply_inverse(
B.apply(I_m))))).data.T
B.as_range_array())))).data.T
else:
I_p = C.range.from_data(sp.eye(self.p))
dtfs = B.apply_adjoint(left_and_right.apply_inverse_adjoint(middle.apply_adjoint(
left_and_right.apply_inverse_adjointi(C.apply_adjoint(I_p))))).data.conj()
left_and_right.apply_inverse_adjointi(C.as_source_array())))).data.conj()
return dtfs


Expand Down
8 changes: 4 additions & 4 deletions src/pymor/reductors/lti.py
View file
Expand Up @@ -177,9 +177,9 @@ def irka(discretization, r, sigma=None, b=None, c=None, tol=1e-4, maxit=100, ver
if sigma is None:
sigma = np.logspace(-1, 1, r)
if b is None:
b = discretization.B.source.from_data(np.ones((r, discretization.m)))
b = discretization.B.source.make_array(np.ones((r, discretization.m)))
if c is None:
c = discretization.C.range.from_data(np.ones((r, discretization.p)))
c = discretization.C.range.make_array(np.ones((r, discretization.p)))

if verbose:
if compute_errors:
Expand Down Expand Up @@ -248,8 +248,8 @@ def irka(discretization, r, sigma=None, b=None, c=None, tol=1e-4, maxit=100, ver
else:
print('{:4d} | {:15.9e}'.format(it + 1, dist[-1]))

Y = rom.B.range.from_data(Y.conj().T)
X = rom.C.source.from_data(X.T)
Y = rom.B.range.make_array(Y.conj().T)
X = rom.C.source.make_array(X.T)
b = rom.B.apply_adjoint(Y)
c = rom.C.apply(X)
R.append(b)
Expand Down

0 comments on commit c331d14

Please sign in to comment.