[petsc-users] schurcomplement of a rectangular system

[Patrick Lacasse]

Hi,

  in my works, I'm computing a formula (involving an L2 projection) of the
form:

(D - M^{-1} N) n

where D,M and B are matrices and n is a vector.

Operator (D - M^{-1} N) appears to be the Schur complement of the system
( M   N )
(  I    D )
and it is convenient for me to use a Mat of type schurcomplement to compute
it.
However, my matrix D is not a square matrix (this system is not invertible,
it is rectangular).

Assertions in src/ksp/ksp/utils/schurm.c prohibit the use of
schurcomplement for rectangular system, but I can't figure where it is
needed by schurcomplement itself.
Could it be possible to remove those two assertions (see my patch)?


Patrick Lacasse

From 65d027de9ccb5e344c206ed0349cd43487ad0c54 Mon Sep 17 00:00:00 2001
From: Patrick Lacasse <patrick.m.laca...@gmail.com>
Date: Tue, 26 Aug 2014 15:59:37 -0400
Subject: [PATCH] Remove an assertion in MatSchurComplement.

Allow to define a matschur for a rectangular system
(A|B)
(-+-)
(C|D)
where D is not square.
---
 src/ksp/ksp/utils/schurm.c | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/ksp/ksp/utils/schurm.c b/src/ksp/ksp/utils/schurm.c
index 80a146c..d015956 100644
--- a/src/ksp/ksp/utils/schurm.c
+++ b/src/ksp/ksp/utils/schurm.c
@@ -262,7 +262,6 @@ PetscErrorCode  MatSchurComplementSetSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,
   if (A11) {
     PetscValidHeaderSpecific(A11,MAT_CLASSID,5);
     PetscCheckSameComm(A00,1,A11,5);
-    if (A11->rmap->n != A11->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A11 %D do not equal local columns %D",A11->rmap->n,A11->cmap->n);
     if (A10->rmap->n != A11->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A10 %D do not equal local rows A11 %D",A10->rmap->n,A11->rmap->n);
   }
 
@@ -401,7 +400,6 @@ PetscErrorCode  MatSchurComplementUpdateSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A
   if (A11) {
     PetscValidHeaderSpecific(A11,MAT_CLASSID,5);
     PetscCheckSameComm(A00,1,A11,5);
-    if (A11->rmap->n != A11->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A11 %D do not equal local columns %D",A11->rmap->n,A11->cmap->n);
     if (A10->rmap->n != A11->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A10 %D do not equal local rows A11 %D",A10->rmap->n,A11->rmap->n);
   }
 
-- 
1.8.4.5