petsc 矩阵测试

来源：互联网发布：linux 播放 4k 视频编辑：程序博客网时间：2024/05/19 19:40

第一个PETSc矩阵测试:

static char help[] = " matrix parallel partition test\n\n"; #include <petscksp.h> int main(int argc, char **args){ Mat A; PetscInt dim, i, j, Ii, J, Istart, Iend, n=6; PetscScalar v ; PetscInitialize(&argc, &args, (char*)0,help); PetscOptionsGetReal(NULL, "-sigma1", &sigma1, NULL); PetscOptionsGetInt(NULL, "-n", &n, NULL); dim = n * n; MatCreate(PETSC_COMM_WORLD, &A); MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim); MatSetFromOptions(A); MatSetUp(A); MatGetOwnershipRange(A,&Istart,&Iend); for (Ii = Istart; Ii < Iend; Ii++) {        v = 1.0; i = Ii/n; j = Ii - i * n;        if(i > 0) {                J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);        }        if(i < n - 1) {                J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);        }        if(j > 0) {                J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);        }        if(j < n - 1) {                J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);        }} MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); PetscPrintf(PETSC_COMM_WORLD, "Matrix output\n"); MatView(A,PETSC_VIEWER_STDOUT_WORLD); MatDestroy(A); PetscFinalize(); return 0;}

PETSC中矩阵按行划分分给各节点，对每个节点中的矩阵子块进行赋值。

输出矩阵：

row 0: (1, 1) (4, 1)
row 1: (0, 1) (2, 1) (5, 1)
row 2: (1, 1) (3, 1) (6, 1)
row 3: (2, 1) (7, 1)
row 4: (0, 1) (5, 1) (8, 1)
row 5: (1, 1) (4, 1) (6, 1) (9, 1)
row 6: (2, 1) (5, 1) (7, 1) (10, 1)
row 7: (3, 1) (6, 1) (11, 1)
row 8: (4, 1) (9, 1) (12, 1)
row 9: (5, 1) (8, 1) (10, 1) (13, 1)
row 10: (6, 1) (9, 1) (11, 1) (14, 1)
row 11: (7, 1) (10, 1) (15, 1)
row 12: (8, 1) (13, 1)
row 13: (9, 1) (12, 1) (14, 1)
row 14: (10, 1) (13, 1) (15, 1)
row 15: (11, 1) (14, 1)

注意：

首先，MatGetOwnershipRange 和 MatSetValue 之后，任何一个cpu节点实际上是可以修改不属于它自身“管辖”范围内矩阵元素的。其中的数据通信，属于petsc封装（省了麻烦）。

本例中给值的矩阵元素 pattern 如下：

（I， J+1）

（I-1 ， J）（I + 1, J）

（I， J-1）

if (i > 0) 语句完成了对当前元素(I,J ) 正下方邻接元素（I， J -1）的赋值;

if (i < n - 1) 语句对 (I, J +1) 赋值。

同理，其他。

这样就形成了“条状”矩阵元素带。显然，对于本程序用了4 cpu节点（按行划分），即整体矩阵的每四行属于一个节点。对当前cpu节点（比如p0)，显然不仅仅对整体“带状矩阵”的前四行相关位置赋值了，而是所有满足if条件的位置，不管该位置是不是在当前cpu节点管辖范围内。