5.2 testlapack3.c

/* testlapack3.c
 *     
 *     格納の仕方の練習をかねて、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)
 */

#include <stdio.h>
#include <matrix.h>

#define _AB(i,j) AB[(j)*nb+(kl+ku+(i)-(j))]

void dgbsv_(int *n, int *kl, int *ku, int *nrhs,
            double *AB, int *ldab, int *ipiv, double *b, int *ldb, int *info);

int max(int i, int j) { return (i > j) ? i : j; }
int min(int i, int j) { return (i < j) ? i : j; }

int main()
{
  int n=8, ku=2, kl=2, nb=2*kl+ku+1, ldab, ldb, nrhs, info;
  ivector ipiv;
  int i, j;
  vector x, b;
  vector AB;
  matrix a;

  /* 変数の確保 */
  AB = new_vector(nb * n);
  ipiv = new_ivector(n);
  x = new_vector(n);
  b = new_vector(n);
  a = new_matrix(n,n);

  /* 連立1次方程式の問題を作る (A と 解 x を準備して、b:= A x を計算) */
  for (i = 0; i < n; i++)
    x[i] = i + 1;
  for (j = 0; j < n; j++)
    for (i = max(0,j-ku); i <= min(n-1,j+kl); i++)
      _AB(i,j) = 10 * (i+1) + (j+1);

  for (i = 0; i < n; i++)
    b[i] = 0;
  for (j = 0; j < n; j++)
    for (i = max(0,j-ku); i <= min(n-1,j+kl); i++)
      b[i] += _AB(i,j) * x[j];

  /* A と b を表示 */
  printf("A=\n");
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++)
      if (i-kl <= j && j <= i+ku)
        printf("%4.0f ", _AB(i,j));
      else
        printf("%4.0f ", 0.0);
    printf("\n");
  }
  printf("b=\n");
  for (i = 0; i < n; i++)
    printf("%g ", b[i]);
  printf("\n");

  /*     NRHS は右辺の個数 (ここでは1)
   *     LDAB は変数 AB の「行の数」 (ここでは =nb=2*kl+ku+1)
   *     LDB は変数 B の「行の数」 (ここでは n)
   */
  nrhs = 1;
  ldab = nb;
  ldb = n;
  printf("\ncall dgbsv\n");
  dgbsv_(&n, &kl, &ku, &nrhs, AB, &ldab, ipiv, b, &ldb, &info);
  if (info == 0)
    printf(" successful\n");
  else if (info > 0)
    printf(" U is singular\n");
  else
    printf("%d-th argument has illegal value.\n", abs(info));
  printf(" return from dgbsv\n");

  /* 計算して得た解 */
  printf("\n計算して得た解=\n");
  for (i = 0; i < n; i++)
    printf("%g ", b[i]);
  printf("\n");

  printf(" 結果が 1,2,3,.. となっていたら成功\n");

  return 0;
}

testlapack3 の実行結果

oyabun% make testlapack3
gcc -W -Wall -O -I/usr/local/include -c testlapack3.c
g77 -O -o testlapack3 testlapack3.o -L/usr/local/lib -llapack -lblas -lmatrix
oyabun% ./testlapack3

A=
  11   12   13    0    0    0    0    0 
  21   22   23   24    0    0    0    0 
  31   32   33   34   35    0    0    0 
   0   42   43   44   45   46    0    0 
   0    0   53   54   55   56   57    0 
   0    0    0   64   65   66   67   68 
   0    0    0    0   75   76   77   78 
   0    0    0    0    0   86   87   88 
b=
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%