[petsc-dev] question on MatGolubKahanComputeExplicitOperator()

Barry Smith bsmith at petsc.dev
Wed Jun 24 00:23:00 CDT 2020

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.



static PetscErrorCode MatGolubKahanComputeExplicitOperator(Mat A,Mat B,Mat C,Mat *H,PetscReal gkbnu)
  PetscErrorCode    ierr;
  Mat               BT,T;
  PetscReal         nrmT,nrmB;

  ierr = MatHermitianTranspose(C,MAT_INITIAL_MATRIX,&T);CHKERRQ(ierr);            /* Test if augmented matrix is symmetric */
  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);

More information about the petsc-dev mailing list