Skip to content

Instantly share code, notes, and snippets.

@gogotanaka
Last active January 18, 2016 15:55
Show Gist options
  • Save gogotanaka/c9acf058a346e5ddfc4c to your computer and use it in GitHub Desktop.
Save gogotanaka/c9acf058a346e5ddfc4c to your computer and use it in GitHub Desktop.
ludecomp.c
#include<stdio.h>
main(){
long i,j,N=4,NRHS=1,IPIV[N],LDA=N,LDB=N,INFO;
double A[N*N],b[N];
/* 変数の説明 N:変数の数 NRHS:右辺の列数 IPIV[N]:Pivot選択の情報
LDA:左辺の行数 LDB:右辺の行数 INFO:ルーチンが成功したかどうかを返す */
/* 配列の引数は列に対して連続になるようにする(例:1行1列の次は2行1列) */
A[0]=1;A[4]=-1;A[8]=0;A[12]=2;
A[1]=0;A[5]=1;A[9]=2;A[13]=3;
A[2]=2;A[6]=0;A[10]=-1;A[14]=3;
A[3]=2;A[7]=-3;A[11]=1;A[15]=4;
b[0]=0;
b[1]=2;
b[2]=1;
b[3]=1;
/* dgetrf_で行列をLU分解する。AにはLU分解後の行列が入る
(Lの対角成分は1なのでそれ以外をUの下三角部分に格納したもの)*/
dgetrf_(&N,&N,A,&N,IPIV,&INFO);
/*
for(i=0;i<N;i++){
for(j=0;j<N;j++)
printf("%lf ",A[i+j*N]);
printf("\n");
}
*/
/* dgetrs_によって、LU分解後の行列Aを用いて行列方程式を解く*/
dgetrs_("N",&N,&NRHS,A,&N,IPIV,b,&LDB,&INFO);
for(i=0;i<N;i++)
printf("%lf\n",b[i]);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment