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