@@ -3923,29 +3923,42 @@ namespace PETSc {
39233923 delete dA;
39243924 }
39253925 delete empty;
3926- if (c == 1 ) {
3927- MatSetType (ptA->_petsc , MATSHELL);
3928- User< LinearSolver< Dmat > > user = nullptr ;
3929- PetscNew (&user);
3930- user->op = nullptr ;
3931- const Polymorphic* op = nargs[0 ] ? dynamic_cast < const Polymorphic* >(nargs[0 ]) : nullptr ;
3932- if (op) {
3933- const OneOperator* codeAt = op->Find (" (" , ArrayOfaType (atype< KN< PetscScalar >* >( ), false ));
3934- if (codeAt) {
3935- user->op = new LinearSolver< Dmat >::MatF_O (ptB->_last - ptB->_first , stack, codeA, -1 , codeAt);
3936- MatShellSetOperation (ptA->_petsc , MATOP_MULT_TRANSPOSE,
3937- (PetscErrorCodeFn *)Op_User< LinearSolver< Mat >, Mat, ' T' >);
3938- }
3939- }
3940- if (!user->op ) user->op = new LinearSolver< Dmat >::MatF_O (ptB->_last - ptB->_first , stack, codeA);
3941- MatShellSetContext (ptA->_petsc , user);
3942- MatShellSetOperation (ptA->_petsc , MATOP_MULT,
3943- (PetscErrorCodeFn *)Op_User< LinearSolver< Mat >, Mat >);
3944- MatShellSetOperation (ptA->_petsc , MATOP_DESTROY, (PetscErrorCodeFn *)ShellDestroy< LinearSolver< Dmat > >);
3945- MatSetUp (ptA->_petsc );
3926+ } else if (ptB->_petsc ) {
3927+ if (c != 1 ) {
3928+ PetscBool assembled;
3929+
3930+ MatAssembled (ptB->_petsc , &assembled);
3931+ MatDuplicate (ptB->_petsc , assembled ? MAT_COPY_VALUES : MAT_DO_NOT_COPY_VALUES, &ptA->_petsc );
3932+ } else {
3933+ PetscInt m, M, n, N;
3934+
3935+ MatGetLocalSize (ptB->_petsc , &m, &n);
3936+ MatGetSize (ptB->_petsc , &M, &N);
3937+ MatCreate (PetscObjectComm ((PetscObject)ptB->_petsc ), &ptA->_petsc );
3938+ MatSetSizes (ptA->_petsc , m, n, M, N);
3939+ }
3940+ }
3941+ if (c == 1 ) {
3942+ ffassert (ptA->_petsc );
3943+ MatSetType (ptA->_petsc , MATSHELL);
3944+ User< LinearSolver< Dmat > > user = nullptr ;
3945+ PetscNew (&user);
3946+ user->op = nullptr ;
3947+ const Polymorphic* op = nargs[0 ] ? dynamic_cast < const Polymorphic* >(nargs[0 ]) : nullptr ;
3948+ if (op) {
3949+ const OneOperator* codeAt = op->Find (" (" , ArrayOfaType (atype< KN< PetscScalar >* >( ), false ));
3950+ if (codeAt) {
3951+ user->op = new LinearSolver< Dmat >::MatF_O (ptB->_last - ptB->_first , stack, codeA, -1 , codeAt);
3952+ MatShellSetOperation (ptA->_petsc , MATOP_MULT_TRANSPOSE,
3953+ (PetscErrorCodeFn *)Op_User< LinearSolver< Mat >, Mat, ' T' >);
3954+ }
39463955 }
3947- } else if (ptB->_petsc ) {
3948- MatDuplicate (ptB->_petsc , MAT_COPY_VALUES, &ptA->_petsc );
3956+ if (!user->op ) user->op = new LinearSolver< Dmat >::MatF_O (ptB->_last - ptB->_first , stack, codeA);
3957+ MatShellSetContext (ptA->_petsc , user);
3958+ MatShellSetOperation (ptA->_petsc , MATOP_MULT,
3959+ (PetscErrorCodeFn *)Op_User< LinearSolver< Mat >, Mat >);
3960+ MatShellSetOperation (ptA->_petsc , MATOP_DESTROY, (PetscErrorCodeFn *)ShellDestroy< LinearSolver< Dmat > >);
3961+ MatSetUp (ptA->_petsc );
39493962 }
39503963 if (c == 0 && nargs[0 ] && GetAny< bool >((*nargs[0 ])(stack))) ptK->destroy ( );
39513964 return ptA;
0 commit comments