[petsc-dev] question on MatGolubKahanComputeExplicitOperator()

[Stefano Zampini]

We usually try to avoid these expensive checks in PETSc (unless they are really needed they are run in debug mode only);


> On Jun 24, 2020, at 8:23 AM, Barry Smith <bsmith at petsc.dev> wrote:
> 
> 
> MatGolubKahanComputeExplicitOperator() first verifies that C is the transpose of B, it then computes the transpose of B to do the MatMatMult().
> 
> 1) Since C is the transpose of B isn't the second transpose not needed? Just use C instead of the BT? 
> 
> 2) If one knows the C = B' one can save the expensive check.  Presumably if the entire large matrix is known to be symmetric one can skip the transpose check? The user
> can provide this information with MatSetOption() when they build the matrix.
> 
>  Thanks
> 
>   Barry
> 
> 
> static PetscErrorCode MatGolubKahanComputeExplicitOperator(Mat A,Mat B,Mat C,Mat *H,PetscReal gkbnu)
> {
>  PetscErrorCode    ierr;
>  Mat               BT,T;
>  PetscReal         nrmT,nrmB;
> 
>  PetscFunctionBegin;
>  ierr = MatHermitianTranspose(C,MAT_INITIAL_MATRIX,&T);CHKERRQ(ierr);            /* Test if augmented matrix is symmetric */
>  ierr = MatAXPY(T,-1.0,B,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
>  ierr = MatNorm(T,NORM_1,&nrmT);CHKERRQ(ierr);
>  ierr = MatNorm(B,NORM_1,&nrmB);CHKERRQ(ierr);
>  if (nrmB > 0) {
>    if (nrmT/nrmB >= PETSC_SMALL) {
>      SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix is not symmetric/hermitian, GKB is not applicable.");
>    }
>  }
>  /* Compute augmented Lagrangian matrix H = A00 + nu*A01*A01'. This corresponds to */
>  /* setting N := 1/nu*I in [Ar13].                                                 */
>  ierr = MatHermitianTranspose(B,MAT_INITIAL_MATRIX,&BT);CHKERRQ(ierr);
>  ierr = MatMatMult(B,BT,MAT_INITIAL_MATRIX,PETSC_DEFAULT,H);CHKERRQ(ierr);       /* H = A01*A01'          */
>  ierr = MatAYPX(*H,gkbnu,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);             /* H = A00 + nu*A01*A01' */
> 
>  ierr = MatDestroy(&BT);CHKERRQ(ierr);
>  ierr = MatDestroy(&T);CHKERRQ(ierr);
>  PetscFunctionReturn(0);
> }
>