Improved block solution with constraints. Take away the annoying poss…

…ibility to use transpose of constraints as an operator.

fem/src/BlockSolve.F90

@@ -82,7 +82,7 @@ FUNCTION CreateSchurApproximation(A,P,Q) RESULT ( S )
         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) )
+          CALL List_AddToMatrixElement(S % ListMatrix, i, k2, -val * Q % Values(k2) )
         END DO
       END DO
     END DO
@@ -93,11 +93,6 @@ FUNCTION CreateSchurApproximation(A,P,Q) RESULT ( S )
     WRITE(Message,'(A,ES12.3)') 'Schur matrix increase factor: ',val
     CALL Info('CreateSchurApproximation',Message)
 
-    val = ListGetCReal( CurrentModel % Simulation,'c test',Found)
-    IF(Found  ) THEN   
-      S % Values = val * S  % Values
-    END IF
-
     !ALLOCATE( S % rhs(S % NumberOfRows))
     !S % rhs = 0.0_dp
 
@@ -648,7 +643,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            
@@ -760,7 +755,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
 
@@ -1055,7 +1051,7 @@ 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.
   !-------------------------------------------------------------------------------------
@@ -1186,7 +1182,8 @@ SUBROUTINE BlockPickMatrixAV( Solver, NoVar )
 
 
   END SUBROUTINE BlockPickMatrixAV
-
+#endif
+
 
   !-------------------------------------------------------------------------------------
   !> Picks vertical and horizontal components of a full matrix.
@@ -2080,8 +2077,6 @@ 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
@@ -3228,7 +3223,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
@@ -3313,6 +3308,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) )
@@ -3321,9 +3318,12 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
     ! Initial guess:
     !-----------------------------------------
     u(1:n) = 0.0_dp
-
+
+    BlockSch = ListCheckPrefix( Params,'Schur Operator' )
+
     IF( BlockGS .OR. BlockSch ) THEN
-      ALLOCATE( vtmp(n), rtmp(n), xtmp(n))
+      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
 
@@ -3339,8 +3339,6 @@ 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)//':')
-
       ! Set pointers to the new linear system
       !-------------------------------------------------------------------
       Var => TotMatrix % SubVector(i) % Var
@@ -3409,8 +3407,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
@@ -3438,7 +3438,6 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
         END IF
       END IF
 
-
       IF( InfoActive( 25 ) ) THEN
         CALL BlockMatrixInfo()
       END IF
@@ -3449,7 +3448,6 @@ 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 )
@@ -3458,7 +3456,6 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
         !-------------------------------------------------------------
         Aij => TotMatrix % SubMatrix(i,k) % Mat 
         IF( Aij % NumberOfRows == 0) CYCLE
-        Isolated = TotMatrix % SubMatrix(i,k) % ParallelIsolatedMatrix 
 
         ! r = P^T b
         Aij => TotMatrix % Submatrix(k,i) % Mat            
@@ -3472,25 +3469,22 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
         ELSE
           CALL CRS_MatrixVectorMultiply(Aij,b,rtmp )
         END IF
-
+        
         ! u = A^-1 r
-        CALL ListPopNamespace()
-        CALL ListPushNameSpace('block '//i2s(k)//i2s(k)//':')
-        Aij => TotMatrix % Submatrix(k,k) % Mat            
+        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
-          CalcLoads = ListGetLogical( ASolver % Values, 'Calculate Loads', Found )
-          CALL ListAddLogical( ASolver % Values, 'Calculate Loads', .FALSE.)
-          CALL SolveLinearSystem( Aij, btmp, xtmp, Var % Norm, Var % DOFs, ASolver )
-          IF (CalcLoads) CALL ListAddLogical( ASolver % Values, 'Calculate Loads', .TRUE.)        
+          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            
@@ -3500,31 +3494,29 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
           CALL CRS_MatrixVectorMultiply(Aij,xtmp,rtmp )
         END IF
 
-        ! 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)) &
-            - rtmp(1:offset(i+1)-offset(i))
-
-
+        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:
+
         CYCLE
       ELSE
-
+        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
-          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.)        
         END IF
-
       END IF
-
-
 
       ! If this was a special preconditioning matrix then update the solution in the scaled system. 
       IF( DoPrecScaling ) THEN
@@ -3570,7 +3562,6 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
           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 
@@ -3630,7 +3621,6 @@ SUBROUTINE BlockMatrixPrec( u,v,ipar )
 
       END IF  ! Gauss-Seidel
 
-
       CALL ListPopNameSpace() ! block ij:
 
     END DO ! j=1,NoVar
@@ -3640,9 +3630,10 @@ 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
@@ -3772,7 +3763,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)
 
@@ -3858,7 +3849,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'
 
 
@@ -3867,7 +3858,6 @@ 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 
@@ -4008,27 +3998,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,&
@@ -4067,7 +4047,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
 
@@ -4418,7 +4398,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
     INTEGER :: ColVar, RowVar, NoVar, BlockDofs, VarDofs
 
     REAL(KIND=dp) :: NonlinearTol, Norm, PrevNorm, Residual, PrevResidual, &
@@ -4466,7 +4446,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)
@@ -4499,7 +4478,7 @@ SUBROUTINE BlockSolveInt(A,x,b,Solver)
       BlockDofs = 0
       CALL BlockPickAV( PSolver, BlockIndex, BlockDofs )  
       SkipVar = .TRUE.
-    ELSE IF( BlockAVOld .OR. BlockNodal .OR. BlockHorVer .OR. BlockHcurl ) THEN
+    ELSE IF( BlockNodal .OR. BlockHorVer .OR. BlockHcurl ) THEN
       BlockDofs = 2
       IF( ListGetLogical( Params,'Block Quadratic Hcurl Faces',Found ) ) BlockDofs = 3
       IF(BlockComplex) THEN
@@ -4561,8 +4540,6 @@ SUBROUTINE BlockSolveInt(A,x,b,Solver)
       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.