この記事は、実際に、PETScの関数を用いて、Ax=bの方程式を解くやり方についての備忘録です。
今回は、大規模疎行列に対して非常に有効な解法であるICCG法( 不完全コレスキー分解付き共役勾配法 )をPETScを用いて簡単に実装したいと思います。

私もそうですが、シミュレーションのコードでPETScのようなライブラリを使用する場合、行列 A やベクトル b の値を保持する配列がすでに存在しており、解くべきxに値する配列(例えば圧力など)もすでに存在しています。

したがって、PETScを用いる関数の外で定義された A と b をPETScが計算できる形にして(VecとMat)計算し、そこで得た解を配列に変換してやるという作業が必要になります。この記事では特にこの部分について、解説します。

また、筆者は数値計算の専門家でもなければ、プログラミングの専門家でもないので、この実装をするととりあえず動くというくらいの認識でお願いします( どうしても英語ドキュメントを読みたくない人向け？？ )。

[ 9/9 追記 ]
src/ksp/ksp/tutorials/ex13.c.htmlでは、Cのarray から PETScのVecに代入せずに計算する方法が記載してあります。(ただし、Vecの方の操作のみ)
Vec側については、こちらの VecCreateSeqWithArray と、VecPlaceArrayとVecResetArrayの３つの関数を用いた方法が良いと思われます。
また、Matについても、MatCreateSeqAIJWithArraysが用意されているので、こちらもご参照ください。

環境

この記事では、公式マニュアルの1.2章でも使用されている$PETSC_DIR/src/ksp/ksp/tutorials/ex1.cを題材にしていきます。

また、今回使用するPETScの環境は"MacでPETScをMPI無しで使う"の記事にある通り、mpi無しで./configureを行っています。(MPIありでも結果は変わらないと思います。)

コンパイルは同様にMacでPETScをMPI無しで使う"の記事で示したのと同様のmakeファイルを用いています。

今回解く方程式は以下の簡単な三重対角行列を解いていきます。

(画像はPC: Preconditionersより引用)

main関数

今回使用するmainプログラムは以下の通りです。( 余計な関数を使わなかったので、長くなって申し訳ないです。)

main.c

#include <petscksp.h>

int main(int argc,char **args){
   double *A,*x,*b;
   int n = 10; 
   int i,j;

   A = malloc(n * n * sizeof(double));
   x = malloc(n * sizeof(double));
   b = malloc(n * sizeof(double));

   for(i=0; i<n; i++){
      x[i] = 0.0;
   }

   printf("---------x before using petsc-------------\n");
   for(i=0; i<n; i++){
      printf("%lf \n",x[i]);
   }
   // set b 
   b[0] = 1.0 ; b[n-1] = 1.0;
   for(i=1; i<n-1; i++){
      b[i] = 0.0;
   }
   printf("--------- b -------------\n");
   for(i=0;i<n;i++){
      printf("%lf \n",b[i]);
   }
   printf("\n");

  //Assemble matrix A [a_ii = 2.0, a_i(i-1) = a_i(i+1) = -1.0]
   A[0] = 2.0 ;A[1] = -1.0; for(i=2;i<n;i++){A[i] = 0.0;} //row 0 
   for(i=1; i<n-1; i++){
     for(j=0; j<n; j++){
        if(j == i-1 || j == i+1){ A[n*i + j] = -1.0; }
        else if(j == i){ A[n*i + j] = 2.0; }
        else { A[n*i + j] = 0.0; } 
      }
   }
   for(i=0; i<n-2;i++){ A[n*(n-1) + i] = 0.0; } A[n*n -2] = -1.0; A[n*n -1] = 2.0; //row n-1

   printf("---------- A --------------\n");
   for(i=0;i<n;i++){
      for(j=0;j<n;j++){
         printf("%lf ",A[n * i + j]);
      }
      printf("\n");
   }

   calcAxb(argc,args,A,x,b,n); //here use PETSc library

   printf("-----------x after using petsc------------\n");
   for(int i=0; i<n; i++){
      printf("%lf \n",x[i]);
   }
   free(x);free(A);free(b);

   return 0;
}

ここでやっているのは、

A と b の定義と x の初期化
それらの確認(print)
calcAxbでPETScを用いて、Ax=bを解く

です。

calcAxb関数

中身はこんな感じです。(ex1.c のプログラムを少しいじって作りました。)

calcAxb.c


PetscErrorCode calcAxb(int argc,char **args,double *A_input, double *x_input, double *b_input,int n)
{
  Vec            x, b, u;      /* approx solution, RHS, exact solution */
  Mat            A;            /* linear system matrix */
  KSP            ksp;          /* linear solver context */
  PC             pc;           /* preconditioner context */
  PetscReal      norm;         /* norm of solution error */
  PetscErrorCode ierr;
  PetscInt       its,*ix;
  PetscMPIInt    size;
  
//ix = PetscInt *indices in VecSetValues(),MatSetValues(),VecGetValues()
  ix = malloc(sizeof(int) * n); 
  for(int k=0; k<n; k++){
     ix[k] = k;
  }

// must do this procedure for uniprocessor 
  ierr = PetscInitialize(&argc,&args,(char*)0,help); if (ierr) return ierr; 
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!");
  ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
         Compute the matrix and right-hand-side vector that define
         the linear system, Ax = b.
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  /* Create vectors. */
  ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) x, "Solution");CHKERRQ(ierr);
  ierr = VecSetSizes(x,PETSC_DECIDE,n);CHKERRQ(ierr);
  ierr = VecSetFromOptions(x);CHKERRQ(ierr); //must, but I don't know this function's mean 
  ierr = VecDuplicate(x,&b);CHKERRQ(ierr); //let b same size of vector as x (not copy)
  ierr = VecDuplicate(x,&u);CHKERRQ(ierr); 

   // set vector values
   VecSetValues(b,n,ix,b_input,INSERT_VALUES); //b_input --> b 

   //If it is finished to insert value in vector, these two functions are necessary.
   ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
   ierr = VecAssemblyEnd(b);CHKERRQ(ierr);

   printf("---- b in petsc -------\n");
   ierr = VecView(b,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

  /* Create matrix.  When using MatCreate(), the matrix format can be specified at runtime. */
 
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = MatSetUp(A);CHKERRQ(ierr);

  // Assemble matrix [a_ii = 2.0, a_i(j-1) = a_i(j+1) = -1.0] 
  ierr = MatSetValues(A,n,ix,n,ix,A_input,INSERT_VALUES); CHKERRQ(ierr);
  ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr);
  
  /* A is symmetric. Set symmetric flag to enable ICC/Cholesky preconditioner */
  ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
  
  /*  These two funtion is needed if assemblying matirix is finished. */
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  printf("----- A in petsc -------\n");
  ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
  
  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
                Create the linear solver and set various options
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  /*
     Set operators. Here the matrix that defines the linear system
     also serves as the matrix that defines the preconditioner.
  */
  ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr); //set A is the matrix to solve 
  
/* Set linear solver defaults for this problem (optional). */
  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
  ierr = PCSetType(pc,PCICC);CHKERRQ(ierr);
  ierr = KSPSetTolerances(ksp,1.e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); //誤差の設定
  ierr = KSPSetType(ksp,KSPCG); CHKERRQ(ierr);

  /* Set runtime options */
 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
                      Solve the linear system
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);

  /* View solver info */
  ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"View matrix and vector \n");CHKERRQ(ierr);
  ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
  
  /* get the x vector  (x --> x_input) */
  ierr = VecGetValues(x,n,ix,x_input);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                      Check the solution and clean up
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = VecSet(u,1.0);CHKERRQ(ierr); //u[i] = 1.0 
  ierr = VecAXPY(x,-1.0,u);CHKERRQ(ierr); //x = x + (-1.0)*u
  ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr); 
  ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr); 
  ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);CHKERRQ(ierr);

  /* Free work space */
  ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&u);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = KSPDestroy(&ksp);CHKERRQ(ierr)
  
  /* finish */
  ierr = PetscFinalize();
  return 0;
}

いじった部分としては

Vec b への値の代入
Mat A への値の代入
解となるVec xを関数外に出す
ソルバーはCG(共役勾配法、CG法)、前処理(PC)はICC(不完全コレスキー分解)

になるでしょうか。

PETScの型とCの型

普通、Cでベクトルや行列を扱う際は配列を使用します。
しかし、PETScではベクトルについてはVec、行列についてはMatの型を定義し計算します。
つまり、上記のmain関数で使用した A や b はそのままのdouble の配列のままでは使用できず、それぞれVecやMatの型に変換する必要があります。

また、Cのデフォルトの型と一致しているものもいくつか存在します。今回のコードでは以下の三つが該当します。

PetscInt = int
PetscReal = PetscScalar = double

おそらくfloatやcharもあると思います。

`Vec`の作成

ベクトル b の作成パートは以下の部分になります。

vector_part.c

 /* Create vectors. */
  ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) x, "Solution");CHKERRQ(ierr);
  ierr = VecSetSizes(x,PETSC_DECIDE,n);CHKERRQ(ierr);
  ierr = VecSetFromOptions(x);CHKERRQ(ierr); //must, but I don't know this function's mean 
  ierr = VecDuplicate(x,&b);CHKERRQ(ierr); //let b same size of vector as x (not copy)
  ierr = VecDuplicate(x,&u);CHKERRQ(ierr); 

   // set vector values
   VecSetValues(b,n,ix,b_input,INSERT_VALUES); //b_input --> b 

   //If it is finished to insert value in vector, these two functions are necessary.
   ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
   ierr = VecAssemblyEnd(b);CHKERRQ(ierr);

   printf("---- b in petsc -------\n");
   ierr = VecView(b,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

実際に、Vec bにdouble b_input[n]を代入しているのはVecSetValues(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)の部分になります。

引数のconst PetscScalar y[]が今回はb_inputに当たります。
PetscScalarがdoubleと等しいのでそのまま代入することができます。

ベクトルの値の設定が終わり次第、VecAssemblyBegin( ) と VecAssemblyEnd( ) を行ってください。

代入結果はVecView( )で見ることができます。

---- b in petsc -------
Vec Object: 1 MPI processes
  type: seq
1.
0.
0.
0.
0.
0.
0.
0.
0.
1.

ちゃんと代入できていますね。

[追記]
すでに、データの入った配列が用意されている場合(上記の場合に相当)にはを使うそうです。

参考 : Petsc for ISPH

`Mat`の作成

マトリックス A の作成パートは以下の部分になります。

Matrix_part.c

  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = MatSetUp(A);CHKERRQ(ierr);

  // Assemble matrix [a_ii = 2.0, a_i(j-1) = a_i(j+1) = -1.0] 
  ierr = MatSetValues(A,n,ix,n,ix,A_input,INSERT_VALUES); CHKERRQ(ierr);
  ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr);
  
  /* A is symmetric. Set symmetric flag to enable ICC/Cholesky preconditioner */
  ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
  
  /*  These two funtion is needed if assemblying matirix is finished. */
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  printf("----- A in petsc -------");
  ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

実際にMat Aにdouble A_input[n*n]の値を代入するのは、MatSetValues(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],const PetscScalar v[],InsertMode addv)になります。

おそらく、sequential matrix (MATSEQAIJやMATSEQBAIJなど) では、A_inputのように1次元配列で値を保持しているため、このように1行で値を代入することができているのだと思います。(未確認なのであまり信じないでください。)

MatSetType( ) では今回はA が　MATSEQAIJ であることを指定しています。( CG法を使用するため )

MatSetOption( ) では、A が MAT_SYMMETRIC（対称行列）であるということを指定しています。

Matについても、値を代入次第、MatAssemblyBegin( ) と MatAssemblyEnd( ) を忘れずに行ってください。

代入結果はMatView( ) で確認することができます。

----- A in petsc -------
Mat Object: 1 MPI processes
  type: seqaij
row 0: (0, 2.)  (1, -1.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)  (8, 0.)  (9, 0.) 
row 1: (0, -1.)  (1, 2.)  (2, -1.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)  (8, 0.)  (9, 0.) 
row 2: (0, 0.)  (1, -1.)  (2, 2.)  (3, -1.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)  (8, 0.)  (9, 0.) 
row 3: (0, 0.)  (1, 0.)  (2, -1.)  (3, 2.)  (4, -1.)  (5, 0.)  (6, 0.)  (7, 0.)  (8, 0.)  (9, 0.) 
row 4: (0, 0.)  (1, 0.)  (2, 0.)  (3, -1.)  (4, 2.)  (5, -1.)  (6, 0.)  (7, 0.)  (8, 0.)  (9, 0.) 
row 5: (0, 0.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, -1.)  (5, 2.)  (6, -1.)  (7, 0.)  (8, 0.)  (9, 0.) 
row 6: (0, 0.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, -1.)  (6, 2.)  (7, -1.)  (8, 0.)  (9, 0.) 
row 7: (0, 0.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, -1.)  (7, 2.)  (8, -1.)  (9, 0.) 
row 8: (0, 0.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, -1.)  (8, 2.)  (9, -1.) 
row 9: (0, 0.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)  (8, -1.)  (9, 2.)

正しく代入できていますね。
元のex1.cでMatView( )を行うと、以下の通りになります。

----- A in petsc -------
Mat Object: 1 MPI processes
  type: seqaij
row 0: (0, 2.)  (1, -1.) 
row 1: (0, -1.)  (1, 2.)  (2, -1.) 
row 2: (1, -1.)  (2, 2.)  (3, -1.) 
row 3: (2, -1.)  (3, 2.)  (4, -1.) 
row 4: (3, -1.)  (4, 2.)  (5, -1.) 
row 5: (4, -1.)  (5, 2.)  (6, -1.) 
row 6: (5, -1.)  (6, 2.)  (7, -1.) 
row 7: (6, -1.)  (7, 2.)  (8, -1.) 
row 8: (7, -1.)  (8, 2.)  (9, -1.) 
row 9: (8, -1.)  (9, 2.)

お気づきの通り、ex1.cでは0以外の値のみを代入したので、0じゃない値のみを保持しています。(こちらが理想)
しかし、今回のコードでは、0もそうでない値も直接全て代入したため 10×10＝100 個の値を保持しています。
ここから分かることは、0の値をMatに代入した場合、その0の値はそのままデータとして保持されるということです。
つまり、PETScでの疎行列の形式ではないということになります。

この疎行列の扱い方については別記事で解説します。

ソルバー`KSP`の設定

ソルバーの設定パートは以下の部分です。

ksp_part.c

 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  /*
     Set operators. Here the matrix that defines the linear system
     also serves as the matrix that defines the preconditioner.
  */
  ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr); //set A is the matrix to solve 
  
/* Set linear solver defaults for this problem (optional). */
  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
  ierr = PCSetType(pc,PCICC);CHKERRQ(ierr);
  ierr = KSPSetTolerances(ksp,1.e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); //誤差の設定
  ierr = KSPSetType(ksp,KSPCG); CHKERRQ(ierr);

  /* Set runtime options */
 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

PCSetType(PC pc,PCType type) で前処理方法の指定、
KSPSetType(KSP ksp, KSPType type)で解法の指定をします。

解ベクトル( x )の受け取り

以下の部分でx_input[]に解であるVec xの値を渡します。

x.c

/* get the x vector  (x --> x_input) */
  ierr = VecGetValues(x,n,ix,x_input);

VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])でxの値をx_inputに渡します。
使い方は、VecSetValues( )と同様です。

計算結果の確認

実行ログは以下の通りです。

$ make 
$ ./a.out
---------x before using petsc-------------
0.000000 
0.000000 
0.000000 
0.000000 
0.000000 
0.000000 
0.000000 
0.000000 
0.000000 
0.000000 
--------- b -------------
1.000000 
0.000000 
0.000000 
0.000000 
0.000000 
0.000000 
0.000000 
0.000000 
0.000000 
1.000000 

---------- A --------------
2.000000 -1.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 
-1.000000 2.000000 -1.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 
0.000000 -1.000000 2.000000 -1.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 
0.000000 0.000000 -1.000000 2.000000 -1.000000 0.000000 0.000000 0.000000 0.000000 0.000000 
0.000000 0.000000 0.000000 -1.000000 2.000000 -1.000000 0.000000 0.000000 0.000000 0.000000 
0.000000 0.000000 0.000000 0.000000 -1.000000 2.000000 -1.000000 0.000000 0.000000 0.000000 
0.000000 0.000000 0.000000 0.000000 0.000000 -1.000000 2.000000 -1.000000 0.000000 0.000000 
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -1.000000 2.000000 -1.000000 0.000000 
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -1.000000 2.000000 -1.000000 
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -1.000000 2.000000 
---- b in petsc -------
Vec Object: 1 MPI processes
  type: seq
1.
0.
0.
0.
0.
0.
0.
0.
0.
1.
----- A in petsc -------
Mat Object: 1 MPI processes
  type: seqaij
row 0: (0, 2.)  (1, -1.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)  (8, 0.)  (9, 0.) 
row 1: (0, -1.)  (1, 2.)  (2, -1.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)  (8, 0.)  (9, 0.) 
row 2: (0, 0.)  (1, -1.)  (2, 2.)  (3, -1.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)  (8, 0.)  (9, 0.) 
row 3: (0, 0.)  (1, 0.)  (2, -1.)  (3, 2.)  (4, -1.)  (5, 0.)  (6, 0.)  (7, 0.)  (8, 0.)  (9, 0.) 
row 4: (0, 0.)  (1, 0.)  (2, 0.)  (3, -1.)  (4, 2.)  (5, -1.)  (6, 0.)  (7, 0.)  (8, 0.)  (9, 0.) 
row 5: (0, 0.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, -1.)  (5, 2.)  (6, -1.)  (7, 0.)  (8, 0.)  (9, 0.) 
row 6: (0, 0.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, -1.)  (6, 2.)  (7, -1.)  (8, 0.)  (9, 0.) 
row 7: (0, 0.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, -1.)  (7, 2.)  (8, -1.)  (9, 0.) 
row 8: (0, 0.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, -1.)  (8, 2.)  (9, -1.) 
row 9: (0, 0.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)  (8, -1.)  (9, 2.) 
KSP Object: 1 MPI processes
  type: cg
  maximum iterations=10000, initial guess is zero
  tolerances:  relative=1e-05, absolute=1e-50, divergence=10000.
  left preconditioning
  using PRECONDITIONED norm type for convergence test
PC Object: 1 MPI processes
  type: icc
    out-of-place factorization
    0 levels of fill
    tolerance for zero pivot 2.22045e-14
    using Manteuffel shift [POSITIVE_DEFINITE]
    matrix ordering: natural
    factor fill ratio given 1., needed 1.
      Factored matrix follows:
        Mat Object: 1 MPI processes
          type: seqsbaij
          rows=10, cols=10
          package used to perform factorization: petsc
          total: nonzeros=55, allocated nonzeros=55
          total number of mallocs used during MatSetValues calls=0
              block size is 1
  linear system matrix = precond matrix:
  Mat Object: 1 MPI processes
    type: seqaij
    rows=10, cols=10
    total: nonzeros=100, allocated nonzeros=200
    total number of mallocs used during MatSetValues calls=10
      using I-node routines: found 2 nodes, limit used is 5
View matrix and vector 
Vec Object: Solution 1 MPI processes
  type: seq
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
Norm of error 2.22045e-16, Iterations 1
-----------x after using petsc------------
1.000000 
1.000000 
1.000000 
1.000000 
1.000000 
1.000000 
1.000000 
1.000000 
1.000000 
1.000000

とりあえず、うまくいきました！　万歳！！

さらに改善するには、、、

疎行列 Mat A を圧縮する
つまり、0以外の値だけをMat Aに代入する必要があります。
こちらについては、公式マニュアルの3.1.1章の手順にしたがう必要があります。
こちらについては次の記事で解説しようと思います。
また、こちらの記事に疎行列やその扱いについてわかりやすく書いてあります。
PCTypeやKSPTypeの変更
今回は前処理についてはICC( 不完全コレスキー分解 )、逆行列を求めるにはCG法( 共役勾配法 )を使用しました。このICCG法は A が正定値対称行列の時のみに使用することができる解法です。このように、Aの形式によって、効果的な解法は異なります。数値計算についての理解を深めるとより良い解法が見つかると思います。(これは戒めでもある。)

ただ、いずれにしろ配列を経由してから値を代入するのではなく、直接MatやVecに値を代入する方法が推奨されているようです。