5.1 testlapack1.f

* testlapack1.f
*
*  格納の仕方の練習をかねて、dgbsv を使ってみる
*    http://www.netlib.org/cgi-bin/netlibfiles.txt?行繋ぐ
*      format=txt&blas=0&filename=lapack%2Fdouble%2Fdgbsv.f
*
*    N    方程式の個数 (=未知数の個数)
*    KL   the number of subdiagonals within the band of 某
*    KU   the number of superdiagonals within the band of 某
*    A    係数行列 (N次正方行列)
*    AB   A の帯の外を省いて詰め込んで作った行列 (2*KL+KU+1行, N列)
*         成分の対応は
*           AB(KL+KU+1+i-j,j)=A(i,j)    max(1,j-KU)≦i≦min(N,j+KL)
*    LD某 (LAPACK用語というよりLINPACK用語) Leading Dimension of 某
*         AB は配列として AB(2*KL+KU+1,N) と定義すれば十分だが
*         それよりも大きな配列として定義することを許すと便利なことがある。
*         AB(LDAB,N) (ただし LDAB≧2*KL+KU+1)
*
      program main
      implicit none
      integer n,ku,kl,nb,ldab,ldb,nrhs
      parameter (n=8,ku=2,kl=2,nb=2*kl+ku+1)
      integer i,j,p,ipiv(n),info
      real*8 a(n,n),ab(nb,n),x(n),b(n),s
      do i=1,n
         do j=1,n
            if ((i .ge. j-ku).and.(i.le.j+kl)) then
               a(i,j)=10*i+j
            else
               a(i,j)=0
            endif
         end do
      end do
      do i=1,n
         write(*,10) (a(i,j),j=1,n)
      end do
 10   format(' ',8(F5.1))
*
      do i=1,nb
         do j=1,n
            ab(i,j)=0
         end do
      end do
*
      write(*,*)
      write(*,*) ' AからABを求める公式を用いる'
      do j=1,n
         do p=max(kl+1,kl+ku+2-j),min(nb,n+kl+ku+1-j)
            ab(p,j)=a(p+j-(kl+ku+1),j)
         end do
      end do
      do i=kl+1,nb
         write(*,10) (ab(i,j),j=1,n)
      end do
*
      do i=1,nb
         do j=1,n
            ab(i,j)=0
         end do
      end do
      write(*,*)
      write(*,*) ' ABからAを求める公式を用いる'
      do j=1,n
         do i=max(1,j-ku),min(n,j+kl)
            ab(kl+ku+1+i-j,j)=a(i,j)
         end do
      end do
      do i=kl+1,nb
         write(*,10) (ab(i,j),j=1,n)
      end do
*
      do i=1,n
         x(i)=i
      end do
*
      write(*,*)
      write(*,*) ' 普通に掛け算'
      do i=1,n
         s=0
         do j=1,n
            s=s+a(i,j)*x(j)
         end do
         b(i)=s
      end do
      write(*,*) (b(i),i=1,n)
*
      write(*,*)
      write(*,*) ' 圧縮してつめた行列で掛け算'
      do i=1,n
         b(i)=0
      end do
      do j=1,n
         do i=max(1,j-ku),min(n,j+kl)
            b(i)=b(i)+ab(kl+ku+1+i-j,j)*x(j)
         end do
      end do
      write(*,*) (b(i),i=1,n)
*
*     NRHS は右辺の個数 (ここでは1)
*     LDAB は変数 AB の「行の数」 (ここでは =nb=2*kl+ku+1)
*     LDB は変数 B の「行の数」 (ここでは n)
      nrhs=1
      ldab=nb
      ldb=n
      write(*,*)
      write(*,*) ' call dgbsv'
      call dgbsv(n,kl,ku,nrhs,ab,ldab,ipiv,b,ldb,info)
      if (info .eq. 0) then
         write(*,*) '  successful'
      else if (info .gt. 0) then
         write(*,*) ' U is singular'
      else
         write(*,*) ' ', abs(i), '-th argument has illegal value.'
      endif
      write(*,*) ' return from dgbsv'
*
      write(*,*) (b(i),i=1,n)
*
      write(*,*) ' 結果が 1,2,3,.. となっていたら成功'
      end

testlapack1 の実行結果

oyabun% make

g77 -O -o testlapack1 testlapack1.f -L/usr/local/lib -llapack -lblas
oyabun% ./testlapack1

  11.0 12.0 13.0  0.0  0.0  0.0  0.0  0.0
  21.0 22.0 23.0 24.0  0.0  0.0  0.0  0.0
  31.0 32.0 33.0 34.0 35.0  0.0  0.0  0.0
   0.0 42.0 43.0 44.0 45.0 46.0  0.0  0.0
   0.0  0.0 53.0 54.0 55.0 56.0 57.0  0.0
   0.0  0.0  0.0 64.0 65.0 66.0 67.0 68.0
   0.0  0.0  0.0  0.0 75.0 76.0 77.0 78.0
   0.0  0.0  0.0  0.0  0.0 86.0 87.0 88.0

  AからABを求める公式を用いる
   0.0  0.0 13.0 24.0 35.0 46.0 57.0 68.0
   0.0 12.0 23.0 34.0 45.0 56.0 67.0 78.0
  11.0 22.0 33.0 44.0 55.0 66.0 77.0 88.0
  21.0 32.0 43.0 54.0 65.0 76.0 87.0  0.0
  31.0 42.0 53.0 64.0 75.0 86.0  0.0  0.0

  ABからAを求める公式を用いる
   0.0  0.0 13.0 24.0 35.0 46.0 57.0 68.0
   0.0 12.0 23.0 34.0 45.0 56.0 67.0 78.0
  11.0 22.0 33.0 44.0 55.0 66.0 77.0 88.0
  21.0 32.0 43.0 54.0 65.0 76.0 87.0  0.0
  31.0 42.0 53.0 64.0 75.0 86.0  0.0  0.0

  普通に掛け算
  74.  230.  505.  890.  1385.  1990.  1994.  1829.

  圧縮してつめた行列で掛け算
  74.  230.  505.  890.  1385.  1990.  1994.  1829.

  call dgbsv
   successful
  return from dgbsv
  1.  2.  3.  4.  5.  6.  7.  8.
  結果が 1,2,3,.. となっていたら成功
oyabun%