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