@@ -10,27 +10,33 @@ struct PlusOperator{T,BW,SZ,O<:Operator{T},BBW,SBBW} <: Operator{T}
10
10
sz:: SZ
11
11
blockbandwidths:: BBW
12
12
subblockbandwidths:: SBBW
13
+ isbandedblockbanded:: Bool
14
+ israggedbelow:: Bool
13
15
14
16
function PlusOperator {T,BW,SZ,O,BBW,SBBW} (opsin:: Vector{O} , bw:: BW ,
15
- sz:: SZ , bbw:: BBW , sbbw:: SBBW ) where {T,O<: Operator{T} ,BW,SZ,BBW,SBBW}
17
+ sz:: SZ , bbw:: BBW , sbbw:: SBBW ,
18
+ ibbb:: Bool , irb:: Bool ) where {T,O<: Operator{T} ,BW,SZ,BBW,SBBW}
16
19
all (x -> size (x) == sz, opsin) || throw (" sizes of operators are incompatible" )
17
- new {T,BW,SZ,O,BBW,SBBW} (opsin, bw, sz, bbw,sbbw)
20
+ new {T,BW,SZ,O,BBW,SBBW} (opsin, bw, sz, bbw,sbbw, ibbb, irb )
18
21
end
19
22
end
20
23
21
24
bandwidthsmax (ops, f= bandwidths) = mapreduce (f, (t1, t2) -> max .(t1, t2), ops, init= (- 720 , - 720 )) #= approximate (-∞,-∞) =#
22
25
23
- function PlusOperator (opsin :: Vector{O} , args... ) where {O<: Operator }
24
- PlusOperator {eltype(O)} (opsin , args... )
26
+ function PlusOperator (ops :: Vector{O} , args... ) where {O<: Operator }
27
+ PlusOperator {eltype(O)} (ops , args... )
25
28
end
26
- function PlusOperator {ET} (opsin:: Vector{O} ,
27
- bw:: Tuple{Any,Any} = bandwidthsmax (opsin),
28
- sz:: Tuple{Any,Any} = size (first (opsin)),
29
- bbw:: Tuple{Any,Any} = bandwidthsmax (opsin, blockbandwidths),
30
- sbbw:: Tuple{Any,Any} = bandwidthsmax (opsin, subblockbandwidths),
29
+ function PlusOperator {ET} (ops:: Vector{O} ,
30
+ bw:: Tuple{Any,Any} = bandwidthsmax (ops),
31
+ sz:: Tuple{Any,Any} = size (first (ops)),
32
+ bbw:: Tuple{Any,Any} = bandwidthsmax (ops, blockbandwidths),
33
+ sbbw:: Tuple{Any,Any} = bandwidthsmax (ops, subblockbandwidths),
34
+ ibbb:: Bool = all (isbandedblockbanded, ops),
35
+ irb:: Bool = all (israggedbelow, ops),
31
36
) where {ET,O<: Operator{ET} }
32
37
33
- PlusOperator {ET,typeof(bw),typeof(sz),O,typeof(bbw),typeof(sbbw)} (opsin, bw, sz, bbw, sbbw)
38
+ PlusOperator {ET,typeof(bw),typeof(sz),O,typeof(bbw),typeof(sbbw)} (
39
+ ops, bw, sz, bbw, sbbw, ibbb, irb)
34
40
end
35
41
36
42
for (OP, mn) in ((:colstart , :min ), (:colstop , :max ), (:rowstart , :min ), (:rowstop , :max ))
@@ -57,7 +63,9 @@ function convert(::Type{Operator{T}}, P::PlusOperator) where {T}
57
63
ops = P. ops
58
64
PlusOperator (eltype (ops) <: Operator{T} ? ops :
59
65
_convertops (Operator{T}, ops),
60
- bandwidths (P), size (P), blockbandwidths (P), subblockbandwidths (P)):: Operator{T}
66
+ bandwidths (P), size (P), blockbandwidths (P),
67
+ subblockbandwidths (P), isbandedblockbanded (P),
68
+ israggedbelow (P)):: Operator{T}
61
69
end
62
70
end
63
71
@@ -221,9 +229,12 @@ struct TimesOperator{T,BW,SZ,O<:Operator{T},BBW,SBBW} <: Operator{T}
221
229
sz:: SZ
222
230
blockbandwidths:: BBW
223
231
subblockbandwidths:: SBBW
232
+ isbandedblockbanded:: Bool
233
+ israggedbelow:: Bool
224
234
225
235
function TimesOperator {T,BW,SZ,O,BBW,SBBW} (ops:: Vector{O} , bw:: BW ,
226
- sz:: SZ , bbw:: BBW , sbbw:: SBBW ) where {T,O<: Operator{T} ,BW,SZ,BBW,SBBW}
236
+ sz:: SZ , bbw:: BBW , sbbw:: SBBW ,
237
+ ibbb:: Bool , irb:: Bool ) where {T,O<: Operator{T} ,BW,SZ,BBW,SBBW}
227
238
# check compatible
228
239
for k = 1 : length (ops)- 1
229
240
size (ops[k], 2 ) == size (ops[k+ 1 ], 1 ) || throw (ArgumentError (" incompatible operator sizes" ))
@@ -241,7 +252,7 @@ struct TimesOperator{T,BW,SZ,O<:Operator{T},BBW,SBBW} <: Operator{T}
241
252
newops = ops
242
253
end
243
254
244
- new {T,BW,SZ,O,BBW,SBBW} (newops, bw, sz, bbw, sbbw)
255
+ new {T,BW,SZ,O,BBW,SBBW} (newops, bw, sz, bbw, sbbw, ibbb, irb )
245
256
end
246
257
end
247
258
@@ -260,18 +271,22 @@ function TimesOperator(ops::AbstractVector{O},
260
271
sz:: Tuple{Any,Any} = _timessize (ops),
261
272
bbw:: Tuple{Any,Any} = bandwidthssum (ops, blockbandwidths),
262
273
sbbw:: Tuple{Any,Any} = bandwidthssum (ops, subblockbandwidths),
263
- ) where {T,O<: Operator{T} }
264
- TimesOperator {T,typeof(bw),typeof(sz),O,typeof(bbw),typeof(sbbw)} (convert_vector (ops),
265
- bw, sz, bbw, sbbw)
274
+ ibbb:: Bool = all (isbandedblockbanded, ops),
275
+ irb:: Bool = all (israggedbelow, ops),
276
+ ) where {O<: Operator }
277
+ TimesOperator {eltype(O),typeof(bw),typeof(sz),O,typeof(bbw),typeof(sbbw)} (
278
+ convert_vector (ops), bw, sz, bbw, sbbw, ibbb, irb)
266
279
end
267
280
268
281
_extractops (A:: TimesOperator , :: typeof (* )) = A. ops
269
282
270
283
function TimesOperator (A:: Operator , B:: Operator )
271
284
v = collateops (* , A, B)
285
+ ibbb = all (isbandedblockbanded, (A, B))
286
+ irb = all (israggedbelow, (A, B))
272
287
TimesOperator (convert_vector (v), _bandwidthssum (A, B), _timessize ((A, B)),
273
288
_bandwidthssum (A, B, blockbandwidths),
274
- _bandwidthssum (A, B, subblockbandwidths))
289
+ _bandwidthssum (A, B, subblockbandwidths), ibbb, irb )
275
290
end
276
291
277
292
@@ -285,7 +300,8 @@ function convert(::Type{Operator{T}}, P::TimesOperator) where {T}
285
300
TimesOperator (eltype (ops) <: Operator{T} ? ops :
286
301
_convertops (Operator{T}, ops),
287
302
bandwidths (P), size (P), blockbandwidths (P),
288
- subblockbandwidths (P)):: Operator{T}
303
+ subblockbandwidths (P), isbandedblockbanded (P),
304
+ israggedbelow (P)):: Operator{T}
289
305
end
290
306
end
291
307
@@ -298,14 +314,13 @@ end
298
314
@inline function _promotetimes (opsin,
299
315
dsp= domainspace (last (opsin)),
300
316
sz= _timessize (opsin),
301
- anytimesop= true ,
302
- )
317
+ anytimesop= true )
303
318
304
319
@assert length (opsin) > 1 " need at least 2 operators"
305
- ops, bw, bbw, sbbw = __promotetimes (opsin, dsp, anytimesop)
306
- TimesOperator (ops, bw, sz, bbw, sbbw)
320
+ ops, bw, bbw, sbbw, ibbb, irb = __promotetimes (opsin, dsp, anytimesop)
321
+ TimesOperator (ops, bw, sz, bbw, sbbw, ibbb, irb )
307
322
end
308
- @inline function __promotetimes (opsin, dsp, anytimesop)
323
+ function __promotetimes (opsin, dsp, anytimesop)
309
324
ops = Vector {Operator{promote_eltypeof(opsin)}} (undef, 0 )
310
325
sizehint! (ops, length (opsin))
311
326
322
337
end
323
338
end
324
339
reverse! (ops), bandwidthssum (ops), bandwidthssum (ops, blockbandwidths),
325
- bandwidthssum (ops, subblockbandwidths)
340
+ bandwidthssum (ops, subblockbandwidths), all (isbandedblockbanded, ops),
341
+ all (israggedbelow, ops)
342
+ end
343
+ @inline function _op_bws (op)
344
+ [op], bandwidths (op), blockbandwidths (op),
345
+ subblockbandwidths (op), isbandedblockbanded (op),
346
+ israggedbelow (op)
326
347
end
327
- _op_bws (op) = [op], bandwidths (op), blockbandwidths (op), subblockbandwidths (op)
328
348
@inline function __promotetimes (opsin:: Tuple{Operator,Operator} , dsp, anytimesop)
329
349
@assert ! any (Base. Fix2 (isa, TimesOperator), opsin) " TimesOperator should have been extracted already"
330
350
@@ -343,9 +363,12 @@ _op_bws(op) = [op], bandwidths(op), blockbandwidths(op), subblockbandwidths(op)
343
363
else
344
364
op2_dsp = op2: dsp
345
365
op1_dsp = op1: rangespace (op2_dsp)
346
- return [op1_dsp, op2_dsp], bandwidthssum ((op1_dsp, op2_dsp)),
347
- bandwidthssum ((op1_dsp, op2_dsp), blockbandwidths),
348
- bandwidthssum ((op1_dsp, op2_dsp), subblockbandwidths)
366
+ ops = (op1_dsp, op2_dsp)
367
+ return [ops... ], bandwidthssum (ops),
368
+ bandwidthssum (ops, blockbandwidths),
369
+ bandwidthssum (ops, subblockbandwidths),
370
+ all (isbandedblockbanded, ops),
371
+ all (israggedbelow, ops)
349
372
end
350
373
end
351
374
@@ -361,9 +384,9 @@ bandwidths(P::PlusOrTimesOp) = P.bandwidths
361
384
blockbandwidths (P:: PlusOrTimesOp ) = P. blockbandwidths
362
385
subblockbandwidths (P:: PlusOrTimesOp ) = P. subblockbandwidths
363
386
364
- isbandedblockbanded (P:: PlusOrTimesOp ) = all (isbandedblockbanded, P . ops)
387
+ isbandedblockbanded (P:: PlusOrTimesOp ) = P . isbandedblockbanded
365
388
366
- israggedbelow (P:: PlusOrTimesOp ) = isbandedbelow (P) || all (israggedbelow, P . ops)
389
+ israggedbelow (P:: PlusOrTimesOp ) = P . israggedbelow
367
390
368
391
Base. stride (P:: TimesOperator ) = mapreduce (stride, gcd, P. ops)
369
392
@@ -552,7 +575,8 @@ for OP in (:(adjoint), :(transpose))
552
575
strictconvert (Vector, reverse! (map ($ OP, A. ops))),
553
576
reverse (bandwidths (A)), reverse (size (A)),
554
577
reverse (blockbandwidths (A)),
555
- reverse (subblockbandwidths (A))
578
+ reverse (subblockbandwidths (A)),
579
+ isbandedblockbanded (A),
556
580
)
557
581
end
558
582
0 commit comments