4−５．行列演算のまとめと演習

例題１．連立方程式 [SGESV]

| 3.1  1.3 -5.7 |   | x₁ |   | -1.3 |
| 1.0 -6.9  5.8 | * | x₂ | = | -0.1 |
| 3.4  7.2 -8.8 |   | x₃ |   |  1.8 |

をLAPACKを用いて解き x₁, x₂, x₃ を求めよ。

回答例：

[Fortran]
A(1,1)=3.1
A(1,2)=1.3
A(1,3)=-5.7
A(2,1)=1.0
A(2,2)=-6.9
A(2,3)=5.8
A(3,1)=3.4
A(3,2)=7.2
A(3,3)=-8.8
と行列を指定し
call SGESV(3, 1, A, 3, pivot, b, 3, ok)
でサブルーチンを呼び出す。
> f77 linear.f -o linear -L. -llapack -lblas
> ./linear
1.00000012
1.
1.

[C]
AT[1]=3.1;
AT[2]=1.0;
AT[3]=3.4;
AT[4]=1.3;
AT[5]=-6.9;
AT[6]=7.2;
AT[7]=-5.7;
AT[8]=5.8;
AT[9]=-8.8;
と列連続な１次元配列で行列を指定し、
sgesv_(&c1, &c2, AT, &c1, pivot, b, &c1, &ok);
でサブルーチンを呼び出す。サブルーチン名の後の"_"を忘れずに。
> gcc linear.c -o linear -L. -llapack -lblas -lm -lg2c
> ./linear
1.000000e+00
1.000000e+00
1.000000e+00

　

例題２．固有値問題 [ZGEEV]

| 3.1-1.8i  1.3+0.2i -5.7-4.3i |   | x₁ |    | x₁ |
| 1.0      -6.9+3.2i  5.8+2.2i | * | x₂ | = k| x₂ |
| 3.4-4.0i  7.2+2.9i -8.8+3.2i |   | x₃ |    | x₃ |

をLAPACKを用いて解き、すべての固有値kを求めよ。

回答例：

[Fortran]
complex*16 A(3,3)
A(1,1)=(3.1, -1.8)
A(1,2)=(1.3, 0.2)
A(1,3)=(-5.7, -4.3)
A(2,1)=(1.0, 0)
A(2,2)=(-6.9, 3.2)
A(2,3)=(5.8, 2.2)
A(3,1)=(3.4, -4.0)
A(3,2)=(7.2, 2.9)
A(3,3)=(-8.8, 3.2)
と複素行列を指定する。
> f77 eigen.f -o eigen -L. -llapack -lblas
> ./eigen
(1.47067795,-3.19468549)
(-0.843010085,7.51121797)
(-13.2276682,0.283467657)

[C]
double AT[2*3*3];
AT[1]=3.1;  AT[2]=-1.8;
AT[3]=1.0;  AT[4]=0.0;
AT[5]=3.4;  AT[6]=-4.0;
AT[7]=1.3;  AT[8]=0.2;
AT[9]=-6.9; AT[10]=3.2;
AT[11]=7.2; AT[12]=2.9;
AT[13]=-5.7;AT[14]=-4.3;
AT[15]=5.8; AT[16]=2.2;
AT[17]=-8.8;AT[18]=3.2;
とまず実虚が連続、次に列が連続な１次元配列で複素行列を指定する。
> gcc eigen.c -o eigen -L. -llapack -lblas -lm -lg2c
> ./eigen
1.470678 -3.194685
-0.843010 7.511218
-13.227668 0.283468
　

《演習》

1. 次の連立１次方程式の解をLAPACKを使って求めよ。ただし実数一般行列A(100x100)、実数ベクトルb(100)はファイルで与える。[ A^T（注：Cで連続して読み込めるように行と列を転置させた形でデータが入っています 。）, b ]　
   
   A * x = b
   
   ヒント１：正しくプログラミング出来ていれば、すべての解の値が誤差内で同じ値になる
   
   ヒント２：ファイルの読み込み方法の具体例
   
   [Fortran]
   Real*4 A(100,100), b(100)
   ...
   open(10,file="a.dat")
   open(11,file="b.dat")
   
   do j=1,100
   read(10,*)(a(i,j),i=1,100)
   enddo
   read(11,*)(b(i),i=1,100)
   
   close(10)
   close(11)
   
   ただし今回のようにデータが列連続で与えられている場合のみ以下でも可！！
   read(10,*)a
   read(11,*)b
   
   [C]
   double AT[N*N];
   double b[N];
   FILE *fpa,*fpb;
   ...
   fpa=fopen("a.dat","r");
   fpb=fopen("b.dat","r");
   for(i=0;i<n*n;i++){
   fscanf(fpa,"%lf",&AT[i]);
   }
   for(i=0;i<n;i++){
   fscanf(fpb,"%lf",&b[i]);
   }
   fclose(fpa);
   fclose(fpb);
   
   　
2. 上の行列Aの固有値をすべて求めよ。
   　
3. ３章の演習問題で作ったGaussの消去法によるCrank-Nicholsonのプログラムをもとにして、LAPACKを使って陰解法（連立方程式）の部分を解くように変更せよ。
   
   ヒント：一般３重対角実行列に対するサブルーチン[単精度SGTSV／倍精度DGTSV]を用いる