基于PETSC的线性方程求解

最新推荐文章于 2023-03-13 22:16:18 发布

Galerkin码农选手

最新推荐文章于 2023-03-13 22:16:18 发布

阅读量2.5k

点赞数 2

分类专栏：高性能计算

本文链接：https://blog.csdn.net/forrestguang/article/details/109312367

版权

高性能计算专栏收录该内容

44 篇文章 29 订阅 ¥239.90 ¥99.00

订阅专栏

使用KSP方法求解小型线性系统

需要注意的是，在PETSC中出现的语法和C语言基本相同，区别大部分都只存在于不同的表达方式
比如在C中，int表示整型，在PETSC中，PetscInt表示整型，这里大家百度用法，或者参考我上次的学习报告，我们这里重点介绍在PETSC中一些简单函数的定义调用，并且使用这种函数来求解实际问题（求解线性系统）。
我们慢慢来解释：第一句static char help[] = “Solve a 4x4 linear system using KSP.\n”;是为了对这个代码进行说明，比如说这个代码在做什么。
然后我们看到下面两句：
extern PetscErrorCode matrix(Mat);// form A
extern PetscErrorCode right(Vec);//form b
这两句在申明函数，函数的定义被我放在了最后，第一个函数的参数是Mat，表示这个函数接受矩阵作为输入，第二个函数接受向量Vec作为输入。

//STARTWHOLE
static char help[] = "Solve a 4x4 linear system using KSP.\n";

#include <petsc.h>

extern PetscErrorCode matrix(Mat);// form A
extern PetscErrorCode right(Vec);//form b

int main(int argc,char **args) {
    PetscErrorCode ierr;
    Vec x,b;
    Mat A;
    KSP ksp;
    

    PetscInitialize(&argc,&args,NULL,help);

    ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr);
    ierr = VecSetSizes(b,PETSC_DECIDE,4);CHKERRQ(ierr);
    ierr = VecSetFromOptions(b);CHKERRQ(ierr);
    ierr = right(b);CHKERRQ(ierr);// set the value of A
    ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
    ierr = VecAssemblyEnd(b);CHKERRQ(ierr);

    ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
    ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,4,4);CHKERRQ(ierr);
    ierr = MatSetFromOptions(A);CHKERRQ(ierr);
    ierr = MatSetUp(A);CHKERRQ(ierr);
    ierr = matrix(A);CHKERRQ(ierr);//set the value of A
    ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

    ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
    ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
    ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
    ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
    ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
    ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);


    KSPDestroy(&ksp);MatDestroy(&A);
    VecDestroy(&x);VecDestroy(&b);
    return PetscFinalize();
}
//ENDWHOLE
PetscErrorCode matrix(Mat A) {
    PetscErrorCode ierr;
    PetscInt i,j[4] = {0,1,2,3};
    PetscReal aA[4][4] = {{1.0,0.0,0.0,0.0},
		    	  {0.0,1.0,3.0,0.0},
		    	  {0.0,0.0,2.0,0.0},
		    	  {0.0,0.0,3.0,4.0}};
    for (i = 0; i < 4; i++) {
	ierr = MatSetValues(A,1,&i,4,j,aA[i],INSERT_VALUES);CHKERRQ(ierr);
    }
    
    return 0;
}
PetscErrorCode right(Vec b) {
    PetscErrorCode ierr;
    PetscInt i;
    PetscReal *ab;
 
    ierr = VecGetArray(b,&ab);CHKERRQ(ierr);
    for (i = 0; i < 4; i++) {
	ab[i] = i;
    }
    ierr = VecRestoreArray(b,&ab);CHKERRQ(ierr);
    return 0;
}

ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr);
ierr = VecSetSizes(b,PETSC_DECIDE,4);CHKERRQ(ierr);
ierr = VecSetFromOptions(b);CHKERRQ(ierr);
ierr = right(b);CHKERRQ(ierr);// set the value of A
ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
ierr = VecAssemblyEnd(b);CHKERRQ(ierr);

这一大段在创建向量，其中除了第4句是在调用函数设置向量里面的元素值以外，其他的都是在创建向量中必不可少的部分。
同样的，后面两段在创建Mat（矩阵）对象和KSP对象，KSP是Krylov Space Method的意思，就是使用这种方法来求解线性系统，至于这种方法是什么，自己百度吧，我的学习文档里有介绍一些。
链接：https://pan.baidu.com/s/1ZoEeGsS6VfciUzU67Qudfw
提取码：svww
如果大家需要更详细的关于PETSC的介绍
link.
里面涉及到矩阵和向量创建赋值的细节，参考官网手册和我的文档。

makefile

这里需要指出重点，我们的这个PETSC代码假设以fun.c存储为C文件，用普通的gcc编译器无法编译，因为这个涉及到PETSC的内容，大家需要参考我上一篇文章下载好PETSC包。同时我们需要makefile文件来帮助编译，makefile的作用大家可以上网查，这里我们针对这个C文件给出makefile的内容。
在这里插入图片描述

include ${PETSC_DIR}/lib/petsc/conf/variables
include ${PETSC_DIR}/lib/petsc/conf/rules
CFLAGS += -pedantic -std=c99
#f和fu都是使用KSP方法求解线性系统

tri: tri.o
	-${CLINKER} -o tri tri.o  ${PETSC_LIB}
	${RM} tri.o
fun: fun.o
	-${CLINKER} -o fun fun.o  ${PETSC_LIB}
	${RM} fun.o
fu: fu.o
	-${CLINKER} -o fu fu.o  ${PETSC_LIB}
	${RM} fu.o
f: f.o
	-${CLINKER} -o f f.o  ${PETSC_LIB}
	${RM} f.o



runtri_1:
	-@../testit.sh tri "-a_mat_view ::ascii_dense" 1 1

runtri_2:
	-@../testit.sh tri "-tri_m 1000 -ksp_rtol 1.0e-4 -ksp_type cg -pc_type bjacobi -sub_pc_type jacobi -ksp_converged_reason" 2 2
runfun_1:
	-@../testit.sh fun "-a_mat_view ::ascii_dense" 1 1

runfun_2:
	-@../testit.sh fun "-fun_m 1000 -ksp_rtol 1.0e-4 -ksp_type cg -pc_type bjacobi -sub_pc_type jacobi -ksp_converged_reason" 2 2

runfu_1:
	-@../testit.sh fu "-a_mat_view ::ascii_dense" 1 1

runfu_2:
	-@../testit.sh fu "-fun_m 1000 -ksp_rtol 1.0e-10 -ksp_type cg -pc_type bjacobi -sub_pc_type jacobi -ksp_converged_reason" 2 2
runf_1:
	-@../testit.sh f "-a_mat_view ::ascii_dense" 1 1

runf_2:
	-@../testit.sh f "-fun_m 1000 -ksp_rtol 1.0e-10 -ksp_type cg -pc_type bjacobi -sub_pc_type jacobi -ksp_converged_reason" 2 2



test_tri: runtri_1 runtri_2
test_fun: runfun_1 runfun_2
test_fu: runfu_1 runfu_2
test_f: runf_1 runf_2


test: test_tri test_fun test_fu test_f

# etc

.PHONY: distclean runtri_1 runtri_2 runfun_1 runfun_2 runfu_1 runfu_2 runf_1 runf_2 test test_tri test_fun test_fu test_f

distclean:
	@rm -f *~ tri fun fu f *tmp

任意大小的线性系统求解

首先我们需要设定好矩阵的大小m，然后开始搭建矩阵和右端项，然后设定误差分析，打印函数值和误差。

//STARTWHOLE
static char help[] = "Solve a mxm linear system using KSP.\n"
" use KSP 方法"
"注意到当m太大的时候，矩阵求解出现异常";

#include <petsc.h>

extern PetscErrorCode matrix(Mat , PetscInt);// form A
extern PetscErrorCode right(Vec ,PetscInt);//form b
extern PetscErrorCode exact(Vec ,PetscInt);//form xexact

int main(int argc,char **args) {
    PetscErrorCode ierr;
    Vec x,b,xexact;
    Mat A;
    KSP ksp;
    PetscInt m = 4;
    PetscReal err,errnorm;

    PetscInitialize(&argc,&args,NULL,help);

    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"f_","options for f",""); CHKERRQ(ierr);
    ierr = PetscOptionsInt("-m","dimension of linear system","f.c",m,&m,NULL); CHKERRQ(ierr);//create an integer
    ierr = PetscOptionsEnd(); CHKERRQ(ierr);

    ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr);
    ierr = VecSetSizes(b,PETSC_DECIDE,m);CHKERRQ(ierr);
    ierr = VecSetFromOptions(b);CHKERRQ(ierr);
    ierr = right(b ,m);CHKERRQ(ierr);// set the value of b
    ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
    ierr = VecAssemblyEnd(b);CHKERRQ(ierr);

    

    ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
    ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,m);CHKERRQ(ierr);
    ierr = MatSetOptionsPrefix(A,"a_"); CHKERRQ(ierr);
    ierr = MatSetFromOptions(A);CHKERRQ(ierr);
    ierr = MatSetUp(A);CHKERRQ(ierr);
    ierr = matrix(A ,m);CHKERRQ(ierr);//set the value of A
    ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

    ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
    ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
    ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
    ierr = VecDuplicate(b,&x);CHKERRQ(ierr);// make best of the location of b to store the x
    ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
/***
    ierr = PetscPrintf(PETSC_COMM_WORLD,
    "error for m = %d system , the numerical solution is \n",m); CHKERRQ(ierr);

    ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);//这行命令会在执行过程中打印x
***/
//必须要把xexact这部分放到x下面
    ierr = VecDuplicate(x,&xexact); CHKERRQ(ierr);// make best of the location of b to store the xexact
    ierr = exact(xexact ,m);CHKERRQ(ierr);// set the value of xexact
    ierr = VecAssemblyBegin(xexact);CHKERRQ(ierr);
    ierr = VecAssemblyEnd(xexact);CHKERRQ(ierr);
/***
    ierr = PetscPrintf(PETSC_COMM_WORLD,
    "error for m = %d system , the exact solution is \n",m); CHKERRQ(ierr);

    ierr = VecView(xexact,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
***/
    ierr = VecAXPY(x,-1.0,xexact); CHKERRQ(ierr);
    ierr = VecNorm(x,NORM_2,&errnorm); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,
    "error for m = %d system is |x-xexact|_2 = %.1e\n",m,errnorm); CHKERRQ(ierr);

    ierr = VecAXPY(xexact,0.0,x); CHKERRQ(ierr);
    ierr = VecNorm(xexact,NORM_2,&err); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,
    "error for m = %d system is |x-xexact|_2/|xexact|_2 = %.1e\n",m,errnorm/err); CHKERRQ(ierr);



    KSPDestroy(&ksp);MatDestroy(&A);
    VecDestroy(&x);VecDestroy(&b);//VecDestroy(&xexact);
    return PetscFinalize();
}
//ENDWHOLE
PetscErrorCode matrix(Mat A , PetscInt m) {
    PetscErrorCode ierr;
    PetscInt i,j[4];
    PetscReal v[4];

    for(i = 0; i < m;i++) {
	if (i == 0) {
	    j[0] = 0,j[1] = 1;
	    v[0] = 6.0,v[1] = 1.0;
	    ierr = MatSetValues(A,1,&i,2,j,v,INSERT_VALUES);CHKERRQ(ierr);
	}
	else {
	    j[0] = i - 1,j[1] = i,j[2] = i + 1;
	    v[0] = 8.0,v[1] = 6.0,v[2] = 1.0;
	    if (i == m - 1) {
		ierr = MatSetValues(A,1,&i,2,j,v,INSERT_VALUES);CHKERRQ(ierr);
	    }
	    else {
		ierr = MatSetValues(A,1,&i,3,j,v,INSERT_VALUES);CHKERRQ(ierr);
	    }
	}
    }
    
    return 0;
}
PetscErrorCode right(Vec b , PetscInt m) {
    PetscErrorCode ierr;
    PetscInt i;
    PetscReal *ab;
 
    ierr = VecGetArray(b,&ab);CHKERRQ(ierr);
    for (i = 0; i < m; i++) {
	if (i == 0) {
	    ab[i] = 7.0;
	}
	else {
	    if (i == m - 1) {
		ab[i] = 14.0;
	    }
	    else {
	        ab[i] = 15.0;
	    }
	}
    }
    ierr = VecRestoreArray(b,&ab);CHKERRQ(ierr);
    return 0;
}

PetscErrorCode exact(Vec xexact , PetscInt m) {
    PetscErrorCode ierr;
    PetscInt i;
    PetscReal *ab;
 
    ierr = VecGetArray(xexact,&ab);CHKERRQ(ierr);
    for (i = 0; i < m; i++) {
	if (i == 0) {
	    ab[i] = 1.0;
	}
	else {
	    if (i == m - 1) {
		ab[i] = 1.0;
	    }
	    else {
	        ab[i] = 1.0;
	    }
	}
    }
    ierr = VecRestoreArray(xexact,&ab);CHKERRQ(ierr);
    return 0;
}

这段代码只有一个部分需要重点注意，那就是如何设定m。

    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"f_","options for f",""); CHKERRQ(ierr);
    ierr = PetscOptionsInt("-m","dimension of linear system","f.c",m,&m,NULL); CHKERRQ(ierr);//create an integer
    ierr = PetscOptionsEnd(); CHKERRQ(ierr);

就是这3行代码，我们假设C文件是f.c，我们执行代码的过程是make f
./f -f_m 10
表示设定矩阵大小为 $10\times10$ ，其他的都一样。

Galerkin码农选手

关注

2
点赞
踩
15

收藏

觉得还不错? 一键收藏
打赏
7
评论
基于PETSC的线性方程求解

使用KSP方法求解小型线性系统需要注意的是，在PETSC中出现的语法和C语言基本相同，区别大部分都只存在于不同的表达方式比如在C中，int表示整型，在PETSC中，PetscInt表示整型，这里大家百度用法，或者参考我上次的学习报告，我们这里重点介绍在PETSC中一些简单函数的定义调用，并且使用这种函数来求解实际问题（求解线性系统）。我们慢慢来解释：第一句static char help[] = “Solve a 4x4 linear system using KSP.\n”;是为了对这个代码进行说明
复制链接

扫一扫