前回の記事で述べた通り、疎行列 A_input をそのままPETScのMat Aに代入しても A_input[n*i+j]=0.0 の値がそのまま非ゼロの値としてA[n*i+j]に代入されてしまいます。
つまり、以下のような状態になります。
----- 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.)
目指すべき形は以下の通りです。
----- A_sparce 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.)
ということで、単純に非ゼロの値のみを代入するやり方と、そもそもA_inputをCSR記法で圧縮した状態で作成しておく方法の2通りをメモしておきます。
非ゼロの値のみを代入
ただ単純にfor文とif文を用いて、double A_input[]からMat A_sparce[]に非ゼロの値だけを代入します。
sparce.c
// Assemble matrix [a_ii = 2.0, a_i(j-1) = a_i(j+1) = -1.0]
for(int i=0; i<n; i++){
for(int j=0; j<n; j++){
if(A_input[n*i+j]==0.0)continue;
else{
ierr = MatSetValues(A_sparce,1,&i,1,&j,&A_input[n*i+j],INSERT_VALUES); CHKERRQ(ierr);
}
}
}
代入する値( A[n*i+j] )はポインタを渡します。
n×nループを回すので、非効率ですが確実です。
PCTypeやKSPTypeによっては、対角項に値が代入されていないとエラーを出すものもあるので、その場合は例え値が0でも対角項の時( i = j )には必ず値を代入するようにしましょう。
そもそもA_input[]をCSR記法で作成する
後日、記述