行列の上三角化によってAx=bの解ベクトルxを求めるC言語の手法について

実現したいこと

行列A,ベクトルx,bに対してAx=bを満たすようなxを求めるプログラムの作成

前提

Mac(Intel Core)のVScodeを使ってC言語のプログラミングを行っています。コンパイルエラーは出ず,表記状もエラーは起きていませんが,実行するとエラーが起きてしまいます。調べてみてもどこから手をつけていい問題なのかわからずに困っています。

発生している問題・エラーメッセージ

./a.out
n= 2

A[0][0]=1
A[0][1]=2
A[1][0]=4
A[1][1]=3
b[0]=2
b[1]=1
A=
1.0 2.0
4.0 3.0

zsh: bus error ./a.out

該当のソースコード

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

//行列表示
void PrintMatrix(double A, int n){
int i,j;
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
printf("%6.1f",A[i
n+j]);
}
printf("\n");
}
}

//ベクトル表示
void PrintVector(double *x, int n){
int i;
for (i = 0; i < n; i++) {
printf("x[%i] = %6.1f\n",i, x[i]);
}
}

//ピボット
int PIBOT(int N, double *A, int k){
double *al;
double D[N];
int p=0;
int l;
*al = A[0];

for(l=0; l<N; l++) D[l] = 0; for(l=k; l<N; l++){ D[l] = *(al+k)* *(al+k);//絶対値で最大値を探索 } //最大値を探して保存 double m = D[0]; for(int a=1; a<N; a++) if(m<D[a]){ m=D[a]; } if(m==0){ printf("zero divison!"); exit(0);//プログラム終了 } for(int b=0; b<N; b++) if(m==D[b]) p = b; return p;

}

//PIBOTで得た行の交換
void CHANGE(int N, double *A, double *b, int k, int p){
int i;
double S;
double T;

T = b[k]; b[k] = b[p]; b[p] = T; for(i=0; i<N; i++){ S = A[k+i]; A[k+i] = A[p*N+i]; A[p*N+i] = S; }

}

//i行 k~N-1列番目の要素を0に
void DASH(int N, double *A, double *b, int k, int i){
int l;
double r = A[i+k]/A[k+k];

for(l=k; l<N; l++){ A[i+l] = A[i+l] - r*A[k+l]; } b[i] = b[i] - r*b[k];

}

//Bの三角化
void DELETE(int N, double *A, double b[N]){
int k;
int p=0;
for(k=0; k<N-1; k++){
p = PIBOT(N,A,k);
CHANGE(N, A, b, k, p);
for(int i=k+1; i<N; i++){
DASH(N, A, b, k, i);
}
}
}

//解を求めてXに格納
void ANS(int N, double *A, double *b, double *x){
int i,j;
x[N-1] = b[N-1]/A[(N-1)+(N-1)];

for(i=N-2; i>=0; i--){ double S=0; for(j=i+1; j<N; j++) S = S + A[i+j] * x[j]; x[i] = (b[i]-S)/A[i+i]; }

}

int main(void) {
int N, i, j;
printf("n= ");
scanf("%d",&N);
printf("\n");

double *A; double *b; double *x; // Allocate A,b,x A = (double *)malloc(N*N*sizeof(double)); if (A==NULL) { printf("Cannot allocate memory.\n"); exit(1); } b = (double *)malloc(N*sizeof(double)); if (b==NULL) { printf("Cannot allocate memory.\n"); exit(1); } x = (double *)malloc(N*sizeof(double)); if (x==NULL) { printf("Cannot allocate memory.\n"); exit(1); } // Input A for (i = 0; i < N; i++) { for (j = 0; j < N; j++) { printf("A[%d][%d]=",i,j); scanf("%lf",&A[i*N+j]); } } // Input b for (i = 0; i < N; i++) { printf("b[%d]=",i); scanf("%lf",&b[i]); } printf("A= \n"); PrintMatrix(A, N); printf("\n"); int p = 0, k; for(k=0; k<N-1; k++){ p = PIBOT(N,A,k); CHANGE(N, A, b, k, p); for(int i=k+1; i<N; i++){ DASH(N, A, b, k, i); } } ANS(N,A,b,x); printf("Answer:\n "); PrintVector(x,N); return 0;

}

試したこと

補足情報(FW/ツールのバージョンなど)

ここにより詳細な情報を記載してください。

コメントを投稿

0 コメント