Saddle Operator

[lrsantos11]

Given operators opA, opM and opN, creates saddle point linear operator.

    SaddleOperator(opA, opM, opN; S = promote_type(storage_type.(opA, opM, opN)))

Creates saddle point linear operator:

    [ opM     opAᵀ ]
    [ opA     opN  ]

I'll contribute on this, using BlockDiagonalOperator as template.

LinearOperators.jl/src/special-operators.jl

Lines 229 to 240 in 9025458

	"""
	BlockDiagonalOperator(M1, M2, ..., Mn; S = promote_type(storage_type.(M1, M2, ..., Mn)))
	Creates a block-diagonal linear operator:
	[ M1 ]
	[ M2 ]
	[ ... ]
	[ Mn ]
	Change `S` to use LinearOperators on GPU.
	"""

I already have implemented outside LO. Its a matter of just including, testing and PRing.

[dpo]

Looking forward to it! Thanks! More generally, SaddleOperator could be a special case of a BlockOperator.

[lrsantos11]

@dpo It was nice to meet you at Seattle.
Regarding BlockOperator, do you think it should this be a super type of SaddleOperator, like
SaddleOperator <: BlockOperator?

[dpo]

It was great to (finally) meet you as well @lrsantos11 !

Now that I think of it, hcat and vcat already let us write K = [M A’ ; A N], so perhaps what you need is already there?!

[dpo]

@lrsantos11 The following already works

m, n = 5, 3
A = rand(m, n)
M = rand(m, m); M = M + M'
N = rand(n, n); N = N + N'
opA = LinearOperator(A)
opM = LinearOperator(M)
opN = LinearOperator(N)
K = [opM  opA ;
     opA' opN]

but it doesn’t let us specify that K is symmetric. We could add some logic to make that happen. But I think dedicated types BlockOperator and SaddleOperator are justified because they could have specific methods to retrieve individual blocks.