c++ - 使用带有 INSERT_VALUES 的 MatSetValuesStencil 替换整行

Question

我正在使用 Petsc Ksp 例程。我使用 MatSetValuesStencil 构造了一个运算符，在此函数的每次调用中，我指定一个长度为 5 的行矩阵值。有时我需要将一行从 5 长度模板完全替换为 3 长度模板。INSERT_VALUES 模式会将两个值保留在未更改的位置上还是会将它们丢弃为零？

Answer 1

idxm未在参数和idxn函数中指定的矩阵元素MatSetValuesStencil(...)保持不变，即使INSERT_VALUES使用。

这是一个从ksp_ex29开始的小代码来测试它：

static char help[] = "Does INSERT_VALUES changes the whole row ? No.\n\n";

#include <petscdm.h>
#include <petscdmda.h>
#include <petscksp.h>

extern PetscErrorCode ComputeMatrix42(DM da,Mat jac);
extern PetscErrorCode ComputeMatrix(DM da,Mat jac);

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **argv)
{

    DM             da;
    PetscErrorCode ierr;
    Mat matrix;

    PetscInitialize(&argc,&argv,(char*)0,help);

    ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,-3,-3,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);CHKERRQ(ierr);

    DMCreateMatrix(da,&matrix);

    ComputeMatrix(da,matrix);

    PetscPrintf(PETSC_COMM_WORLD,"A matrix of negative terms : \n");
    MatView(matrix, PETSC_VIEWER_STDOUT_WORLD );

    ComputeMatrix42(da,matrix);

    PetscPrintf(PETSC_COMM_WORLD,"The diagonal, i-1 and i+1 are set to 42 : \n");
    MatView(matrix, PETSC_VIEWER_STDOUT_WORLD );

    ierr = DMDestroy(&da);CHKERRQ(ierr);
    ierr = MatDestroy(&matrix);CHKERRQ(ierr);
    ierr = PetscFinalize();

    return 0;
}

#undef __FUNCT__
#define __FUNCT__ "ComputeMatrix"
PetscErrorCode ComputeMatrix(DM da,Mat jac)
{
    PetscErrorCode ierr;
    PetscInt       i,j,mx,my,xm,ym,xs,ys;
    PetscScalar    v[5];
    MatStencil     row, col[5];

    PetscFunctionBeginUser;
    ierr      = DMDAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
    ierr      = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
    for (j=ys; j<ys+ym; j++) {
        for (i=xs; i<xs+xm; i++) {
            row.i = i; row.j = j;
            v[0] = -1;              col[0].i = i;   col[0].j = j-1;
            v[1] = -1;              col[1].i = i-1; col[1].j = j;
            v[2] = -13; col[2].i = i;   col[2].j = j;
            v[3] = -1;              col[3].i = i+1; col[3].j = j;
            v[4] = -1;              col[4].i = i;   col[4].j = j+1;
            ierr = MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);CHKERRQ(ierr);
        }
    }

    ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

    PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "ComputeMatrix42"
PetscErrorCode ComputeMatrix42(DM da,Mat jac)
{
    PetscErrorCode ierr;
    PetscInt       i,j,mx,my,xm,ym,xs,ys;
    PetscScalar    v[3];
    MatStencil     row, col[3];

    PetscFunctionBeginUser;
    ierr      = DMDAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
    ierr      = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
    for (j=ys; j<ys+ym; j++) {
        for (i=xs; i<xs+xm; i++) {
            row.i = i; row.j = j;
            v[0] = 42;              col[0].i = i-1; col[0].j = j;
            v[1] = 42; col[1].i = i;   col[1].j = j;
            v[2] = 42;              col[2].i = i+1; col[2].j = j;
            ierr = MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);CHKERRQ(ierr);
        }
    }

    ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

    PetscFunctionReturn(0);
}

使用以下 makefile 编译它：

include ${PETSC_DIR}/conf/variables
include ${PETSC_DIR}/conf/rules

main: main.o  chkopts
    -${CLINKER} -o main main.o  ${PETSC_LIB}
    ${RM} main.o

输出 ：

A matrix of negative terms : 
Mat Object: 1 MPI processes
  type: seqaij
row 0: (0, -13)  (1, -1)  (2, -1)  (3, -1)  (6, -1) 
row 1: (0, -1)  (1, -13)  (2, -1)  (4, -1)  (7, -1) 
row 2: (0, -1)  (1, -1)  (2, -13)  (5, -1)  (8, -1) 
row 3: (0, -1)  (3, -13)  (4, -1)  (5, -1)  (6, -1) 
row 4: (1, -1)  (3, -1)  (4, -13)  (5, -1)  (7, -1) 
row 5: (2, -1)  (3, -1)  (4, -1)  (5, -13)  (8, -1) 
row 6: (0, -1)  (3, -1)  (6, -13)  (7, -1)  (8, -1) 
row 7: (1, -1)  (4, -1)  (6, -1)  (7, -13)  (8, -1) 
row 8: (2, -1)  (5, -1)  (6, -1)  (7, -1)  (8, -13) 
The diagonal, i-1 and i+1 are set to 42 : 
Mat Object: 1 MPI processes
  type: seqaij
row 0: (0, 42)  (1, 42)  (2, 42)  (3, -1)  (6, -1) 
row 1: (0, 42)  (1, 42)  (2, 42)  (4, -1)  (7, -1) 
row 2: (0, 42)  (1, 42)  (2, 42)  (5, -1)  (8, -1) 
row 3: (0, -1)  (3, 42)  (4, 42)  (5, 42)  (6, -1) 
row 4: (1, -1)  (3, 42)  (4, 42)  (5, 42)  (7, -1) 
row 5: (2, -1)  (3, 42)  (4, 42)  (5, 42)  (8, -1) 
row 6: (0, -1)  (3, -1)  (6, 42)  (7, 42)  (8, 42) 
row 7: (1, -1)  (4, -1)  (6, 42)  (7, 42)  (8, 42) 
row 8: (2, -1)  (5, -1)  (6, 42)  (7, 42)  (8, 42)

我正在使用 PETSC 3.5.2。