[petsc-dev] question on MatGolubKahanComputeExplicitOperator()

[Barry Smith]

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);
}