cygwinのlapackで連立方程式をとく。

C

cygwinのlapackで連立方程式をとく。
まず、準備として、setup.exeからlapackっぽいのをいれておく。
とりあえずのプログラム：
http://auewe.hatenablog.com/entry/2013/05/04/075717を参考にさせていただきました。
http://www.jsces.org/Issue/Journal/pdf/nakata-0411.pdf


ポイントは行列は横にみる。 N×M行列のAijは

a[(i-1)+M*(j-1)]

ただし,ladのときは違う。


簡単な例：
a(1,1)=1
a(1,2)=4
a(2,1)=5
a(2,2)=6
b(1)=3
b(2)=4

#include <stdio.h>
#include<stdlib.h>
#define N 2
int main(void)
{
  long i,j,n=N,nrhs=1,ipiv[N],lda=N,ldb=N,info;
   double* A;
   double* b;
   A=(double*)malloc(sizeof(double)*N*N);
   b=(double*)malloc(sizeof(double)*N);

   i=1;j=1;
   A[(i-1)+N*(j-1)]=1;
   i=1;j=2;
   A[(i-1)+N*(j-1)]=4;
   i=2;j=1;
   A[(i-1)+N*(j-1)]=5;
   i=2;j=2;
   A[(i-1)+N*(j-1)]=6;

    i=1;
    b[i-1]=3;
    i=2;
    b[i-1]=4;



  dgesv_(&n,&nrhs,A,&lda,ipiv,b,&ldb,&info);
  for(i=0;i<N;i++)
    printf("%lf\n",b[i]);
  return 0;
}


コンパイル

gcc -O3 solve.c -llapack

ランダムな例：
サンプルプログラム
$ cat lapack_test.c

#include <stdio.h>
#include<stdlib.h>
#define N 2000   // 2000x2000行列
int main(void)
{
  long i,j,n=N,nrhs=1,ipiv[N],lda=N,ldb=N,info;
   double* A;
   double* b;
   A=(double*)malloc(sizeof(double)*N*N);

   b=(double*)malloc(sizeof(double)*N);

 
  for(i=1;i<N*N;i++){
        A[i]=rand();
  }

  for(i=1;i<N;i++){

        b[i]=rand();
  }


  dgesv_(&n,&nrhs,A,&lda,ipiv,b,&ldb,&info);
  for(i=0;i<N;i++)
    printf("%lf\n",b[i]);
  return 0;
}

コンパイル

gcc -O3 lapack_test.c -llapack

実行

./a.exe

結果

0.346819

  • 0.144207

0.165927

  • 1.834181
  • 0.967015
  • 0.040578

1.918667

  • 0.092941
  • 2.201499

...