diff --git a/fem/src/BlockSolve.F90 b/fem/src/BlockSolve.F90
index 50501fa096..59c9959653 100644
--- a/fem/src/BlockSolve.F90
+++ b/fem/src/BlockSolve.F90
@@ -31,7 +31,8 @@
 MODULE BlockSolve
 
  USE DefUtils
-
+ USE ListMatrix
+ 
  IMPLICIT NONE
 
   TYPE(BlockMatrix_t), POINTER, SAVE :: TotMatrix
@@ -44,7 +45,64 @@ MODULE BlockSolve
 
 CONTAINS
 
-#define USEPERM 0
+
+  ! Create matrix S=P((diag(A))^-1)Q
+  !------------------------------------------------------------------------  
+  FUNCTION CreateSchurApproximation(A,P,Q) RESULT ( S ) 
+
+    TYPE(Matrix_t), POINTER :: A, P, Q
+    TYPE(Matrix_t), POINTER :: S
+
+    INTEGER :: n, i, j, k, j2, k2
+    REAL(KIND=dp) :: val
+    LOGICAL :: Found
+    
+    CALL Info('CreateSchurApproximation','Creating Shcur complement for preconditioning!',Level=20)
+
+    NULLIFY(S)
+    S => AllocateMatrix()
+    S % FORMAT = MATRIX_LIST
+
+    IF(.NOT. ASSOCIATED(P)) THEN
+      CALL Info('CreateSchurApproximation','Constraint matrix not associated!')
+    END IF
+    n = P % NumberOfRows
+    IF(n == 0) THEN
+      CALL Info('CreateSchurApproximation','No rows in Constraint matrix!')
+      RETURN
+    END IF
+          
+    ! Add the corner entry to give the max size for list.  
+    CALL List_AddToMatrixElement(S % ListMatrix, n, n, 0.0_dp ) 
+
+    DO i=1,n
+      DO j=P % Rows(i),P % Rows(i+1)-1
+        k = P % Cols(j)
+        val = P % Values(k) / A % Values(A % Diag(k))        
+        DO j2= Q % Rows(k),Q % Rows(k+1)-1
+          k2 = Q % Cols(j2)
+          CALL List_AddToMatrixElement(S % ListMatrix, i, k2, -val * Q % Values(k2) )
+        END DO
+      END DO
+    END DO
+
+    CALL List_toCRSMatrix(S)
+    
+    val = 1.0_dp * SIZE(S % Values) / SIZE(P % Values)
+    WRITE(Message,'(A,ES12.3)') 'Schur matrix increase factor: ',val
+    CALL Info('CreateSchurApproximation',Message)
+
+    !ALLOCATE( S % rhs(S % NumberOfRows))
+    !S % rhs = 0.0_dp
+ 
+    IF( InfoActive(20) ) THEN
+      CALL VectorValuesRange(P % Values,SIZE(P % Values),'Constraint')
+      CALL VectorValuesRange(S % Values,SIZE(S % Values),'Schur')
+    END IF          
+    
+  END FUNCTION CreateSchurApproximation
+
+  
   
   !-----------------------------------------------------------------------------------
   !> If a block variable does not exist it will be created. 
@@ -230,15 +288,6 @@ SUBROUTINE BlockInitMatrix( Solver, BlockMatrix, BlockDofs, FieldDofs, SkipVar )
     CHARACTER(*), PARAMETER :: Caller = 'BlockInitMatrix'
 
             
-    Params => Solver % Values
-
-    IsComplex = ListGetLogical( Params,'Linear System Complex',Found)
-    IF( IsComplex ) THEN
-      CALL Info(Caller,'Assuming block matrix to be complex!',Level=8)
-    ELSE
-      CALL Info(Caller,'Assuming block matrix to be real!',Level=20)
-    END IF
-    
     BlockMatrix => Solver % BlockMatrix
     IF (ASSOCIATED(BlockMatrix)) THEN
       CALL Info(Caller,'Using existing block matrix',Level=10)
@@ -246,7 +295,15 @@ SUBROUTINE BlockInitMatrix( Solver, BlockMatrix, BlockDofs, FieldDofs, SkipVar )
     END IF
 
     CALL Info(Caller,'Initializing block matrix',Level=10)
-    
+
+    Params => Solver % Values
+    IsComplex = ListGetLogical( Params,'Linear System Complex',Found)
+    IF( IsComplex ) THEN
+      CALL Info(Caller,'Assuming block matrix to be complex!',Level=8)
+    ELSE
+      CALL Info(Caller,'Assuming block matrix to be real!',Level=20)
+    END IF
+        
     ALLOCATE(Solver % BlockMatrix)
     BlockMatrix => Solver % BlockMatrix
  
@@ -278,11 +335,12 @@ SUBROUTINE BlockInitMatrix( Solver, BlockMatrix, BlockDofs, FieldDofs, SkipVar )
     ELSE
       CALL Info(Caller,'Inheriting blocks from variable dofs',Level=8)
       NoVar = BlockDofs
+      CALL Info(Caller,'Inheriting block count '//I2S(NoVar)//' from variable dofs',Level=8)
     END IF    
 
 
     IF( BlockMatrix % NoVar == NoVar ) THEN
-      CALL Info(Caller,'Reusing existing blockmatrix',Level=6)
+      CALL Info(Caller,'Reusing existing blockmatrix structures!',Level=6)
       RETURN
     ELSE IF( BlockMatrix % Novar /= 0 ) THEN
       CALL Fatal(Caller,'Previous blockmatrix was of different size?')
@@ -295,6 +353,7 @@ SUBROUTINE BlockInitMatrix( Solver, BlockMatrix, BlockDofs, FieldDofs, SkipVar )
 
     
     ALLOCATE( BlockMatrix % SubMatrix(NoVar,NoVar) )
+    CALL Info(Caller,'Allocating block matrix of size '//I2S(NoVar),Level=10)
     DO i=1,NoVar
       DO j=1,NoVar
         DO k=1,2
@@ -318,9 +377,6 @@ SUBROUTINE BlockInitMatrix( Solver, BlockMatrix, BlockDofs, FieldDofs, SkipVar )
     ALLOCATE( BlockMatrix % SubMatrixActive(NoVar,NoVar) )
     BlockMatrix % SubMatrixActive = .FALSE.
 
-    ALLOCATE( BlockMatrix % SubMatrixTranspose(NoVar,NoVar) )
-    BlockMatrix % SubMatrixTranspose = .FALSE.
-
     ALLOCATE( BlockMatrix % SubVector(NoVar))
     DO i=1,NoVar
       BlockMatrix % Subvector(i) % Var => NULL()
@@ -333,7 +389,7 @@ SUBROUTINE BlockInitMatrix( Solver, BlockMatrix, BlockDofs, FieldDofs, SkipVar )
     
     IF( PRESENT( SkipVar ) ) THEN
       IF( SkipVar ) THEN
-        CALL Info(Caller,'Skipping creation of block variables',Level=10)
+        CALL Info(Caller,'Skipping creation of block variables for now',Level=10)
         BlockMatrix % ParentMatrix => Solver % Matrix
         RETURN
       END IF
@@ -362,6 +418,8 @@ SUBROUTINE BlockInitMatrix( Solver, BlockMatrix, BlockDofs, FieldDofs, SkipVar )
     END IF        
     IF( UseSolverMatrix ) THEN
       BlockMatrix % SubVector(1) % Var => Solver % Variable
+!      BlockMatrix % SubMatrix(1,1) % Mat => Solver % Matrix        
+      
       n = SIZE( Solver % Variable % Values )
       BlockMatrix % Offset(1) = 0
       BlockMatrix % Offset(2) = n      
@@ -414,9 +472,9 @@ SUBROUTINE BlockInitMatrix( Solver, BlockMatrix, BlockDofs, FieldDofs, SkipVar )
       IF(ASSOCIATED( Var ) ) THEN
         CALL Info(Caller,'Using existing variable > '//VarName//' <')		
       ELSE		
-        CALL Info(Caller,'Variable > '//VarName//' < does not exist, creating')
         PSolver => Solver
         IF( BlockMatrix % GotBlockStruct ) THEN
+          CALL Info(Caller,'Variable > '//VarName//' < does not exist, creating from existing Perm')
           j = COUNT( BlockMatrix % BlockStruct == i ) 
           IF( j == 0 ) THEN
             CALL Fatal(Caller,'Invalid > Block Structure < given!')
@@ -425,6 +483,7 @@ SUBROUTINE BlockInitMatrix( Solver, BlockMatrix, BlockDofs, FieldDofs, SkipVar )
         ELSE
           Var => CreateBlockVariable(PSolver, i, VarName )
         END IF
+        CALL Info(Caller,'Variable > '//VarName//' < does not exist, creating')
       END IF
       
       BlockMatrix % SubVector(i) % Var => Var
@@ -462,10 +521,12 @@ SUBROUTINE BlockInitVar( Solver, BlockMatrix, BlockIndex )
     TYPE(Mesh_t), POINTER :: Mesh
     REAL(KIND=dp), POINTER :: Vals(:)
     INTEGER, POINTER :: VarPerm(:)
+    TYPE(Matrix_t), POINTER :: B  
     
     Params => Solver % Values
     Mesh => Solver % Mesh
     NoVar = BlockMatrix % NoVar
+    BlockMatrix % Offset = 0
     
     DO i=1,NoVar
       Amat => BlockMatrix % Submatrix(i,i) % Mat 
@@ -485,34 +546,49 @@ SUBROUTINE BlockInitVar( Solver, BlockMatrix, BlockIndex )
         Vals = 0.0_dp
         
         IF( PRESENT(BlockIndex) ) THEN
+          CALL Info('BlockInitVar','Using BlockIndex to pick variable and perm',Level=20)
           NULLIFY(VarPerm)
           m = SIZE(Solver % Variable % Perm) 
           ALLOCATE(VarPerm(m))
           VarPerm = 0
+#if 0
+          ! This works only when dofs follow each other systematically. 
           WHERE( BlockIndex == i )
             VarPerm = Solver % Variable % Perm - BlockMatrix % Offset(i)
           END WHERE          
+#else
+          k = 0
+          DO j=1,SIZE(Solver % Variable % Perm)
+            IF(BlockIndex(j) == i) THEN
+              k = k+1
+              VarPerm(j) = k
+            END IF
+          END DO
+#endif
+          IF( ANY( VarPerm < 0) ) CALL Fatal('BlockInitVar','Negative Perm!')
+          
           CALL VariableAdd( Mesh % Variables,Mesh,PSolver,VarName,1,Vals,&
               Output = .FALSE., Perm = VarPerm )                    
         ELSE
-        CALL VariableAdd( Mesh % Variables,Mesh,PSolver,VarName,1,Vals,&
-            Output = .FALSE. )
+          CALL VariableAdd( Mesh % Variables,Mesh,PSolver,VarName,1,Vals,&
+              Output = .FALSE. )
         END IF
         Var => VariableGet( Mesh % Variables, VarName )
+      ELSE
+        Var % Values = 0.0_dp
       END IF
       BlockMatrix % SubVector(i) % Var => Var
 
-      ! Take the monolithic solution as initial guess
-      !DO j=1,n
-      !  k = Amat % InvPerm(j)
-      !  Var % Values(j) = Solver % Variable % Values(k)
-      !END DO
+      IF(PRESENT(BlockIndex)) THEN
+        B => BlockMatrix % SubMatrix(i,i) % Mat
+        Var % Values = Solver % Variable % Values(B % InvPerm)
+      END IF
 
     END DO
         
     BlockMatrix % TotSize = BlockMatrix % Offset( NoVar + 1 )
 
-    CALL Info('BlockInitVar','All done',Level=12)
+    CALL Info('BlockInitVar','Block variables initialized!',Level=12)
       
   END SUBROUTINE BlockInitVar
 
@@ -587,7 +663,7 @@ SUBROUTINE BlockPickMatrix( Solver, NoVar )
         
     ReuseMatrix = ListGetLogical( Params,'Block Matrix Reuse',Found)
     EliminateZero = ListGetLogical( Params, &
-                         'Block Eliminate Zero Submatrices', Found )
+        'Block Eliminate Zero Submatrices', Found )
 
     DO RowVar=1,NoVar
       DO ColVar=1,NoVar            
@@ -699,7 +775,8 @@ SUBROUTINE BlockPickDofsPhysical( Solver, BlockIndex, NoVar )
     INTEGER, POINTER :: Perm(:)
     
     
-    CALL Info('BlockPickDofsPhysical','Picking block matrix of size '//I2S(NoVar)//' from monolithic one',Level=10)
+    CALL Info('BlockPickDofsPhysical','Picking block matrix of size '&
+        //I2S(NoVar)//' from monolithic one',Level=10)
 
     n_bf = CurrentModel % NumberOfBodyForces
 
@@ -909,7 +986,9 @@ SUBROUTINE BlockPickAV( Solver, BlockIndex, NoVar )
   END SUBROUTINE BlockPickAV
 
   
-  !-------------------------------------------------------------------------------------
+
+
+    !-------------------------------------------------------------------------------------
   !> Picks the components of a full matrix when blockindex table is given.
   !-------------------------------------------------------------------------------------
   SUBROUTINE BlockPickMatrixPerm( Solver, BlockIndex, NoVar )
@@ -993,8 +1072,12 @@ SUBROUTINE BlockPickMatrixPerm( Solver, BlockIndex, NoVar )
               
   END SUBROUTINE BlockPickMatrixPerm
 
-    
+
+
+
+
   
+#if 0   
   !-------------------------------------------------------------------------------------
   !> Picks the components of a full matrix to the submatrices of a block matrix assuming AV solver.
   !-------------------------------------------------------------------------------------
@@ -1125,22 +1208,24 @@ SUBROUTINE BlockPickMatrixAV( Solver, NoVar )
 
     
   END SUBROUTINE BlockPickMatrixAV
-
+#endif
+  
 
   !-------------------------------------------------------------------------------------
   !> Picks vertical and horizontal components of a full matrix.
   !-------------------------------------------------------------------------------------
-  SUBROUTINE BlockPickMatrixHorVer( Solver, NoVar, Cart )
+  SUBROUTINE BlockPickMatrixHorVer( Solver, NoVar, Cart, DTag)
 
     TYPE(Solver_t) :: Solver
     INTEGER :: Novar
     LOGICAL :: Cart
+    INTEGER :: DTag(:)
 
     INTEGER::i,j,k,n,t,ne,dofs,nd,ni,nn,ndir(3),ic,kc
     TYPE(Matrix_t), POINTER :: A,B
     TYPE(Nodes_t), SAVE :: Nodes, EdgeNodes
     INTEGER :: ActiveCoordinate
-    INTEGER, ALLOCATABLE :: DTag(:), DPerm(:)
+    INTEGER, ALLOCATABLE :: DPerm(:)
     INTEGER, POINTER :: Indexes(:)
     TYPE(Mesh_t), POINTER :: Mesh
     REAL(KIND=dp) :: Wlen, Wproj, Wtol, u, v, w, DetJ, Normal(3), MaxCoord, MinCoord
@@ -1158,18 +1243,15 @@ SUBROUTINE BlockPickMatrixHorVer( Solver, NoVar, Cart )
     
     CALL Info('BlockPickMatrixHorVer','Dividing matrix into vertical and horizontal dofs',Level=10)
 
-
+    NoVar = 2
     n = MAXVAL(Solver % Variable % Perm)
     Mesh => Solver % Mesh 
     
     A => Solver % Matrix
     dofs = Solver % Variable % Dofs
     
-    n = A % NumberOfRows / dofs
-    
-    ALLOCATE( DTag(n), DPerm(n*dofs)  ) 
+    n = A % NumberOfRows / dofs    
     DTag = 0
-    DPerm = 0
         
     PiolaVersion = ListGetLogical( Solver % Values,'Use Piola Transform', Found )
     ActiveCoordinate = ListGetInteger( Solver % Values,'Active Coordinate',Found )
@@ -1178,7 +1260,7 @@ SUBROUTINE BlockPickMatrixHorVer( Solver, NoVar, Cart )
     Normal(ActiveCoordinate) = 1.0_dp
     
     Wtol = 1.0e-3
-      
+
     
     DO t=1,Solver % NumberOfActiveElements
       Element => Mesh % Elements( Solver % ActiveElements(t) )
@@ -1275,9 +1357,19 @@ SUBROUTINE BlockPickMatrixHorVer( Solver, NoVar, Cart )
           Wproj = ( MaxCoord - MinCoord ) / Wlen
 
           IF( WProj > 1.0_dp - Wtol ) THEN  
-            DTag(j) = 1  
+            IF( dofs == 1 ) THEN
+              DTag(j) = 1  
+            ELSE
+              DTag(2*j-1) = 1
+              DTag(2*j) = 1
+            END IF
           ELSE IF( Wproj < Wtol ) THEN
-            DTag(j) = 2  
+            IF( dofs == 1 ) THEN
+              DTag(j) = 2  
+            ELSE
+              DTag(2*j-1) = 2
+              DTag(2*j) = 2
+            END IF
           END IF
         END DO
 
@@ -1307,11 +1399,21 @@ SUBROUTINE BlockPickMatrixHorVer( Solver, NoVar, Cart )
           Wproj = ABS( SUM( WBasis(i,:) * Normal ) ) / Wlen 
 
           IF( WProj > 1.0_dp - Wtol ) THEN  
-            IF( DTag(j) == 2 ) PRINT *,'Vertical edge '//I2S(j)//' is also horizontal?'
-            DTag(j) = 1  ! set to be vertical
+            !IF( DTag(j) == 2 ) PRINT *,'Vertical edge '//I2S(j)//' is also horizontal?'
+            IF( dofs == 1 ) THEN
+              DTag(j) = 1  ! set to be vertical
+            ELSE
+              DTag(2*j-1) = 1  ! set to be vertical
+              DTag(2*j) = 1  ! set to be vertical
+            END IF
           ELSE IF( Wproj < Wtol ) THEN
-            IF( DTag(j) == 1 ) PRINT *,'Horizontal edge '//I2S(j)//' is also vertical?'
-            DTag(j) = 2  ! set to be horizontal
+            !IF( DTag(j) == 1 ) PRINT *,'Horizontal edge '//I2S(j)//' is also vertical?'
+            IF( dofs == 1 ) THEN
+              DTag(j) = 2  ! set to be horizontal
+            ELSE
+              DTag(2*j-1) = 2  ! set to be horizontal
+              DTag(2*j) = 2  ! set to be horizontal
+            END IF
           ELSE
             PRINT *,'Edge '//I2S(j)//' direction undefined: ',Wproj
           END IF
@@ -1319,82 +1421,12 @@ SUBROUTINE BlockPickMatrixHorVer( Solver, NoVar, Cart )
 
       END IF
     END DO
-    
-
-    ! Number vertical and horizontal (or all cartesian) dofs separately.
-    ndir = 0
-    DO i=1,n
-      DO j=1,dofs
-        k = dofs*(i-1)+j
-        ndir(DTag(i)) = ndir(DTag(i)) + 1
-        DPerm(k) = ndir(DTag(i))
-      END DO
-    END DO
-
-    !PRINT *,'Cartesian dofs:',ndir(1:NoVar)
-
-    i = n - SUM( ndir ) 
-    IF( i > 0 ) THEN      
-      CALL Fatal('BlockPickMatrixHorVer','Could not determine all nodes: '&
-          //I2S(i))
-    END IF
-
-
-    ! Allocate vectors if not present
-    DO i=1,NoVar
-      DO j=1,NoVar
-        B => TotMatrix % SubMatrix(i,j) % Mat
-        IF( ASSOCIATED( B % Values ) ) B % Values = 0.0_dp
-      END DO
-      B => TotMatrix % SubMatrix(i,i) % Mat      
-      IF(.NOT. ASSOCIATED( B % InvPerm ) ) ALLOCATE( B % InvPerm(ndir(i)) )
-      IF(.NOT. ASSOCIATED( B % Rhs) ) ALLOCATE(B % Rhs(ndir(i)) )
-      !PRINT *,'a complex', a % complex
-      !B % COMPLEX = A % COMPLEX
-    END DO
-    
-
-    DO i=1,A % NumberOfRows
-      ic = (i-1)/dofs+1
-      
-      DO j=A % Rows(i+1)-1,A % Rows(i),-1
-        k = A % Cols(j)
-        kc = (k-1)/dofs+1
-        
-        IF( DTag(ic) < 1 .OR. DTag(ic) > NoVar ) THEN
-          PRINT *,'i:',i,ic,Dtag(ic)
-        END IF
         
-        IF( DTag(kc) < 1 .OR. DTag(kc) > NoVar ) THEN
-          PRINT *,'k:',k,kc,Dtag(kc)
-        END IF
-        
-        B => TotMatrix % SubMatrix(DTag(ic),DTag(kc)) % Mat
-        
-        IF( Dperm(i) < 1 .OR. DPerm(k) < 1 ) THEN
-          PRINT *,'ik',Dperm(i),Dperm(k)
-          STOP EXIT_ERROR
-        END IF
-        CALL AddToMatrixElement(B,Dperm(i),DPerm(k),A % Values(j))
-      END DO
-
-      B => TotMatrix % SubMatrix(DTag(ic),DTag(ic)) % Mat      
-      B % Rhs(Dperm(i)) = A % Rhs(i)          
-      B % InvPerm(Dperm(i)) = i          
-    END DO
-    
-    DO i=1,NoVar
-      DO j=1,NoVar
-        B => TotMatrix % SubMatrix(i,j) % Mat        
-        IF (B % FORMAT == MATRIX_LIST) THEN
-          CALL List_toCRSMatrix(B)
-        END IF
+    IF( InfoActive(20) ) THEN
+      DO i=1,NoVar
+        j = COUNT(DTag==i)
+        CALL Info('BlockPickMatrixHorVer','There are '//I2S(j)//' dofs group '//I2S(i))
       END DO
-    END DO
-
-    
-    IF( ASSOCIATED( A % ConstraintMatrix ) ) THEN
-      CALL Warn('BlockPickMatrixHorVer','Cannot deal with constraints')
     END IF
     
   END SUBROUTINE BlockPickMatrixHorVer
@@ -1403,18 +1435,19 @@ END SUBROUTINE BlockPickMatrixHorVer
   !-------------------------------------------------------------------------------------
   !> Makes a quadratic H(curl) approximation to have a block structure
   !-------------------------------------------------------------------------------------
-  SUBROUTINE BlockPickMatrixHcurl( Solver, NoVar, DoCmplx )
+  SUBROUTINE BlockPickMatrixHcurl( Solver, NoVar, DoCmplx, BlockTag )
 
     TYPE(Solver_t) :: Solver
     INTEGER :: Novar
     LOGICAL :: DoCmplx
+    INTEGER :: BlockTag(:)
 
     INTEGER :: i,j,k,n,m,n0,dofs,ic,kc
     TYPE(Matrix_t), POINTER :: A,B
-    INTEGER, ALLOCATABLE :: BlockTag(:), BlockPerm(:), nblock(:)
     TYPE(Mesh_t), POINTER :: Mesh
     TYPE(Element_t), POINTER :: Element, Edge
-
+    LOGICAL :: Found
+    
     Mesh => Solver % Mesh
 
     IF(.NOT. ASSOCIATED( Mesh % Edges ) ) THEN
@@ -1422,18 +1455,22 @@ SUBROUTINE BlockPickMatrixHcurl( Solver, NoVar, DoCmplx )
     END IF
     IF(.NOT. ASSOCIATED( Mesh % Faces ) ) THEN
       CALL Fatal('BlockPickMatrixHcurl','This subroutine needs Faces!')
-    END IF
-        
+    END IF        
     CALL Info('BlockPickMatrixHcurl','Arranging a quadratic H(curl) approximation into blocks',Level=10)
 
+    NoVar = 2
+    IF( ListGetLogical( Solver % Values,'Block Quadratic Hcurl Faces',Found ) ) NoVar = 3
+    IF(DoCmplx) THEN
+      NoVar = 2 * NoVar 
+      IF( ListGetLogical( Solver % Values,'Block Quadratic Hcurl semicomplex',Found ) ) NoVar = 3
+    END IF
+
+    
     A => Solver % Matrix
     dofs = Solver % Variable % Dofs
 
     m = A % NumberOfRows    
     n = m
-    IF(.NOT. DoCmplx) n = n / dofs
-      
-    ALLOCATE( BlockTag(n), BlockPerm(m), nblock(NoVar) )
 
     ! Set the default blocks
     IF( DoCmplx ) THEN
@@ -1462,14 +1499,19 @@ SUBROUTINE BlockPickMatrixHcurl( Solver, NoVar, DoCmplx )
       j = n0 + 2*i-1
       k = Solver % Variable % Perm(j)
       IF(k==0) CYCLE
-      IF( DoCmplx ) THEN
-        BlockTag(2*k-1) = 1
-        BlockTag(2*k) = 2 
-      ELSE
-        ! If BlockTag array were created for all DOFs with DOFs>1,
-        ! it would have a repeated entries occuring in clusters of
-        ! size DOFs. Therefore a smaller array can be used to tag DOFs.
+
+      ! If BlockTag array were created for all DOFs with DOFs>1,
+      ! it would have a repeated entries occuring in clusters of
+      ! size DOFs. Therefore a smaller array can be used to tag DOFs.
+      IF( dofs == 1 ) THEN
         BlockTag(k) = 1
+      ELSE IF(dofs == 2 ) THEN
+        BlockTag(2*k-1) = 1
+        IF( DoCmplx ) THEN
+          BlockTag(2*k) = 2 
+        ELSE
+          BlockTag(2*k) = 1
+        END IF
       END IF
     END DO
 
@@ -1481,208 +1523,96 @@ SUBROUTINE BlockPickMatrixHcurl( Solver, NoVar, DoCmplx )
       ELSE
         DO j = n0 + 2*Mesh % NumberOfEdges + 1, SIZE(Solver % Variable % Perm)
           k = Solver % Variable % Perm(j)
-          IF(k>0) BlockTag(k) = 3
+          IF(k==0) CYCLE
+
+          IF(dofs == 1 ) THEN
+            BlockTag(k) = 3
+          ELSE IF( dofs == 2 ) THEN
+            BlockTag(2*k-1) = 3
+            BlockTag(2*k) = 3
+          END IF
         END DO
       END IF
     END IF
-    
-    ! Number each group of DOFs separately
-    BlockPerm = 0
-    nblock = 0
-    IF( DoCmplx ) THEN
-      DO i=1,n
-        nblock(BlockTag(i)) = nblock(BlockTag(i)) + 1
-        BlockPerm(i) = nblock(BlockTag(i))
-      END DO
-    ELSE
-      DO i=1,n
-        n0 = dofs*(i-1)
-        DO j=1,dofs
-          nblock(BlockTag(i)) = nblock(BlockTag(i)) + 1
-          k = n0 + j
-          BlockPerm(k) = nblock(BlockTag(i))
-        END DO
-      END DO
-    END IF
-
-    DO i=1,NoVar
-      CALL Info('BlockPickMatrixHcurl','Block vector '//I2S(i)//' size: '//I2S(nblock(i)),Level=6)
-    END DO
-
-    ! Allocate vectors if not present
-    DO i=1,NoVar
-      DO j=1,NoVar
-        B => TotMatrix % SubMatrix(i,j) % Mat
-        IF( ASSOCIATED( B % Values ) ) B % Values = 0.0_dp
-      END DO
-      B => TotMatrix % SubMatrix(i,i) % Mat      
-      IF(.NOT. ASSOCIATED( B % InvPerm ) ) ALLOCATE( B % InvPerm(nblock(i)) )
-      IF(.NOT. ASSOCIATED( B % Rhs) ) ALLOCATE(B % Rhs(nblock(i)) )
-      !PRINT *,'a complex', a % complex
-      !B % COMPLEX = A % COMPLEX
-    END DO
-    
-
-    DO i=1,A % NumberOfRows
-      IF( DoCmplx ) THEN
-        ic = i
-      ELSE       
-        ic = (i-1)/dofs+1
-      END IF
-        
-      DO j=A % Rows(i+1)-1,A % Rows(i),-1
-        k = A % Cols(j)
-
-        IF( DoCmplx ) THEN
-          kc = k
-        ELSE
-          kc = (k-1)/dofs+1
-        END IF
-          
-        IF( BlockTag(ic) < 1 .OR. BlockTag(ic) > NoVar ) THEN
-          PRINT *,'i:',i,ic,BlockTag(ic)
-        END IF
-        
-        IF( BlockTag(kc) < 1 .OR. BlockTag(kc) > NoVar ) THEN
-          PRINT *,'k:',k,kc,BlockTag(kc)
-        END IF
-        
-        B => TotMatrix % SubMatrix(BlockTag(ic),BlockTag(kc)) % Mat
-        
-        IF( BlockPerm(i) < 1 .OR. BlockPerm(k) < 1 ) THEN
-          PRINT *,'ik',BlockPerm(i),BlockPerm(k)
-          STOP EXIT_ERROR
-        END IF
-        CALL AddToMatrixElement(B,BlockPerm(i),BlockPerm(k),A % Values(j))
-      END DO
 
-      B => TotMatrix % SubMatrix(BlockTag(ic),BlockTag(ic)) % Mat      
-      B % Rhs(BlockPerm(i)) = A % Rhs(i)          
-      B % InvPerm(BlockPerm(i)) = i          
-    END DO
-    
-    DO i=1,NoVar
-      DO j=1,NoVar
-        B => TotMatrix % SubMatrix(i,j) % Mat        
-        IF (B % FORMAT == MATRIX_LIST) THEN
-          CALL List_toCRSMatrix(B)
-        END IF
-        CALL Info('BlockPickMatrixHcurl','Matrix '//I2S(10*i+j)//' nonzeros: '//I2S(SIZE(B % Values)),Level=6)
+    IF( InfoActive(20) ) THEN
+      DO i=1,NoVar
+        j = COUNT(BlockTag==i)
+        CALL Info('BlockMatrixHCurl','There are '//I2S(j)//' dofs group '//I2S(i))
       END DO
-    END DO
-    
-    IF( ASSOCIATED( A % ConstraintMatrix ) ) THEN
-      CALL Warn('BlockPickMatrixHcurl','Cannot deal with constraints')
     END IF
-    
+        
   END SUBROUTINE BlockPickMatrixHcurl
   
 
+
   !-------------------------------------------------------------------------------------
-  !> Picks the components of a full matrix to the submatrices of a block matrix assuming AV solver.
+  !> Picks apart nodal and non-nodal dofs.
   !-------------------------------------------------------------------------------------
-  SUBROUTINE BlockPickMatrixNodal( Solver, NoVar )
+  SUBROUTINE BlockPickMatrixNodal( Solver, NoVar, DoCmplx, BlockTag )
 
     TYPE(Solver_t) :: Solver
     INTEGER :: Novar
+    LOGICAL :: DoCmplx
+    INTEGER :: BlockTag(:)
 
-    INTEGER :: RowVar, ColVar, Dofs
-    TYPE(Matrix_t), POINTER :: SolverMatrix
-    TYPE(Matrix_t), POINTER :: Amat
-    INTEGER::i,j,k,ii,l,ll,i_aa,i_vv,i_av,i_va,n_a,n_v,ntot,rdof,cdof
-    TYPE(Matrix_t), POINTER :: B_aa,B_av,B_va,B_vv,C_aa,C_vv,A,CM
-    REAL(KIND=DP) :: SumAbsMat, val
-    
-    CALL Info('BlockPickMatrixNodal','Picking nodal and non-nodal block matrices from monolithic one',Level=10)
-
-    SolverMatrix => Solver % Matrix 
-    
-    A => SolverMatrix
-    i_aa=0; i_vv=0; i_av=0; i_va=0;
-    i = Solver % Mesh % NumberOfNodes
-    n_v = MAXVAL( Solver % Variable % Perm(1:i) )
-    ntot = MAXVAL( Solver % Variable % Perm ) 
-    n_a = ntot - n_v    
-    
-    dofs = Solver % Variable % Dofs
-
-    !PRINT *,'Dofs, n_v, n_a: ',dofs,n_v,n_a
+    INTEGER :: i,j,k,n,dofs
+    TYPE(Matrix_t), POINTER :: A
+    TYPE(Mesh_t), POINTER :: Mesh
+    LOGICAL :: Found
     
-    B_vv => TotMatrix % SubMatrix(1,1) % Mat
-    B_va => TotMatrix % SubMatrix(1,2) % Mat
-    B_av => TotMatrix % SubMatrix(2,1) % Mat
-    B_aa => TotMatrix % SubMatrix(2,2) % Mat
-
-    IF(ASSOCIATED(B_aa % Values)) B_aa % Values=0._dp
-    IF(ASSOCIATED(B_av % Values)) B_av % Values=0._dp
-    IF(ASSOCIATED(B_va % Values)) B_va % Values=0._dp
-    IF(ASSOCIATED(B_vv % Values)) B_vv % Values=0._dp
+    Mesh => Solver % Mesh
 
-    IF( .NOT. ASSOCIATED( B_aa % Rhs ) ) ALLOCATE(B_aa % Rhs(dofs*n_a))
-    IF( .NOT. ASSOCIATED( B_vv % Rhs ) ) ALLOCATE(B_vv % Rhs(dofs*n_v))
+    CALL Info('BlockPickMatrixNodal','Separates nodal and non-nodal dofs from each other',Level=10)
 
+    NoVar = 2
+    IF(DoCmplx) NoVar = 2 * NoVar 
     
-    DO i=1,A % NumberOfRows
-      
-      ii = (i-1)/dofs+1
-      rdof = (i-1)/ntot+1
+    A => Solver % Matrix
+    dofs = Solver % Variable % Dofs
 
-      IF( ii < n_v ) THEN
-        B_vv % Rhs(i-(rdof-1)*n_a) = A % Rhs(i)
-      ELSE
-        B_aa % Rhs(i-rdof*n_v) = A % Rhs(i)
-      END IF
-      
-      DO k=A % Rows(i+1)-1,A % Rows(i),-1
-        l = A % Cols(k)
+    IF( DoCmplx ) THEN
+      BlockTag(1::2) = 3
+      BlockTag(2::2) = 4
+    ELSE
+      BlockTag = 2
+    END IF
         
-        ll = (l-1)/dofs+1
-        cdof = (l-1)/ntot+1
-        val = A % Values(k)
-                    
-        IF( ii <= n_v ) THEN
-          IF( ll <= n_v ) THEN
-            CALL AddToMatrixElement(B_vv,i-(rdof-1)*n_a,l-(cdof-1)*n_a,val)
-          ELSE
-            IF( ABS( val ) > TINY( val ) ) THEN
-              CALL AddToMatrixElement(B_va,i-(rdof-1)*n_a,l-cdof*n_v,val)
-            END IF
-          END IF
+    DO i=1, Mesh % NumberOfNodes
+      k = Solver % Variable % Perm(j)
+      IF(k==0) CYCLE
+
+      IF( dofs == 1 ) THEN
+        BlockTag(k) = 1
+      ELSE IF(dofs == 2 ) THEN        
+        BlockTag(2*k-1) = 1
+        IF( DoCmplx ) THEN
+          BlockTag(2*k) = 2 
         ELSE
-          IF( ll <= n_v ) THEN
-            IF( ABS( val ) > TINY( val ) ) THEN
-              CALL AddToMatrixElement(B_av,i-rdof*n_v,l-(cdof-1)*n_a,val)
-            END IF
-          ELSE
-            CALL AddToMatrixElement(B_aa,i-rdof*n_v,l-cdof*n_v,val)
-          END IF
+          BlockTag(2*k) = 1
         END IF
-      END DO
+      END IF
     END DO
 
-    IF (B_aa % Format == MATRIX_LIST) THEN
-      CALL List_toCRSMatrix(B_aa)
-      CALL List_toCRSMatrix(B_av)
-      CALL List_toCRSMatrix(B_va)
-      CALL List_toCRSMatrix(B_vv)
-    END IF
-   
-    CM => A % ConstraintMatrix
-    IF( ASSOCIATED( CM ) ) THEN
-      CALL Fatal('BlockPickMatrixNodal','There would be some constraints to pick too!')
+    IF( InfoActive(20) ) THEN
+      DO i=1,NoVar
+        j = COUNT(BlockTag==i)
+        CALL Info('BlockMatrixNodal','There are '//I2S(j)//' dofs group '//I2S(i))
+      END DO
     END IF
-    
+        
   END SUBROUTINE BlockPickMatrixNodal
-
+  
 
 
   !-------------------------------------------------------------------------------------
   !> Picks the components of the constraint matrix.
   !-------------------------------------------------------------------------------------
-  SUBROUTINE BlockPickConstraint( Solver, NoVar )
+  SUBROUTINE BlockPickConstraint( Solver, NoVar, SkipPrec )
     
     TYPE(Solver_t), TARGET :: Solver
     INTEGER :: Novar
+    LOGICAL :: SkipPrec
 
     INTEGER :: RowVar, ColVar
     TYPE(Matrix_t), POINTER :: SolverMatrix
@@ -1729,16 +1659,23 @@ SUBROUTINE BlockPickConstraint( Solver, NoVar )
     CALL Info('BlockPickConstraint','Number of constraint matrices: '//I2S(NoCon),Level=10)
     
     InheritCM = (NoVar == 1 ) .AND. (NoCon == 1 )
-    
+    IF( InheritCM ) THEN
+      CALL info('BlockPickConstraint','Inheriting constraint matrix',Level=20)
+    END IF
+      
     
     IF( NoVar == 1 ) THEN
       IF( InheritCM ) THEN
+        C1 => A % ConstraintMatrix 
         TotMatrix % Submatrix(NoVar+1,1) % Mat => A % ConstraintMatrix
-        CALL Info('BlockPickConstraint','Using constraint matrix as is!',Level=10)
+        CALL Info('BlockPickConstraint',&
+            'Using constraint matrix ('//I2S(NoVar+1)//',1) as is!',Level=10)
+        i =  A % ConstraintMatrix % NumberOfRows
+        CALL Info('BlockPickConstraint','Number of rows in matrix: '//I2S(i),Level=20)
         IF( PrecTrue ) THEN
           C1prec => TotMatrix % Submatrix(NoVar+1,NoVar+1) % Mat                  
         ELSE
-          C1prec => TotMatrix % Submatrix(NoVar+1,NoVar+1) % PrecMat        
+          C1prec => TotMatrix % Submatrix(NoVar+1,NoVar+1) % PrecMat
         END IF
       ELSE      
         C1 => TotMatrix % Submatrix(NoVar+1,1) % Mat
@@ -1748,8 +1685,7 @@ SUBROUTINE BlockPickConstraint( Solver, NoVar )
           C1prec => TotMatrix % Submatrix(NoVar+1,NoVar+1) % PrecMat        
         END IF
       END IF
-      TotMatrix % SubMatrixTranspose(NoVar+1,1) = .TRUE.
-      
+
     ELSE IF(BlockAV) THEN          
       C1 => TotMatrix % SubMatrix(NoVar+1,1) % Mat
       C2 => TotMatrix % Submatrix(NoVar+2,2) % Mat
@@ -1760,8 +1696,6 @@ SUBROUTINE BlockPickConstraint( Solver, NoVar )
         C1prec => TotMatrix % Submatrix(NoVar+1,NoVar+1) % PrecMat
         C2prec => TotMatrix % Submatrix(NoVar+2,NoVar+2) % PrecMat
       END IF
-      TotMatrix % SubMatrixTranspose(NoVar+1,1) = .TRUE.
-      TotMatrix % SubMatrixTranspose(NoVar+2,2) = .TRUE.
     ELSE
       CALL Fatal('BlockPickConstraint','Not done for vectors!')
     END IF
@@ -1779,6 +1713,7 @@ SUBROUTINE BlockPickConstraint( Solver, NoVar )
     
 
     DO DoPrec = 0, 1
+      IF( DoPrec == 1 .AND. SkipPrec ) CYCLE
       
       i1 = 0
       i2 = 0
@@ -1845,47 +1780,66 @@ SUBROUTINE BlockPickConstraint( Solver, NoVar )
       END DO
 
       ! It is more efficient to set the last entry of the list matrix first      
-      IF( DoPrec == 0 ) THEN
+#if 0 
+      IF( DoPrec == 0 .AND. .NOT. SkipPrec ) THEN
         CALL AddToMatrixElement(C1prec,i1,i1,0.0_dp)
       END IF
-               
+#endif
+      
     END DO
       
     CALL Info('BlockPickConstraint','Setting format of constraint blocks to CRS',Level=20)
     IF(.NOT. InheritCM ) THEN
       CALL List_toCRSMatrix(C1)
     END IF
-    CALL List_toCRSMatrix(C1prec)
-
+    IF(.NOT. SkipPrec ) CALL List_toCRSMatrix(C1prec)
+      
     IF( BlockAV ) THEN
       CALL List_toCRSMatrix(C2)    
-      CALL List_toCRSMatrix(C2prec)
+      IF(.NOT. SkipPrec) CALL List_toCRSMatrix(C2prec)
     END IF
 
     
     IF( ListGetLogical( Solver % Values,'Save Prec Matrix', Found ) ) THEN   
       CALL SaveProjector(C1prec,.TRUE.,"CM")
     END IF
-    
-    
+        
     VarName = "lambda"
     Var => VariableGet( Solver % Mesh % Variables, VarName )
     IF(ASSOCIATED( Var ) ) THEN
-      CALL Info('BlockPickConstraint','Using existing variable > '//VarName//' <')		
+      n = SIZE( Var % Values ) 
+      CALL Info('BlockPickConstraint','Using existing variable > '//VarName//' <')		      
     ELSE		
       CALL Info('BlockPickConstraint','Variable > '//VarName//' < does not exist, creating')
-      PSolver => Solver
-      
+      PSolver => Solver      
       n = i1
       Var => CreateBlockVariable(PSolver, NoVar+1, VarName, 1, ConsPerm )      
-
-      TotMatrix % SubVector(NoVar+1) % Var => Var      
-      TotMatrix % Offset(NoVar+2) = TotMatrix % Offset(NoVar+1) + n
-      TotMatrix % MaxSize = MAX( TotMatrix % MaxSize, n )
-      TotMatrix % TotSize = TotMatrix % TotSize + n
     END IF
     
+    TotMatrix % SubVector(NoVar+1) % Var => Var      
+    TotMatrix % Offset(NoVar+2) = TotMatrix % Offset(NoVar+1) + n
+    TotMatrix % MaxSize = MAX( TotMatrix % MaxSize, n )
+    TotMatrix % TotSize = TotMatrix % TotSize + n
 
+    IF(ASSOCIATED(TotMatrix % SubMatrix(1,NoVar+1) % Mat)) & 
+        CALL FreeMatrix(TotMatrix % SubMatrix(1,NoVar+1) % Mat)
+    IF( BlockAV ) THEN
+      IF(ASSOCIATED(TotMatrix % SubMatrix(1,NoVar+2) % Mat)) & 
+          CALL FreeMatrix(TotMatrix % SubMatrix(1,NoVar+2) % Mat)
+    END IF
+        
+    TotMatrix % SubMatrix(1,NoVar+1) % Mat => CRS_Transpose( C1 )
+    IF( BlockAV ) THEN
+      TotMatrix % SubMatrix(1,NoVar+2) % Mat => CRS_Transpose( C2 )
+      TotMatrix % SubMatrixActive(1,NoVar+2) = .TRUE.
+    END IF
+    TotMatrix % SubMatrix(1,NoVar+1) % ParallelIsolatedMatrix = &
+        TotMatrix % SubMatrix(NoVar+1,1) % ParallelIsolatedMatrix
+    IF( BlockAV ) THEN
+      TotMatrix % SubMatrix(1,NoVar+2) % ParallelIsolatedMatrix = &
+          TotMatrix % SubMatrix(NoVar+2,1) % ParallelIsolatedMatrix 
+    END IF
+    
   END SUBROUTINE BlockPickConstraint
 
 
@@ -1910,9 +1864,9 @@ SUBROUTINE BlockPrecMatrix( Solver, NoVar )
     TYPE(Variable_t), POINTER :: AVar
     
     CALL Info('BlockPrecMatrix','Checking for tailored preconditioning matrices',Level=6)
-    
+
     Params => Solver % Values
-    
+
     ! The user may give a user defined preconditioner matrix
     !-----------------------------------------------------------
     DO RowVar=1,NoVar
@@ -1995,14 +1949,26 @@ SUBROUTINE BlockPrecMatrix( Solver, NoVar )
         Amat => TotMatrix % Submatrix(RowVar,RowVar) % PrecMat 
         CALL Info('BlockPrecMatrix','Creating preconditioner from matrix ('&
             //I2S(RowVar)//','//I2S(CopyVar)//')',Level=6)
-        PRINT *,'proj matrix sum:',SUM( ABS( PMat % Values ) ) 
-        PRINT *,'orig diag matrix sum:',SUM( ABS( TotMatrix % Submatrix(RowVar,RowVar) % Mat % Values ) )
         
         CALL DiagonalMatrixSumming( Solver, PMat, Amat )
         Amat % Values = Coeff * Amat % Values
       END IF
     END DO
     
+    IF( ListGetLogical( Params,'Create Schur Matrix Approximation',GotIt ) ) THEN
+      CALL Info('BlockPrecMatrix','Generating block '//I2S(11*NoVar),Level=7)
+      IF(NoVar == 1) THEN
+        CALL Fatal('BlockPrecMatrix','We should have more than one block')
+      END IF
+      Pmat => CreateSchurApproximation( &
+          TotMatrix % Submatrix(1,1) % Mat, &
+          TotMatrix % Submatrix(NoVar,1) % Mat, &
+          TotMatrix % Submatrix(1,NoVar) % Mat )
+      TotMatrix % Submatrix(NoVar,NoVar) % PrecMat => Pmat
+      NULLIFY(Pmat)
+    END IF
+
+    
     str = ListGetString( Params,'Block Matrix Schur Variable', GotIt)      
     IF( GotIt ) THEN
       AVAr => VariableGet( Solver % Mesh % Variables, str )
@@ -2017,8 +1983,9 @@ SUBROUTINE BlockPrecMatrix( Solver, NoVar )
       END IF
       CALL Info('BlockPrecMatrix','Using Schur matrix to precondition block '//I2S(NoVar))
       TotMatrix % Submatrix(NoVar,NoVar) % PrecMat => AVar % Solver % Matrix
-    END IF  
-
+    END IF
+    
+    
     ! When we have an inner-outer iteration, we could well have a different matrix
     ! assembled for the purpose of preconditioning. Use it here, if available.
     IF(ListGetLogical( Params,'Block Nested System',GotIt ) ) THEN
@@ -2378,7 +2345,6 @@ SUBROUTINE BlockUpdateRhs( BlockMatrix, ThisRow )
     REAL(KIND=dp), POINTER :: x(:),rtmp(:),rhs(:)
     REAL(KIND=dp) :: bnorm
     TYPE(Variable_t), POINTER :: Var
-    LOGICAL :: GotRhs, Trans
     
     CALL Info('BlockUpdateRhs','Computing block r.h.s',Level=8)
 
@@ -2403,21 +2369,15 @@ SUBROUTINE BlockUpdateRhs( BlockMatrix, ThisRow )
       !-----------------------------------------------------------
       IF(.NOT. ALLOCATED( BlockMatrix % SubVector(NoRow) % rhs )) THEN
         ALLOCATE( BlockMatrix % SubVector(NoRow) % rhs(n) )
-        BlockMatrix % SubVector(NoRow) % rhs = 0.0_dp
         CALL Info('BlockUpdateRhs','Creating rhs for component: '//I2S(NoRow),Level=12)
       END IF
       rhs => BlockMatrix % SubVector(NoRow) % rhs
+      rhs = 0.0_dp
       
       A => BlockMatrix % SubMatrix( NoRow, NoRow ) % Mat
-      GotRhs = .FALSE.
       IF( ASSOCIATED( A ) ) THEN
-        IF( ASSOCIATED( A % rhs ) ) THEN
-          GotRhs = .TRUE.
-          rhs = A % rhs
-        END IF
+        IF( ASSOCIATED( A % rhs ) ) rhs = A % rhs
       END IF
-
-      IF( .NOT. GotRhs ) rhs = 0.0_dp
       
       DO NoCol = 1,NoVar           
         ! This ensures that the diagonal itself is not subtracted
@@ -2427,22 +2387,10 @@ SUBROUTINE BlockUpdateRhs( BlockMatrix, ThisRow )
         Var => BlockMatrix % SubVector(NoCol) % Var
         x => Var % Values
 
-        Trans = .FALSE.        
         A => BlockMatrix % SubMatrix( NoRow, NoCol ) % Mat
-        IF( A % NumberOfRows == 0 ) THEN
-          IF( BlockMatrix % SubMatrixTranspose(NoCol,NoRow) ) THEN         
-            A => TotMatrix % SubMatrix(NoCol,NoRow) % Mat
-            Trans = .TRUE.
-          END IF
-        END IF
-                       
         IF( A % NumberOfRows == 0 ) CYCLE
 
-        IF(.NOT. Trans) THEN
-          CALL CRS_MatrixVectorMultiply( A, x, rtmp)              
-        ELSE 
-          CALL CRS_TransposeMatrixVectorMultiply( A, x, rtmp )
-        END IF
+        CALL CRS_MatrixVectorMultiply( A, x, rtmp)              
         rhs(1:n) = rhs(1:n) - rtmp(1:n) 
       END DO
 
@@ -2486,7 +2434,7 @@ SUBROUTINE BlockMatrixVectorProd( u,v,ipar )
     TYPE(Matrix_t), POINTER :: A
     REAL(KIND=dp), POINTER :: b(:)
 
-    LOGICAL :: Trans, Isolated
+    LOGICAL :: Isolated
     LOGICAL :: DoSum , DoAMGXMV, Found
 
     DoAMGXMV = ListGetLogical( SolverRef % Values, 'Block AMGX M-V', Found)
@@ -2513,29 +2461,10 @@ SUBROUTINE BlockMatrixVectorProd( u,v,ipar )
         j1 = offset(j)+1
         j2 = offset(j+1)
 
-        Trans = .FALSE.
-        Isolated = .FALSE.
         A => TotMatrix % SubMatrix(i,j) % Mat
-        
-        IF( A % NumberOfRows == 0 ) THEN
-          IF( TotMatrix % SubMatrixTranspose(j,i) ) THEN
-            A => TotMatrix % SubMatrix(j,i) % Mat
-            IF( A % NumberOfRows >  0 ) THEN
-              Trans = .TRUE.
-              CALL Info('BlockMatrixVectorProd','Multiplying with transpose of submatrix ('&
-                  //I2S(j)//','//I2S(i)//')',Level=10)
-              IF( isParallel ) THEN
-                Isolated = TotMatrix % SubMatrix(j,i) % ParallelIsolatedMatrix
-                IF( .NOT. Isolated ) THEN
-                  CALL Fatal('BlockMatrixVectorProd','Only isolated matrices may be transposed in parallel!')
-                END IF
-              END IF
-            END IF
-          END IF
-        ELSE
-          Isolated = TotMatrix % SubMatrix(i,j) % ParallelIsolatedMatrix 
-        END IF
+                
         IF( A % NumberOfRows == 0) CYCLE
+        Isolated = TotMatrix % SubMatrix(i,j) % ParallelIsolatedMatrix 
         
         CALL Info('BlockMatrixVectorProd','Multiplying with submatrix ('&
             //I2S(i)//','//I2S(j)//')',Level=15)          
@@ -2544,23 +2473,15 @@ SUBROUTINE BlockMatrixVectorProd( u,v,ipar )
           IF( ASSOCIATED( A % ParMatrix ) ) THEN
             CALL ParallelMatrixVector( A, u(j1:j2), s  )
           ELSE IF( Isolated ) THEN
-            IF(.NOT. Trans ) THEN
-              CALL CRS_MatrixVectorMultiply( A, u(j1:j2), s )
-            ELSE
-              CALL CRS_TransposeMatrixVectorMultiply( A, u(j1:j2), s )              
-            END IF
+            CALL CRS_MatrixVectorMultiply( A, u(j1:j2), s )
           ELSE
             CALL Fatal('BlockMatrixVectorProd','Cannot make the matric-vector product in parallel!')
           END IF
         ELSE
-          IF( .NOT. Trans ) THEN
-            IF ( DoAMGXMV ) THEN
-              CALL AMGXMatrixVectorMultiply(A, u(j1:j2), s, SolverRef )
-            ELSE
-              CALL CRS_MatrixVectorMultiply( A, u(j1:j2), s )
-            END IF
+          IF ( DoAMGXMV ) THEN
+            CALL AMGXMatrixVectorMultiply(A, u(j1:j2), s, SolverRef )
           ELSE
-            CALL CRS_TransposeMatrixVectorMultiply( A, u(j1:j2), s )
+            CALL CRS_MatrixVectorMultiply( A, u(j1:j2), s )
           END IF
         END IF
           
@@ -2607,11 +2528,13 @@ END SUBROUTINE BlockMatrixVectorProd
   SUBROUTINE ParallelShrinkPerm()
     TYPE(Matrix_t), POINTER :: A, Adiag    
     INTEGER, POINTER :: BlockPerm(:), ParBlockPerm(:), ParPerm(:)    
-    INTEGER :: i,j,k,l,n
+    INTEGER :: i,j,k,l,n,m
     INTEGER, ALLOCATABLE :: ShrinkPerm(:),RenumPerm(:)
-
+    LOGICAL :: Halt
+    
     n = TotMatrix % TotSize
-           
+    Halt = .FALSE.
+    
     ! Example of blockperm: block solution u to monolithic solution v
     ! v(1:n) = u(BlockPerm(1:n)) 
 
@@ -2632,14 +2555,29 @@ SUBROUTINE ParallelShrinkPerm()
       IF(.NOT. ASSOCIATED( Adiag % ParallelInfo ) ) CYCLE
       k = 0
       l = 0
+      m = 0
       ShrinkPerm(1:Adiag % NumberOfRows) = 0 
       DO i=1,Adiag % NumberOfRows
         k = k + 1
-        IF (Parenv % MyPE /= Adiag % ParallelInfo % NeighbourList(i) % Neighbours(1)) CYCLE
+        IF (Parenv % MyPE /= Adiag % ParallelInfo % NeighbourList(i) % Neighbours(1)) THEN
+          m = m+1
+          CYCLE          
+        END IF
         l = l+1
         ShrinkPerm(k) = l
       END DO
 
+      IF(ASSOCIATED(TotMatrix % Submatrix(j,j) % ParPerm )) THEN
+        k = SIZE(TotMatrix % SubMatrix(j,j) % ParPerm)
+        IF(k /= l) THEN
+          CALL Warn('ParallelShrinkPerm','Different size for ParPerm '&
+              //I2S(l)//': '//I2S(k)//' vs. '//I2S(l))
+          PRINT *,'ParBlockPerm skipped1',ParEnv % MyPe, m, j
+          DEALLOCATE(TotMatrix % SubMatrix(j,j) % ParPerm )
+          Halt = .TRUE.
+        END IF
+      END IF
+
       IF(.NOT. ASSOCIATED(TotMatrix % Submatrix(j,j) % ParPerm ) ) THEN
         ALLOCATE( TotMatrix % Submatrix(j,j) % ParPerm(l) )
       END IF
@@ -2652,12 +2590,16 @@ SUBROUTINE ParallelShrinkPerm()
       END DO      
     END DO
 
+    m = 0 
     
     ShrinkPerm = 0    
     IF( ASSOCIATED( A ) ) THEN
       l = 0
       DO i=1,A % NumberOfRows
-        IF (Parenv % MyPE /= A % ParallelInfo % NeighbourList(i) % Neighbours(1)) CYCLE
+        IF (Parenv % MyPE /= A % ParallelInfo % NeighbourList(i) % Neighbours(1)) THEN
+          m = m+1
+          CYCLE
+        END IF
         l = l+1
         ShrinkPerm(i) = l
       END DO
@@ -2682,10 +2624,14 @@ SUBROUTINE ParallelShrinkPerm()
     k = 0
     l = 0      
     DO j=1,TotMatrix % Novar
+      m = 0
       A => TotMatrix % Submatrix(j,j) % Mat
       DO i=1,A % NumberOfRows
         k = k + 1
-        IF (Parenv % MyPE /= A % ParallelInfo % NeighbourList(i) % Neighbours(1)) CYCLE
+        IF (Parenv % MyPE /= A % ParallelInfo % NeighbourList(i) % Neighbours(1)) THEN
+          m = m+1
+          CYCLE
+        END IF
         l = l+1
       END DO
       TotMatrix % ParOffset(j+1) = l 
@@ -2697,8 +2643,11 @@ SUBROUTINE ParallelShrinkPerm()
     
     ! We can only make the ParBlockPerm if also BlockPerm exists!
     BlockPerm => TotMatrix % BlockPerm 
-    IF(.NOT. ASSOCIATED(BlockPerm) ) RETURN
-
+    IF(.NOT. ASSOCIATED(BlockPerm) ) THEN
+      IF(Halt) STOP
+      RETURN
+    END IF
+      
     ! Dense numbering for the block system dofs
     ALLOCATE(RenumPerm(n))
     RenumPerm = 1    
@@ -2715,6 +2664,15 @@ SUBROUTINE ParallelShrinkPerm()
     END DO
 
     ! Allocate the parallel permutation, if not already done
+    IF(ASSOCIATED(TotMatrix % ParBlockPerm)) THEN
+      k = SIZE(TotMatrix % ParBlockPerm)
+      IF(k /= l) THEN
+        CALL Warn('ParallelShrinkPerm','Different size for ParBlockPerm: '//I2S(k)//' vs. '//I2S(l))
+        PRINT *,'ParBlockPerm skipped2',ParEnv % MyPe, m
+        DEALLOCATE(TotMatrix % ParBlockPerm)
+        Halt = .TRUE.
+      END IF
+    END IF
     IF(.NOT. ASSOCIATED(TotMatrix % ParBlockPerm) ) THEN
       ALLOCATE(TotMatrix % ParBlockPerm(l))
     END IF
@@ -2727,6 +2685,8 @@ SUBROUTINE ParallelShrinkPerm()
       IF(RenumPerm(i)==0) CYCLE
       ParBlockPerm(RenumPerm(i)) = ShrinkPerm(BlockPerm(i))
     END DO
+
+    IF(Halt) STOP 
     
   END SUBROUTINE ParallelShrinkPerm
 
@@ -2788,8 +2748,10 @@ SUBROUTINE BlockMatrixVectorProdMono( u,v,ipar )
         i = MAXVAL( BlockPerm )
         j = MINVAL( BlockPerm )        
         IF( j <= 0 .OR. i > n ) THEN
-          CALL Fatal('BlockMatrixVectorProdMono','Invalid sizes!')
+          CALL Fatal('BlockMatrixVectorProdMono','Invalid sizes: '&
+              //I2S(i)//' '//I2S(j)//' '//I2S(n)//' '//I2S(SIZE(BlockPerm)))
         END IF
+        CALL Info('BlockMatrixVectorProdMono','BlockPermSizes: '//I2S(i)//' '//I2S(j)//' '//I2S(n))
       END IF
     END IF
 
@@ -2838,7 +2800,7 @@ SUBROUTINE CreateBlockMatrixScaling( )
 
     INTEGER :: i,j,k,l,n,m,NoVar,istat
     REAL(KIND=dp) :: nrm, tmp, blocknrm
-    TYPE(Matrix_t), POINTER :: A, Atrans
+    TYPE(Matrix_t), POINTER :: A
     REAL(KIND=dp), POINTER :: b(:), Diag(:), Values(:)
     LOGICAL :: ComplexMatrix, GotIt, DiagOnly, PrecScale
     INTEGER, POINTER :: Rows(:), Cols(:)
@@ -2917,15 +2879,6 @@ SUBROUTINE CreateBlockMatrixScaling( )
 
         IF(.NOT.Found) THEN
           A => TotMatrix % Submatrix(k,l) % Mat          
-          IF( A % NumberOfRows == 0 ) THEN
-            IF( TotMatrix % SubMatrixTranspose(l,k) ) THEN                    
-              CALL Info(Caller,'Creating the transpose for real!',Level=12)
-              Atrans => CRS_Transpose( TotMatrix % SubMatrix(l,k) % Mat ) 
-              TotMatrix % Submatrix(k,l) % Mat => Atrans
-              TotMatrix % SubMatrixTranspose(l,k) = .FALSE.
-              A => Atrans
-            END IF
-          END IF
         END IF
         
         IF(.NOT. ASSOCIATED( A  ) ) CYCLE
@@ -2995,43 +2948,118 @@ SUBROUTINE BlockMatrixInfo()
 !------------------------------------------------------------------------------
     IMPLICIT NONE
     INTEGER :: i,j,k,l,n,m,NoVar
+    INTEGER(KIND=8) :: ll
     REAL(KIND=dp) :: nrm, tmp, blocknrm
-    TYPE(Matrix_t), POINTER :: A, Atrans
+    TYPE(Matrix_t), POINTER :: A
     REAL(KIND=dp), POINTER :: b(:), Diag(:), Values(:)
     LOGICAL :: ComplexMatrix, GotIt, DiagOnly
-    INTEGER, POINTER :: Rows(:), Cols(:)
+    INTEGER, POINTER :: Rows(:), Cols(:)    
     LOGICAL :: Found
-    
-    
+    CHARACTER(:), ALLOCATABLE :: pre
+
+
+    CALL Info('BlockMatrixInfo','')
     CALL Info('BlockMatrixInfo','Showing some ranges of block matrix stuff',Level=10)
     
     NoVar = TotMatrix % NoVar
     m = 0
+
+    pre = 'BlockInfo'//I2S(ParEnv % MyPe)//': '
     
-    PRINT *,'BlockInfo:',NoVar   
-    
+    PRINT *,pre,NoVar, ParEnv % Mype
     DO k=1,NoVar
       DO l=1,NoVar
-        
+
         A => TotMatrix % SubMatrix(k,l) % Mat
         IF( .NOT. ASSOCIATED( A ) ) CYCLE
         IF( A % NumberOfRows == 0 ) CYCLE
-        
-        n = TotMatrix % offset(k+1) - TotMatrix % offset(k)
+
+        n = 0
+        m = 0
+
+        IF( ASSOCIATED( TotMatrix % Offset ) ) THEN
+          n = TotMatrix % offset(k+1) - TotMatrix % offset(k)
+        END IF
         IF(isParallel ) THEN
-          m = TotMatrix % ParOffset(k+1) -TotMatrix % ParOffset(k)
+          IF( ASSOCIATED( TotMatrix % ParOffset ) ) THEN
+            m = TotMatrix % ParOffset(k+1) -TotMatrix % ParOffset(k)
+          END IF
         END IF
-          
-        PRINT *,'BlockInfo:',k,l,A % NumberOfRows, n, m, A % COMPLEX
+
+        PRINT *,TRIM(pre)//' size',k,l,A % NumberOfRows, n, m, A % COMPLEX
 
         Rows   => A % Rows
         Cols   => A % Cols
         Values => A % Values
 
-        PRINT *,'BlockInfo: A range',SUM( Values ), MINVAL( Values ), MAXVAL( Values ) 
+        PRINT *,pre,'A'//I2S(10*k+l)//' range',SUM( Values ), SUM( ABS( Values ) ), &
+            MINVAL( Values ), MAXVAL( Values ) 
+
+        A => TotMatrix % SubMatrix(k,k) % PrecMat
+        IF( .NOT. ASSOCIATED( A ) ) CYCLE
+        IF( A % NumberOfRows == 0 ) CYCLE
+
+        Values => A % Values
+        PRINT *,TRIM(pre),'B'//I2S(11*k)//' range',SUM( Values ), SUM( ABS( Values ) ), &
+            MINVAL( Values ), MAXVAL( Values ) 
       END DO
     END DO
-    
+
+
+    DO k=1,NoVar      
+      PRINT *,'BlockInfo rhs:',k
+
+      A => TotMatrix % SubMatrix(k,k) % Mat
+      IF( ASSOCIATED( A ) ) THEN
+        b => A % rhs
+        IF(ASSOCIATED(b)) THEN
+          PRINT *,pre, 'A rhs range',k,SUM( b ), SUM( ABS( b ) ), &
+              MINVAL( b ), MAXVAL( b ) 
+        END IF
+      END IF
+
+      b => TotMatrix % SubVector(k) % rhs
+      IF(ASSOCIATED(b)) THEN
+        PRINT *,pre,'b range',k,SUM( b ), SUM( ABS( b ) ), &
+            MINVAL( b ), MAXVAL( b ) 
+      END IF
+    END DO
+
+    IF(isParallel) THEN
+      DO k=1,NoVar      
+        PRINT *,'BlockInfo ParallelInfo'//I2S(k)
+
+        A => TotMatrix % SubMatrix(k,k) % Mat
+        IF(.NOT. ASSOCIATED(A % ParallelInfo) ) THEN
+          PRINT *,pre,'No ParallelInfo associated', k
+          CYCLE
+        END IF
+
+        n = A % NumberOfRows
+        j = 0
+        ll = 0
+        IF( ASSOCIATED(A % ParallelInfo % NeighbourList ) ) THEN
+          i = SIZE( A % ParallelInfo % NeighbourList)
+          DO l=1,i
+            j = j+SIZE(A % ParallelInfo % NeighbourList(l) % Neighbours) 
+            ll = ll+SUM(A % ParallelInfo % NeighbourList(l) % Neighbours) 
+          END DO
+          PRINT *,pre,'BlockInfo NeighbourList:',n,i,j,ll
+        END IF
+        IF( ASSOCIATED(A % ParallelInfo % GInterface) ) THEN
+          i = SIZE( A % ParallelInfo % Ginterface)
+          j = COUNT(A % ParallelInfo % GInterface)
+          PRINT *,pre,'BlockInfo GInterface', n,i,j
+        END IF
+        IF( ASSOCIATED(A % ParallelInfo % GlobalDOFs) ) THEN
+          i = SIZE( A % ParallelInfo % GlobalDOFs)
+          ll = SUM(A % ParallelInfo % GlobalDOFs)
+          PRINT *,pre,'BlockInfo GlobalDofs', n,i,ll
+        END IF
+
+      END DO
+    END IF
+
   END SUBROUTINE BlockMatrixInfo
 !------------------------------------------------------------------------------
 
@@ -3176,7 +3204,7 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
 !---------------------------------------------------------------------------------
     REAL(KIND=dp), POINTER :: rtmp(:),vtmp(:),xtmp(:),btmp(:),diagtmp(:),b(:),x(:), a_rhs_save(:)
     REAL(KIND=dp), POINTER CONTIG :: rhs_save(:)
-    INTEGER :: i,j,k,l,n,m,NoVar,nc,kk,ll
+    INTEGER :: i,j,k,l,n,m,NoVar,nc,kk,ll,istat
     TYPE(Solver_t), POINTER :: Solver, Solver_save, ASolver
     INTEGER, POINTER :: Offset(:), BlockPerm(:),ParPerm(:)
     TYPE(ValueList_t), POINTER :: Params
@@ -3185,8 +3213,8 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
     TYPE(Variable_t), POINTER :: Var, Var_save
     REAL(KIND=dp) :: nrm
     LOGICAL :: GotOrder, BlockGS, Found, NS, ScaleSystem, DoSum, &
-        IsComplex, BlockScaling, DoDiagScaling, DoPrecScaling, UsePrecMat, Trans, &
-        Isolated, NoNestedScaling, DoAMGXmv, CalcLoads
+        IsComplex, BlockScaling, DoDiagScaling, DoPrecScaling, UsePrecMat, &
+        Isolated, NoNestedScaling, DoAMGXmv, CalcLoads, BlockSch
     CHARACTER(:), ALLOCATABLE :: str
     INTEGER(KIND=AddrInt) :: AddrFunc
     EXTERNAL :: AddrFunc
@@ -3210,6 +3238,8 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
     !---------------------------------------------------------------
     BlockOrder => ListGetIntegerArray( Params,'Block Order',GotOrder)
     BlockGS = ListGetLogical( Params,'Block Gauss-Seidel',Found)
+    BlockSch = ListGetLogical( Params,'Block Schur',Found)
+    
     
     NoVar = TotMatrix % NoVar
     Solver => TotMatrix % Solver
@@ -3223,7 +3253,8 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
     END IF
       
     IF( n /= Offset(NoVar+1) ) THEN
-      CALL Fatal('BlockMatrixPrec','There is a mismatch between sizes!')
+      CALL Warn('BlockMatrixPrec','There is a mismatch between sizes: '&
+          //I2S(n)//' vs. '//I2S(Offset(NoVar+1)))
     END IF
     
     ! Save the initial solver stuff
@@ -3258,7 +3289,8 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
 
     ! Always treat the inner iterations as truly complex if they are
     CALL ListAddLogical( Params,'Linear System Skip Complex',.FALSE.) 
-    
+    CALL ListAddLogical( Params,'Linear System Skip Loads',.TRUE.)
+
     IF (isParallel) THEN
       ALLOCATE( x(TotMatrix % MaxSize), b(TotMatrix % MaxSize) )
     END IF
@@ -3266,9 +3298,12 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
     ! Initial guess:
     !-----------------------------------------
     u(1:n) = 0.0_dp
-    
-    IF( BlockGS ) THEN
-      ALLOCATE( vtmp(n), rtmp(n), xtmp(n))
+
+    BlockSch = ListCheckPrefix( Params,'Schur Operator' )
+
+    IF( BlockGS .OR. BlockSch ) THEN
+      ALLOCATE( vtmp(n), rtmp(n), xtmp(n),STAT=istat)
+      IF ( istat /= 0 ) CALL Fatal('BlockMatrixPrec','Memory allocation error for wrk space!')
       vtmp(1:n) = v(1:n)
     END IF
     
@@ -3284,7 +3319,10 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
       WRITE(Message,'(A,I0)') 'Solving block: ',i
       CALL Info('BlockMatrixPrec',Message,Level=8)
 
-      CALL ListPushNameSpace('block '//i2s(i)//i2s(i)//':')
+      ! We do probably not want to compute the change within each iteration
+      CALL ListAddLogical( Params,'Skip Advance Nonlinear iter',.TRUE.)         
+      CALL ListAddLogical( Params,'Skip Compute Nonlinear Change',.TRUE.)         
+
       
       ! Set pointers to the new linear system
       !-------------------------------------------------------------------
@@ -3341,11 +3379,7 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
           IF( CRS_CopyMatrixPrec( TotMatrix % Submatrix(k,k) % Mat, A ) ) EXIT
         END DO
       END IF
-      
-      ! We do probably not want to compute the change within each iteration
-      CALL ListAddLogical( Asolver % Values,'Skip Advance Nonlinear iter',.TRUE.)         
-      CALL ListAddLogical( Asolver % Values,'Skip Compute Nonlinear Change',.TRUE.)         
-       
+             
       IF( BlockScaling ) CALL BlockMatrixScaling(.TRUE.,i,i,b,UsePrecMat)
 
       ! The special preconditioning matrices have not been scaled with the monolithic system.
@@ -3354,8 +3388,10 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
       DoPrecScaling = DoDiagScaling .AND. UsePrecMat
       IF( DoPrecScaling ) THEN
         n = A % NumberOfRows
-        ALLOCATE( btmp(n), diagtmp(n) )
-
+        ALLOCATE( btmp(n), diagtmp(n), STAT=istat)
+        IF ( istat /= 0 ) CALL Fatal('BlockMatrixPrec',&
+            'Memory allocation error for scaling wrk space!')
+        
         IF( TotMatrix % GotBlockStruct ) THEN
           k = TotMatrix % InvBlockStruct(i)
           IF( k <= 0 ) THEN
@@ -3383,7 +3419,6 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
         END IF
       END IF
               
-
       IF( InfoActive( 25 ) ) THEN
         CALL BlockMatrixInfo()
       END IF
@@ -3394,20 +3429,76 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
         nc = 1
       END IF
       
+      k = ListGetInteger( Params,'Schur Operator '//I2S(i),Found )
+      IF( k > 0 ) THEN
+        CALL Info('BlockMatrixPrec','Perform Schur complement operation for block '//I2S(i), Level=7 )
+
+        ! The residual is used only as a temporary vector
+        !-------------------------------------------------------------
+        Aij => TotMatrix % SubMatrix(i,k) % Mat 
+        IF( Aij % NumberOfRows == 0) CYCLE
+
+        ! r = P^T b
+        Aij => TotMatrix % Submatrix(k,i) % Mat            
+        IF( BlockGS ) THEN
+          b => vtmp(offset(i)+1:offset(i+1))
+        ELSE      
+          b => v(offset(i)+1:offset(i+1))
+        END IF
+        IF(DoAMGXMV) THEN
+          CALL AMGXMatrixVectorMultiply(Aij,b,rtmp,SolverRef )
+        ELSE
+          CALL CRS_MatrixVectorMultiply(Aij,b,rtmp )
+        END IF
+        
+        ! u = A^-1 r
+        CALL ListPushNameSpace('block '//i2s(11*k)//':')
+
+        Aij => TotMatrix % Submatrix(k,k) % Mat            
+        Isolated = TotMatrix % SubMatrix(k,k) % ParallelIsolatedMatrix 
+        IF(DoAMGXMv) THEN
+          ScaleSystem = ListGetLogical( Params,'Linear System Scaling', Found )
+          IF(.NOT. Found) ScaleSystem = .TRUE.
+          IF ( ScaleSystem ) CALL ScaleLinearSystem(ASolver, Aij,btmp,xtmp )
+          CALL AMGXSolver( Aij, xtmp, btmp, ASolver )
+          IF( ScaleSystem ) CALL BackScaleLinearSystem(ASolver,Aij,btmp,xtmp)
+        ELSE
+          CALL SolveLinearSystem( Aij, rtmp, xtmp, Var % Norm, Var % DOFs, ASolver )
+        END IF
+        
+        ! x = -P u 
+        rtmp = 0.0_dp
+        Aij => TotMatrix % Submatrix(i,k) % Mat            
+        IF(DoAMGXMV) THEN
+          CALL AMGXMatrixVectorMultiply(Aij,xtmp,rtmp,SolverRef )
+        ELSE
+          CALL CRS_MatrixVectorMultiply(Aij,xtmp,rtmp )
+        END IF
+
+        BLOCK 
+          REAL(KIND=dp) :: cmult
+          cmult = ListGetCReal( Params,'Schur Multiplier '//I2S(i),Found, DefValue=1.0_dp )          
+          ! Up-date the off-diagonal entries to the r.h.s. 
+          u(offset(i)+1:offset(i+1)) = u(offset(i)+1:offset(i+1)) &
+              - cmult * rtmp(1:offset(i+1)-offset(i))
+        END BLOCK
+        CALL ListPopNameSpace() ! block ij:
 
-      IF(DoAMGXMv) THEN
-        ScaleSystem = ListGetLogical( Params,'Linear System Scaling', Found )
-        IF(.NOT. Found) ScaleSystem = .TRUE.
-        IF ( ScaleSystem ) CALL ScaleLinearSystem(ASolver, A,btmp,x )
-        CALL AMGXSolver( A, x, btmp, ASolver )
-        IF( ScaleSystem ) CALL BackScaleLinearSystem(ASolver,A,btmp,x)
+        CYCLE
       ELSE
-        CalcLoads = ListGetLogical( ASolver % Values, 'Calculate Loads', Found )
-        CALL ListAddLogical( ASolver % Values, 'Calculate Loads', .FALSE.)
-        CALL SolveLinearSystem( A, btmp, x, Var % Norm, Var % DOFs, ASolver )
-        IF (CalcLoads) CALL ListAddLogical( ASolver % Values, 'Calculate Loads', .TRUE.)        
+        CALL ListPushNameSpace('block '//i2s(11*i)//':')
+      
+        IF(DoAMGXMv) THEN
+          ScaleSystem = ListGetLogical( Params,'Linear System Scaling', Found )
+          IF(.NOT. Found) ScaleSystem = .TRUE.
+          IF ( ScaleSystem ) CALL ScaleLinearSystem(ASolver, A,btmp,x )
+          CALL AMGXSolver( A, x, btmp, ASolver )
+          IF( ScaleSystem ) CALL BackScaleLinearSystem(ASolver,A,btmp,x)
+        ELSE
+          CALL SolveLinearSystem( A, btmp, x, Var % Norm, Var % DOFs, ASolver )
+        END IF
       END IF
-
+      
       ! If this was a special preconditioning matrix then update the solution in the scaled system. 
       IF( DoPrecScaling ) THEN
         x(1:n) = x(1:n) / diagtmp(1:n)
@@ -3447,33 +3538,11 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
           
           ! The residual is used only as a temporary vector
           !-------------------------------------------------------------
-          Trans = .FALSE.
-          Isolated = .FALSE.          
           Aij => TotMatrix % SubMatrix(k,i) % Mat 
 
-          ! We may use transpose that is actually never created as only its action is used!
-          IF( Aij % NumberOfRows == 0 ) THEN
-            IF( TotMatrix % SubMatrixTranspose(i,k) ) THEN
-              Aij => TotMatrix % SubMatrix(i,k) % Mat
-              IF( Aij % NumberOfRows >  0 ) THEN
-                Trans = .TRUE.
-                CALL Info('BlockMatrixPrec','Multiplying with transpose of submatrix ('&
-                    //I2S(i)//','//I2S(k)//')',Level=8)
-                IF(isParallel ) THEN
-                  IF( TotMatrix % SubMatrix(i,k) % ParallelIsolatedMatrix ) THEN
-                    Isolated = .TRUE.
-                  ELSE
-                    CALL Fatal('BlockMatrixPrec','Only isolated matrices may be transposed in parallel!')
-                  END IF
-                END IF
-              END IF
-            END IF
-          ELSE
-            Isolated = TotMatrix % SubMatrix(k,i) % ParallelIsolatedMatrix 
-          END IF
           IF( Aij % NumberOfRows == 0) CYCLE
+          Isolated = TotMatrix % SubMatrix(k,i) % ParallelIsolatedMatrix 
           
-
           IF (isParallel) THEN
             IF(ASSOCIATED(Aij % ParMatrix)) THEN
               ! x is packed, r is full 
@@ -3482,12 +3551,7 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
             ELSE IF( Isolated ) THEN
               ! If our matrix is not active on shared nodes we may apply serial Mv
               ! and pack the results to include only the dofs owned by the partition.
-              IF( .NOT. Trans ) THEN
-                CALL CRS_MatrixVectorMultiply(Aij,x,rtmp)
-              ELSE
-                CALL CRS_TransposeMatrixVectorMultiply( Aij, x, rtmp )
-              END IF
-
+              CALL CRS_MatrixVectorMultiply(Aij,x,rtmp)
               ParPerm => TotMatrix % Submatrix(i,i) % ParPerm
 
 #if 0
@@ -3523,14 +3587,10 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
             END IF
           ELSE
             Aij => TotMatrix % Submatrix(k,i) % Mat            
-            IF( .NOT. Trans ) THEN
-                IF(DoAMGXMV) THEN
-                  CALL AMGXMatrixVectorMultiply(Aij, x, rtmp, SolverRef )
-                ELSE
-                  CALL CRS_MatrixVectorMultiply(Aij,x,rtmp )
-                END IF
+            IF(DoAMGXMV) THEN
+              CALL AMGXMatrixVectorMultiply(Aij, x, rtmp, SolverRef )
             ELSE
-              CALL CRS_TransposeMatrixVectorMultiply( Aij, x, rtmp )
+              CALL CRS_MatrixVectorMultiply(Aij,x,rtmp )
             END IF
           END IF
 
@@ -3540,7 +3600,7 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
 
         END DO ! l=j+1,NoVar
         
-      END IF
+      END IF  ! Gauss-Seidel
 
       CALL ListPopNameSpace() ! block ij:
       
@@ -3551,15 +3611,16 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
     CALL ListPopNameSpace('block:') ! block:
 
     CALL ListAddLogical( Params,'Linear System Refactorize',.FALSE. )
-    CALL ListAddLogical( Asolver % Values,'Skip Advance Nonlinear iter',.FALSE.)
-    CALL ListAddLogical( Asolver % Values,'Skip Compute Nonlinear Change',.FALSE.)
-    
+    CALL ListAddLogical( Params,'Skip Advance Nonlinear iter',.FALSE.)
+    CALL ListAddLogical( Params,'Skip Compute Nonlinear Change',.FALSE.)
+    CALL ListAddLogical( Params,'Linear System Skip Loads',.FALSE.)
+
     Solver => Solver_save
     Solver % Matrix => mat_save
     Solver % Matrix % RHS => rhs_save
     Solver % Variable => Var_save
 
-    IF( BlockGS ) THEN
+    IF( BlockGS .OR. BlockSch ) THEN
       DEALLOCATE( vtmp, rtmp, xtmp ) 
     END IF
 
@@ -3683,7 +3744,7 @@ SUBROUTINE BlockStandardIter( Solver, MaxChange )
         ALLOCATE( dx( SIZE( Var % Values ) ) )
         dx = 0.0_dp
 
-        CALL ListPushNamespace('block '//i2s(RowVar)//i2s(RowVar)//':')          
+        CALL ListPushNamespace('block '//i2s(11*RowVar)//':')          
 
         IF( BlockScaling ) CALL BlockMatrixScaling(.TRUE.,i,i,b)
               
@@ -3769,7 +3830,7 @@ SUBROUTINE BlockKrylovIter( Solver, MaxChange )
     INTEGER :: Rounds, OutputInterval, PolynomialDegree
     INTEGER, POINTER :: Offset(:),poffset(:),BlockStruct(:),ParPerm(:)
     INTEGER :: i,j,k,l,ia,ib,istat
-    LOGICAL :: LS, BlockAV,BlockAVOld,Found, UseMono
+    LOGICAL :: LS, BlockAV,Found, UseMono
     CHARACTER(*), PARAMETER :: Caller = 'BlockKrylovIter'
  
     
@@ -3778,14 +3839,14 @@ SUBROUTINE BlockKrylovIter( Solver, MaxChange )
     !CALL ListPushNameSpace('outer:')
     Params => Solver % Values
     
-    BlockAVOld = ListGetLogical(Params,'Block A-V System Old', Found)
     BlockAV = ListGetLogical(Params,'Block A-V System', Found)
 
     ndim = TotMatrix % TotSize 
     NoVar = TotMatrix % NoVar
 
     TotMatrix % NoIters = 0
-    
+
+    isParallel = (ParEnv % PEs > 1)
     offset => TotMatrix % Offset
     IF(isParallel) THEN
       poffset => TotMatrix % ParOffset
@@ -3845,6 +3906,7 @@ SUBROUTINE BlockKrylovIter( Solver, MaxChange )
         x = 0.0_dp
         b = 0.0_dp
       END IF
+      CALL Info(Caller,'Parallel initializations done!',Level=25)
     END IF
       
     CALL Info(Caller,'Initializing monolithic system vectors',Level=18)
@@ -3892,8 +3954,13 @@ SUBROUTINE BlockKrylovIter( Solver, MaxChange )
       
     prevXnorm = CompNorm(b,ndim)
     WRITE( Message,'(A,ES12.5)') 'Rhs norm at start: ',PrevXnorm
-    CALL Info(Caller,Message,Level=10)
+    CALL Info(Caller,Message,Level=10)    
+    IF( PrevXNorm < EPSILON( PrevXNorm ) ) THEN
+      CALL Warn(Caller,"With zero'ish r.h.s. it does not make sense to call the solver!")
+      RETURN
+    END IF
 
+    
     prevXnorm = CompNorm(x,ndim)
     WRITE( Message,'(A,ES12.5)') 'Solution norm at start: ',PrevXnorm
     CALL Info(Caller,Message,Level=10)
@@ -3919,27 +3986,17 @@ SUBROUTINE BlockKrylovIter( Solver, MaxChange )
           Solver,ndim=ndim,MatvecF=mvProc,PrecF=precProc,&
           DotF=AddrFunc(SParDotProd), NormF=AddrFunc(SParNorm))
 
-      IF (BlockAVOld) THEN        
-        IF(ASSOCIATED(SolverMatrix)) THEN
-          ! Communicate the packed solution among all partitions on the shared nodes
-          SolverMatrix % ParMatrix % SplittedMatrix % TmpXvec = x(1:poffset(NoVar+1))
-          CALL ParallelUpdateResult(SolverMatrix,x,r)
-        ELSE
-          CALL Fatal(Caller,'No matrix to apply.')
-        END IF
-      ELSE
-        DO i=1,NoVar
-          TotMatrix % SubMatrix(i,i) % Mat % ParMatrix % SplittedMatrix % &
-              TmpXvec = x(poffset(i)+1:poffset(i+1))
-        END DO
+      DO i=1,NoVar
+        TotMatrix % SubMatrix(i,i) % Mat % ParMatrix % SplittedMatrix % &
+            TmpXvec = x(poffset(i)+1:poffset(i+1))
+      END DO
 
-        ! Communicate blockwise information.
-        ! After this the packed x is again a full vector with redundant shared dofs
-        DO i=1,NoVar
-          CALL ParallelUpdateResult(TotMatrix % SubMatrix(i,i) % Mat, &
-              x(offset(i)+1:offset(i+1)), r )
-        END DO
-      END IF
+      ! Communicate blockwise information.
+      ! After this the packed x is again a full vector with redundant shared dofs
+      DO i=1,NoVar
+        CALL ParallelUpdateResult(TotMatrix % SubMatrix(i,i) % Mat, &
+            x(offset(i)+1:offset(i+1)), r )
+      END DO
     ELSE
       IF( ListGetLogical( Params,'Linear System test',Found ) ) THEN
         CALL IterSolver( A,x,b,&
@@ -3978,7 +4035,7 @@ SUBROUTINE BlockKrylovIter( Solver, MaxChange )
     IF( TotMatrix % GotBlockStruct ) THEN
       SVar => CurrentModel % Solver % Variable
       Svar % Values(TotMatrix % BlockPerm) = x
-    ELSE IF (BlockAV .OR. BlockAVOld ) THEN
+    ELSE IF (BlockAV ) THEN
       SolverVar % Values = x
     END IF
 
@@ -4329,7 +4386,7 @@ SUBROUTINE BlockSolveInt(A,x,b,Solver)
     INTEGER :: i,j,k,l,n,nd,NonLinIter,tests,NoTests,iter
     LOGICAL :: GotIt, GotIt2, BlockPrec, BlockAV, &
         BlockHdiv, BlockHcurl, BlockHorVer, BlockCart, BlockNodal, BlockDomain, &
-        BlockDummy, BlockComplex, BlockAVOld
+        BlockDummy, BlockComplex, SkipPrec
     INTEGER :: ColVar, RowVar, NoVar, BlockDofs, VarDofs
     
     REAL(KIND=dp) :: NonlinearTol, Norm, PrevNorm, Residual, PrevResidual, &
@@ -4343,7 +4400,7 @@ SUBROUTINE BlockSolveInt(A,x,b,Solver)
     INTEGER, POINTER :: VarPerm(:)
     INTEGER, POINTER :: BlockPerm(:)
     INTEGER, POINTER :: SlaveSolvers(:)
-    LOGICAL :: GotSlaveSolvers, SkipVar
+    LOGICAL :: GotSlaveSolvers, DoMyOwnVars,OldMatrix
     
     
     TYPE(Matrix_t), POINTER :: Amat, SaveCM
@@ -4363,6 +4420,11 @@ SUBROUTINE BlockSolveInt(A,x,b,Solver)
 
     isParallel = ParEnv % PEs > 1
     
+    OldMatrix = ASSOCIATED( Solver % BlockMatrix )
+    IF( OldMatrix ) THEN
+      CALL Info('BlockSolveInit','Using old block matrix structures!',Level=12)
+    END IF
+    
     
     ! Determine some parameters related to the block strategy
     !------------------------------------------------------------------------------
@@ -4377,7 +4439,6 @@ SUBROUTINE BlockSolveInt(A,x,b,Solver)
     BlockScaling = ListGetLogical( Params,'Block Scaling',GotIt)
 
     ! Different strategies on how to split the initial monolithic matrix into blocks
-    BlockAVOld = ListGetLogical( Params,'Block A-V System Old', GotIt)
     BlockAV = ListGetLogical( Params,'Block A-V System', GotIt)
     BlockHcurl = ListGetLogical( Params,'Block Quadratic Hcurl System', GotIt)   
     BlockHdiv = ListGetLogical( Params,'Block Hdiv system',GotIt)
@@ -4391,61 +4452,72 @@ SUBROUTINE BlockSolveInt(A,x,b,Solver)
     SlaveSolvers =>  ListGetIntegerArray( Params, &
          'Block Solvers', GotSlaveSolvers )
 
-    SkipVar = .FALSE.
-    IF( BlockDomain ) THEN
-      n = MAXVAL( Solver % Variable % Perm )
-      ALLOCATE( BlockIndex(n) )
-      BlockDofs = 0
-      CALL BlockPickDofsPhysical( PSolver, BlockIndex, BlockDofs )  
-      SkipVar = .TRUE.
-    ELSE IF( BlockHdiv ) THEN
-      n = MAXVAL( Solver % Variable % Perm )
-      ALLOCATE( BlockIndex(n) )
-      BlockDofs = 0
-      CALL BlockPickHdiv( PSolver, BlockIndex, BlockDofs )  
-      SkipVar = .TRUE.
-    ELSE IF( BlockAV ) THEN
-      n = MAXVAL( Solver % Variable % Perm )
+    DoMyOwnVars = .FALSE.
+    
+    IF( BlockDomain .OR. BlockHdiv .OR. BlockHcurl .OR. BlockAV .OR. BlockHorVer .OR. &
+        BlockCart .OR. BlockNodal ) THEN
+      ! These take monolithic splitting and only return the "BlockIndex" which tells to which
+      ! block each dof belongs to. Then the same routine can split the matrices for all. 
+      !-----------------------------------------------------------------------------------------------
+      n = SIZE( Solver % Variable % Values )      
       ALLOCATE( BlockIndex(n) )
+      BlockIndex = 0
       BlockDofs = 0
-      CALL BlockPickAV( PSolver, BlockIndex, BlockDofs )  
-      SkipVar = .TRUE.
-    ELSE IF( BlockAVOld .OR. BlockNodal .OR. BlockHorVer .OR. BlockHcurl ) THEN
-      BlockDofs = 2
-      IF( ListGetLogical( Params,'Block Quadratic Hcurl Faces',Found ) ) BlockDofs = 3
-      IF(BlockComplex) THEN
-        BlockDofs = 2 * BlockDofs 
-        IF( ListGetLogical( Params,'Block Quadratic Hcurl semicomplex',Found ) ) BlockDofs = 3
+      DoMyOwnVars = .TRUE.
+      
+      IF( BlockDomain ) THEN
+        CALL BlockPickDofsPhysical( PSolver, BlockIndex, BlockDofs )  
+      ELSE IF( BlockHdiv ) THEN
+        CALL BlockPickHdiv( PSolver, BlockIndex, BlockDofs )  
+      ELSE IF( BlockHCurl ) THEN
+        CALL BlockPickMatrixHcurl( PSolver, BlockDofs, BlockComplex, BlockIndex )
+      ELSE IF( BlockAV ) THEN
+        CALL BlockPickAV( PSolver, BlockIndex, BlockDofs )  
+      ELSE IF( BlockNodal ) THEN
+        CALL BlockPickMatrixNodal( PSolver, BlockDofs, BlockComplex, BlockIndex )
+      ELSE IF(BlockHorVer .OR. BlockCart) THEN
+        CALL BlockPickMatrixHorVer( PSolver, BlockDofs, BlockCart, BlockIndex )
       END IF
-      SkipVar = .TRUE.
-    ELSE IF( BlockCart ) THEN
-      BlockDofs = 3
-      SkipVar = .TRUE.
-    ELSE IF( GotSlaveSolvers ) THEN
-      BlockDofs = SIZE( SlaveSolvers )
-    ELSE IF( BlockDummy ) THEN
-      BlockDofs = 1
+
     ELSE
-      BlockDofs = Solver % Variable % Dofs      
+      IF( BlockCart ) THEN
+        BlockDofs = 3
+      ELSE IF( GotSlaveSolvers ) THEN
+        BlockDofs = SIZE( SlaveSolvers )
+      ELSE IF( BlockDummy ) THEN
+        BlockDofs = 1
+      ELSE
+        BlockDofs = Solver % Variable % Dofs      
+      END IF
     END IF
-    VarDofs = BlockDofs
 
+    VarDofs = BlockDofs
+      
     
     HaveConstraint = 0
     IF ( ASSOCIATED(A % ConstraintMatrix) )  HaveConstraint = 1
     HaveConstraint = ParallelReduction(HaveConstraint)
-     
+    
     HaveAdd = 0
     IF ( ASSOCIATED(A % AddMatrix) )  THEN
       IF ( A % AddMatrix % NumberOFRows > 0 ) HaveAdd = 1
     END IF
     HaveAdd = ParallelReduction(HaveAdd)
 
-    IF( HaveConstraint > 0 ) BlockDofs = BlockDofs + 1
-    IF( HaveAdd > 0 ) BlockDofs = BlockDofs + 1    
-   
-    CALL BlockInitMatrix( Solver, TotMatrix, BlockDofs, VarDofs, SkipVar )
-    
+    IF( HaveConstraint > 0 ) THEN
+      CALL Info('BlockSolveInt','Block system has ConstraintMatrix!',Level=10)
+      BlockDofs = BlockDofs + 1
+    END IF
+
+    IF( HaveAdd > 0 ) THEN
+      CALL Info('BlockSolveInt','Block system has AddMatrix!',Level=10)
+      BlockDofs = BlockDofs + 1    
+    END IF
+
+    IF(.NOT. OldMatrix ) THEN
+      CALL BlockInitMatrix( Solver, TotMatrix, BlockDofs, VarDofs, DoMyOwnVars )
+    END IF
+      
     NoVar = TotMatrix % NoVar
     TotMatrix % Solver => Solver
     
@@ -4459,42 +4531,30 @@ SUBROUTINE BlockSolveInt(A,x,b,Solver)
 
     SaveRHS => SolverMatrix % RHS
     SolverMatrix % RHS => b
-    
+
+
     IF( .NOT. GotSlaveSolvers ) THEN
       CALL Info('BlockSolveInt','Splitting monolithic matrix into pieces',Level=10)
-      IF( BlockDomain .OR. BlockHdiv .OR. BlockAV ) THEN
-        CALL BlockPickMatrixPerm( Solver, BlockIndex, VarDofs )
-        !DEALLOCATE( BlockIndex ) 
-      ELSE IF( BlockAVOld ) THEN
-        CALL BlockPickMatrixAV( Solver, VarDofs )
-      ELSE IF( BlockHcurl ) THEN
-        CALL BlockPickMatrixHcurl( Solver, VarDofs, BlockComplex )
-        IF(VarDofs == 3 ) BlockComplex = .FALSE.
-      ELSE IF( BlockHorVer .OR. BlockCart ) THEN
-        CALL BlockPickMatrixHorVer( Solver, VarDofs, BlockCart )       
-      ELSE IF( BlockNodal ) THEN
-        CALL BlockPickMatrixNodal( Solver, VarDofs )        
-      ELSE IF( BlockDummy .OR. VarDofs == 1 ) THEN
-        CALL Info('BlockSolveInt','Using the original matrix as the (1,1) block!',Level=10)
-        TotMatrix % SubMatrix(1,1) % Mat => SolverMatrix        
-        TotMatrix % SubMatrix(1,1) % Mat % Complex = ListGetLogical(Params,'Linear System Complex',Found)        
-      ELSE
-        CALL BlockPickMatrix( Solver, NoVar ) !VarDofs )
-        VarDofs = NoVar
-      END IF
-
-      IF( SkipVar ) THEN
+      IF( DoMyOwnVars ) THEN
+        CALL BlockPickMatrixPerm( Solver, BlockIndex, VarDofs )      
         IF(ASSOCIATED(BlockIndex)) THEN
           CALL BlockInitVar( Solver, TotMatrix, BlockIndex )
           DEALLOCATE(BlockIndex)
         ELSE
           CALL BlockInitVar( Solver, TotMatrix )
         END IF
+      ELSE IF( BlockDummy .OR. VarDofs == 1 ) THEN
+        CALL Info('BlockSolveInt','Using the original matrix as the (1,1) block!',Level=10)
+        TotMatrix % SubMatrix(1,1) % Mat => SolverMatrix        
+        TotMatrix % SubMatrix(1,1) % Mat % Complex = ListGetLogical(Params,'Linear System Complex',Found)        
+      ELSE
+        ! Default splitting of matrices. 
+        CALL BlockPickMatrix( Solver, NoVar ) 
+        VarDofs = NoVar
       END IF
-
-      CALL BlockPrecMatrix( Solver, VarDofs ) 
       CALL Info('BlockSolveInt','Block matrix system created',Level=12)
     END IF
+
     
     ! Currently we cannot have both structure-structure and fluid-structure couplings!
     IF( ListGetLogical( Solver % Values,'Structure-Structure Coupling',Found ) ) THEN
@@ -4504,12 +4564,17 @@ SUBROUTINE BlockSolveInt(A,x,b,Solver)
     END IF
         
     IF( HaveConstraint > 0 ) THEN
-      CALL BlockPickConstraint( Solver, VarDofs )
+      ! Do not try to create preconditioner if we have Schur operator requested. 
+      SkipPrec =  ListCheckPrefix( Params,'Schur Operator' )
+      CALL BlockPickConstraint( Solver, VarDofs, SkipPrec )
       ! Storing pointer to CM so that the SolverMatrix won't be treated with the constraints
       SaveCM => Solver % Matrix % ConstraintMatrix
       Solver % Matrix % ConstraintMatrix => NULL()
     END IF
-
+    
+    ! Create preconditioners for block matrices.
+    CALL BlockPrecMatrix( Solver, NoVar )
+    
     Found = .FALSE.
     DO RowVar=1,NoVar
       Amat => TotMatrix % SubMatrix(RowVar,RowVar) % Mat
@@ -4529,9 +4594,8 @@ SUBROUTINE BlockSolveInt(A,x,b,Solver)
         CALL Warn('BlockSolveInt','Could not continue without Perm in parallel!!')
       END IF
     END IF
-      
 
-    IF (isParallel) THEN
+    IF (isParallel .AND. .NOT. OldMatrix ) THEN
       CALL Info('BlockSolveInt','Initializing parallel block matrices',Level=12)
       DO RowVar=1,NoVar
         DO ColVar=1,NoVar
@@ -4566,8 +4630,8 @@ SUBROUTINE BlockSolveInt(A,x,b,Solver)
             END IF
             !CALL SParIterActiveBarrier()
           ELSE
-            CALL Info('BlockSolverInt',&
-                'Block matrix is not square matrix - cannot initialize parallel structures!',Level=20)
+            CALL Info('BlockSolverInt','Block matrix '//I2S(10*RowVar+ColVar)//&
+                ' is not square matrix - cannot initialize parallel structures!',Level=20)
           END IF
         END DO
       END DO
@@ -4577,6 +4641,14 @@ SUBROUTINE BlockSolveInt(A,x,b,Solver)
       Solver % Variable  => SolverVar
       CALL Info('BlockSolveInt','Initialization of block matrix finished',Level=20)
     END IF
+
+
+    IF(InfoActive(25)) THEN
+      CALL Info('BlockSolveInt','Initial Block matrix system information:')
+      CALL  BlockMatrixInfo()
+    END IF
+   
+
     
     !------------------------------------------------------------------------------
     ! Finally solve the system using 'outer: ' as the optional namespace
@@ -4642,7 +4714,7 @@ SUBROUTINE BlockSolveInt(A,x,b,Solver)
       Solver % Matrix % ConstraintMatrix => SaveCM 
     END IF
 
-    IF( BlockHorVer .OR. BlockCart .OR. BlockDomain .OR. BlockHcurl ) THEN
+    IF( DoMyOwnVars ) THEN
       CALL BlockBackCopyVar( Solver, TotMatrix )
     END IF
 
@@ -4654,6 +4726,7 @@ SUBROUTINE BlockSolveInt(A,x,b,Solver)
     CALL Info('BlockSolveInt','-------------------------------------------------',Level=5)
 
   END SUBROUTINE BlockSolveInt
+  
 END MODULE BlockSolve